diff --git a/src/polynomials.jl b/src/polynomials.jl index a7d2670..93ee257 100644 --- a/src/polynomials.jl +++ b/src/polynomials.jl @@ -165,25 +165,23 @@ function chebyshev_filter(k) Φ1 = zeros(k, k) ϕ, ψ1, ψ2 = chebyshev_ϕ_ψ(k) - # ---------------------------------------------------------- - - # x = Symbol('x') - # kUse = 2*k - # roots = Poly(chebyshevt(kUse, 2*x-1)).all_roots() - # x_m = np.array([rt.evalf(20) for rt in roots]).astype(np.float64) - # # x_m[x_m==0.5] = 0.5 + 1e-8 # add small noise to avoid the case of 0.5 belonging to both phi(2x) and phi(2x-1) - # # not needed for our purpose here, we use even k always to avoid - # wm = np.pi / kUse / 2 - - # for ki in range(k): - # for kpi in range(k): - # H0[ki, kpi] = 1/np.sqrt(2) * (wm * phi[ki](x_m/2) * phi[kpi](x_m)).sum() - # G0[ki, kpi] = 1/np.sqrt(2) * (wm * psi(psi1, psi2, ki, x_m/2) * phi[kpi](x_m)).sum() - # H1[ki, kpi] = 1/np.sqrt(2) * (wm * phi[ki]((x_m+1)/2) * phi[kpi](x_m)).sum() - # G1[ki, kpi] = 1/np.sqrt(2) * (wm * psi(psi1, psi2, ki, (x_m+1)/2) * phi[kpi](x_m)).sum() - - # PHI0[ki, kpi] = (wm * phi[ki](2*x_m) * phi[kpi](2*x_m)).sum() * 2 - # PHI1[ki, kpi] = (wm * phi[ki](2*x_m-1) * phi[kpi](2*x_m-1)).sum() * 2 + k_use = 2k + c = convert(Polynomial, gen_poly(Chebyshev, k_use)) + x_m = roots(c(Polynomial([-1, 2]))) # 2x-1 + # x_m[x_m==0.5] = 0.5 + 1e-8 # add small noise to avoid the case of 0.5 belonging to both phi(2x) and phi(2x-1) + # not needed for our purpose here, we use even k always to avoid + wm = π / k_use / 2 + + for ki in 0:(k-1) + for kpi in 0:(k-1) + H0[ki+1, kpi+1] = 1/sqrt(2) * sum(wm .* ϕ[ki+1].(x_m/2) .* ϕ[kpi+1].(x_m)) + H1[ki+1, kpi+1] = 1/sqrt(2) * sum(wm .* ϕ[ki+1].((x_m.+1)/2) .* ϕ[kpi+1].(x_m)) + G0[ki+1, kpi+1] = 1/sqrt(2) * sum(wm .* ψ(ψ1, ψ2, ki, x_m/2) .* ϕ[kpi+1].(x_m)) + G1[ki+1, kpi+1] = 1/sqrt(2) * sum(wm .* ψ(ψ1, ψ2, ki, (x_m.+1)/2) .* ϕ[kpi+1].(x_m)) + Φ0[ki+1, kpi+1] = 2*sum(wm .* ϕ[ki+1].(2x_m) .* ϕ[kpi+1].(2x_m)) + Φ1[ki+1, kpi+1] = 2*sum(wm .* ϕ[ki+1].(2 .* x_m .- 1) .* ϕ[kpi+1].(2 .* x_m .- 1)) + end + end zero_out!(H0) zero_out!(H1) diff --git a/test/polynomials.jl b/test/polynomials.jl index 1163b7e..fdb53df 100644 --- a/test/polynomials.jl +++ b/test/polynomials.jl @@ -85,25 +85,25 @@ end @testset "chebyshev_filter" begin - # H0, H1, G0, G1, Φ1, Φ2 = NeuralOperators.chebyshev_filter(3) + H0, H1, G0, G1, Φ0, Φ1 = NeuralOperators.chebyshev_filter(3) - # @test H0 ≈ [0.70710678 0. 0. ; - # -0.5 0.35355339 0. ; - # -0.25 -0.70710678 0.1767767] - # @test H1 ≈ [0.70710678 0. 0. ; - # 0.5 0.35355339 0. ; - # -0.25 0.70710678 0.1767767] - # @test G0 ≈ [0.60944614 0.77940383 0. ; - # 0.66325172 1.02726613 1.14270252; - # 0.61723435 0.90708619 1.1562954 ] - # @test G1 ≈ [-0.60944614 0.77940383 0. ; - # 0.66325172 -1.02726613 1.14270252; - # -0.61723435 0.90708619 -1.1562954 ] - # @test Φ1 ≈ [1. -0.40715364 -0.21440101; - # -0.40715364 0.84839559 -0.44820615; - # -0.21440101 -0.44820615 0.84002127] - # @test Φ2 ≈ [1. 0.40715364 -0.21440101; - # 0.40715364 0.84839559 0.44820615; - # -0.21440101 0.44820615 0.84002127] + @test H0 ≈ [0.70710678 0. 0. ; + -0.5 0.35355339 0. ; + -0.25 -0.70710678 0.1767767] + @test H1 ≈ [0.70710678 0. 0. ; + 0.5 0.35355339 0. ; + -0.25 0.70710678 0.1767767] + @test G0 ≈ [0.60944614 0.77940383 0. ; + 0.66325172 1.02726613 1.14270252; + 0.61723435 0.90708619 1.1562954 ] + @test G1 ≈ [-0.60944614 0.77940383 0. ; + 0.66325172 -1.02726613 1.14270252; + -0.61723435 0.90708619 -1.1562954 ] + @test Φ0 ≈ [1. -0.40715364 -0.21440101; + -0.40715364 0.84839559 -0.44820615; + -0.21440101 -0.44820615 0.84002127] + @test Φ1 ≈ [1. 0.40715364 -0.21440101; + 0.40715364 0.84839559 0.44820615; + -0.21440101 0.44820615 0.84002127] end end