Type safety: curve.py

splipy/curve.py

@@ -1,27 +1,39 @@
 # -*- coding: utf-8 -*-
+from __future__ import annotations
 
 from itertools import chain
 from bisect import bisect_left, bisect_right
+from typing import Optional, Sequence, Union, Tuple
 
 import numpy as np
+from numpy.typing import ArrayLike, NDArray
+from numpy import floating
 import scipy.sparse.linalg as splinalg
 
 from .basis import BSplineBasis
 from .splineobject import SplineObject
 from .utils import ensure_listlike, is_singleton
+from .types import OneOrMore, Float
 
 __all__ = ['Curve']
 
 
 class Curve(SplineObject):
     """Curve()
 
-    Represents a curve: an object with a one-dimensional parameter space."""
+    Represents a curve: an object with a one-dimensional parameter space.
+    """
 
     _intended_pardim = 1
 
-    def __init__(self, basis=None, controlpoints=None, rational=False, **kwargs):
-        """  Construct a curve with the given basis and control points.
+    def __init__(
+        self,
+        basis: Optional[BSplineBasis] = None,
+        controlpoints: Optional[ArrayLike] = None,
+        rational: bool = False,
+        raw: bool = False,
+    ):
+        """Construct a curve with the given basis and control points.
 
         The default is to create a linear one-element mapping from (0,1) to the
         unit interval.
@@ -32,48 +44,16 @@ def __init__(self, basis=None, controlpoints=None, rational=False, **kwargs):
             control points are interpreted as pre-multiplied with the weight,
             which is the last coordinate)
         """
-        super(Curve, self).__init__([basis], controlpoints, rational, **kwargs)
+        super().__init__([basis], controlpoints, rational, raw)
 
-    def evaluate(self, *params):
-        """  Evaluate the object at given parametric values.
-
-        This function returns an *n1* × *n2* × ... × *dim* array, where *ni* is
-        the number of evaluation points in direction *i*, and *dim* is the
-        physical dimension of the object.
-
-        If there is only one evaluation point, a vector of length *dim* is
-        returned instead.
-
-        :param u,v,...: Parametric coordinates in which to evaluate
-        :type u,v,...: float or [float]
-        :return: Geometry coordinates
-        :rtype: numpy.array
-        """
-        squeeze = is_singleton(params[0])
-        params = [ensure_listlike(p) for p in params]
-
-        self._validate_domain(*params)
-
-        # Evaluate the derivatives of the corresponding bases at the corresponding points
-        # and build the result array
-        N = self.bases[0].evaluate(params[0], sparse=True)
-        result = N @ self.controlpoints
-
-        # For rational objects, we divide out the weights, which are stored in the
-        # last coordinate
-        if self.rational:
-            for i in range(self.dimension):
-                result[..., i] /= result[..., -1]
-            result = np.delete(result, self.dimension, -1)
-
-        # Squeeze the singleton dimensions if we only have one point
-        if squeeze:
-            result = result.reshape(self.dimension)
-
-        return result
-
-    def derivative(self, t, d=1, above=True, tensor=True):
-        """  Evaluate the derivative of the curve at the given parametric values.
+    def derivative(
+        self,
+        *params: OneOrMore[Float],
+        d: OneOrMore[int] = 1,
+        above: OneOrMore[bool] = True,
+        tensor: bool = True,
+    ) -> NDArray[floating]:
+        """Evaluate the derivative of the curve at the given parametric values.
 
         This function returns an *n* × *dim* array, where *n* is the number of
         evaluation points, and *dim* is the physical dimension of the curve.
