SpM-lab · shinaoka · Nov 17, 2024 · Nov 12, 2024 · Nov 12, 2024 · Nov 12, 2024
diff --git a/test/Project.toml b/test/Project.toml
@@ -4,4 +4,5 @@ name = "SparseIR Tests"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/poly.jl b/test/poly.jl
@@ -1,10 +1,301 @@
 using Test
 using SparseIR
 using SparseIR.LinearAlgebra
+using StableRNGs: StableRNG
 
 isdefined(Main, :sve_logistic) || include("_conftest.jl")
 
 @testset "poly.jl" begin
+    @testset "StableRNG" begin
+        # Useful for porting poly.jl to C++
+        # https://github.com/SpM-lab/SparseIR.jl/issues/51
+        rng = StableRNG(2024)
+        data = rand(rng, 3, 3)
+        knots = rand(rng, size(data, 2) + 1) |> sort
+        @test data == [
+            0.8177021060277301 0.7085670484724618 0.5033588232863977; 
+            0.3804323567786363 0.7911959541742282 0.8268504271915096; 
+            0.5425813266814807 0.38397463704084633 0.21626598379927042
+        ]
+        @test knots == [
+            0.507134318967235, 0.5766150365607372, 0.7126662232433161, 0.7357313003784003
+        ]
+        @assert issorted(knots)
+
+        drng = StableRNG(999)
+        randsymm = rand(drng, 1:10)
+        randsymm == 9
+        ddata = rand(drng, 3, 3)
+        ddata == [
+            0.5328437345518631 0.8443074122979211 0.6722336389122814; 
+            0.1799506228788046 0.6805545318460489 0.17641780726469292; 
+            0.13124858727993338 0.2193663343416914 0.7756615110113394
+        ]
+    end
+
+    @testset "PiecewiseLegendrePoly(data::Matrix, knots::Vector, l::Int)" begin
+        rng = StableRNG(2024)
+        data = rand(rng, 3, 3)
+        knots = rand(rng, size(data, 2) + 1) |> sort
+        l = 3
+
+        pwlp = SparseIR.PiecewiseLegendrePoly(data, knots, l)
+        @test pwlp.data == data
+        @test pwlp.xmin == first(knots)
+        @test pwlp.xmax == last(knots)
+        @test pwlp.knots == knots
+        @test pwlp.polyorder == size(data, 1)
+        @test pwlp.symm == 0
+    end
+
+    @testset "PiecewiseLegendrePoly(data, p::PiecewiseLegendrePoly; symm=symm(p))" begin
+        rng = StableRNG(2024)
+        data = rand(rng, 3, 3)
+        knots = rand(rng, size(data, 2) + 1) |> sort
+        l = 3
+
+        pwlp = SparseIR.PiecewiseLegendrePoly(data, knots, l)
+
+        drng = StableRNG(999)
+        randsymm = rand(drng, Int)
+        ddata = rand(drng, 3, 3)
+        ddata_pwlp = SparseIR.PiecewiseLegendrePoly(ddata, pwlp; symm=randsymm)
+
+        @test ddata_pwlp.data == ddata
+        @test ddata_pwlp.symm == randsymm
+        for n in fieldnames(SparseIR.PiecewiseLegendrePoly)
+            n === :data && continue
+            n === :symm && continue
+            @test getfield(pwlp, n) == getfield(ddata_pwlp, n)
+        end
+    end
+
+    @testset "PiecewiseLegendrePolyVector" begin
+#=
+julia> # The following data and knots are generated by
+julia> using SparseIR
+julia> sve_result = SparseIR.SVEResult(LogisticKernel(1.0))
+julia> data1 = sve_result.u[1].data
+julia> data2 = sve_result.u[2].data
+julia> data3 = sve_result.u[3].data
+julia> knots1 = sve_result.u[1].knots
+julia> knots2 = sve_result.u[2].knots
+julia> knots3 = sve_result.u[3].knots
+julia> l1 = sve_result.u[1].l
+julia> l2 = sve_result.u[2].l
+julia> l3 = sve_result.u[3].l
+=#
+        begin
+            data1 = reshape(
+                [
+                    0.49996553669802485
+                    -0.009838135710548356
+                    0.003315915376286483
+                    -2.4035906967802686e-5
+                    3.4824832610792906e-6
+                    -1.6818592059096e-8
+                    1.5530850593697272e-9
+                    -5.67191158452736e-12
+                    3.8438802553084145e-13
+                    -1.12861464373688e-15
+                    -1.4028528586225198e-16
+                    5.199431653846204e-18
+                    -3.490774002228127e-16
+                    4.339342349553959e-18
+                    -8.247505551908268e-17
+                    7.379549188001237e-19
+                    0.49996553669802485
+                    0.009838135710548356
+                    0.003315915376286483
+                    2.4035906967802686e-5
+                    3.4824832610792906e-6
+                    1.6818592059096e-8
+                    1.5530850593697272e-9
+                    5.67191158452736e-12
+                    3.8438802553084145e-13
+                    1.12861464373688e-15
+                    -1.4028528586225198e-16
+                    -5.199431653846204e-18
+                    -3.490774002228127e-16
+                    -4.339342349553959e-18
+                    -8.247505551908268e-17
+                    -7.379549188001237e-19
+                ],
+                16,
+                2,
+            )
+
+            knots1 = [-1.0, 0.0, 1.0]
+            l1 = 0
+        end
+
+        begin
+            data2 = reshape(
+                [
+                    -0.43195475509329695
+                    0.436151579050162
+                    -0.005257007544885257
+                    0.0010660519696441624
+                    -6.611545612452212e-6
+                    7.461310619506964e-7
+                    -3.2179499894475862e-9
+                    2.5166526274315926e-10
+                    -8.387341925898803e-13
+                    5.008268649326024e-14
+                    3.7750894390998034e-17
+                    -2.304983535459561e-16
+                    3.0252856483620636e-16
+                    -1.923751082183687e-16
+                    7.201014354168769e-17
+                    -3.2715804561902326e-17
+                    0.43195475509329695
+                    0.436151579050162
+                    0.005257007544885257
+                    0.0010660519696441624
+                    6.611545612452212e-6
+                    7.461310619506964e-7
+                    3.2179499894475862e-9
+                    2.5166526274315926e-10
+                    8.387341925898803e-13
+                    5.008268649326024e-14
+                    -3.7750894390998034e-17
+                    -2.304983535459561e-16
+                    -3.0252856483620636e-16
+                    -1.923751082183687e-16
+                    -7.201014354168769e-17
+                    -3.2715804561902326e-17
+                ],
+                16,
+                2,
+            )
+
+            knots2 = [-1.0, 0.0, 1.0]
+            l2 = 1
+        end
+
+        begin
+            data3 = reshape(
+                [
+                    -0.005870438661638806
+                    -0.8376202388555938
+                    0.28368166184926036
+                    -0.0029450618222246236
+                    0.0004277118923277169
+                    -2.4101642603229184e-6
+                    2.2287962786878678e-7
+                    -8.875091544426018e-10
+                    6.021488924175155e-11
+                    -1.8705305570705647e-13
+                    9.924398482443944e-15
+                    4.299521053905097e-16
+                    -1.0697019178666955e-16
+                    3.6972269778329906e-16
+                    -8.848885164903329e-17
+                    6.327687614609368e-17
+                    -0.005870438661638806
+                    0.8376202388555938
+                    0.28368166184926036
+                    0.0029450618222246236
+                    0.0004277118923277169
+                    2.4101642603229184e-6
+                    2.2287962786878678e-7
+                    8.875091544426018e-10
+                    6.021488924175155e-11
+                    1.8705305570705647e-13
+                    9.924398482443944e-15
+                    -4.299521053905097e-16
+                    -1.0697019178666955e-16
+                    -3.6972269778329906e-16
+                    -8.848885164903329e-17
+                    -6.327687614609368e-17
+                ],
+                16,
+                2,
+            )
+
+            knots3 = [-1.0, 0.0, 1.0]
+            l3 = 2
+        end
+
+        pwlp1 = SparseIR.PiecewiseLegendrePoly(data1, knots1, l1)
+        pwlp2 = SparseIR.PiecewiseLegendrePoly(data2, knots2, l2)
+        pwlp3 = SparseIR.PiecewiseLegendrePoly(data3, knots3, l3)
+
+        polys = SparseIR.PiecewiseLegendrePolyVector([pwlp1, pwlp2, pwlp3])
+
+        @test length(polys) == 3
+
+        polys = SparseIR.PiecewiseLegendrePolyVector(
+            [
+                SparseIR.PiecewiseLegendrePoly(data1, knots1, l1)
+                SparseIR.PiecewiseLegendrePoly(data2, knots2, l2)
+                SparseIR.PiecewiseLegendrePoly(data3, knots3, l3)
+            ],
+        )
+
+        @test SparseIR.xmin(polys) == SparseIR.xmin(pwlp1)
+        @test SparseIR.xmax(polys) == SparseIR.xmax(pwlp1)
+        @test SparseIR.knots(polys) == SparseIR.knots(pwlp1)
+        @test SparseIR.Δx(polys) == SparseIR.Δx(pwlp1)
+        @test SparseIR.polyorder(polys) == SparseIR.polyorder(pwlp1)
+        @test SparseIR.norms(polys) == SparseIR.norms(pwlp1)
+        @test SparseIR.symm(polys) == SparseIR.symm.([pwlp1, pwlp2, pwlp3])
+
+        @testset "polys(x::Float64)" begin
+            x = rand(StableRNG(42))
+            @test polys(x) == [pwlp1(x), pwlp2(x), pwlp3(x)]
+            @test SparseIR.data(polys) ==
+                cat(SparseIR.data.([pwlp1, pwlp2, pwlp3])..., dims = 3)
+        end
+
+        @testset "polys(x::Array)" begin
+            x = rand(StableRNG(42), 2, 1, 4)
+            tar = polys(x)
+            @test size(tar) == (3, 2, 1, 4)
+            ref = reshape(vcat([[pwlp1(e), pwlp2(e), pwlp3(e)] for e in x]...), 3, 2, 1, 4)
+            @test tar == ref
+        end
+
+        @testset "PiecewiseLegendrePolyVector(polys::PiecewiseLegendrePolyVector, knots::AbstractVector)" begin
+            new_knots = [-1.0, 0.0, 1.0]
+
+            new_polys = SparseIR.PiecewiseLegendrePolyVector(
+                polys, new_knots,
+                symm=zeros(Int, size(SparseIR.data(polys), 3)),
+            )
+
+            @test SparseIR.data(new_polys) == SparseIR.data(polys)
+            @test SparseIR.knots(new_polys) == SparseIR.knots(new_polys)
+            @test SparseIR.Δx(new_polys) == diff(new_knots)
+        end
+    end
+
+    @testset "deriv" begin
+        # independent from sve.jl
+        # https://github.com/SpM-lab/SparseIR.jl/issues/51
+        rng = StableRNG(2024)
+
+        data = rand(rng, 3, 3)
+        knots = rand(rng, size(data, 2) + 1) |> sort
+        l = 3
+        pwlp = SparseIR.PiecewiseLegendrePoly(data, knots, l)
+
+        n = 1
+        ddata = SparseIR.legder(pwlp.data, n)
+        ddata .*= pwlp.inv_xs'
+
+        deriv_pwlp = SparseIR.deriv(pwlp)
+
+        @test deriv_pwlp.data == ddata
+        @test deriv_pwlp.symm == 0
+
+        for n in fieldnames(SparseIR.PiecewiseLegendrePoly)
+            n === :data && continue
+            n === :symm && continue
+            @test getfield(pwlp, n) == getfield(deriv_pwlp, n)
+        end
+    end
+
     @testset "shape" begin
         u, s, v = SparseIR.part(sve_logistic[42])
         l = length(s)
