Improve performance of multivector arithmetic by avoiding copies

clifford/_multivector.py

@@ -15,10 +15,22 @@ class MultiVector(object):
     Parameters
     -------------
     layout: instance of :class:`clifford.Layout`
-        the layout of the algebra
+        The layout of the algebra
 
     value : sequence of length ``layout.gaDims``
-        the coefficients of the base blades
+        The coefficients of the base blades
+
+    dtype : numpy.dtype
+        The datatype to use for the multivector, if no
+        value was passed.
+
+        .. versionadded:: 1.1.0
+
+    copy : bool
+        If True, the default, create a copy of `value` to store in this
+        multivector.
+
+        .. versionadded:: 1.3.0
 
     Notes
     ------
@@ -33,7 +45,7 @@ class MultiVector(object):
     * ``M[N]`` : blade projection
     """
 
-    def __init__(self, layout, value=None, string=None, *, dtype: np.dtype = np.float64) -> None:
+    def __init__(self, layout, value=None, string=None, *, dtype: np.dtype = np.float64, copy: bool = True) -> None:
         """Constructor."""
 
         self.layout = layout
@@ -45,7 +57,7 @@ def __init__(self, layout, value=None, string=None, *, dtype: np.dtype = np.floa
             else:
                 self.value = layout.parse_multivector(string).value
         else:
-            self.value = np.array(value)
+            self.value = np.array(value, copy=copy)
             if self.value.shape != (self.layout.gaDims,):
                 raise ValueError(
                     "value must be a sequence of length %s" %
@@ -85,7 +97,7 @@ def _newMV(self, newValue=None, *, dtype: np.dtype = None) -> 'MultiVector':
         if newValue is None and dtype is None:
             raise TypeError("Must specify either a type or value")
 
-        return self.__class__(self.layout, newValue, dtype=dtype)
+        return self.__class__(self.layout, newValue, dtype=dtype, copy=False)
 
     # numeric special methods
     # binary
@@ -110,7 +122,7 @@ def vee(self, other) -> 'MultiVector':
         functions instead, as these work in degenerate metrics like PGA too,
         and are equivalent but faster in other metrics.
         """
-        return self.layout.MultiVector(value=self.layout.vee_func(self.value, other.value))
+        return self.layout.MultiVector(value=self.layout.vee_func(self.value, other.value), copy=False)
 
     def __and__(self, other) -> 'MultiVector':
         """ Alias for :meth:`~MultiVector.vee` """
@@ -243,10 +255,10 @@ def __rsub__(self, other) -> 'MultiVector':
         return self._newMV(newValue)
 
     def right_complement(self) -> 'MultiVector':
-        return self.layout.MultiVector(value=self.layout.right_complement_func(self.value))
+        return self.layout.MultiVector(value=self.layout.right_complement_func(self.value), copy=False)
 
     def left_complement(self) -> 'MultiVector':
-        return self.layout.MultiVector(value=self.layout.left_complement_func(self.value))
+        return self.layout.MultiVector(value=self.layout.left_complement_func(self.value), copy=False)
 
     def __truediv__(self, other) -> 'MultiVector':
         """Division, :math:`M N^{-1}`"""
@@ -757,7 +769,7 @@ def dual(self, I=None) -> 'MultiVector':
         I defaults to the pseudoscalar.
         """
         if I is None:
-            return self.layout.MultiVector(value=self.layout.dual_func(self.value))
+            return self.layout.MultiVector(value=self.layout.dual_func(self.value), copy=False)
         else:
             Iinv = I.inv()
 

clifford/tools/g3/__init__.py

@@ -191,7 +191,7 @@ def np_to_euc_mv(np_in):
     output[1] = np_in[0]
     output[2] = np_in[1]
     output[3] = np_in[2]
-    return layout.MultiVector(output)
+    return layout.MultiVector(output, copy=False)
 
 
 def euc_mv_to_np(euc_point):

clifford/tools/g3c/__init__.py

@@ -598,7 +598,7 @@ def get_line_intersection(L3, Ldd):
     Pd = 0.5*(Xdd+Xddd)
     P = -(Pd*ninf*Pd)(1)/(2*(Pd|einf)**2)[0]
     """
-    return layout.MultiVector(val_get_line_intersection(L3.value, Ldd.value))
+    return layout.MultiVector(val_get_line_intersection(L3.value, Ldd.value), copy=False)
 
 
 @numba.njit