@@ -89,26 +69,29 @@ def derivative(self, t, d=1, above=True, tensor=True):
         :return: Derivative array
         :rtype: numpy.array
         """
-        if not is_singleton(d):
-            d = d[0]
-        if not self.rational or d < 2 or d > 3:
-            return super(Curve, self).derivative(t, d=d, above=above, tensor=tensor)
 
-        t = ensure_listlike(t)
-        result = np.zeros((len(t), self.dimension))
+        t, = params
+        dd = ensure_listlike(d, 1)[0]
+
+        if not self.rational or dd < 2 or dd > 3:
+            return super().derivative(t, d=dd, above=above, tensor=tensor)
+
+        tt = ensure_listlike(t)
+        eval_above, = ensure_listlike(above, dups=1)
+        result = np.zeros((len(tt), self.dimension))
 
-        d2 = np.array(self.bases[0].evaluate(t, 2, above) @ self.controlpoints)
-        d1 = np.array(self.bases[0].evaluate(t, 1, above) @ self.controlpoints)
-        d0 = np.array(self.bases[0].evaluate(t) @ self.controlpoints)
+        d2 = np.array(self.bases[0].evaluate(tt, 2, eval_above) @ self.controlpoints)
+        d1 = np.array(self.bases[0].evaluate(tt, 1, eval_above) @ self.controlpoints)
+        d0 = np.array(self.bases[0].evaluate(tt) @ self.controlpoints)
         W  = d0[:, -1]  # W(t)
         W1 = d1[:, -1]  # W'(t)
         W2 = d2[:, -1]  # W''(t)
-        if d == 2:
+        if dd == 2:
             for i in range(self.dimension):
                 result[:, i] = (d2[:, i] * W * W - 2 * W1 *
                                (d1[:, i] * W - d0[:, i] * W1) - d0[:, i] * W2 * W) / W / W / W
-        if d == 3:
-            d3 = np.array(self.bases[0].evaluate(t, 3, above) @ self.controlpoints)
+        if dd == 3:
+            d3 = np.array(self.bases[0].evaluate(tt, 3, eval_above) @ self.controlpoints)
             W3 = d3[:,-1]    # W'''(t)
             W6 = W*W*W*W*W*W # W^6
             for i in range(self.dimension):
@@ -124,8 +107,8 @@ def derivative(self, t, d=1, above=True, tensor=True):
 
         return result
 
-    def binormal(self, t, above=True):
-        """  Evaluate the normalized binormal of the curve at the given parametric value(s).
+    def binormal(self, t: OneOrMore[Float], above: bool = True) -> NDArray[floating]:
+        """Evaluate the normalized binormal of the curve at the given parametric value(s).
 
         This function returns an *n* × 3 array, where *n* is the number of
         evaluation points.
@@ -175,8 +158,8 @@ def binormal(self, t, above=True):
 
         return result / magnitude
 
-    def normal(self, t, above=True):
-        """  Evaluate the normal of the curve at the given parametric value(s).
+    def normal(self, t: OneOrMore[Float], above: bool = True) -> NDArray[floating]:
+        """Evaluate the normal of the curve at the given parametric value(s).
 
         This function returns an *n* × 3 array, where *n* is the number of
         evaluation points.
@@ -200,8 +183,8 @@ def normal(self, t, above=True):
 
         return np.cross(B,T)
 
-    def curvature(self, t, above=True):
-        """  Evaluate the curvaure at specified point(s). The curvature is defined as
+    def curvature(self, t: OneOrMore[Float], above: bool = True) -> NDArray[floating]:
+        """Evaluate the curvaure at specified point(s). The curvature is defined as
 
         .. math:: \\frac{|\\boldsymbol{v}\\times \\boldsymbol{a}|}{|\\boldsymbol{v}|^3}
 
@@ -217,23 +200,22 @@ def curvature(self, t, above=True):
         w = np.cross(v,a)
 
         if len(v.shape) == 1: # single evaluation point
-            magnitude = np.linalg.norm(w)
-            speed     = np.linalg.norm(v)
-        else:                 # multiple evaluation points
-            if self.dimension == 2:
-                # for 2D-cases np.cross() outputs scalars
-                # (the z-component of the cross product)
-                magnitude = np.abs(w)
-            else:
-                # for 3D, it is vectors
-                magnitude = np.linalg.norm( w,  axis= -1)
+            return np.linalg.norm(w) / np.power(np.linalg.norm(v), 2)
+
+        if self.dimension == 2:
+            # for 2D-cases np.cross() outputs scalars
+            # (the z-component of the cross product)
+            magnitude = np.abs(w)
+        else:
+            # for 3D, it is vectors
+            magnitude = np.linalg.norm(w, axis= -1)
 
-            speed = np.linalg.norm( v, axis=-1)
+        speed = np.linalg.norm(v, axis=-1)
 
         return magnitude / np.power(speed,3)
 
-    def torsion(self, t, above=True):
-        """  Evaluate the torsion for a 3D curve at specified point(s). The torsion is defined as
+    def torsion(self, t: OneOrMore[Float], above: bool = True) -> NDArray[floating]:
+        """Evaluate the torsion for a 3D curve at specified point(s). The torsion is defined as
 
         .. math:: \\frac{(\\boldsymbol{v}\\times \\boldsymbol{a})\\cdot (d\\boldsymbol{a}/dt)}{|\\boldsymbol{v}\\times \\boldsymbol{a}|^2}
 
@@ -245,12 +227,10 @@ def torsion(self, t, above=True):
         """
         if self.dimension == 2:
             # no torsion for 2D curves
