Skip to content
This repository was archived by the owner on Sep 28, 2024. It is now read-only.

Commit 60f1faa

Browse files
yuehhuayuehhua
committed
complete chebyshev polynomials
1 parent 66670c3 commit 60f1faa

File tree

3 files changed

+102
-72
lines changed

3 files changed

+102
-72
lines changed

‎src/Transform/polynomials.jl

Lines changed: 62 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -61,73 +61,71 @@ function legendre_ϕ_ψ(k)
6161
return ϕ, ψ1, ψ2
6262
end
6363

64-
# function chebyshev_ϕ_ψ(k)
65-
# ϕ_coefs = zeros(k, k)
66-
# ϕ_2x_coefs = zeros(k, k)
67-
68-
# p = Polynomial([-1, 2]) # 2x-1
69-
# p2 = Polynomial([-1, 4]) # 4x-1
70-
71-
# for ki in 0:(k-1)
72-
# if ki == 0
73-
# ϕ_coefs[ki+1, 1:(ki+1)] .= sqrt(2/π)
74-
# ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(4/π)
75-
# else
76-
# c = convert(Polynomial, gen_poly(Chebyshev, ki)) # Chebyshev of n=ki
77-
# ϕ_coefs[ki+1, 1:(ki+1)] .= 2/sqrt(π) .* coeffs(c(p))
78-
# ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(2) * 2/sqrt(π) .* coeffs(c(p2))
79-
# end
80-
# end
81-
82-
# ϕ = [ϕ_(ϕ_coefs[i, :]) for i in 1:k]
83-
84-
# k_use = 2k
85-
86-
# # phi = [partial(phi_, phi_coeff[i,:]) for i in range(k)]
64+
function chebyshev_ϕ_ψ(k)
65+
ϕ_coefs = zeros(k, k)
66+
ϕ_2x_coefs = zeros(k, k)
67+
68+
p = Polynomial([-1, 2]) # 2x-1
69+
p2 = Polynomial([-1, 4]) # 4x-1
70+
71+
for ki in 0:(k-1)
72+
if ki == 0
73+
ϕ_coefs[ki+1, 1:(ki+1)] .= sqrt(2/π)
74+
ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(4/π)
75+
else
76+
c = convert(Polynomial, gen_poly(Chebyshev, ki)) # Chebyshev of n=ki
77+
ϕ_coefs[ki+1, 1:(ki+1)] .= 2/sqrt(π) .* coeffs(c(p))
78+
ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(2) * 2/sqrt(π) .* coeffs(c(p2))
79+
end
80+
end
81+
82+
ϕ = [ϕ_(ϕ_coefs[i, :]) for i in 1:k]
83+
84+
k_use = 2k
85+
c = convert(Polynomial, gen_poly(Chebyshev, k_use))
86+
x_m = roots(c(p))
87+
# 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)
88+
# not needed for our purpose here, we use even k always to avoid
89+
wm = π / k_use / 2
90+
91+
ψ1_coefs = zeros(k, k)
92+
ψ2_coefs = zeros(k, k)
93+
94+
ψ1 = Array{Any,1}(undef, k)
95+
ψ2 = Array{Any,1}(undef, k)
96+
97+
for ki in 0:(k-1)
98+
ψ1_coefs[ki+1, :] .= ϕ_2x_coefs[ki+1, :]
99+
for i in 0:(k-1)
100+
proj_ = sum(wm .* ϕ[i+1].(x_m) .* sqrt(2) .* ϕ[ki+1].(2*x_m))
101+
ψ1_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :)
102+
ψ2_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :)
103+
end
104+
105+
for j in 0:(ki-1)
106+
proj_ = sum(wm .* ψ1[j+1].(x_m) .* sqrt(2) .* ϕ[ki+1].(2*x_m))
107+
ψ1_coefs[ki+1, :] .-= proj_ .* view(ψ1_coefs, j+1, :)
108+
ψ2_coefs[ki+1, :] .-= proj_ .* view(ψ2_coefs, j+1, :)
109+
end
110+
111+
ψ1[ki+1] = ϕ_(ψ1_coefs[ki+1,:]; lb=0., ub=0.5)
112+
ψ2[ki+1] = ϕ_(ψ2_coefs[ki+1,:]; lb=0.5, ub=1.)
87113

88-
# # x = Symbol('x')
89-
# # kUse = 2*k
90-
# # roots = Poly(chebyshevt(kUse, 2*x-1)).all_roots()
91-
# # x_m = np.array([rt.evalf(20) for rt in roots]).astype(np.float64)
92-
# # # 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)
93-
# # # not needed for our purpose here, we use even k always to avoid
94-
# # wm = np.pi / kUse / 2
114+
norm1 = sum(wm .* ψ1[ki+1].(x_m) .* ψ1[ki+1].(x_m))
115+
norm2 = sum(wm .* ψ2[ki+1].(x_m) .* ψ2[ki+1].(x_m))
95116

