Mercurial > repos > public > sbplib_julia
diff src/LazyTensors/lazy_tensor_operations.jl @ 562:8f7919a9b398 feature/boundary_ops
Merge with default
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Mon, 30 Nov 2020 18:30:24 +0100 |
parents | a5caa934b35f 53828d3ed132 |
children | 1c512e796c6d |
line wrap: on
line diff
--- a/src/LazyTensors/lazy_tensor_operations.jl Thu Nov 26 09:03:54 2020 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Mon Nov 30 18:30:24 2020 +0100 @@ -14,10 +14,7 @@ # TODO: Do boundschecking on creation! export LazyTensorMappingApplication -# TODO: Go through and remove unneccerary type parameters on functions -Base.getindex(ta::LazyTensorMappingApplication{T,0}, I::Index) where T = apply(ta.t, ta.o, I) -Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Index,R}) where {T,R} = apply(ta.t, ta.o, I...) -Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Int,R}) where {T,R} = apply(ta.t, ta.o, Index{Unknown}.(I)...) +Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Any,R}) where {T,R} = apply(ta.t, ta.o, I...) Base.size(ta::LazyTensorMappingApplication) = range_size(ta.t) # TODO: What else is needed to implement the AbstractArray interface? @@ -50,8 +47,8 @@ Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm) Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm -apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Index,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) -apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmt.tm, v, I...) +apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) +apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmt.tm, v, I...) range_size(tmt::LazyTensorMappingTranspose) = domain_size(tmt.tm) domain_size(tmt::LazyTensorMappingTranspose) = range_size(tmt.tm) @@ -67,11 +64,11 @@ end # TODO: Boundschecking in constructor. -apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...) -apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...) +apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) + apply(tmBinOp.tm2, v, I...) +apply(tmBinOp::LazyTensorMappingBinaryOperation{:-,T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...) -range_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = range_size(tmBinOp.tm1) -domain_size(tmBinOp::LazyTensorMappingBinaryOperation{Op,T,R,D}) where {Op,T,R,D} = domain_size(tmBinOp.tm1) +range_size(tmBinOp::LazyTensorMappingBinaryOperation) = range_size(tmBinOp.tm1) +domain_size(tmBinOp::LazyTensorMappingBinaryOperation) = domain_size(tmBinOp.tm1) Base.:+(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:+,T,R,D}(tm1,tm2) Base.:-(tm1::TensorMapping{T,R,D}, tm2::TensorMapping{T,R,D}) where {T,R,D} = LazyTensorMappingBinaryOperation{:-,T,R,D}(tm1,tm2) @@ -95,11 +92,11 @@ range_size(tm::TensorMappingComposition) = range_size(tm.t1) domain_size(tm::TensorMappingComposition) = domain_size(tm.t2) -function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{S,R} where S) where {T,R,K,D} +function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,K,D} apply(c.t1, c.t2*v, I...) end -function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{S,D} where S) where {T,R,K,D} +function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,K,D} apply_transpose(c.t2, c.t1'*v, I...) end @@ -132,7 +129,7 @@ range_size(llm::LazyLinearMap) = size(llm.A)[[llm.range_indicies...]] domain_size(llm::LazyLinearMap) = size(llm.A)[[llm.domain_indicies...]] -function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Index,R}) where {T,R,D} +function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} view_index = ntuple(i->:,ndims(llm.A)) for i ∈ 1:R view_index = Base.setindex(view_index, Int(I[i]), llm.range_indicies[i]) @@ -141,7 +138,7 @@ return sum(A_view.*v) end -function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Index,D}) where {T,R,D} +function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} apply(LazyLinearMap(llm.A, llm.domain_indicies, llm.range_indicies), v, I...) end @@ -240,8 +237,6 @@ # Resolve ambiguity between the two previous methods InflatedTensorMapping(I1::IdentityMapping{T}, I2::IdentityMapping{T}) where T = InflatedTensorMapping(I1,I2,IdentityMapping{T}()) -# TODO: Implement syntax and constructors for products of different combinations of InflatedTensorMapping and IdentityMapping - # TODO: Implement some pretty printing in terms of ⊗. E.g InflatedTensorMapping(I(3),B,I(2)) -> I(3)⊗B⊗I(2) function range_size(itm::InflatedTensorMapping) @@ -261,30 +256,56 @@ end function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} - view_index, inner_index = split_index(itm, I...) + dim_before = range_dim(itm.before) + dim_domain = domain_dim(itm.tm) + dim_range = range_dim(itm.tm) + dim_after = range_dim(itm.after) + + view_index, inner_index = split_index(Val(dim_before), Val(dim_domain), Val(dim_range), Val(dim_after), I...) v_inner = view(v, view_index...) return apply(itm.tm, v_inner, inner_index...) end +function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} + dim_before = range_dim(itm.before) + dim_domain = domain_dim(itm.tm) + dim_range = range_dim(itm.tm) + dim_after = range_dim(itm.after) + + view_index, inner_index = split_index(Val(dim_before), Val(dim_range), Val(dim_domain), Val(dim_after), I...) + + v_inner = view(v, view_index...) + return apply_transpose(itm.tm, v_inner, inner_index...) +end + """ - split_index(...) + split_index(::Val{dim_before}, ::Val{dim_view}, ::Val{dim_index}, ::Val{dim_after}, I...) -Splits the multi-index into two parts. One part for the view that the inner TensorMapping acts on, and one part for indexing the result +Splits the multi-index `I` into two parts. One part which is expected to be +used as a view, and one which is expected to be used as an index. Eg. ``` -(1,2,3,4) -> (1,:,:,4), (2,3) +split_index(Val(1),Val(3),Val(2),Val(1),(1,2,3,4)) -> (1,:,:,:,4), (2,3) ``` + +`dim_view` controls how many colons are in the view, and `dim_index` controls +how many elements are extracted from the middle. +`dim_before` and `dim_after` decides the length of the index parts before and after the colons in the view index. + +Arguments should satisfy `length(I) == dim_before+B_domain+dim_after`. + +The returned values satisfy + * `length(view_index) == dim_before + dim_view + dim_after` + * `length(I_middle) == dim_index` """ -function split_index(itm::InflatedTensorMapping{T,R,D}, I::Vararg{Any,R}) where {T,R,D} - I_before = slice_tuple(I, Val(1), Val(range_dim(itm.before))) - I_after = slice_tuple(I, Val(R-range_dim(itm.after)+1), Val(R)) +function split_index(::Val{dim_before}, ::Val{dim_view}, ::Val{dim_index}, ::Val{dim_after}, I...) where {dim_before,dim_view, dim_index,dim_after} + I_before, I_middle, I_after = split_tuple(I, Val(dim_before), Val(dim_index)) - view_index = (I_before..., ntuple((i)->:,domain_dim(itm.tm))..., I_after...) - inner_index = slice_tuple(I, Val(range_dim(itm.before)+1), Val(R-range_dim(itm.after))) + view_index = (I_before..., ntuple((i)->:, dim_view)..., I_after...) - return (view_index, inner_index) + return view_index, I_middle end # TODO: Can this be replaced by something more elegant while still being type stable? 2020-10-21 @@ -302,6 +323,32 @@ end """ + split_tuple(t::Tuple{...}, ::Val{M}) where {N,M} + +Split the tuple `t` into two parts. the first part is `M` long. +E.g +``` +split_tuple((1,2,3,4),Val(3)) -> (1,2,3), (4,) +``` +""" +function split_tuple(t::NTuple{N},::Val{M}) where {N,M} + return slice_tuple(t,Val(1), Val(M)), slice_tuple(t,Val(M+1), Val(N)) +end + +""" + split_tuple(t::Tuple{...},::Val{M},::Val{K}) where {N,M,K} + +Same as `split_tuple(t::NTuple{N},::Val{M})` but splits the tuple in three parts. With the first +two parts having lenght `M` and `K`. +""" +function split_tuple(t::NTuple{N},::Val{M},::Val{K}) where {N,M,K} + p1, tail = split_tuple(t, Val(M)) + p2, p3 = split_tuple(tail, Val(K)) + return p1,p2,p3 +end + + +""" flatten_tuple(t) Takes a nested tuple and flattens the whole structure