diff --git a/clifford/_multivector.py b/clifford/_multivector.py index 506cb680..067d5d93 100644 --- a/clifford/_multivector.py +++ b/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() diff --git a/clifford/tools/g3/__init__.py b/clifford/tools/g3/__init__.py index 67fc1238..a60ae959 100644 --- a/clifford/tools/g3/__init__.py +++ b/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): diff --git a/clifford/tools/g3c/__init__.py b/clifford/tools/g3c/__init__.py index 3db4b10e..5c6d6c45 100644 --- a/clifford/tools/g3c/__init__.py +++ b/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,7 +618,7 @@ 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): @@ -626,7 +626,7 @@ 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]) diff --git a/clifford/tools/g3c/cost_functions.py b/clifford/tools/g3c/cost_functions.py index 4ab1a0e0..582b7222 100644 --- a/clifford/tools/g3c/cost_functions.py +++ b/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): diff --git a/clifford/tools/g3c/object_fitting.py b/clifford/tools/g3c/object_fitting.py index 291f3864..39bc8d62 100644 --- a/clifford/tools/g3c/object_fitting.py +++ b/clifford/tools/g3c/object_fitting.py @@ -49,7 +49,7 @@ def fit_circle(point_list): """ Performs Leo Dorsts circle fitting technique """ - return layout.MultiVector(val_fit_circle(np.array([p.value for p in point_list]))) + return layout.MultiVector(val_fit_circle(np.array([p.value for p in point_list])), copy=False) @numba.njit @@ -85,7 +85,7 @@ def fit_line(point_list): """ Does line fitting with combo J.Lasenbys method and L. Dorsts """ - return layout.MultiVector(val_fit_line(np.array([p.value for p in point_list]))) + return layout.MultiVector(val_fit_line(np.array([p.value for p in point_list])), copy=False) @numba.njit @@ -124,7 +124,7 @@ def fit_sphere(point_list): """ Performs Leo Dorsts sphere fitting technique """ - return layout.MultiVector(val_fit_sphere(np.array([p.value for p in point_list]))) + return layout.MultiVector(val_fit_sphere(np.array([p.value for p in point_list])), copy=False) @numba.njit @@ -158,4 +158,4 @@ def fit_plane(point_list): """ Does plane fitting with combo J.Lasenbys method and L. Dorsts """ - return layout.MultiVector(val_fit_plane(np.array([p.value for p in point_list]))) + return layout.MultiVector(val_fit_plane(np.array([p.value for p in point_list])), copy=False) diff --git a/clifford/tools/g3c/rotor_estimation.py b/clifford/tools/g3c/rotor_estimation.py index 2d7f4844..c11e8942 100644 --- a/clifford/tools/g3c/rotor_estimation.py +++ b/clifford/tools/g3c/rotor_estimation.py @@ -84,7 +84,7 @@ def extract_rotor_from_TRS_mat_est(mat_est): Given a matrix of the form [TRS_left@~TRS_right] returns TRS """ sph = (up(e1)^up(-e1)^up(e2)^up(e3)).normal()*I5 - sph2 = layout.MultiVector(mat_est@sph.value).normal() + sph2 = layout.MultiVector(mat_est@sph.value, copy=False).normal() t = down((sph2*einf*sph2)(1)) T = generate_translation_rotor(t) S = generate_dilation_rotor(get_radius_from_sphere(sph2*I5)/get_radius_from_sphere(sph*I5)) @@ -192,7 +192,7 @@ def de_keninck_twist(Y, X, guess=None): """ if guess is None: guess = (1.0 + 0*e1) - return layout.MultiVector(val_de_keninck_twist(Y.value, X.value, guess.value)) + return layout.MultiVector(val_de_keninck_twist(Y.value, X.value, guess.value), copy=False) def average_estimator(reference_model, query_model): diff --git a/clifford/tools/g3c/rotor_parameterisation.py b/clifford/tools/g3c/rotor_parameterisation.py index 940bce93..705e3c28 100644 --- a/clifford/tools/g3c/rotor_parameterisation.py +++ b/clifford/tools/g3c/rotor_parameterisation.py @@ -195,8 +195,8 @@ def ga_exp(B): Fast implementation of the translation and rotation specific exp function """ if np.sum(np.abs(B.value)) < np.finfo(float).eps: - return layout.MultiVector(unit_scalar_mv.value) - return layout.MultiVector(val_exp(B.value)) + return layout.MultiVector(unit_scalar_mv.value, copy=False) + return layout.MultiVector(val_exp(B.value), copy=False) def interpolate_TR_rotors(R_n_plus_1, R_n, interpolation_fraction): @@ -295,7 +295,7 @@ def TR_biv_params_to_rotor(x): """ Converts between the parameters of a TR bivector and the rotor that it is generating """ - return layout.MultiVector(val_TR_biv_params_to_rotor(x)) + return layout.MultiVector(val_TR_biv_params_to_rotor(x), copy=False) rotorconversion = TR_biv_params_to_rotor @@ -318,4 +318,4 @@ def R_biv_params_to_rotor(x): """ Converts between the parameters of a R bivector and the rotor that it is generating """ - return layout.MultiVector(val_R_biv_params_to_rotor(x)) + return layout.MultiVector(val_R_biv_params_to_rotor(x), copy=False) diff --git a/docs/changelog.rst b/docs/changelog.rst index 7da11288..7c61b442 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -18,6 +18,8 @@ Changes in 1.3.x * Faster algebra construction. ``Cl(3)`` is now 100\ |times| faster, and ``Cl(6)`` is 20\ |times| faster. This is achieved by deferring product JIT compilation until the product is used for the first time. + * The :class:`MultiVector` constructor gained a ``copy`` argument, which can be + set to ``False`` to allow construction without copying. Bugs fixed ---------- diff --git a/docs/tutorials/PerformanceCliffordTutorial.ipynb b/docs/tutorials/PerformanceCliffordTutorial.ipynb index f78aa6d0..d5a37cb4 100644 --- a/docs/tutorials/PerformanceCliffordTutorial.ipynb +++ b/docs/tutorials/PerformanceCliffordTutorial.ipynb @@ -188,7 +188,7 @@ "outputs": [], "source": [ "def apply_rotor_faster(R,mv):\n", - " return layout.MultiVector(layout.gmt_func(R.value,layout.gmt_func(mv.value,layout.adjoint_func(R.value))) )" + " return layout.MultiVector(layout.gmt_func(R.value,layout.gmt_func(mv.value,layout.adjoint_func(R.value))), copy=False)" ] }, { @@ -253,7 +253,7 @@ " return gmt_func(R_val,gmt_func(mv_val,adjoint_func(R_val)))\n", "\n", "def apply_rotor_wrapped(R,mv):\n", - " return cf.MultiVector(layout,apply_rotor_val(R.value,mv.value))" + " return layout.MultiVector(apply_rotor_val(R.value, mv.value), copy=False)" ] }, { @@ -312,7 +312,7 @@ " return gmt_func(R_val,gmt_func(mv_val,adjoint_func(R_val)))\n", "\n", "def apply_rotor_wrapped_numba(R,mv):\n", - " return cf.MultiVector(layout,apply_rotor_val_numba(R.value,mv.value))" + " return layout.MultiVector(apply_rotor_val_numba(R.value,mv.value), copy=False)" ] }, { @@ -374,7 +374,7 @@ " return -dual_val( omt_func( dual_val(mv_a_val) , dual_val(mv_b_val)) )\n", "\n", "def meet_wrapped(mv_a,mv_b):\n", - " return cf.layout.MultiVector(meet_val(mv_a.value, mv_b.value))\n", + " return layout.MultiVector(meet_val(mv_a.value, mv_b.value), copy=False)\n", "\n", "sphere = (up(0)^up(e1)^up(e2)^up(e3)).normal()\n", "print(sphere.meet(line_one).normal().normal())\n", @@ -555,7 +555,7 @@ " return left_rotor_mult(R_val,right_rotor_mult(mv_val,adjoint_func(R_val)))\n", "\n", "def sparse_apply_rotor(R,mv):\n", - " return cf.MultiVector(layout,sparse_apply_rotor_val(R.value,mv.value))" + " return layout.MultiVector(sparse_apply_rotor_val(R.value, mv.value), copy=False)" ] }, { @@ -619,7 +619,7 @@ " return -dual_sparse_val( sparse_omt_2_1( dual_sparse_val(mv_a_val) , dual_sparse_val(mv_b_val)) )\n", "\n", "def meet_sparse_3_4(mv_a,mv_b):\n", - " return cf.layout.MultiVector(meet_sparse_3_4_val(mv_a.value, mv_b.value))\n", + " return layout.MultiVector(meet_sparse_3_4_val(mv_a.value, mv_b.value), copy=False)\n", "\n", "print(sphere.meet(line_one).normal().normal())\n", "print(meet_sparse_3_4(line_one,sphere).normal())"