Adds integral functions and fixes type hinting for np.array

pybezier/bezier_curve.py

@@ -5,7 +5,7 @@
 
 class BezierCurve(object):
 
-    def __init__(self, points : np.array, initial_time : float = 0, final_time : float = 1):
+    def __init__(self, points : np.ndarray, initial_time : float = 0, final_time : float = 1):
         if initial_time >= final_time:
             raise ValueError("Initial time must be smaller than final time.")
         self.points = points
@@ -15,10 +15,10 @@ def __init__(self, points : np.array, initial_time : float = 0, final_time : flo
         self.final_time = final_time
         self.duration = final_time - initial_time
 
-    def initial_point(self) -> np.array:
+    def initial_point(self) -> np.ndarray:
         return self.points[0]
 
-    def final_point(self) -> np.array:
+    def final_point(self) -> np.ndarray:
         return self.points[-1]
 
     def _berstein(self, time : float | List[float], n : int) -> float:
@@ -96,8 +96,16 @@ def  __neg__(self) -> Self:
     def derivative(self) -> Self:
         points = (self.points[1:] - self.points[:-1]) * (self.degree / self.duration)
         return BezierCurve(points, self.initial_time, self.final_time)
+
+    def integral(self, initial_condition : np.ndarray | None = None) -> Self:
+        points = self.points * self.duration / self.degree
+        points = np.vstack([np.zeros(self.dimension), points])
+        points = np.cumsum(points, axis=0)
+        if initial_condition is not None:
+            points += initial_condition
+        return BezierCurve(points, self.initial_time, self.final_time)
 
-    def split(self, time : float) -> Tuple[Self, Self]:
+    def domain_split(self, time : float) -> Tuple[Self, Self]:
         if time < self.initial_time:
             raise ValueError("Split time must be larger than initial time.")
         if time > self.final_time:
@@ -132,7 +140,7 @@ def l2_squared(self) -> float:
                 a += b * self.points[i].dot(self.points[j])
         return self.duration * a / (2 * self.degree + 1)
 
-    def integral_of_convex(self, f : Callable) -> float:
+    def integral_of_convex_function(self, f : Callable) -> float:
         c = self.duration / (self.degree + 1)
         return c * sum(f(point) for point in self.points)
 

pybezier/composite_bezier_curve.py

@@ -1,5 +1,5 @@
 import numpy as np
-from typing import List, Self
+from typing import List, Callable, Self
 from collections.abc import Iterable
 from numbers import Number
 from pybezier.bezier_curve import BezierCurve
@@ -27,14 +27,14 @@ def curve_segment(self, time : float) -> int:
             segment += 1
         return segment
 
-    def __call__(self, time : float) -> np.array:
+    def __call__(self, time : float) -> np.ndarray:
         segment = self.curve_segment(time)
         return self[segment](time)
 
-    def initial_point(self) -> np.array:
+    def initial_point(self) -> np.ndarray:
         return self[0].initial_point()
 
-    def final_point(self) -> np.array:
+    def final_point(self) -> np.ndarray:
         return self[-1].final_point()
 
     def __iter__(self) -> Iterable[BezierCurve]:
@@ -81,13 +81,13 @@ def __rsub__(self, composite_curve : Self | Number) -> Self:
     def  __neg__(self) -> Self:
         return 0 - self
 
-    def knot_points(self) -> np.array:
+    def knot_points(self) -> np.ndarray:
         # assumes that the curve is continuous
         knots = [curve.initial_point() for curve in self]
         knots.append(self.final_point())
         return np.array(knots)
 
-    def durations(self) -> np.array:
+    def durations(self) -> np.ndarray:
         return np.array([curve.duration for curve in self])
 
     def concatenate(self, composite_curve : Self) -> Self:
@@ -101,10 +101,20 @@ def concatenate(self, composite_curve : Self) -> Self:
 
     def derivative(self) -> Self:
         return CompositeBezierCurve([curve.derivative() for curve in self])
+
+    def integral(self, initial_condition : np.ndarray | None = None) -> Self:
+        curves = []
+        for curve in self:
+            curves.append(curve.integral(initial_condition))
+            initial_condition = curves[-1].final_point()
+        return CompositeBezierCurve(curves)
 
     def l2_squared(self) -> float:
         return sum(curve.l2_squared() for curve in self)
 
+    def integral_of_convex_function(self, f : Callable) -> float:
+        return sum(curve.integral_of_convex_function(f) for curve in self)
+
     def plot_components(self, n : int = 51, legend : bool = True, **kwargs):
         for i, curve in enumerate(self):
             curve.plot_components(n, legend, **kwargs)

tests/test_bezier_curve.py

@@ -109,15 +109,26 @@ def test_neg(self):
             np.testing.assert_array_almost_equal(neg(time), -self.curve(time))
 
     def test_derivative(self):
-        der = self.curve.derivative()
-        time_step = 1e-9
+        derivative = self.curve.derivative()
+        time_step = 1e-6
         for time in np.linspace(self.initial_time, self.final_time - time_step):
-            value = (self.curve(time + time_step) - self.curve(time)) / time_step
-            np.testing.assert_array_almost_equal(der(time), value)
-
-    def test_split(self):
+            numerical_derivative = (self.curve(time + time_step) - self.curve(time)) / time_step
+            np.testing.assert_array_almost_equal(derivative(time), numerical_derivative)
+
+    def test_integral(self):
+        initial_conditions = [None, np.ones(self.curve.dimension)]
+        for initial_condition in initial_conditions:
+            integral = self.curve.integral(initial_condition)
+            value = 0 if initial_condition is None else initial_condition
+            np.testing.assert_array_almost_equal(integral(self.initial_time), value)
+            time_step = 1e-6
+            for time in np.linspace(self.initial_time, self.final_time - time_step):
+                numerical_integral = integral(time) + time_step * self.curve(time)
+                np.testing.assert_array_almost_equal(integral(time + time_step), numerical_integral)
+
+    def test_domain_split(self):
         split_time = (self.initial_time + self.final_time) / 2
-        curve_1, curve_2 = self.curve.split(split_time)
+        curve_1, curve_2 = self.curve.domain_split(split_time)
         for time in self.time_samples:
             if time < split_time:
                 np.testing.assert_array_almost_equal(self.curve(time), curve_1(time))
@@ -131,14 +142,17 @@ def test_l2_squared(self):
         integral = np.trapezoid(values, times)
         self.assertAlmostEqual(self.curve.l2_squared(), integral)
 
-    def test_integral_of_convex(self):
+    def test_integral_of_convex_function(self):
         f = lambda point: point.dot(point)
-        self.assertTrue(self.curve.integral_of_convex(f) >= self.curve.l2_squared())
+        value = self.curve.l2_squared()
+        upper_bound = self.curve.integral_of_convex_function(f)
+        self.assertTrue(value <= upper_bound)
         # upper bound for curve lenght is equal to distance of control points
+        diff_points = self.points[:-1] - self.points[1:]
+        value = sum(np.linalg.norm(diff_points, axis=1))
         derivative = self.curve.derivative()
-        value = sum(np.linalg.norm(y - x) for x, y in zip(self.points[:-1], self.points[1:]))
-
-        self.assertAlmostEqual(derivative.integral_of_convex(np.linalg.norm), value)
+        upper_bound = derivative.integral_of_convex_function(np.linalg.norm)
+        self.assertAlmostEqual(value, upper_bound)
 
 if __name__ == '__main__':
     unittest.main()

tests/test_composite_bezier_curve.py

@@ -120,11 +120,22 @@ def test_neg(self):
             np.testing.assert_array_almost_equal(neg(time), -self.composite_curve(time))
 
     def test_derivative(self):
-        der = self.composite_curve.derivative()
-        time_step = 1e-9
+        derivative = self.composite_curve.derivative()
+        time_step = 1e-7
         for time in np.linspace(self.initial_time, self.final_time - time_step):
-            value = (self.composite_curve(time + time_step) - self.composite_curve(time)) / time_step
-            np.testing.assert_array_almost_equal(der(time), value)
+            numerical_derivative = (self.composite_curve(time + time_step) - self.composite_curve(time)) / time_step
+            np.testing.assert_array_almost_equal(derivative(time), numerical_derivative)
+
+    def test_integral(self):
+        initial_conditions = [None, np.ones(self.composite_curve.dimension)]
+        for initial_condition in initial_conditions:
+            integral = self.composite_curve.integral(initial_condition)
+            value = 0 if initial_condition is None else initial_condition
+            np.testing.assert_array_almost_equal(integral(self.initial_time), value)
+            time_step = 1e-6
+            for time in np.linspace(self.initial_time, self.final_time - time_step):
+                numerical_integral = integral(time) + time_step * self.composite_curve(time)
+                np.testing.assert_array_almost_equal(integral(time + time_step), numerical_integral)
 
     def test_knot_points(self):
         for i, point in enumerate(self.composite_curve.knot_points()):
@@ -152,5 +163,17 @@ def test_l2_squared(self):
         integral = np.trapezoid(values, times)
         self.assertAlmostEqual(self.composite_curve.l2_squared(), integral, places=4)
 
+    def test_integral_of_convex_function(self):
+        f = lambda point: point.dot(point)
+        value = self.composite_curve.l2_squared()
+        upper_bound = self.composite_curve.integral_of_convex_function(f)
+        self.assertTrue(value <= upper_bound)
+        # upper bound for curve lenght is equal to distance of control points
+        diff_points = np.vstack([curve.points[:-1] - curve.points[1:] for curve in self.composite_curve])
+        value = sum(np.linalg.norm(diff_points, axis=1))
+        derivative = self.composite_curve.derivative()
+        upper_bound = derivative.integral_of_convex_function(np.linalg.norm)
+        self.assertAlmostEqual(value, upper_bound)
+
 if __name__ == '__main__':
     unittest.main()