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

Commit 66670c3

Browse files
yuehhuayuehhua
committed
complete legendre polynomials
1 parent c868f71 commit 66670c3

File tree

6 files changed

+71
-49
lines changed

6 files changed

+71
-49
lines changed

‎src/Transform/polynomials.jl

Lines changed: 30 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -22,23 +22,24 @@ function legendre_ϕ_ψ(k)
2222
ϕ_2x_coefs[ki+1, 1:(ki+1)] .= sqrt(2*(2*ki+1)) .* coeffs(l(p2))
2323
end
2424

25-
ψ1_coefs .= ϕ_2x_coefs
25+
ψ1_coefs = zeros(k, k)
2626
ψ2_coefs = zeros(k, k)
2727
for ki in 0:(k-1)
28+
ψ1_coefs[ki+1, :] .= ϕ_2x_coefs[ki+1, :]
2829
for i in 0:(k-1)
2930
a = ϕ_2x_coefs[ki+1, 1:(ki+1)]
3031
b = ϕ_coefs[i+1, 1:(i+1)]
3132
proj_ = proj_factor(a, b)
32-
view(ψ1_coefs, ki+1, :) .-= proj_ .* view(ϕ_coefs, i+1, :)
33-
view(ψ2_coefs, ki+1, :) .-= proj_ .* view(ϕ_coefs, i+1, :)
33+
ψ1_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :)
34+
ψ2_coefs[ki+1, :] .-= proj_ .* view(ϕ_coefs, i+1, :)
3435
end
3536

3637
for j in 0:(k-1)
3738
a = ϕ_2x_coefs[ki+1, 1:(ki+1)]
3839
b = ψ1_coefs[j+1, :]
3940
proj_ = proj_factor(a, b)
40-
view(ψ1_coefs, ki+1, :) .-= proj_ .* view(ψ1_coefs, j+1, :)
41-
view(ψ2_coefs, ki+1, :) .-= proj_ .* view(ψ2_coefs, j+1, :)
41+
ψ1_coefs[ki+1, :] .-= proj_ .* view(ψ1_coefs, j+1, :)
42+
ψ2_coefs[ki+1, :] .-= proj_ .* view(ψ2_coefs, j+1, :)
4243
end
4344

4445
a = ψ1_coefs[ki+1, :]
@@ -129,16 +130,11 @@ end
129130
# end
130131

131132
function legendre_filter(k)
132-
# x = Symbol('x')
133-
# H0 = np.zeros((k,k))
134-
# H1 = np.zeros((k,k))
135-
# G0 = np.zeros((k,k))
136-
# G1 = np.zeros((k,k))
137-
# PHI0 = np.zeros((k,k))
138-
# PHI1 = np.zeros((k,k))
139-
# phi, psi1, psi2 = get_phi_psi(k, base)
140-
141-
# ----------------------------------------------------------
133+
H0 = zeros(k, k)legendre
134+
H1 = zeros(k, k)
135+
G0 = zeros(k, k)
136+
G1 = zeros(k, k)
137+
ϕ, ψ1, ψ2 = legendre_ϕ_ψ(k)
142138

143139
# roots = Poly(legendre(k, 2*x-1)).all_roots()
144140
# x_m = np.array([rt.evalf(20) for rt in roots]).astype(np.float64)
@@ -150,29 +146,23 @@ function legendre_filter(k)
150146
# G0[ki, kpi] = 1/np.sqrt(2) * (wm * psi(psi1, psi2, ki, x_m/2) * phi[kpi](x_m)).sum()
151147
# H1[ki, kpi] = 1/np.sqrt(2) * (wm * phi[ki]((x_m+1)/2) * phi[kpi](x_m)).sum()
152148
# G1[ki, kpi] = 1/np.sqrt(2) * (wm * psi(psi1, psi2, ki, (x_m+1)/2) * phi[kpi](x_m)).sum()
153-
154-
# PHI0 = np.eye(k)
155-
# PHI1 = np.eye(k)
156-
157-
# ----------------------------------------------------------
158149