-            t = ensure_listlike(t)
-            return np.zeros(len(t))
-        elif self.dimension == 3:
+            tt = ensure_listlike(t)
+            return np.zeros(len(tt))
+        elif self.dimension != 3:
             # only allow 3D curves
-            pass
-        else:
             raise ValueError('dimension must be 2 or 3')
 
         # compute derivative
@@ -268,34 +248,8 @@ def torsion(self, t, above=True):
 
         return nominator / np.power(magnitude, 2)
 
-    def raise_order(self, amount, direction=None):
-        """  Raise the polynomial order of the curve.
-
-        :param int amount: Number of times to raise the order
-        :return: self
-        """
-        if amount < 0:
-            raise ValueError('Raise order requires a non-negative parameter')
-        elif amount == 0:
-            return
-
-        # create the new basis
-        newBasis = self.bases[0].raise_order(amount)
-
-        # set up an interpolation problem. This is in projective space, so no problems for rational cases
-        interpolation_pts_t = newBasis.greville()  # parametric interpolation points (t)
-        N_old = self.bases[0].evaluate(interpolation_pts_t)
-        N_new = newBasis.evaluate(interpolation_pts_t, sparse=True)
-        interpolation_pts_x = N_old @ self.controlpoints  # projective interpolation points (x,y,z,w)
-
-        # solve the interpolation problem
-        self.controlpoints = np.array(splinalg.spsolve(N_new, interpolation_pts_x))
-        self.bases = [newBasis]
-
-        return self
-
-    def append(self, curve):
-        """  Extend the curve by merging another curve to the end of it.
+    def append(self, curve: Curve) -> Curve:
+        """Extend the curve by merging another curve to the end of it.
 
         The curves are glued together in a C0 fashion with enough repeated
         knots. The function assumes that the end of this curve perfectly
@@ -349,21 +303,20 @@ def append(self, curve):
 
         return self
 
-    def continuity(self, knot):
-        """  Get the parametric continuity of the curve at a given point. Will
+    def continuity(self, knot: float) -> Union[int, float]:
+        """Get the parametric continuity of the curve at a given point. Will
         return p-1-m, where m is the knot multiplicity and inf between knots"""
         return self.bases[0].continuity(knot)
 
-    def get_kinks(self):
-        """  Get the parametric coordinates at all points which have C0-
+    def get_kinks(self) -> NDArray[floating]:
+        """Get the parametric coordinates at all points which have C0-
         continuity"""
-        return [k for k in self.knots(0) if self.continuity(k)<1]
+        return np.array([k for k in self.knots(0) if self.continuity(k)<1])
 
-    def length(self, t0=None, t1=None):
+    def length(self, t0: Optional[float] = None, t1: Optional[float] = None) -> float:
         """ Computes the euclidian length of the curve in geometric space
 
         .. math:: \\int_{t_0}^{t_1}\\sqrt{x(t)^2 + y(t)^2 + z(t)^2} dt
-
         """
         (x,w) = np.polynomial.legendre.leggauss(self.order(0)+1)
         knots = self.knots(0)
