Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Polynomial 2d basis Implementation and Test Suite #71

Merged
merged 3 commits into from
Mar 11, 2024
Merged

Polynomial 2d basis Implementation and Test Suite #71

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4e467a3
poly basis file 2d
paulvxx Mar 10, 2024
15c7b4e
added test case
paulvxx Mar 10, 2024
a0cb747
converted test case to pytest format
paulvxx Mar 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions src/pymgm_test/utils/polynomialBasis2D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
import numpy as np

# Polynomial basis for 2d array
# Pts is a set of points to obtain a polynomial basis for (numpy array)
# Pts should be an n by 2 array (since a 2d point is a 2-tuple) where n is the number of points
# l is the maximum degree of the polynomial basis
# Returns a tuple (Pbasis, lapP) of the polynomial basis (Pbasis) and Laplacian matrix (lapP) [both are numpy arrays]
def poly_basis(Pts, l):
# invalid degree (must be non-negative)
if l < 0: return (None, None)

terms = (l+1)*(l+2) // 2
Pbasis = np.zeros((len(Pts), terms))
lapP = np.zeros((len(Pts), terms))

# fill the PolyBasis array (P) with the different terms (evaluated at the points)
# e.g. for l = 1 would simply be a linear basis [1, x, y]
# for l = 2 would include terms of up to degree 2 [1, x, y, x^2, x*y, y^2]

for ind in range(len(Pts)):
t = 0
for d in range(l+1):
for i in range(d+1):
Pbasis[ind][t] = Pts[ind][0]**(d-i) * Pts[ind][1]**(i)
if (d-i >= 2):
lapP[ind][t] += (d-i)*(d-i-1) * Pts[ind][0]**(d-i-2) * Pts[ind][1]**(i)
if (i >= 2):
lapP[ind][t] += (i)*(i-1) * Pts[ind][0]**(d-i) * Pts[ind][1]**(i-2)
t += 1

return (Pbasis, lapP)
116 changes: 116 additions & 0 deletions test/test_polyBasis2D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
import pytest
import numpy as np
from src.pymgm_test.utils.polynomialBasis2D import poly_basis