96-
# # psi1_coeff = np.zeros((k, k))
97-
# # psi2_coeff = np.zeros((k, k))
98-
99-
# # psi1 = [[] for _ in range(k)]
100-
# # psi2 = [[] for _ in range(k)]
101-
102-
# # for ki in range(k):
103-
# # psi1_coeff[ki,:] = phi_2x_coeff[ki,:]
104-
# # for i in range(k):
105-
# # proj_ = (wm * phi[i](x_m) * np.sqrt(2)* phi[ki](2*x_m)).sum()
106-
# # psi1_coeff[ki,:] -= proj_ * phi_coeff[i,:]
107-
# # psi2_coeff[ki,:] -= proj_ * phi_coeff[i,:]
108-
109-
# # for j in range(ki):
110-
# # proj_ = (wm * psi1[j](x_m) * np.sqrt(2) * phi[ki](2*x_m)).sum()
111-
# # psi1_coeff[ki,:] -= proj_ * psi1_coeff[j,:]
112-
# # psi2_coeff[ki,:] -= proj_ * psi2_coeff[j,:]
113-
114-
# # psi1[ki] = partial(phi_, psi1_coeff[ki,:], lb = 0, ub = 0.5)
115-
# # psi2[ki] = partial(phi_, psi2_coeff[ki,:], lb = 0.5, ub = 1)
116-
117-
# # norm1 = (wm * psi1[ki](x_m) * psi1[ki](x_m)).sum()
118-
# # norm2 = (wm * psi2[ki](x_m) * psi2[ki](x_m)).sum()
119-
120-
# # norm_ = np.sqrt(norm1 + norm2)
121-
# # psi1_coeff[ki,:] /= norm_
122-
# # psi2_coeff[ki,:] /= norm_
123-
# # psi1_coeff[np.abs(psi1_coeff)<1e-8] = 0
124-
# # psi2_coeff[np.abs(psi2_coeff)<1e-8] = 0
125-
126-
# # psi1[ki] = partial(phi_, psi1_coeff[ki,:], lb = 0, ub = 0.5+1e-16)
127-
# # psi2[ki] = partial(phi_, psi2_coeff[ki,:], lb = 0.5+1e-16, ub = 1)
117+
norm_ = sqrt(norm1 + norm2)
118+
ψ1_coefs[ki+1, :] ./= norm_
119+
ψ2_coefs[ki+1, :] ./= norm_
120+
zero_out!(ψ1_coefs)
121+
zero_out!(ψ2_coefs)
128122

129-
# # return phi, psi1, psi2
130-
# end
123+
ψ1[ki+1] = ϕ_(ψ1_coefs[ki+1,:]; lb=0., ub=0.5+1e-16)
124+
ψ2[ki+1] = ϕ_(ψ2_coefs[ki+1,:]; lb=0.5+1e-16, ub=1.)
125+
end
126+
127+
return ϕ, ψ1, ψ2
128+
end
131129

132130
function legendre_filter(k)
133131
H0 = zeros(k, k)legendre

‎src/Transform/utils.jl

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,10 @@
1-
# function ϕ_(ϕ_coefs; lb::Real=0., ub::Real=1.)
2-
# mask =
3-
# return Polynomial(ϕ_coefs)
4-
# end
5-
6-
# def phi_(phi_c, x, lb = 0, ub = 1):
7-
# mask = np.logical_or(x<lb, x>ub) * 1.0
8-
# return np.polynomial.polynomial.Polynomial(phi_c)(x) * (1-mask)
1+
function ϕ_(ϕ_coefs; lb::Real=0., ub::Real=1.)
2+
function partial(x)
3+
mask = (lb ≤ x ≤ ub) * 1.
4+
return Polynomial(ϕ_coefs)(x) * mask
5+
end
6+
return partial
7+
end
98

109
function ψ(ψ1, ψ2, i, inp)
1110
mask = (inp ≤ 0.5) * 1.0

‎test/polynomials.jl

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,4 +31,37 @@
3131
@test coeffs(ψ2[2]) ≈ [29.44486372867492, -79.67433714817852, 51.96152422707261]
3232
@test coeffs(ψ2[3]) ≈ [-29.068883707507286, 80.49844719001908, -53.665631460012115]
3333
end
34+
35+
@testset "chebyshev_ϕ_ψ" begin
36+
ϕ, ψ1, ψ2 = NeuralOperators.chebyshev_ϕ_ψ(3)
37+
@test ϕ[1](0) ≈ 0.7978845608028654
38+
@test ϕ[1](1) ≈ 0.7978845608028654
39+
@test ϕ[1](2) ≈ 0.
40+
@test ϕ[2](0) ≈ -1.1283791670955126
41+
@test ϕ[2](1) ≈ 1.1283791670955126
42+
@test ϕ[2](2) ≈ 0.
43+
@test ϕ[3](0) ≈ 1.1283791670955126
44+
@test ϕ[3](1) ≈ 1.1283791670955126
45+
@test ϕ[3](2) ≈ 0.
46+
47+
@test ψ1[1](0) ≈ -0.5560622352843183
48+
@test ψ1[1](1) ≈ 0.
49+
@test ψ1[1](2) ≈ 0.
50+
@test ψ1[2](0) ≈ 0.932609257876051
51+
@test ψ1[2](1) ≈ 0.
52+
@test ψ1[2](2) ≈ 0.
53+
@test ψ1[3](0) ≈ 1.0941547380212637
54+
@test ψ1[3](1) ≈ 0.
55+
@test ψ1[3](2) ≈ 0.
56+
57+
@test ψ2[1](0) ≈ -0.
58+
@test ψ2[1](1) ≈ 0.5560622352843181
59+
@test ψ2[1](2) ≈ 0.
60+
@test ψ2[2](0) ≈ 0.
61+
@test ψ2[2](1) ≈ 0.9326092578760665
62+
@test ψ2[2](2) ≈ 0.
63+
@test ψ2[3](0) ≈ 0.
64+
@test ψ2[3](1) ≈ -1.0941547380212384
65+
@test ψ2[3](2) ≈ 0.
66+
end
3467
end

0 commit comments

Comments
 (0)