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

Commit

Permalink
complete chebyshev polynomials
Browse files Browse the repository at this point in the history
  • Loading branch information
@yuehhua
yuehhua committed Aug 31, 2022
1 parent 1a71626 commit e5afafb
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 72 deletions.
126 changes: 62 additions & 64 deletions src/Transform/polynomials.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
15 changes: 7 additions & 8 deletions src/Transform/utils.jl
View file
Original file line number Diff line number Diff line change
@@ -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(x<lb, x>ub) * 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
Expand Down
33 changes: 33 additions & 0 deletions test/polynomials.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit e5afafb

Please sign in to comment.