@@ -618,15 +618,15 @@ def midpoint_between_lines(L1, L2):
     Gets the point that is maximally close to both lines
     Hadfield and Lasenby AGACSE2018
     """
-    return layout.MultiVector(val_midpoint_between_lines(L1.value, L2.value))
+    return layout.MultiVector(val_midpoint_between_lines(L1.value, L2.value), copy=False)
 
 
 def midpoint_of_line_cluster(line_cluster):
     """
     Gets a center point of a line cluster
     Hadfield and Lasenby AGACSE2018
     """
-    return layout.MultiVector(val_midpoint_of_line_cluster(MVArray(line_cluster).value))
+    return layout.MultiVector(val_midpoint_of_line_cluster(MVArray(line_cluster).value), copy=False)
 
 
 @numba.njit
@@ -809,7 +809,7 @@ def generate_dilation_rotor(scale):
     if abs(scale - 1.0) < 0.00001:
         u = np.zeros(32)
         u[0] = 1.0
-        return layout.MultiVector(u)
+        return layout.MultiVector(u, copy=False)
     gamma = math.log(scale)
     return math.cosh(gamma/2) + math.sinh(gamma/2)*(ninf^no)
 
@@ -818,7 +818,7 @@ def generate_translation_rotor(euc_vector_a):
     """
     Generates a rotor that translates objects along the euclidean vector euc_vector_a
     """
-    return layout.MultiVector(val_generate_translation_rotor(euc_vector_a.value))
+    return layout.MultiVector(val_generate_translation_rotor(euc_vector_a.value), copy=False)
 
 
 @numba.njit
@@ -845,7 +845,7 @@ def meet(A, B):
     The meet algorithm as described in "A Covariant Approach to Geometry"
     I5*((I5*A) ^ (I5*B))
     """
-    return layout.MultiVector(meet_val(A.value, B.value))
+    return layout.MultiVector(meet_val(A.value, B.value), copy=False)
 
 
 @numba.njit
@@ -878,7 +878,7 @@ def intersect_line_and_plane_to_point(line, plane):
     ans = val_intersect_line_and_plane_to_point(line.value, plane.value)
     if ans[0] == -1.:
         return None
-    return layout.MultiVector(ans)
+    return layout.MultiVector(ans, copy=False)
 
 
 @numba.njit
@@ -897,7 +897,7 @@ def normalise_n_minus_1(mv):
     """
     Normalises a conformal point so that it has an inner product of -1 with einf
     """
-    return layout.MultiVector(val_normalise_n_minus_1(mv.value))
+    return layout.MultiVector(val_normalise_n_minus_1(mv.value), copy=False)
 
 
 def quaternion_and_vector_to_rotor(quaternion, vector):
@@ -952,8 +952,8 @@ def point_pair_to_end_points(T):
     Extracts the end points of a point pair bivector
     """
     output = val_point_pair_to_end_points(T.value)
-    A = layout.MultiVector(output[0, :])
-    B = layout.MultiVector(output[1, :])
+    A = layout.MultiVector(output[0, :], copy=False)
+    B = layout.MultiVector(output[1, :], copy=False)
     return A, B
 
 
@@ -1037,13 +1037,13 @@ def positive_root(sigma):
     Square Root of Rotors - Evaluates the positive root
     """
     res_val = positive_root_val(sigma.value)
-    return layout.MultiVector(res_val)
+    return layout.MultiVector(res_val, copy=False)
 
 
 def negative_root(sigma):
     """ Square Root of Rotors - Evaluates the negative root """
     res_val = negative_root_val(sigma.value)
-    return layout.MultiVector(res_val)
+    return layout.MultiVector(res_val, copy=False)
 
 
 @numba.njit
@@ -1090,7 +1090,7 @@ def val_annihilate_k(K_val, C_val):
 
 def annihilate_k(K, C):
     """ Removes K from C = KX via (K[0] - K[4])*C """
-    return layout.MultiVector(val_annihilate_k(K.value, C.value))
+    return layout.MultiVector(val_annihilate_k(K.value, C.value), copy=False)
 
 
 @numba.njit
