searchsortedfirst and sample (#208)

* searchsortedfirst and sample

* Update Project.toml

* test with ∞ domains

* BiInfMap

Author: dlfivefifty

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"]

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

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)

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
+
+

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

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