159-
# H0[np.abs(H0)<1e-8] = 0
160-
# H1[np.abs(H1)<1e-8] = 0
161-
# G0[np.abs(G0)<1e-8] = 0
162-
# G1[np.abs(G1)<1e-8] = 0
150+
zero_out!(H0)
151+
zero_out!(H1)
152+
zero_out!(G0)
153+
zero_out!(G1)
163154

164-
# return H0, H1, G0, G1, PHI0, PHI1
155+
return H0, H1, G0, G1, I(k), I(k)
165156
end
166157

167158
function chebyshev_filter(k)
168-
# x = Symbol('x')
169-
# H0 = np.zeros((k,k))
170-
# H1 = np.zeros((k,k))
171-
# G0 = np.zeros((k,k))
172-
# G1 = np.zeros((k,k))
173-
# PHI0 = np.zeros((k,k))
174-
# PHI1 = np.zeros((k,k))
175-
# phi, psi1, psi2 = get_phi_psi(k, base)
159+
H0 = zeros(k, k)
160+
H1 = zeros(k, k)
161+
G0 = zeros(k, k)
162+
G1 = zeros(k, k)
163+
Φ0 = zeros(k, k)
164+
Φ1 = zeros(k, k)
165+
ϕ, ψ1, ψ2 = chebyshev_ϕ_ψ(k)
176166

177167
# ----------------------------------------------------------
178168

@@ -193,16 +183,13 @@ function chebyshev_filter(k)
193183

194184
# PHI0[ki, kpi] = (wm * phi[ki](2*x_m) * phi[kpi](2*x_m)).sum() * 2
195185
# PHI1[ki, kpi] = (wm * phi[ki](2*x_m-1) * phi[kpi](2*x_m-1)).sum() * 2
196-
197-
# PHI0[np.abs(PHI0)<1e-8] = 0
198-
# PHI1[np.abs(PHI1)<1e-8] = 0
199-
200-
# ----------------------------------------------------------
201186

202-
# H0[np.abs(H0)<1e-8] = 0
203-
# H1[np.abs(H1)<1e-8] = 0
204-
# G0[np.abs(G0)<1e-8] = 0
205-
# G1[np.abs(G1)<1e-8] = 0
187+
zero_out!(H0)
188+
zero_out!(H1)
189+
zero_out!(G0)
190+
zero_out!(G1)
191+
zero_out!(Φ0)
192+
zero_out!(Φ1)
206193

207-
# return H0, H1, G0, G1, PHI0, PHI1
194+
return H0, H1, G0, G1, Φ0, Φ1
208195
end

‎src/Transform/utils.jl

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,6 @@ end
3131

3232
function proj_factor(a, b; complement::Bool=false)
3333
prod_ = convolve(a, b)
34-
zero_out!(prod_)
3534
r = collect(1:length(prod_))
3635
s = complement ? (1 .- 0.5 .^ r) : (0.5 .^ r)
3736
proj_ = sum(prod_ ./ r .* s)

‎src/Transform/wavelet_transform.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -65,11 +65,11 @@ struct MWT_CZ1d{T,S,R,Q,P}
6565
end
6666

6767
function MWT_CZ1d(k::Int=3, α::Int=5, L::Int=0, c::Int=1; base::Symbol=:legendre, init=Flux.glorot_uniform)
68-
H0, H1, G0, G1, ϕ0, ϕ1 = get_filter(base, k)
69-
H0r = zero_out!(H0 * ϕ0)
70-
G0r = zero_out!(G0 * ϕ0)
71-
H1r = zero_out!(H1 * ϕ1)
72-
G1r = zero_out!(G1 * ϕ1)
68+
H0, H1, G0, G1, Φ0, Φ1 = get_filter(base, k)
69+
H0r = zero_out!(H0 * Φ0)
70+
G0r = zero_out!(G0 * Φ0)
71+
H1r = zero_out!(H1 * Φ1)
72+
G1r = zero_out!(G1 * Φ1)
7373

