From f40542dfb6ed83a761bb0a137133587f3eec2db2 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 19 Oct 2025 20:51:27 +0100 Subject: [PATCH 1/4] searchsortedfirst and sample --- Project.toml | 9 ++++++--- src/ContinuumArrays.jl | 4 +++- src/sort.jl | 26 ++++++++++++++++++++++++++ test/test_splines.jl | 10 +++++++++- 4 files changed, 44 insertions(+), 5 deletions(-) create mode 100644 src/sort.jl diff --git a/Project.toml b/Project.toml index 404ab51..89e8721 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,7 +38,8 @@ 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" julia = "1.10" @@ -47,8 +48,10 @@ julia = "1.10" 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/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_splines.jl b/test/test_splines.jl index a3b0369..1ec2f7b 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,10 @@ import ContinuumArrays: basis, AdjointBasisLayout, ExpansionLayout, BasisLayout, @test F[:,1:2][0.1,:] ≈ F[0.1,1:2] end 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 From 49e33567df744afb0be22e4dc30dde6c47ea031b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 20 Oct 2025 00:16:55 +0100 Subject: [PATCH 2/4] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 89e8721..3aa8c0f 100644 --- a/Project.toml +++ b/Project.toml @@ -42,6 +42,7 @@ QuasiArrays = "0.13" Random = "1.0" RecipesBase = "1.0" StaticArrays = "1.0" +StatsBase = "0.34" julia = "1.10" [extras] From 17611dbe89d0f1f866ded5706a06b1aeab3389ef Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 20 Oct 2025 12:02:14 +0100 Subject: [PATCH 3/4] =?UTF-8?q?test=20with=20=E2=88=9E=20domains?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/maps.jl | 4 +-- test/test_chebyshev.jl | 59 ++++++++++++++++++++++++++++++++++++++++-- test/test_splines.jl | 6 +++++ 3 files changed, 65 insertions(+), 4 deletions(-) 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/test/test_chebyshev.jl b/test/test_chebyshev.jl index f92feea..064696d 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -58,6 +58,23 @@ 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 FooDomain end struct FooBasis <: Basis{Float64} end @@ -146,6 +163,44 @@ 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 + end end @testset "Broadcasted" begin @@ -173,7 +228,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 +264,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 1ec2f7b..2fd3e7d 100644 --- a/test/test_splines.jl +++ b/test/test_splines.jl @@ -674,6 +674,12 @@ Random.seed!(24543) 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) From d77d90eaf4fac82a9517e427f822c5dfc55a6d63 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 20 Oct 2025 12:26:38 +0100 Subject: [PATCH 4/4] BiInfMap --- test/test_chebyshev.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index 064696d..be455cb 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -75,6 +75,20 @@ 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 @@ -200,6 +214,14 @@ Base.:(==)(::FooBasis, ::FooBasis) = true 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