@pytest.mark.parametrize(
# Test cases
"points, degree, expected",
[
#expected = (Pbasis, lapP)
#Pbasis is the expected polynomial basis
#lapP is the expected Laplacian matrix

# Points (0,1), (1,0)
# max degree = -1 (invalid)
# Should return (None, None)
(np.array([[0,1],[1,0]]), -1, (None, None) ),

# Points (0,1), (1,0)
# max degree = 0 (constant)
# Should return ([[1.], [1.]], [[0.], [0.]])
(np.array([[0,1],[1,0]]), 0, ([[1.], [1.]], [[0.], [0.]]) ),

# Points (0,1), (1,0), (1,1)
# max degree = 1 (linear)
# Should return ([[1., 0., 1.], [1., 1., 0.], [1., 1., 1.]], [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]])
(np.array([[0,1],[1,0],[1,1]]), 1,
([[1., 0., 1.],
[1., 1., 0.],
[1., 1., 1.]],

[[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]]) ),
(np.array([[0,1],[1,0],[1,1]]), 2,
([[1., 0., 1., 0., 0., 1.],
[1., 1., 0., 1., 0., 0.],
[1., 1., 1., 1., 1., 1.]],

[[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.]]) ),
(np.array([[2,3],[-2,4],[3,5]]), 2,
([[1., 2., 3., 4., 6., 9.],
[1., -2., 4., 4., -8., 16.],
[1., 3., 5., 9., 15., 25.]],

[[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.]]) ),
(np.array([[4,4],[2,6],[0,-3],[1,7]]), 2,
([[1., 4., 4., 16., 16., 16.],
[1., 2., 6., 4., 12., 36.],
[1., 0., -3., 0., 0., 9.],
[1., 1., 7., 1., 7., 49.]],

[[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.],
[0., 0., 0., 2., 0., 2.]]) ),
(np.array([[2,1],[3,4],[-1,1],[0,4]]), 3,
([[1.0, 2.0, 1.0, 4.0, 2.0, 1.0, 8.0, 4.0, 2.0, 1.0],
[1.0, 3.0, 4.0, 9.0, 12.0, 16.0, 27.0, 36.0, 48.0, 64.0],
[1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0],
[1.0, 0.0, 4.0, 0.0, 0.0, 16.0, 0.0, 0.0, 0.0, 64.0]],

[[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 12.0, 2.0, 4.0, 6.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 18.0, 8.0, 6.0, 24.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -6.0, 2.0, -2.0, 6.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 0.0, 8.0, 0.0, 24.0]]) ),
(np.array([[1,0],[-1,3],[6,-1],[-2,5],[1,1]]), 3,
([[1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
[1.0, -1.0, 3.0, 1.0, -3.0, 9.0, -1.0, 3.0, -9.0, 27.0],
[1.0, 6.0, -1.0, 36.0, -6.0, 1.0, 216.0, -36.0, 6.0, -1.0],
[1.0, -2.0, 5.0, 4.0, -10.0, 25.0, -8.0, 20.0, -50.0, 125.0],
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]],

[[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 6.0, 0.0, 2.0, 0.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -6.0, 6.0, -2.0, 18.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 36.0, -2.0, 12.0, -6.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -12.0, 10.0, -4.0, 30.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 6.0, 2.0, 2.0, 6.0]]) ),
(np.array([[1, 2], [-1, 1]]), 4,
([[1, 1, 2, 1, 2, 4, 1, 2, 4, 8, 1, 2, 4, 8, 16],
[1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1]],

[[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 6.0, 4.0, 2.0, 12.0, 12.0, 12.0, 10.0, 12.0, 48.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -6.0, 2.0, -2.0, 6.0, 12.0, -6.0, 4.0, -6.0, 12.0]]
) ),
(np.array([[6, 7], [-5, -3], [8, -5], [3, 6], [4, 9]]), 4,
([[1.0, 6.0, 7.0, 36.0, 42.0, 49.0, 216.0, 252.0, 294.0, 343.0, 1296.0, 1512.0, 1764.0, 2058.0, 2401.0],
[1.0, -5.0, -3.0, 25.0, 15.0, 9.0, -125.0, -75.0, -45.0, -27.0, 625.0, 375.0, 225.0, 135.0, 81.0],
[1.0, 8.0, -5.0, 64.0, -40.0, 25.0, 512.0, -320.0, 200.0, -125.0, 4096.0, -2560.0, 1600.0, -1000.0, 625.0],
[1.0, 3.0, 6.0, 9.0, 18.0, 36.0, 27.0, 54.0, 108.0, 216.0, 81.0, 162.0, 324.0, 648.0, 1296.0],
[1.0, 4.0, 9.0, 16.0, 36.0, 81.0, 64.0, 144.0, 324.0, 729.0, 256.0, 576.0, 1296.0, 2916.0, 6561.0]],
[[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 36.0, 14.0, 12.0, 42.0, 432.0, 252.0, 170.0, 252.0, 588.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, -30.0, -6.0, -10.0, -18.0, 300.0, 90.0, 68.0, 90.0, 108.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 48.0, -10.0, 16.0, -30.0, 768.0, -240.0, 178.0, -240.0, 300.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 18.0, 12.0, 6.0, 36.0, 108.0, 108.0, 90.0, 108.0, 432.0],
[0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 24.0, 18.0, 8.0, 54.0, 192.0, 216.0, 194.0, 216.0, 972.0]]
) )
]
)


# Test the poly_basis function
def test_poly_basis(points, degree, expected):
num_points = len(points)
basis = poly_basis(points, degree)
if (degree < 0):
assert basis[0] == None and basis[1] == None
return

for i in range(num_points):
assert np.all(basis[0][i] == expected[0][i])
assert np.all(basis[1][i] == expected[1][i])
Loading