7474
dim = c*k
7575
A = SpectralConv(dim=>dim, (α,); init=init)

‎test/Transform/Transform.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
@testset "Transform" begin
2+
include("polynomials.jl")
23
include("fourier_transform.jl")
34
include("chebyshev_transform.jl")
45
include("wavelet_transform.jl")

‎test/polynomials.jl

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
@testset "polynomials" begin
2+
@testset "legendre_ϕ_ψ" begin
3+
ϕ, ψ1, ψ2 = NeuralOperators.legendre_ϕ_ψ(10)
4+
5+
@test all(coeffs(ϕ[1]) .≈ [1.])
6+
@test all(coeffs(ϕ[2]) .≈ [-1.7320508075688772, 3.4641016151377544])
7+
@test all(coeffs(ϕ[3]) .≈ [2.23606797749979, -13.416407864998739, 13.416407864998739])
8+
@test all(coeffs(ϕ[4]) .≈ [-2.6457513110645907, 31.74901573277509, -79.37253933193772, 52.91502622129181])
9+
@test all(coeffs(ϕ[5]) .≈ [3.0, -60.0, 270.0, -420.0, 210.0])
10+
@test all(coeffs(ϕ[6]) .≈ [-3.3166247903554, 99.498743710662, -696.491205974634, 1857.309882599024,
11+
-2089.4736179239017, 835.7894471695607])
12+
@test all(coeffs(ϕ[7]) .≈ [3.605551275463989, -151.43315356948753, 1514.3315356948754, -6057.326142779501,
13+
11357.486517711566, -9994.588135586178, 3331.529378528726])
14+
@test all(coeffs(ϕ[8]) .≈ [-3.872983346207417, 216.88706738761536, -2927.9754097328073, 16266.530054071152,
15+
-44732.957648695665, 64415.45901412176, -46522.27595464349, 13292.078844183856])
16+
@test all(coeffs(ϕ[9]) .≈ [4.123105625617661, -296.86360504447157, 5195.113088278253, -38097.49598070719,
17+
142865.60992765194, -297160.46864951606, 346687.21342443535, -212257.47760679715,
18+
53064.36940169929])
19+
@test all(coeffs(ϕ[10]) .≈ [-4.358898943540674, 392.30090491866065, -8630.619908210534, 80552.45247663166,
20+
-392693.20582357934, 1099540.9763060221, -1832568.2938433702, 1795168.9409077913,
21+
-953683.4998572641, 211929.66663494756])
22+
23+
ϕ, ψ1, ψ2 = NeuralOperators.legendre_ϕ_ψ(3)
24+
@test coeffs(ϕ[1]) ≈ [1.]
25+
@test coeffs(ϕ[2]) ≈ [-1.7320508075688772, 3.4641016151377544]
26+
@test coeffs(ϕ[3]) ≈ [2.23606797749979, -13.416407864998739, 13.416407864998739]
27+
@test coeffs(ψ1[1]) ≈ [-1.0000000000000122, 6.000000000000073]
28+
@test coeffs(ψ1[2]) ≈ [1.7320508075691732, -24.248711305967735, 51.96152422707261]
29+
@test coeffs(ψ1[3]) ≈ [2.2360679774995615, -26.832815729994504, 53.665631459989214]
30+
@test coeffs(ψ2[1]) ≈ [-5.000000000000066, 6.000000000000073]
31+
@test coeffs(ψ2[2]) ≈ [29.44486372867492, -79.67433714817852, 51.96152422707261]
32+
@test coeffs(ψ2[3]) ≈ [-29.068883707507286, 80.49844719001908, -53.665631460012115]
33+
end
34+
end

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ using CUDA
33
using Flux
44
using GeometricFlux
55
using Graphs
6+
using Polynomials
67
using Zygote
78
using Test
89

0 commit comments

Comments
 (0)