diff --git a/Project.toml b/Project.toml index 404ab51..3aa8c0f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ContinuumArrays" uuid = "7ae1f121-cc2c-504b-ac30-9b923412ae5c" -version = "0.19.6" +version = "0.20.0" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" @@ -38,17 +38,21 @@ Infinities = "0.1" IntervalSets = "0.7" LazyArrays = "2" Makie = "0.20, 0.21, 0.22, 0.24" -QuasiArrays = "0.12.2" +QuasiArrays = "0.13" +Random = "1.0" RecipesBase = "1.0" StaticArrays = "1.0" +StatsBase = "0.34" julia = "1.10" [extras] Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Base64", "FastTransforms", "RecipesBase", "Makie", "Test"] +test = ["Base64", "FastTransforms", "Makie", "Random", "RecipesBase", "StatsBase", "Test"] diff --git a/src/ContinuumArrays.jl b/src/ContinuumArrays.jl index 2f67c24..b0dce96 100644 --- a/src/ContinuumArrays.jl +++ b/src/ContinuumArrays.jl @@ -19,7 +19,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout, LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize, _cutdim, AbstractQuasiFill, UnionDomain, sum_size, sum_layout, _cumsum, cumsum_layout, applylayout, equals_layout, layout_broadcasted, PolynomialLayout, dot_size, - diff_layout, diff_size, AbstractQuasiVecOrMat, vec_layout + diff_layout, diff_size, AbstractQuasiVecOrMat, vec_layout, searchsortedfirst_layout import InfiniteArrays: Infinity, InfAxes import AbstractFFTs: Plan @@ -122,4 +122,6 @@ function copy(d::Dot{<:ExpansionLayout,<:ExpansionLayout,<:AbstractQuasiArray,<: c' * (P'Q) * d end +include("sort.jl") + end diff --git a/src/maps.jl b/src/maps.jl index 08625a8..2ca2efe 100644 --- a/src/maps.jl +++ b/src/maps.jl @@ -22,13 +22,13 @@ Base.union(d::Map) = axes(invmap(d),1) for find in (:findfirst, :findlast) @eval function $find(f::Base.Fix2{typeof(isequal)}, d::Map) f.x in d || return nothing - $find(isequal(invmap(d)[f.x]), union(d)) + $find(isequal(invmap(d)[f.x]), axes(d,1)) end end @eval function findall(f::Base.Fix2{typeof(isequal)}, d::Map) f.x in d || return eltype(axes(d,1))[] - findall(isequal(invmap(d)[f.x]), union(d)) + findall(isequal(invmap(d)[f.x]), axes(d,1)) end function Base.getindex(d::Map, x::Inclusion) diff --git a/src/sort.jl b/src/sort.jl new file mode 100644 index 0000000..2f9e1d8 --- /dev/null +++ b/src/sort.jl @@ -0,0 +1,26 @@ +# gives a generalization of midpoint for when `a` or `b` is infinite +function genmidpoint(a::T, b::T) where T + if isinf(a) && isinf(b) + zero(T) + elseif isinf(a) + b - 100 + elseif isinf(b) + a + 100 + else + (a+b)/2 + end +end + + +function searchsortedfirst_layout(::ExpansionLayout, f, x; iterations=47) + d = axes(f,1) + a,b = first(d), last(d) + + for k=1:iterations #TODO: decide 47 + m= genmidpoint(a,b) + (f[m] ≤ x) ? (a = m) : (b = m) + end + (a+b)/2 +end + + diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index f92feea..be455cb 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -58,6 +58,37 @@ Base.getindex(d::InvQuadraticMap, x::Number) = sqrt((x+1)/2) ContinuumArrays.invmap(::QuadraticMap{T}) where T = InvQuadraticMap{T}() ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}() +struct InfMap{T} <: Map{T} + s +end +struct InvInfMap{T} <: Map{T} + s +end + +InfMap(s=1) = InfMap{Float64}(s) +InvInfMap(s=1) = InvInfMap{Float64}(s) + +Base.getindex(m::InfMap, r::Number) = 1-2/(m.s * r+1) +Base.axes(m::InfMap{T}) where T = (Inclusion(m.s * (0..Inf)),) +Base.axes(::InvInfMap{T}) where T = (Inclusion(-1..1),) +Base.getindex(m::InvInfMap, x::Number) = m.s*( 2/(1-x) - 1) +ContinuumArrays.invmap(m::InfMap{T}) where T = InvInfMap{T}(m.s) +ContinuumArrays.invmap(m::InvInfMap{T}) where T = InfMap{T}(m.s) + +struct BiInfMap{T} <: Map{T} end +struct InvBiInfMap{T} <: Map{T} end + +BiInfMap() = BiInfMap{Float64}() +InvBiInfMap() = InvBiInfMap{Float64}() + +Base.getindex(m::BiInfMap, y::Number) = iszero(y) ? y : (-1 + sqrt(1 + 4y^2))/(2y) +Base.axes(m::BiInfMap{T}) where T = (Inclusion(-Inf..Inf),) +Base.axes(::InvBiInfMap{T}) where T = (Inclusion(-1..1),) +Base.getindex(m::InvBiInfMap, x::Number) = x/(1-x^2) +ContinuumArrays.invmap(m::BiInfMap{T}) where T = InvBiInfMap{T}() +ContinuumArrays.invmap(m::InvBiInfMap{T}) where T = BiInfMap{T}() + + struct FooDomain end struct FooBasis <: Basis{Float64} end @@ -146,6 +177,52 @@ Base.:(==)(::FooBasis, ::FooBasis) = true @test x == Inclusion(0..1) @test M \ exp.(x) ≈ T \ exp.(sqrt.((axes(T,1) .+ 1)/2)) end + + @testset "InvMap" begin + m = InfMap() + mi = InvInfMap() + @test 0.1 ∈ m + @test -0.1 ∈ m + @test 2 ∉ m + @test 2 ∈ mi + @test 0.1 ∈ mi + @test -0.1 ∉ mi + + @test m[findfirst(isequal(0.1), m)] ≈ 0.1 + @test m[findlast(isequal(0.1), m)] ≈ 0.1 + @test m[findall(isequal(0.1), m)] ≈ [0.1] + + @test m[Inclusion(0..Inf)] ≡ m + @test_throws BoundsError m[Inclusion(-1..1)] + T = Chebyshev(5) + M = T[m,:] + @test MemoryLayout(M) isa MappedBasisLayout + @test MemoryLayout(M[:,1:3]) isa MappedBasisLayout + @test M[0.1,:] ≈ T[1-2/(0.1+1),:] + x = axes(M,1) + @test x == Inclusion(0..Inf) + @test M \ exp.(-x) ≈ T \ exp.(-(2 ./ (1 .- axes(T,1)) .- 1)) + + f = M/M\(1 .- exp.(-x)) + @test f[0.1] ≈ 1 - exp(-0.1) atol=1E-2 + + @test f[searchsortedfirst(f, 0.5)] ≈ 0.5 + + M = T[InfMap(-1),:] + @test axes(M,1) == Inclusion(-Inf .. 0) + x = axes(M,1) + f = M/M\(exp.(x)) + @test f[-0.1] ≈ exp(-0.1) atol=1E-2 + @test f[searchsortedfirst(f, 0.5)] ≈ 0.5 + + M = T[BiInfMap(),:] + @test axes(M,1) == Inclusion(-Inf .. Inf) + x = axes(M,1) + f = M/M\(atan.(x)) + @test f[-0.1] ≈ atan(-0.1) atol=1E-2 + @test f[0] ≈ 0 atol=1E-10 + @test f[searchsortedfirst(f, 0.5)] ≈ 0.5 + end end @testset "Broadcasted" begin @@ -173,7 +250,7 @@ Base.:(==)(::FooBasis, ::FooBasis) = true @test (2T)'*(T*(1:5)) ≈ T'*(2T*(1:5)) ≈ T'BroadcastQuasiMatrix(*, 2, T*(1:5)) @test T' * (a .* (T * (1:5))) ≈ T' * ((a .* T) * (1:5)) @test T'BroadcastQuasiMatrix(*, 2, 2T) == 4*(T'T) - + @test LazyArrays.simplifiable(*, T', T*(1:5)) == Val(true) @test LazyArrays.simplifiable(*, T', (a .* (T * (1:5)))) == Val(true) @test LazyArrays.simplifiable(*, T', a .* T) == Val(true) @@ -209,6 +286,6 @@ Base.:(==)(::FooBasis, ::FooBasis) = true @testset "Adjoint*Basis not defined" begin @test_throws ErrorException Chebyshev(5)'LinearSpline([-1,1]) - @test_throws ErrorException FooBasis()'FooBasis() + @test_throws ErrorException FooBasis()'FooBasis() end end diff --git a/test/test_splines.jl b/test/test_splines.jl index a3b0369..2fd3e7d 100644 --- a/test/test_splines.jl +++ b/test/test_splines.jl @@ -1,8 +1,10 @@ -using ContinuumArrays, LinearAlgebra, Base64, FillArrays, QuasiArrays, BandedMatrices, BlockArrays, Test +using ContinuumArrays, LinearAlgebra, Base64, FillArrays, QuasiArrays, BandedMatrices, BlockArrays, StatsBase, Random, Test using QuasiArrays: ApplyQuasiArray, ApplyStyle, MemoryLayout, mul, MulQuasiMatrix, Vec import LazyArrays: MulStyle, LdivStyle, arguments, applied, apply, simplifiable import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, SubBasisLayout, AdjointMappedBasisLayouts, MappedBasisLayout, plan_grid_transform, weaklaplacian +Random.seed!(24543) + @testset "Splines" begin @testset "HeavisideSpline" begin H = HeavisideSpline([1,2,3]) @@ -671,4 +673,16 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, @test F[:,1:2][0.1,:] ≈ F[0.1,1:2] end end + + @testset "searchsortedfirst" begin + L = LinearSpline(range(-1,1,1000)) + f = expand(L, exp) + @test searchsortedfirst(f, 1) ≈ 0 atol=1E-6 + end + + @testset "sample" begin + H = HeavisideSpline(range(0,1,1000)) + f = expand(H, exp) + @test mean(sample(f, 1000)) ≈ 1/(ℯ-1) atol=1E-1 + end end