Adding speedup for _curve_helpers.reduce_pseudo_inverse.

src/bezier/_curve.pxd

@@ -13,6 +13,9 @@
 """Cython wrapper for ``curve.f90``."""
 
 
+from libcpp cimport bool as bool_t
+
+
 cdef extern from "bezier/curve.h":
     void evaluate_curve_barycentric(
         int *degree, int *dimension, double *nodes, int *num_vals,
@@ -47,3 +50,6 @@ cdef extern from "bezier/curve.h":
     void get_curvature(
         int *num_nodes, int *dimension, double *nodes, double *tangent_vec,
         double *s, double *curvature)
+    void reduce_pseudo_inverse(
+        int *num_nodes, int *dimension, double *nodes, double *reduced,
+        bool_t *not_implemented)

src/bezier/_curve_helpers.py

@@ -807,7 +807,7 @@ def _locate_point(nodes, degree, point):
     return newton_refine(nodes, degree, point, s_approx)
 
 
-def reduce_pseudo_inverse(nodes, degree):
+def _reduce_pseudo_inverse(nodes, degree):
     """Performs degree-reduction for a B |eacute| zier curve.
 
     Does so by using the pseudo-inverse of the degree elevation
@@ -946,6 +946,7 @@ def full_reduce(nodes):
     get_curvature = _get_curvature
     newton_refine = _newton_refine
     locate_point = _locate_point
+    reduce_pseudo_inverse = _reduce_pseudo_inverse
 else:
     subdivide_nodes = _curve_speedup.subdivide_nodes
     evaluate_multi = _curve_speedup.evaluate_multi
@@ -956,4 +957,5 @@ def full_reduce(nodes):
     get_curvature = _curve_speedup.get_curvature
     newton_refine = _curve_speedup.newton_refine
     locate_point = _curve_speedup.locate_point
+    reduce_pseudo_inverse = _curve_speedup.reduce_pseudo_inverse
 # pylint: enable=invalid-name

src/bezier/_curve_speedup.pyx

@@ -14,6 +14,7 @@
 """Cython "wrapped" interface for `_curve_helpers`."""
 
 
+from libcpp cimport bool as bool_t
 import numpy as np
 from numpy cimport ndarray as ndarray_t
 
@@ -214,3 +215,27 @@ def get_curvature(
     )
 
     return curvature
+
+
+def reduce_pseudo_inverse(double[::1, :] nodes, int unused_degree):
+    cdef int num_nodes, dimension
+    cdef bool_t not_implemented
+    cdef ndarray_t[double, ndim=2, mode='fortran'] reduced
+
+    # NOTE: We don't check that there are ``degree + 1`` rows.
+    num_nodes, dimension = np.shape(nodes)
+
+    reduced = np.empty((num_nodes - 1, dimension), order='F')
+
+    bezier._curve.reduce_pseudo_inverse(
+        &num_nodes,
+        &dimension,
+        &nodes[0, 0],
+        &reduced[0, 0],
+        &not_implemented,
+    )
+
+    if not_implemented:
+        raise NotImplementedError(num_nodes - 1)
+
+    return reduced

src/bezier/curve.f90

@@ -21,7 +21,7 @@ module curve
        evaluate_curve_barycentric, evaluate_multi, specialize_curve_generic, &
        specialize_curve_quadratic, specialize_curve, evaluate_hodograph, &
        subdivide_nodes_generic, subdivide_nodes, newton_refine, locate_point, &
-       elevate_nodes, get_curvature
+       elevate_nodes, get_curvature, reduce_pseudo_inverse
 
   ! For ``locate_point``.
   type :: LocateCandidate
@@ -477,4 +477,48 @@ subroutine get_curvature( &
 
   end subroutine get_curvature
 
+  subroutine reduce_pseudo_inverse( &
+       num_nodes, dimension_, nodes, reduced, not_implemented) &
+       bind(c, name='reduce_pseudo_inverse')
+
+    integer(c_int), intent(in) :: num_nodes, dimension_
+    real(c_double), intent(in) :: nodes(num_nodes, dimension_)
+    real(c_double), intent(out) :: reduced(num_nodes - 1, dimension_)
+    logical(c_bool), intent(out) :: not_implemented
+
+    not_implemented = .FALSE.
+    if (num_nodes == 2) then
+       reduced(1, :) = 0.5_dp * (nodes(1, :) + nodes(2, :))
+    else if (num_nodes == 3) then
+       reduced(1, :) = (5 * nodes(1, :) + 2 * nodes(2, :) - nodes(3, :)) / 6
+       reduced(2, :) = (-nodes(1, :) + 2 * nodes(2, :) + 5 * nodes(3, :)) / 6
+    else if (num_nodes == 4) then
+       reduced(1, :) = ( &
+            19 * nodes(1, :) + 3 * nodes(2, :) - &
+            3 * nodes(3, :) + nodes(4, :)) / 20
+       reduced(2, :) = 0.25_dp * ( &
+            -nodes(1, :) + 3 * nodes(2, :) + &
+            3 * nodes(3, :) - nodes(4, :))
+       reduced(3, :) = ( &
+            nodes(1, :) - 3 * nodes(2, :) + &
+            3 * nodes(3, :) + 19 * nodes(4, :)) / 20
+    else if (num_nodes == 5) then
+       reduced(1, :) = ( &
+            69 * nodes(1, :) + 4 * nodes(2, :) - 6 * nodes(3, :) + &
+            4 * nodes(4, :) - nodes(5, :)) / 70
+       reduced(2, :) = ( &
+            -53 * nodes(1, :) + 212 * nodes(2, :) + 102 * nodes(3, :) - &
+            68 * nodes(4, :) + 17 * nodes(5, :)) / 210
+       reduced(3, :) = ( &
+            17 * nodes(1, :) - 68 * nodes(2, :) + 102 * nodes(3, :) + &
+            212 * nodes(4, :) - 53 * nodes(5, :)) / 210
+       reduced(4, :) = ( &
+            -nodes(1, :) + 4 * nodes(2, :) - 6 * nodes(3, :) + &
+            4 * nodes(4, :) + 69 * nodes(5, :)) / 70
+    else
+       not_implemented = .TRUE.
+    end if
+
+  end subroutine reduce_pseudo_inverse
+
 end module curve

src/bezier/include/bezier/curve.h

@@ -15,6 +15,24 @@
 
 #if defined (__cplusplus)
 extern "C" {
+#elif !defined (_MSC_VER)
+#include <stdbool.h>
+#endif
+
+#if !defined (__cplusplus) && defined (_MSC_VER)
+#  if !defined (FALSE)
+#    define FALSE (0)
+#  endif
+#  if !defined (TRUE)
+#    define TRUE (!FALSE)
+#  endif
+#  if _MSC_VER >= 1800
+#    include <stdbool.h>
+#  else
+#    define bool int
+#    define true TRUE
+#    define false FALSE
+#  endif
 #endif
 
 void evaluate_curve_barycentric(
@@ -49,6 +67,9 @@ void elevate_nodes(
 void get_curvature(
     int *num_nodes, int *dimension, double *nodes, double *tangent_vec,
     double *s, double *curvature);
+void reduce_pseudo_inverse(
+    int *num_nodes, int *dimension, double *nodes, double *reduced,
+    bool *not_implemented);
 
 #if defined (__cplusplus)
 }

tests/test__curve_helpers.py

@@ -757,15 +757,15 @@ def _call_function_under_test(nodes, degree, point):
         return _curve_speedup.locate_point(nodes, degree, point)
 
 
-class Test_reduce_pseudo_inverse(utils.NumPyTestCase):
+class Test__reduce_pseudo_inverse(utils.NumPyTestCase):
 
     EPS = 0.5**52
 
     @staticmethod
     def _call_function_under_test(nodes, degree):
         from bezier import _curve_helpers
 
-        return _curve_helpers.reduce_pseudo_inverse(nodes, degree)
+        return _curve_helpers._reduce_pseudo_inverse(nodes, degree)
 
     def test_to_constant(self):
         nodes = np.asfortranarray([
@@ -876,6 +876,16 @@ def test_unsupported_degree(self):
             self._call_function_under_test(nodes, degree)
 
 
+@unittest.skipIf(utils.WITHOUT_SPEEDUPS, 'No speedups available')
+class Test_speedup__reduce_pseudo_inverse(Test__reduce_pseudo_inverse):
+
+    @staticmethod
+    def _call_function_under_test(nodes, degree):
+        from bezier import _curve_speedup
+
+        return _curve_speedup.reduce_pseudo_inverse(nodes, degree)
+
+
 class Test_projection_error(unittest.TestCase):
 
     @staticmethod