@@ -1275,7 +1275,7 @@ def motor_between_rounds(X1, X2):
     Calculate the motor between any pair of rounds of the same grade
     Line up the carriers, then line up the centers
     """
-    return layout.MultiVector(val_motor_between_rounds(X1.value, X2.value))
+    return layout.MultiVector(val_motor_between_rounds(X1.value, X2.value), copy=False)
 
 
 @numba.njit
@@ -1333,7 +1333,7 @@ def motor_between_objects(X1, X2):
     """
     Calculates a motor that takes X1 to X2
     """
-    return layout.MultiVector(val_motor_between_objects(X1.value, X2.value))
+    return layout.MultiVector(val_motor_between_objects(X1.value, X2.value), copy=False)
 
 
 def calculate_S_over_mu(X1, X2):
@@ -1492,7 +1492,7 @@ def val_normalised(mv_val):
 
 def normalised(mv):
     """ fast version of the normal() function """
-    return layout.MultiVector(val_normalised(mv.value))
+    return layout.MultiVector(val_normalised(mv.value), copy=False)
 
 
 @numba.njit
@@ -1517,12 +1517,12 @@ def val_rotor_between_lines(L1_val, L2_val):
 
 def rotor_between_lines(L1, L2):
     """ return the rotor between two lines """
-    return layout.MultiVector(val_rotor_between_lines(L1.value, L2.value))
+    return layout.MultiVector(val_rotor_between_lines(L1.value, L2.value), copy=False)
 
 
 def rotor_between_planes(P1, P2):
     """ return the rotor between two planes """
-    return layout.MultiVector(val_rotor_rotor_between_planes(P1.value, P2.value))
+    return layout.MultiVector(val_rotor_rotor_between_planes(P1.value, P2.value), copy=False)
 
 
 @numba.njit
@@ -1614,7 +1614,7 @@ def random_circle():
     mv_a = val_random_euc_mv()
     mv_b = val_random_euc_mv()
     mv_c = val_random_euc_mv()
-    return layout.MultiVector(val_normalised(omt_func(omt_func(val_up(mv_a), val_up(mv_b)), val_up(mv_c))))
+    return layout.MultiVector(val_normalised(omt_func(omt_func(val_up(mv_a), val_up(mv_b)), val_up(mv_c))), copy=False)
 
 
 def random_sphere_at_origin():
@@ -1664,7 +1664,7 @@ def val_apply_rotor(mv_val, rotor_val):
 
 def apply_rotor(mv_in, rotor):
     """ Applies rotor to multivector in a fast way """
-    return layout.MultiVector(val_apply_rotor(mv_in.value, rotor.value))
+    return layout.MultiVector(val_apply_rotor(mv_in.value, rotor.value), copy=False)
 
 
 @numba.njit
@@ -1675,7 +1675,7 @@ def val_apply_rotor_inv(mv_val, rotor_val, rotor_val_inv):
 
 def apply_rotor_inv(mv_in, rotor, rotor_inv):
     """ Applies rotor to multivector in a fast way takes pre computed adjoint"""
-    return layout.MultiVector(val_apply_rotor_inv(mv_in.value, rotor.value, rotor_inv.value))
+    return layout.MultiVector(val_apply_rotor_inv(mv_in.value, rotor.value, rotor_inv.value), copy=False)
 
 
 @numba.njit
@@ -1698,14 +1698,14 @@ def val_convert_2D_polar_line_to_conformal_line(rho, theta):
     p1_val = val_convert_2D_point_to_conformal(x1, y1)
     p2_val = val_convert_2D_point_to_conformal(x2, y2)
     line_val = omt_func(omt_func(p1_val, p2_val), ninf_val)
-    line_val = line_val/abs(layout.MultiVector(line_val))
+    line_val = line_val/abs(layout.MultiVector(line_val), copy=False)
     return line_val
 
 
 def convert_2D_polar_line_to_conformal_line(rho, theta):
     """ Converts a 2D polar line to a conformal line """
     line_val = val_convert_2D_polar_line_to_conformal_line(rho, theta)
-    return layout.MultiVector(line_val)
+    return layout.MultiVector(line_val, copy=False)
 
 
 @numba.njit
@@ -1718,7 +1718,7 @@ def val_up(mv_val):
 
 def fast_up(mv):
     """ Fast up mapping """
-    return layout.MultiVector(val_up(mv.value))
+    return layout.MultiVector(val_up(mv.value), copy=False)
 
 
 @numba.njit
@@ -1743,14 +1743,14 @@ def val_down(mv_val):
 
 def fast_down(mv):
     """ A fast version of down() """
-    return layout.MultiVector(val_down(mv.value))
+    return layout.MultiVector(val_down(mv.value), copy=False)
 
 
 def val_distance_point_to_line(point, line):
     """
     Returns the euclidean distance between a point and a line
     """
-    return float(abs(layout.MultiVector(omt_func(point, line))))
+    return float(abs(layout.MultiVector(omt_func(point, line))), copy=False)
 
 
 @numba.njit
@@ -1764,7 +1764,7 @@ def val_convert_2D_point_to_conformal(x, y):
 
 def convert_2D_point_to_conformal(x, y):
     """ Convert a 2D point to conformal """
-    return layout.MultiVector(val_convert_2D_point_to_conformal(x, y))
+    return layout.MultiVector(val_convert_2D_point_to_conformal(x, y), copy=False)
 
 
 def distance_polar_line_to_euc_point_2d(rho, theta, x, y):
@@ -1789,7 +1789,7 @@ def fast_dual(a):
     """
     Fast dual
     """
-    return layout.MultiVector(dual_func(a.value))
+    return layout.MultiVector(dual_func(a.value), copy=False)
 
 
 class ConformalMVArray(cf.MVArray):
@@ -1854,7 +1854,7 @@ def from_value_array(value_array):
 
 
 v_dual = np.vectorize(fast_dual, otypes=[ConformalMVArray])
-v_new_mv = np.vectorize(lambda v: layout.MultiVector(v), otypes=[ConformalMVArray], signature='(n)->()')
+v_new_mv = np.vectorize(lambda v: layout.MultiVector(v), otypes=[ConformalMVArray], signature='(n)->()', copy=False)
 v_up = np.vectorize(fast_up, otypes=[ConformalMVArray])
 v_down = np.vectorize(fast_down, otypes=[ConformalMVArray])
 v_apply_rotor_inv = np.vectorize(apply_rotor_inv, otypes=[ConformalMVArray])

clifford/tools/g3c/cost_functions.py

@@ -62,7 +62,7 @@ def midpoint_and_error_of_line_cluster_eig(line_cluster):
 
     point_val = np.zeros(32)
     point_val[1:6] = np.matmul(mat2solve, start)
-    new_mv = layout.MultiVector(point_val)
+    new_mv = layout.MultiVector(point_val, copy=False)
     new_mv = normalise_n_minus_1((new_mv * einf * new_mv)(1))
     return new_mv, val_point_to_line_cluster_distance(new_mv.value, line_cluster_array)
 
@@ -84,7 +84,7 @@ def midpoint_and_error_of_line_cluster_svd(line_cluster):
 
     point_val = np.zeros(32)
     point_val[np.array(layout.gradeList) == grade_val] = v[:, 1]
-    new_mv = layout.MultiVector(point_val)
+    new_mv = layout.MultiVector(point_val, copy=False)
     # new_mv = normalise_n_minus_1(new_mv * einf * new_mv)
     new_point = normalise_n_minus_1(new_mv)  # up(down(new_mv) / 2)
     return new_point, val_point_to_line_cluster_distance(new_point.value, line_cluster_array)
@@ -98,7 +98,7 @@ def midpoint_and_error_of_line_cluster(line_cluster):
     """
     line_cluster_array = np.array([l.value for l in line_cluster])
     cp_val = val_midpoint_of_line_cluster(line_cluster_array)
-    return layout.MultiVector(cp_val), val_point_to_line_cluster_distance(cp_val, line_cluster_array)
+    return layout.MultiVector(cp_val), val_point_to_line_cluster_distance(cp_val, line_cluster_array, copy=False)
 
 
 def midpoint_and_error_of_line_cluster_grad(line_cluster):
@@ -109,7 +109,7 @@ def midpoint_and_error_of_line_cluster_grad(line_cluster):
     """
     line_cluster_array = np.array([l.value for l in line_cluster])
     cp_val = val_midpoint_of_line_cluster_grad(line_cluster_array)
-    return layout.MultiVector(cp_val), val_point_to_line_cluster_distance(cp_val, line_cluster_array)
+    return layout.MultiVector(cp_val), val_point_to_line_cluster_distance(cp_val, line_cluster_array, copy=False)
 
 
 def line_plane_cost(line, plane):