@@ -64,6 +355,91 @@ isdefined(Main, :sve_logistic) || include("_conftest.jl")
         @test all(isapprox.(overlap(u[1], u), ref; rtol=0, atol))
     end
 
+    @testset "overlap(poly::PiecewiseLegendrePoly, f::F)" begin
+        # independent from sve.jl
+        # https://github.com/SpM-lab/SparseIR.jl/issues/51
+        rng = StableRNG(2024)
+
+        data = rand(rng, 3, 3)
+        knots = rand(rng, size(data, 2) + 1) |> sort
+        l = 3
+        pwlp = SparseIR.PiecewiseLegendrePoly(data, knots, l)
+
+        if Sys.isapple() && Sys.ARCH === :aarch64
+            # On macOS (arm64-apple-darwin22.4.0), we get
+            ∫pwlp, ∫pwlp_err = (0.4934184996836404, 8.326672684688674e-17)
+        else
+            ∫pwlp, ∫pwlp_err = (0.4934184996836403, 2.7755575615628914e-17)
+        end
+
+        @test overlap(pwlp, identity) ≈ ∫pwlp
+        @test all(overlap(pwlp, identity, return_error=true) .≈ (∫pwlp, ∫pwlp_err))
+    end
+
+    @testset "roots(poly::PiecewiseLegendrePoly; tol=1e-10, alpha=Val(2))" begin
+        # https://github.com/SpM-lab/SparseIR.jl/issues/51
+        #=
+        The following data and knots are generated by
+        julia> using SparseIR
+        julia> using Test
+        julia> Λ = 1.0
+        julia> sve_result = SparseIR.SVEResult(SparseIR.LogisticKernel(Λ))
+        julia> basis = SparseIR.FiniteTempBasis{SparseIR.Fermionic}(1, Λ; sve_result)
+        =#
+
+        data = reshape([
+            0.16774734206553019
+            0.49223680914312595
+            -0.8276728567928646
+            0.16912891046582143
+            -0.0016231275318572044
+            0.00018381683946452256
+            -9.699355027805034e-7
+            7.60144228530804e-8
+            -2.8518324490258146e-10
+            1.7090590205708293e-11
+            -5.0081401126025e-14
+            2.1244236198427895e-15
+            2.0478095258000225e-16
+            -2.676573801530628e-16
+            2.338165820094204e-16
+            -1.2050663212312096e-16
+            -0.16774734206553019
+            0.49223680914312595
+            0.8276728567928646
+            0.16912891046582143
+            0.0016231275318572044
+            0.00018381683946452256
+            9.699355027805034e-7
+            7.60144228530804e-8
+            2.8518324490258146e-10
+            1.7090590205708293e-11
+            5.0081401126025e-14
+            2.1244236198427895e-15
+            -2.0478095258000225e-16
+            -2.676573801530628e-16
+            -2.338165820094204e-16
+            -1.2050663212312096e-16
+        ], 16, 2)
+
+        knots = [0.0, 0.5, 1.0]
+        l = 3
+
+        # expected behavior
+        #=
+        julia> @test basis.u[4].data == data
+        julia> @test basis.u[4].knots == knots
+        julia> @test basis.u[4].l == l
+        =#
+
+        pwlp = SparseIR.PiecewiseLegendrePoly(data, knots, l)
+        @test SparseIR.roots(pwlp) == [
+            0.1118633448586015
+            0.4999999999999998
+            0.8881366551413985
+        ]
+    end
+
     @testset "eval unique" begin
         u, s, v = SparseIR.part(sve_logistic[42])
         û = SparseIR.PiecewiseLegendreFTVector(u, Fermionic())