Improve and fix spline (#610)

* vectorize Cox - de Boor recursion

Co-authored-by: Filippo Olivo <folivo@filippoolivo.com>
Co-authored-by: ajacoby9 <a99jacoby@gmail.com>

* fix logic and extend tests

* add b-spline surface

---------

Co-authored-by: Filippo Olivo <folivo@filippoolivo.com>
Co-authored-by: GiovanniCanali <giovanni.canali98@yahoo.it>
This commit is contained in:
ajacoby9
2025-10-15 06:31:43 -04:00
committed by GitHub
parent 9d465fc667
commit 0844282335
7 changed files with 839 additions and 177 deletions

View File

@@ -1,81 +1,175 @@
import torch
import pytest
from scipy.interpolate import BSpline
from pina.model import Spline
from pina import LabelTensor
data = torch.rand((20, 3))
input_vars = 3
output_vars = 4
valid_args = [
{
"knots": torch.tensor([0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0]),
"control_points": torch.tensor([0.0, 0.0, 1.0, 0.0, 0.0]),
"order": 3,
},
{
"knots": torch.tensor(
[-2.0, -2.0, -2.0, -2.0, -1.0, 0.0, 1.0, 2.0, 2.0, 2.0, 2.0]
),
"control_points": torch.tensor([0.0, 0.0, 0.0, 6.0, 0.0, 0.0, 0.0]),
"order": 4,
},
# {'control_points': {'n': 5, 'dim': 1}, 'order': 2},
# {'control_points': {'n': 7, 'dim': 1}, 'order': 3}
# Utility quantities for testing
order = torch.randint(1, 8, (1,)).item()
n_ctrl_pts = torch.randint(order, order + 5, (1,)).item()
n_knots = order + n_ctrl_pts
# Input tensor
points = [
LabelTensor(torch.rand(100, 1), ["x"]),
LabelTensor(torch.rand(2, 100, 1), ["x"]),
]
def scipy_check(model, x, y):
from scipy.interpolate._bsplines import BSpline
import numpy as np
# Function to compare with scipy implementation
def check_scipy_spline(model, x, output_):
spline = BSpline(
# Define scipy spline
scipy_spline = BSpline(
t=model.knots.detach().numpy(),
c=model.control_points.detach().numpy(),
k=model.order - 1,
)
y_scipy = spline(x).flatten()
y = y.detach().numpy()
np.testing.assert_allclose(y, y_scipy, atol=1e-5)
# Compare outputs
torch.allclose(
output_,
torch.tensor(scipy_spline(x), dtype=output_.dtype),
atol=1e-5,
rtol=1e-5,
)
# Define all possible combinations of valid arguments for Spline class
valid_args = [
{
"order": order,
"control_points": torch.rand(n_ctrl_pts),
"knots": torch.linspace(0, 1, n_knots),
},
{
"order": order,
"control_points": torch.rand(n_ctrl_pts),
"knots": {"n": n_knots, "min": 0, "max": 1, "mode": "auto"},
},
{
"order": order,
"control_points": torch.rand(n_ctrl_pts),
"knots": {"n": n_knots, "min": 0, "max": 1, "mode": "uniform"},
},
{
"order": order,
"control_points": None,
"knots": torch.linspace(0, 1, n_knots),
},
{
"order": order,
"control_points": None,
"knots": {"n": n_knots, "min": 0, "max": 1, "mode": "auto"},
},
{
"order": order,
"control_points": None,
"knots": {"n": n_knots, "min": 0, "max": 1, "mode": "uniform"},
},
{
"order": order,
"control_points": torch.rand(n_ctrl_pts),
"knots": None,
},
]
@pytest.mark.parametrize("args", valid_args)
def test_constructor(args):
Spline(**args)
# Should fail if order is not a positive integer
with pytest.raises(AssertionError):
Spline(
order=-1, control_points=args["control_points"], knots=args["knots"]
)
def test_constructor_wrong():
# Should fail if control_points is not None or a torch.Tensor
with pytest.raises(ValueError):
Spline()
Spline(
order=args["order"], control_points=[1, 2, 3], knots=args["knots"]
)
# Should fail if knots is not None, a torch.Tensor, or a dict
with pytest.raises(ValueError):
Spline(
order=args["order"], control_points=args["control_points"], knots=5
)
# Should fail if both knots and control_points are None
with pytest.raises(ValueError):
Spline(order=args["order"], control_points=None, knots=None)
# Should fail if knots is not one-dimensional
with pytest.raises(ValueError):
Spline(
order=args["order"],
control_points=args["control_points"],
knots=torch.rand(n_knots, 4),
)
# Should fail if control_points is not one-dimensional
with pytest.raises(ValueError):
Spline(
order=args["order"],
control_points=torch.rand(n_ctrl_pts, 4),
knots=args["knots"],
)
# Should fail if the number of knots != order + number of control points
# If control points are None, they are initialized to fulfill this condition
if args["control_points"] is not None:
with pytest.raises(ValueError):
Spline(
order=args["order"],
control_points=args["control_points"],
knots=torch.linspace(0, 1, n_knots + 1),
)
# Should fail if the knot dict is missing required keys
with pytest.raises(ValueError):
Spline(
order=args["order"],
control_points=args["control_points"],
knots={"n": n_knots, "min": 0, "max": 1},
)
# Should fail if the knot dict has invalid 'mode' key
with pytest.raises(ValueError):
Spline(
order=args["order"],
control_points=args["control_points"],
knots={"n": n_knots, "min": 0, "max": 1, "mode": "invalid"},
)
@pytest.mark.parametrize("args", valid_args)
def test_forward(args):
min_x = args["knots"][0]
max_x = args["knots"][-1]
xi = torch.linspace(min_x, max_x, 1000)
@pytest.mark.parametrize("pts", points)
def test_forward(args, pts):
# Define the model
model = Spline(**args)
yi = model(xi).squeeze()
scipy_check(model, xi, yi)
return
# Evaluate the model
output_ = model(pts)
assert output_.shape == pts.shape
# Compare with scipy implementation only for interpolant knots (mode: auto)
if isinstance(args["knots"], dict) and args["knots"]["mode"] == "auto":
check_scipy_spline(model, pts, output_)
@pytest.mark.parametrize("args", valid_args)
def test_backward(args):
min_x = args["knots"][0]
max_x = args["knots"][-1]
xi = torch.linspace(min_x, max_x, 100)
model = Spline(**args)
yi = model(xi)
fake_loss = torch.sum(yi)
assert model.control_points.grad is None
fake_loss.backward()
assert model.control_points.grad is not None
@pytest.mark.parametrize("pts", points)
def test_backward(args, pts):
# dim_in, dim_out = 3, 2
# fnn = FeedForward(dim_in, dim_out)
# data.requires_grad = True
# output_ = fnn(data)
# l=torch.mean(output_)
# l.backward()
# assert data._grad.shape == torch.Size([20,3])
# Define the model
model = Spline(**args)
# Evaluate the model
output_ = model(pts)
loss = torch.mean(output_)
loss.backward()
assert model.control_points.grad.shape == model.control_points.shape

View File

@@ -0,0 +1,180 @@
import torch
import random
import pytest
from pina.model import SplineSurface
from pina import LabelTensor
# Utility quantities for testing
orders = [random.randint(1, 8) for _ in range(2)]
n_ctrl_pts = random.randint(max(orders), max(orders) + 5)
n_knots = [orders[i] + n_ctrl_pts for i in range(2)]
# Input tensor
points = [
LabelTensor(torch.rand(100, 2), ["x", "y"]),
LabelTensor(torch.rand(2, 100, 2), ["x", "y"]),
]
@pytest.mark.parametrize(
"knots_u",
[
torch.rand(n_knots[0]),
{"n": n_knots[0], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[0], "min": 0, "max": 1, "mode": "uniform"},
None,
],
)
@pytest.mark.parametrize(
"knots_v",
[
torch.rand(n_knots[1]),
{"n": n_knots[1], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[1], "min": 0, "max": 1, "mode": "uniform"},
None,
],
)
@pytest.mark.parametrize(
"control_points", [torch.rand(n_ctrl_pts, n_ctrl_pts), None]
)
def test_constructor(knots_u, knots_v, control_points):
# Skip if knots_u, knots_v, and control_points are all None
if (knots_u is None or knots_v is None) and control_points is None:
return
SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=knots_v,
control_points=control_points,
)
# Should fail if orders is not list of two elements
with pytest.raises(ValueError):
SplineSurface(
orders=[orders[0]],
knots_u=knots_u,
knots_v=knots_v,
control_points=control_points,
)
# Should fail if both knots and control_points are None
with pytest.raises(ValueError):
SplineSurface(
orders=orders,
knots_u=None,
knots_v=None,
control_points=None,
)
# Should fail if control_points is not a torch.Tensor when provided
with pytest.raises(ValueError):
SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=knots_v,
control_points=[[0.0] * n_ctrl_pts] * n_ctrl_pts,
)
# Should fail if control_points is not of the correct shape when provided
# It assumes that at least one among knots_u and knots_v is not None
if knots_u is not None or knots_v is not None:
with pytest.raises(ValueError):
SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=knots_v,
control_points=torch.rand(n_ctrl_pts + 1, n_ctrl_pts + 1),
)
# Should fail if there are not enough knots_u to define the control points
with pytest.raises(ValueError):
SplineSurface(
orders=orders,
knots_u=torch.linspace(0, 1, orders[0]),
knots_v=knots_v,
control_points=None,
)
# Should fail if there are not enough knots_v to define the control points
with pytest.raises(ValueError):
SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=torch.linspace(0, 1, orders[1]),
control_points=None,
)
@pytest.mark.parametrize(
"knots_u",
[
torch.rand(n_knots[0]),
{"n": n_knots[0], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[0], "min": 0, "max": 1, "mode": "uniform"},
],
)
@pytest.mark.parametrize(
"knots_v",
[
torch.rand(n_knots[1]),
{"n": n_knots[1], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[1], "min": 0, "max": 1, "mode": "uniform"},
],
)
@pytest.mark.parametrize(
"control_points", [torch.rand(n_ctrl_pts, n_ctrl_pts), None]
)
@pytest.mark.parametrize("pts", points)
def test_forward(knots_u, knots_v, control_points, pts):
# Define the model
model = SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=knots_v,
control_points=control_points,
)
# Evaluate the model
output_ = model(pts)
assert output_.shape == (*pts.shape[:-1], 1)
@pytest.mark.parametrize(
"knots_u",
[
torch.rand(n_knots[0]),
{"n": n_knots[0], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[0], "min": 0, "max": 1, "mode": "uniform"},
],
)
@pytest.mark.parametrize(
"knots_v",
[
torch.rand(n_knots[1]),
{"n": n_knots[1], "min": 0, "max": 1, "mode": "auto"},
{"n": n_knots[1], "min": 0, "max": 1, "mode": "uniform"},
],
)
@pytest.mark.parametrize(
"control_points", [torch.rand(n_ctrl_pts, n_ctrl_pts), None]
)
@pytest.mark.parametrize("pts", points)
def test_backward(knots_u, knots_v, control_points, pts):
# Define the model
model = SplineSurface(
orders=orders,
knots_u=knots_u,
knots_v=knots_v,
control_points=control_points,
)
# Evaluate the model
output_ = model(pts)
loss = torch.mean(output_)
loss.backward()
assert model.control_points.grad.shape == model.control_points.shape