Mercurial > repos > public > sbplib_julia
changeset 960:e79debd10f7d feature/variable_derivatives
Merge feature/tensormapping_application_promotion
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Mon, 14 Mar 2022 10:14:38 +0100 |
parents | e9752c1e92f8 (current diff) 86889fc5b63f (diff) |
children | bc12be1b1ae5 |
files | src/LazyTensors/lazy_tensor_operations.jl src/SbpOperators/boundaryops/boundary_operator.jl src/SbpOperators/volumeops/constant_interior_scaling_operator.jl src/SbpOperators/volumeops/volume_operator.jl test/SbpOperators/boundaryops/boundary_operator_test.jl |
diffstat | 9 files changed, 106 insertions(+), 45 deletions(-) [+] |
line wrap: on
line diff
diff -r e9752c1e92f8 -r e79debd10f7d src/LazyTensors/lazy_tensor_operations.jl --- a/src/LazyTensors/lazy_tensor_operations.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Mon Mar 14 10:14:38 2022 +0100 @@ -1,5 +1,15 @@ using Sbplib.RegionIndices +export LazyTensorMappingApplication +export LazyTensorMappingTranspose +export TensorMappingComposition +export LazyLinearMap +export IdentityMapping +export InflatedTensorMapping +export LazyOuterProduct +export ⊗ +export SizeMismatch + """ LazyTensorMappingApplication{T,R,D} <: LazyArray{T,R} @@ -9,12 +19,16 @@ With a mapping `m` and a vector `v` the LazyTensorMappingApplication object can be created by `m*v`. The actual result will be calcualted when indexing into `m*v`. """ -struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{T,R,D}, AA<:AbstractArray{T,D}} <: LazyArray{T,R} +struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R} t::TM o::AA + + function LazyTensorMappingApplication(t::TensorMapping{<:Any,R,D}, o::AbstractArray{<:Any,D}) where {R,D} + T = promote_type(eltype(t), eltype(o)) + return new{T,R,D,typeof(t), typeof(o)}(t,o) + end end # TODO: Do boundschecking on creation! -export LazyTensorMappingApplication Base.getindex(ta::LazyTensorMappingApplication{T,R}, I::Vararg{Any,R}) where {T,R} = apply(ta.t, ta.o, I...) Base.getindex(ta::LazyTensorMappingApplication{T,1}, I::CartesianIndex{1}) where {T} = apply(ta.t, ta.o, I.I...) # Would otherwise be caught in the previous method. @@ -43,15 +57,14 @@ struct LazyTensorMappingTranspose{T,R,D, TM<:TensorMapping{T,R,D}} <: TensorMapping{T,D,R} tm::TM end -export LazyTensorMappingTranspose # # TBD: Should this be implemented on a type by type basis or through a trait to provide earlier errors? # Jonatan 2020-09-25: Is the problem that you can take the transpose of any TensorMapping even if it doesn't implement `apply_transpose`? Base.adjoint(tm::TensorMapping) = LazyTensorMappingTranspose(tm) Base.adjoint(tmt::LazyTensorMappingTranspose) = tmt.tm -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...) +apply(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} = apply_transpose(tmt.tm, v, I...) +apply_transpose(tmt::LazyTensorMappingTranspose{T,R,D}, v::AbstractArray{<:Any,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,8 +80,8 @@ end # TODO: Boundschecking in constructor. -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...) +apply(tmBinOp::LazyTensorMappingBinaryOperation{:+,T,R,D}, v::AbstractArray{<:Any,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{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} = apply(tmBinOp.tm1, v, I...) - apply(tmBinOp.tm2, v, I...) range_size(tmBinOp::LazyTensorMappingBinaryOperation) = range_size(tmBinOp.tm1) domain_size(tmBinOp::LazyTensorMappingBinaryOperation) = domain_size(tmBinOp.tm1) @@ -90,16 +103,15 @@ return new{T,R,K,D, typeof(t1), typeof(t2)}(t1,t2) end end -export TensorMappingComposition 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{Any,R}) where {T,R,K,D} +function apply(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,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{Any,D}) where {T,R,K,D} +function apply_transpose(c::TensorMappingComposition{T,R,K,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,K,D} apply_transpose(c.t2, c.t1'*v, I...) end @@ -127,12 +139,11 @@ return new{T,R,D,RD,AA}(A,range_indicies,domain_indicies) end end -export LazyLinearMap 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{Any,R}) where {T,R,D} +function apply(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,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 +152,7 @@ return sum(A_view.*v) end -function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{T,R}, I::Vararg{Any,D}) where {T,R,D} +function apply_transpose(llm::LazyLinearMap{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg{Any,D}) where {T,R,D} apply(LazyLinearMap(llm.A, llm.domain_indicies, llm.range_indicies), v, I...) end @@ -155,7 +166,6 @@ struct IdentityMapping{T,D} <: TensorMapping{T,D,D} size::NTuple{D,Int} end -export IdentityMapping IdentityMapping{T}(size::NTuple{D,Int}) where {T,D} = IdentityMapping{T,D}(size) IdentityMapping{T}(size::Vararg{Int,D}) where {T,D} = IdentityMapping{T,D}(size) @@ -164,8 +174,8 @@ range_size(tmi::IdentityMapping) = tmi.size domain_size(tmi::IdentityMapping) = tmi.size -apply(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...] -apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{T,D}, I::Vararg{Any,D}) where {T,D} = v[I...] +apply(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...] +apply_transpose(tmi::IdentityMapping{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = v[I...] """ Base.:∘(tm, tmi) @@ -212,7 +222,6 @@ return new{T,R,D,D_before,R_middle,D_middle,D_after, typeof(tm)}(before, tm, after) end end -export InflatedTensorMapping """ InflatedTensorMapping(before, tm, after) InflatedTensorMapping(before,tm) @@ -258,7 +267,7 @@ ) end -function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg{Any,R}) where {T,R,D} +function apply(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,R}) where {T,R,D} dim_before = range_dim(itm.before) dim_domain = domain_dim(itm.tm) dim_range = range_dim(itm.tm) @@ -270,7 +279,7 @@ 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} +function apply_transpose(itm::InflatedTensorMapping{T,R,D}, v::AbstractArray{<:Any,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) @@ -398,7 +407,6 @@ ``` """ function LazyOuterProduct end -export LazyOuterProduct function LazyOuterProduct(tm1::TensorMapping{T}, tm2::TensorMapping{T}) where T itm1 = InflatedTensorMapping(tm1, IdentityMapping{T}(range_size(tm2))) @@ -414,7 +422,6 @@ LazyOuterProduct(tms::Vararg{TensorMapping}) = foldl(LazyOuterProduct, tms) ⊗(a::TensorMapping, b::TensorMapping) = LazyOuterProduct(a,b) -export ⊗ function check_domain_size(tm::TensorMapping, sz) @@ -427,7 +434,6 @@ tm::TensorMapping sz end -export SizeMismatch function Base.showerror(io::IO, err::SizeMismatch) print(io, "SizeMismatch: ")
diff -r e9752c1e92f8 -r e79debd10f7d src/LazyTensors/tensor_mapping.jl --- a/src/LazyTensors/tensor_mapping.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/src/LazyTensors/tensor_mapping.jl Mon Mar 14 10:14:38 2022 +0100 @@ -1,10 +1,16 @@ +export TensorMapping +export apply +export apply_transpose +export range_dim, domain_dim +export range_size, domain_size + """ TensorMapping{T,R,D} Describes a mapping of a `D` dimension tensor to an `R` dimension tensor. The action of the mapping is implemented through the method ```julia - apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T} + apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T} ``` The size of the range and domain that the operator works with should be returned by @@ -21,23 +27,20 @@ ``` """ abstract type TensorMapping{T,R,D} end -export TensorMapping """ - apply(t::TensorMapping{T,R,D}, v::AbstractArray{T,D}, I::Vararg) where {R,D,T} + apply(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,D}, I::Vararg) where {R,D,T} Return the result of the mapping for a given index. """ function apply end -export apply """ - apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{T,R}, I::Vararg) where {R,D,T} + apply_transpose(t::TensorMapping{T,R,D}, v::AbstractArray{<:Any,R}, I::Vararg) where {R,D,T} Return the result of the transposed mapping for a given index. """ function apply_transpose end -export apply_transpose """ range_dim(::TensorMapping) @@ -51,7 +54,6 @@ """ domain_dim(::TensorMapping{T,R,D}) where {T,R,D} = D -export range_dim, domain_dim """ range_size(M::TensorMapping) @@ -67,7 +69,6 @@ """ function domain_size end -export range_size, domain_size """ eltype(::TensorMapping{T})
diff -r e9752c1e92f8 -r e79debd10f7d src/SbpOperators/boundaryops/boundary_operator.jl --- a/src/SbpOperators/boundaryops/boundary_operator.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Mon Mar 14 10:14:38 2022 +0100 @@ -62,27 +62,27 @@ LazyTensors.range_size(op::BoundaryOperator) = () LazyTensors.domain_size(op::BoundaryOperator) = (op.size,) -function LazyTensors.apply(op::BoundaryOperator{T,Lower}, v::AbstractVector{T}) where T +function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector) apply_stencil(op.stencil,v,1) end -function LazyTensors.apply(op::BoundaryOperator{T,Upper}, v::AbstractVector{T}) where T +function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector) apply_stencil_backwards(op.stencil,v,op.size) end -function LazyTensors.apply_transpose(op::BoundaryOperator{T,Lower}, v::AbstractArray{T,0}, i::Index{Lower}) where T +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower}) return op.stencil[Int(i)-1]*v[] end -function LazyTensors.apply_transpose(op::BoundaryOperator{T,Upper}, v::AbstractArray{T,0}, i::Index{Upper}) where T +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper}) return op.stencil[op.size[1] - Int(i)]*v[] end # Catch all combinations of Lower, Upper and Interior not caught by the two previous methods. -function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i::Index) where T - return zero(T) +function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i::Index) + return zero(eltype(v)) end -function LazyTensors.apply_transpose(op::BoundaryOperator{T}, v::AbstractArray{T,0}, i) where T +function LazyTensors.apply_transpose(op::BoundaryOperator, v::AbstractArray{<:Any,0}, i) return LazyTensors.apply_transpose_with_region(op, v, closure_size(op), op.size[1], i) end
diff -r e9752c1e92f8 -r e79debd10f7d src/SbpOperators/volumeops/constant_interior_scaling_operator.jl --- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Mon Mar 14 10:14:38 2022 +0100 @@ -28,19 +28,19 @@ LazyTensors.domain_size(op::ConstantInteriorScalingOperator) = (op.size,) # TBD: @inbounds in apply methods? -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Lower}) return op.closure_weights[Int(i)]*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Interior}) return op.interior_weight*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i::Index{Upper}) return op.closure_weights[op.size[1]-Int(i)+1]*v[Int(i)] end -function LazyTensors.apply(op::ConstantInteriorScalingOperator{T}, v::AbstractVector{T}, i) where T +function LazyTensors.apply(op::ConstantInteriorScalingOperator, v::AbstractVector, i) return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i) end
diff -r e9752c1e92f8 -r e79debd10f7d src/SbpOperators/volumeops/volume_operator.jl --- a/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Mon Mar 14 10:14:38 2022 +0100 @@ -43,18 +43,18 @@ LazyTensors.range_size(op::VolumeOperator) = op.size LazyTensors.domain_size(op::VolumeOperator) = op.size -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Lower}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Lower}) return @inbounds apply_stencil(op.closure_stencils[Int(i)], v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Interior}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Interior}) return apply_stencil(op.inner_stencil, v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i::Index{Upper}) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i::Index{Upper}) return @inbounds Int(op.parity)*apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], v, Int(i)) end -function LazyTensors.apply(op::VolumeOperator{T}, v::AbstractVector{T}, i) where T +function LazyTensors.apply(op::VolumeOperator, v::AbstractVector, i) return LazyTensors.apply_with_region(op, v, closure_size(op), op.size[1], i) end
diff -r e9752c1e92f8 -r e79debd10f7d test/LazyTensors/lazy_tensor_operations_test.jl --- a/test/LazyTensors/lazy_tensor_operations_test.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Mon Mar 14 10:14:38 2022 +0100 @@ -80,6 +80,28 @@ v = [[1 2];[3 4]] @test m*v == [[2 4];[6 8]] @test (m*v)[2,1] == 6 + + @testset "Promotion" begin + m = ScalingOperator{Int,1}(2,(3,)) + v = [1.,2.,3.] + @test m*v isa AbstractVector{Float64} + @test m*v == [2.,4.,6.] + + m = ScalingOperator{Int,2}(2,(2,2)) + v = [[1. 2.];[3. 4.]] + @test m*v == [[2. 4.];[6. 8.]] + @test (m*v)[2,1] == 6. + + m = ScalingOperator{ComplexF64,1}(2. +2. *im,(3,)) + v = [1.,2.,3.] + @test m*v isa AbstractVector{ComplexF64} + @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + + m = ScalingOperator{ComplexF64,1}(1,(3,)) + v = [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + @test m*v isa AbstractVector{ComplexF64} + @test m*v == [2. + 2. *im, 4. + 4. *im, 6. + 6. *im] + end end @testset "TensorMapping binary operations" begin @@ -107,6 +129,8 @@ @test range_size(A+B) == range_size(A) == range_size(B) @test domain_size(A+B) == domain_size(A) == domain_size(B) + + @test ((A+B)*ComplexF64[1.1,1.2,1.3])[3] isa ComplexF64 end @@ -129,6 +153,9 @@ v = rand(2) @test (Ã∘B̃)'*v ≈ B'*A'*v rtol=1e-14 + + @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64 + @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64 end @testset "LazyLinearMap" begin @@ -194,6 +221,10 @@ @test I*v == v @test I'*v == v + v = rand(ComplexF64,sz...) + @test I*v == v + @test I'*v == v + @test range_size(I) == sz @test domain_size(I) == sz end @@ -322,6 +353,16 @@ end end + @testset "application to other type" begin + tm = InflatedTensorMapping(I(3,2), A, I(4)) + + v = rand(ComplexF64, domain_size(tm)...) + @test (tm*v)[1,2,3,1] isa ComplexF64 + + v = rand(ComplexF64, domain_size(tm')...) + @test (tm'*v)[1,2,2,1] isa ComplexF64 + end + @testset "Inference of application" begin struct ScalingOperator{T,D} <: TensorMapping{T,D,D} λ::T
diff -r e9752c1e92f8 -r e79debd10f7d test/SbpOperators/boundaryops/boundary_operator_test.jl --- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon Mar 14 10:14:38 2022 +0100 @@ -72,6 +72,14 @@ @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] + + v = evalOn(g_1D, x->1. +x*im) + @test (op_l*v)[] isa ComplexF64 + + u = fill(1. +im) + @test (op_l'*u)[1] isa ComplexF64 + @test (op_l'*u)[5] isa ComplexF64 + @test (op_l'*u)[11] isa ComplexF64 end @testset "2D" begin
diff -r e9752c1e92f8 -r e79debd10f7d test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl --- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Mon Mar 14 10:14:38 2022 +0100 @@ -24,6 +24,9 @@ @test a*v == [.1,.2,.5,.5,.5,.2,.1] @test a'*v == [.1,.2,.5,.5,.5,.2,.1] + @test (a*rand(ComplexF64, domain_size(a)... ))[1] isa ComplexF64 + @test (a'*rand(ComplexF64, domain_size(a')...))[1] isa ComplexF64 + @test range_size(a) == (7,) @test domain_size(a) == (7,)
diff -r e9752c1e92f8 -r e79debd10f7d test/SbpOperators/volumeops/volume_operator_test.jl --- a/test/SbpOperators/volumeops/volume_operator_test.jl Mon Mar 14 09:46:22 2022 +0100 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Mon Mar 14 10:14:38 2022 +0100 @@ -89,6 +89,8 @@ @testset "Application" begin @test op_x*v ≈ rx rtol = 1e-14 @test op_y*v ≈ ry rtol = 1e-14 + + @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64 end @testset "Regions" begin