From a6108136f12a0486c7f7ea385f464c71d9dc3590 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Sun, 17 Nov 2024 23:20:05 +0900 Subject: [PATCH] Add tests for PiecewiseLegendrePolyVector (#55) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add StableRNGs * Increase tests for poly.jl * Use ≈ (isapprox) * Format matrix data * Use all * Update reference data * Add tests for PiecewiseLegendrePolyVector --- test/Project.toml | 1 + test/poly.jl | 376 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 377 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index d680726..a2478de 100644 --- 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 index a5b1655..29bf354 100644 --- 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]) û = SparseIR.PiecewiseLegendreFTVector(u, Fermionic())