@@ -385,8 +338,8 @@ def length(self, t0=None, t1=None):
         detJ = np.sqrt(np.sum(dx**2, axis=1))
         return np.dot(detJ, w)
 
-    def rebuild(self, p, n):
-        """  Creates an approximation to this curve by resampling it using a
+    def rebuild(self, p: int, n: int) -> Curve:
+        """Creates an approximation to this curve by resampling it using a
         uniform knot vector of order *p* with *n* control points.
 
         :param int p: Polynomial discretization order
@@ -411,8 +364,8 @@ def rebuild(self, p, n):
         # return new resampled curve
         return Curve(basis, controlpoints)
 
-    def error(self, target):
-        """  Computes the L2 (squared and per knot span) and max error between
+    def error(self, target: Curve) -> Tuple[NDArray[floating], float]:
+        """Computes the L2 (squared and per knot span) and max error between
         this curve and a target curve
 
         .. math:: ||\\boldsymbol{x_h}(t)-\\boldsymbol{x}(t)||_{L^2(t_1,t_2)}^2 = \\int_{t_1}^{t_2}
@@ -451,10 +404,10 @@ def arclength_circle(t):
             error = np.sum(error**2, axis=1) # |x-xh|^2
             err2.append(np.dot(error, wg))   # integrate over domain
             err_inf = max(np.max(np.sqrt(error)), err_inf)
-        return (err2, err_inf)
+        return (np.array(err2), err_inf)
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return str(self.bases[0]) + '\n' + str(self.controlpoints)
 
-    def get_derivative_curve(self):
-        return super(Curve, self).get_derivative_spline(0)
+    def get_derivative_curve(self) -> Curve:
+        return super().get_derivative_spline(0)

test/curve_test.py

@@ -388,17 +388,17 @@ def expect_derivative_3(x):
         # insert a few more knots to spice things up
         crv.insert_knot([.2, .71])
 
-        self.assertAlmostEqual(crv.derivative(0.32, 1)[0], expect_derivative(0.32))
-        self.assertAlmostEqual(crv.derivative(0.32, 1)[1], 0)
-        self.assertAlmostEqual(crv.derivative(0.71, 1)[0], expect_derivative(0.71))
+        self.assertAlmostEqual(crv.derivative(0.32, d=1)[0], expect_derivative(0.32))
+        self.assertAlmostEqual(crv.derivative(0.32, d=1)[1], 0)
+        self.assertAlmostEqual(crv.derivative(0.71, d=1)[0], expect_derivative(0.71))
 
-        self.assertAlmostEqual(crv.derivative(0.22, 2)[0], expect_derivative_2(0.22))
-        self.assertAlmostEqual(crv.derivative(0.22, 2)[1], 0)
-        self.assertAlmostEqual(crv.derivative(0.86, 2)[0], expect_derivative_2(0.86))
+        self.assertAlmostEqual(crv.derivative(0.22, d=2)[0], expect_derivative_2(0.22))
+        self.assertAlmostEqual(crv.derivative(0.22, d=2)[1], 0)
+        self.assertAlmostEqual(crv.derivative(0.86, d=2)[0], expect_derivative_2(0.86))
 
-        self.assertAlmostEqual(crv.derivative(0.22, 3)[0], expect_derivative_3(0.22))
-        self.assertAlmostEqual(crv.derivative(0.22, 3)[1], 0)
-        self.assertAlmostEqual(crv.derivative(0.86, 3)[0], expect_derivative_3(0.86))
+        self.assertAlmostEqual(crv.derivative(0.22, d=3)[0], expect_derivative_3(0.22))
+        self.assertAlmostEqual(crv.derivative(0.22, d=3)[1], 0)
+        self.assertAlmostEqual(crv.derivative(0.86, d=3)[0], expect_derivative_3(0.86))
 
     def test_tangent_and_normal(self):
         crv = cf.circle()