Mercurial > repos > public > sbplib_julia
changeset 402:1936e38fe51e
Merge feature/lazy_linear_map
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Mon, 05 Oct 2020 10:45:30 +0200 |
parents | 5c10cd0ed1fe (diff) ab0a57bedaa8 (current diff) |
children | 618b7ee73b25 48d57f185f86 16dc5b19843d |
files | src/LazyTensors/lazy_tensor_operations.jl |
diffstat | 7 files changed, 50 insertions(+), 39 deletions(-) [+] |
line wrap: on
line diff
--- a/Notes.md Mon Oct 05 10:43:00 2020 +0200 +++ b/Notes.md Mon Oct 05 10:45:30 2020 +0200 @@ -1,13 +1,14 @@ # Notes ## Known size of range and domain? -It might be a good idea let tensormappings know the size of their range and domain as a constant. This probably can't be enforced on the abstract type but maybe we should write our difference operators this way. Having this as default should clean up the thinking around adjoints of boundary operators. It could also simplify getting high performance out of repeated application of regioned TensorMappings. Is there any reason to use a trait to differentiate between fixed size and unknown size? -## Test setup -Once we figure out how to organize the subpackages we should update test folders to Project. As of writing this there seems to be and issue with this approach combined with dev'ed packages so we can't do it yet. It seems that Pkg might fix this in the future. +When do we need to know the size of the range and domain? + * When indexing to provide boundschecking? + * When doing specialised computations for different parts of the range/domain? + * More? -Some steps to imporve the situation right now is to combine everything to one package and use the `@includetests` macro from [TestSetExtensions](https://github.com/ssfrr/TestSetExtensions.jl) together with `Pkg.test(test_args="...")` to selectively run tests. + Maybe if we should have dynamic sizing it could be only for the range. `domain_size` would not be implemented. And the `range_size` would be a function of a vector that the TensorMapping is applied to. ## Reasearch and thinking - [ ] Use a trait to indicate if a TensorMapping uses indices with regions. @@ -18,9 +19,11 @@ - [ ] Figure out repeated application of regioned TensorMappings. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work? - [ ] Check how the native julia doc generator works - [ ] Check if Vidars design docs fit in there - - [ ] Formalize how range_size() and domain_size() are supposed to work in TensorMappings where dim(domain) != dim(range) (add tests or document) - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters. - [ ] Specificera operatorer i TOML eller något liknande? H.. H_gamma etc.) - [ ] Dispatch in Lower() instead of the type Lower so `::Lower` instead of `::Type{Lower}` ??? Seems better unless there is some specific reason to use the type instead of the value. + - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent. + - [ ] Can we have a trait to tell if a TensorMapping is transposable? + - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
--- a/TODO.md Mon Oct 05 10:43:00 2020 +0200 +++ b/TODO.md Mon Oct 05 10:45:30 2020 +0200 @@ -9,9 +9,13 @@ - [ ] Add 1D operators (D1, D2, e, d ... ) as TensorOperators - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation. - [ ] Add a quick and simple way of running all tests for all subpackages. + - [ ] Replace getindex hack for flatteing tuples with flatten_tuple. + - [ ] Fix indexing signatures. We should make sure we are not too specific. For the "inbetween" layers we don't know what type of index is coming so we should use `I...` instead of `I::Vararg{Int,R}` + - [ ] Use `@inferred` in a lot of tests. -## Other +## Repo - [ ] Add Vidar to the authors list + - [ ] Rename repo to Sbplib.jl # Wrap up tasks - [ ] Kolla att vi har @inbounds och @propagate_inbounds på rätt ställen
--- a/src/Grids/EquidistantGrid.jl Mon Oct 05 10:43:00 2020 +0200 +++ b/src/Grids/EquidistantGrid.jl Mon Oct 05 10:45:30 2020 +0200 @@ -7,17 +7,15 @@ export EquidistantGrid struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid - size::NTuple{Dim, Int} # First coordinate direction stored first + size::NTuple{Dim, Int} limit_lower::NTuple{Dim, T} limit_upper::NTuple{Dim, T} - inverse_spacing::NTuple{Dim, T} # Reciprocal of grid spacing # General constructor function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where T @assert all(size.>0) @assert all(limit_upper.-limit_lower .!= 0) - inverse_spacing = (size.-1)./ abs.(limit_upper.-limit_lower) - return new{Dim,T}(size, limit_lower, limit_upper, inverse_spacing) + return new{Dim,T}(size, limit_lower, limit_upper) end end @@ -39,21 +37,24 @@ return length(grid.size) end -# Returns the reciprocal of the spacing of the grid -# -function inverse_spacing(grid::EquidistantGrid) - return grid.inverse_spacing -end -export inverse_spacing + +""" + spacing(grid::EquidistantGrid) -# Returns the reciprocal of the spacing of the grid -# +The spacing between the grid points of the grid. +""" +spacing(grid::EquidistantGrid) = abs.(grid.limit_upper.-grid.limit_lower)./(grid.size.-1) # TODO: Evaluate if divisions affect performance -function spacing(grid::EquidistantGrid) - return 1.0./grid.inverse_spacing -end export spacing +""" + spacing(grid::EquidistantGrid) + +The reciprocal of the spacing between the grid points of the grid. +""" +inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid) +export inverse_spacing + # Computes the points of an EquidistantGrid as an array of tuples with # the same dimension as the grid. # @@ -79,9 +80,3 @@ return EquidistantGrid(size, limit_lower, limit_upper) end export restrict - -function pointsalongdim(grid::EquidistantGrid, dim::Integer) - @assert dim<=dimension(grid) - @assert dim>0 - points = collect(range(grid.limit_lower[dim],stop=grid.limit_upper[dim],length=grid.size[dim])) -end
--- a/src/Grids/Grids.jl Mon Oct 05 10:43:00 2020 +0200 +++ b/src/Grids/Grids.jl Mon Oct 05 10:45:30 2020 +0200 @@ -14,4 +14,6 @@ include("AbstractGrid.jl") include("EquidistantGrid.jl") +# TODO: Rename AbstractGrid to Grid and move definition here. + end # module
--- a/src/LazyTensors/lazy_tensor_operations.jl Mon Oct 05 10:43:00 2020 +0200 +++ b/src/LazyTensors/lazy_tensor_operations.jl Mon Oct 05 10:45:30 2020 +0200 @@ -7,13 +7,15 @@ 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} <: LazyArray{T,R} - t::TensorMapping{T,R,D} - o::AbstractArray{T,D} +struct LazyTensorMappingApplication{T,R,D, TM<:TensorMapping{T,R,D}, AA<:AbstractArray{T,D}} <: LazyArray{T,R} + t::TM + o::AA end # TODO: Do boundschecking on creation! export LazyTensorMappingApplication +# TODO: Go through and remove unneccerary type parameters on functions + Base.:*(tm::TensorMapping{T,R,D}, o::AbstractArray{T,D}) where {T,R,D} = LazyTensorMappingApplication(tm,o) Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Index,R}) where {T,R,D} = apply(ta.t, ta.o, I...) Base.getindex(ta::LazyTensorMappingApplication{T,R,D}, I::Vararg{Int,R}) where {T,R,D} = apply(ta.t, ta.o, Index{Unknown}.(I)...)
--- a/src/SbpOperators/laplace/secondderivative.jl Mon Oct 05 10:43:00 2020 +0200 +++ b/src/SbpOperators/laplace/secondderivative.jl Mon Oct 05 10:45:30 2020 +0200 @@ -8,14 +8,13 @@ h_inv::T # The grid spacing could be included in the stencil already. Preferable? innerStencil::Stencil{T,N} closureStencils::NTuple{M,Stencil{T,K}} - parity::Parity size::NTuple{1,Int} end export SecondDerivative function SecondDerivative(grid::EquidistantGrid{1}, innerStencil, closureStencils) - h_inv = grid.inverse_spacing[1] - return SecondDerivative(h_inv, innerStencil, closureStencils, even, size(grid)) + h_inv = inverse_spacing(grid)[1] + return SecondDerivative(h_inv, innerStencil, closureStencils, size(grid)) end LazyTensors.range_size(D2::SecondDerivative) = D2.size @@ -36,7 +35,7 @@ function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, I::Index{Upper}) where T N = length(v) # TODO: Use domain_size here instead? N = domain_size(D2,size(v)) - return @inbounds D2.h_inv*D2.h_inv*Int(D2.parity)*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I)) + return @inbounds D2.h_inv*D2.h_inv*apply_stencil_backwards(D2.closureStencils[N-Int(I)+1], v, Int(I)) end function LazyTensors.apply(D2::SecondDerivative{T}, v::AbstractVector{T}, index::Index{Unknown}) where T
--- a/test/testSbpOperators.jl Mon Oct 05 10:43:00 2020 +0200 +++ b/test/testSbpOperators.jl Mon Oct 05 10:45:30 2020 +0200 @@ -4,6 +4,8 @@ using Sbplib.RegionIndices using Sbplib.LazyTensors +# TODO: Remove collects for all the tests with TensorApplications + @testset "SbpOperators" begin # @testset "apply_quadrature" begin @@ -141,12 +143,16 @@ Q = Quadrature(g, op.quadratureClosure) - v = ones(Float64, size(g)) - @test Q isa TensorMapping{T,2,2} where T @test Q' isa TensorMapping{T,2,2} where T - @test sum(collect(Q*v)) ≈ (Lx*Ly) - @test collect(Q*v) == collect(Q'*v) + + v = ones(Float64, size(g)) + @test sum(Q*v) ≈ Lx*Ly + + v = 2*ones(Float64, size(g)) + @test_broken sum(Q*v) ≈ 2*Lx*Ly + + @test Q*v == Q'*v end @testset "InverseDiagonalInnerProduct" begin @@ -175,7 +181,7 @@ @test Qinv isa TensorMapping{T,2,2} where T @test Qinv' isa TensorMapping{T,2,2} where T - @test collect(Qinv*(Q*v)) ≈ v + @test_broken collect(Qinv*(Q*v)) ≈ v @test collect(Qinv*v) == collect(Qinv'*v) end #