diff --git a/src/Transform/polynomials.jl b/src/Transform/polynomials.jl index e4e0658f..8bee986b 100644 --- a/src/Transform/polynomials.jl +++ b/src/Transform/polynomials.jl @@ -61,73 +61,71 @@ function legendre_ϕ_ψ(k) return ϕ, ψ1, ψ2 end -# function chebyshev_ϕ_ψ(k) -# ϕ_coefs = zeros(k, k) -# ϕ_2x_coefs = zeros(k, k) - -# p = Polynomial([-1, 2]) # 2x-1 -# p2 = Polynomial([-1, 4]) # 4x-1 - -# for ki in 0:(k-1) -# if ki == 0 -# ϕ_coefs[ki+1, 1:(ki+1)] .= sqrt(2/π) -# ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(4/π) -# else -# c = convert(Polynomial, gen_poly(Chebyshev, ki)) # Chebyshev of n=ki -# ϕ_coefs[ki+1, 1:(ki+1)] .= 2/sqrt(π) .* coeffs(c(p)) -# ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(2) * 2/sqrt(π) .* coeffs(c(p2)) -# end -# end - -# ϕ = [ϕ_(ϕ_coefs[i, :]) for i in 1:k] - -# k_use = 2k - -# # phi = [partial(phi_, phi_coeff[i,:]) for i in range(k)] +function chebyshev_ϕ_ψ(k) + ϕ_coefs = zeros(k, k) + ϕ_2x_coefs = zeros(k, k) + + p = Polynomial([-1, 2]) # 2x-1 + p2 = Polynomial([-1, 4]) # 4x-1 + + for ki in 0:(k-1) + if ki == 0 + ϕ_coefs[ki+1, 1:(ki+1)] .= sqrt(2/π) + ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(4/π) + else + c = convert(Polynomial, gen_poly(Chebyshev, ki)) # Chebyshev of n=ki + ϕ_coefs[ki+1, 1:(ki+1)] .= 2/sqrt(π) .* coeffs(c(p)) + ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(2) * 2/sqrt(π) .* coeffs(c(p2)) + end + end + + ϕ = [ϕ_(ϕ_coefs[i, :]) for i in 1:k] + + k_use = 2k + c = convert(Polynomial, gen_poly(Chebyshev, k_use)) + x_m = roots(c(p)) + # 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 + + ψ1_coefs = zeros(k, k) + ψ2_coefs = zeros(k, k) + + ψ1 = Array{Any,1}(undef, k) + ψ2 = Array{Any,1}(undef, k) + + for ki in 0:(k-1) + ψ1_coefs[ki+1, :] .= ϕ_2x_coefs[ki+1, :] + for i in 0:(k-1) + proj_ = sum(wm .* ϕ[i+1].(x_m) .* sqrt(2) .* ϕ[ki+1].(2*x_m)) + ψ1_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :) + ψ2_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :) + end + + for j in 0:(ki-1) + proj_ = sum(wm .* ψ1[j+1].(x_m) .* sqrt(2) .* ϕ[ki+1].(2*x_m)) + ψ1_coefs[ki+1, :] .-= proj_ .* view(ψ1_coefs, j+1, :) + ψ2_coefs[ki+1, :] .-= proj_ .* view(ψ2_coefs, j+1, :) + end + + ψ1[ki+1] = ϕ_(ψ1_coefs[ki+1,:]; lb=0., ub=0.5) + ψ2[ki+1] = ϕ_(ψ2_coefs[ki+1,:]; lb=0.5, ub=1.) -# # 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 + norm1 = sum(wm .* ψ1[ki+1].(x_m) .* ψ1[ki+1].(x_m)) + norm2 = sum(wm .* ψ2[ki+1].(x_m) .* ψ2[ki+1].(x_m)) -# # psi1_coeff = np.zeros((k, k)) -# # psi2_coeff = np.zeros((k, k)) - -# # psi1 = [[] for _ in range(k)] -# # psi2 = [[] for _ in range(k)] - -# # for ki in range(k): -# # psi1_coeff[ki,:] = phi_2x_coeff[ki,:] -# # for i in range(k): -# # proj_ = (wm * phi[i](x_m) * np.sqrt(2)* phi[ki](2*x_m)).sum() -# # psi1_coeff[ki,:] -= proj_ * phi_coeff[i,:] -# # psi2_coeff[ki,:] -= proj_ * phi_coeff[i,:] - -# # for j in range(ki): -# # proj_ = (wm * psi1[j](x_m) * np.sqrt(2) * phi[ki](2*x_m)).sum() -# # psi1_coeff[ki,:] -= proj_ * psi1_coeff[j,:] -# # psi2_coeff[ki,:] -= proj_ * psi2_coeff[j,:] - -# # psi1[ki] = partial(phi_, psi1_coeff[ki,:], lb = 0, ub = 0.5) -# # psi2[ki] = partial(phi_, psi2_coeff[ki,:], lb = 0.5, ub = 1) - -# # norm1 = (wm * psi1[ki](x_m) * psi1[ki](x_m)).sum() -# # norm2 = (wm * psi2[ki](x_m) * psi2[ki](x_m)).sum() - -# # norm_ = np.sqrt(norm1 + norm2) -# # psi1_coeff[ki,:] /= norm_ -# # psi2_coeff[ki,:] /= norm_ -# # psi1_coeff[np.abs(psi1_coeff)<1e-8] = 0 -# # psi2_coeff[np.abs(psi2_coeff)<1e-8] = 0 - -# # psi1[ki] = partial(phi_, psi1_coeff[ki,:], lb = 0, ub = 0.5+1e-16) -# # psi2[ki] = partial(phi_, psi2_coeff[ki,:], lb = 0.5+1e-16, ub = 1) + norm_ = sqrt(norm1 + norm2) + ψ1_coefs[ki+1, :] ./= norm_ + ψ2_coefs[ki+1, :] ./= norm_ + zero_out!(ψ1_coefs) + zero_out!(ψ2_coefs) -# # return phi, psi1, psi2 -# end + ψ1[ki+1] = ϕ_(ψ1_coefs[ki+1,:]; lb=0., ub=0.5+1e-16) + ψ2[ki+1] = ϕ_(ψ2_coefs[ki+1,:]; lb=0.5+1e-16, ub=1.) + end + + return ϕ, ψ1, ψ2 +end function legendre_filter(k) H0 = zeros(k, k)legendre diff --git a/src/Transform/utils.jl b/src/Transform/utils.jl index f33f894f..ae64aa18 100644 --- a/src/Transform/utils.jl +++ b/src/Transform/utils.jl @@ -1,11 +1,10 @@ -# function ϕ_(ϕ_coefs; lb::Real=0., ub::Real=1.) -# mask = -# return Polynomial(ϕ_coefs) -# end - -# def phi_(phi_c, x, lb = 0, ub = 1): -# mask = np.logical_or(xub) * 1.0 -# return np.polynomial.polynomial.Polynomial(phi_c)(x) * (1-mask) +function ϕ_(ϕ_coefs; lb::Real=0., ub::Real=1.) + function partial(x) + mask = (lb ≤ x ≤ ub) * 1. + return Polynomial(ϕ_coefs)(x) * mask + end + return partial +end function ψ(ψ1, ψ2, i, inp) mask = (inp ≤ 0.5) * 1.0 diff --git a/test/polynomials.jl b/test/polynomials.jl index 31fb81d6..09fb11d5 100644 --- a/test/polynomials.jl +++ b/test/polynomials.jl @@ -31,4 +31,37 @@ @test coeffs(ψ2[2]) ≈ [29.44486372867492, -79.67433714817852, 51.96152422707261] @test coeffs(ψ2[3]) ≈ [-29.068883707507286, 80.49844719001908, -53.665631460012115] end + + @testset "chebyshev_ϕ_ψ" begin + ϕ, ψ1, ψ2 = NeuralOperators.chebyshev_ϕ_ψ(3) + @test ϕ[1](0) ≈ 0.7978845608028654 + @test ϕ[1](1) ≈ 0.7978845608028654 + @test ϕ[1](2) ≈ 0. + @test ϕ[2](0) ≈ -1.1283791670955126 + @test ϕ[2](1) ≈ 1.1283791670955126 + @test ϕ[2](2) ≈ 0. + @test ϕ[3](0) ≈ 1.1283791670955126 + @test ϕ[3](1) ≈ 1.1283791670955126 + @test ϕ[3](2) ≈ 0. + + @test ψ1[1](0) ≈ -0.5560622352843183 + @test ψ1[1](1) ≈ 0. + @test ψ1[1](2) ≈ 0. + @test ψ1[2](0) ≈ 0.932609257876051 + @test ψ1[2](1) ≈ 0. + @test ψ1[2](2) ≈ 0. + @test ψ1[3](0) ≈ 1.0941547380212637 + @test ψ1[3](1) ≈ 0. + @test ψ1[3](2) ≈ 0. + + @test ψ2[1](0) ≈ -0. + @test ψ2[1](1) ≈ 0.5560622352843181 + @test ψ2[1](2) ≈ 0. + @test ψ2[2](0) ≈ 0. + @test ψ2[2](1) ≈ 0.9326092578760665 + @test ψ2[2](2) ≈ 0. + @test ψ2[3](0) ≈ 0. + @test ψ2[3](1) ≈ -1.0941547380212384 + @test ψ2[3](2) ≈ 0. + end end