Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve performance of multivector arithmetic by avoiding copies #280

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions clifford/_multivector.py
Original file line number Diff line number Diff line change
Expand Up @@ -122,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` """
Expand Down Expand Up @@ -255,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}`"""
Expand Down Expand Up @@ -769,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()

Expand Down
2 changes: 1 addition & 1 deletion clifford/tools/g3/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
58 changes: 29 additions & 29 deletions clifford/tools/g3c/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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):
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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():
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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):
Expand All @@ -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):
Expand Down Expand Up @@ -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])
Expand Down
8 changes: 4 additions & 4 deletions clifford/tools/g3c/cost_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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)
Expand All @@ -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):
Expand All @@ -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):
Expand Down
8 changes: 4 additions & 4 deletions clifford/tools/g3c/object_fitting.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
4 changes: 2 additions & 2 deletions clifford/tools/g3c/rotor_estimation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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):
Expand Down
8 changes: 4 additions & 4 deletions clifford/tools/g3c/rotor_parameterisation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This might cause problems as the unit_scalar is a constant scoped to the whole file

Copy link
Member Author

@eric-wieser eric-wieser Mar 12, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed - however, I don't think this is any worse than the fact that clifford.g3c.e3 is global, and I could do something careless like

my_vector = e3
my_vector.value[:] = e2.value  # whoops, now e3 == e2

Perhaps I'll put my patch in to fix that first - the blades that are part of the predefined algebras should definitely be readonly. I'm not sure if layout.blades should be though.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

omg I hate it

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah I have no strong feeling about layout.blades but the predefined blades definitely have to be readonly

Copy link
Member Author

@eric-wieser eric-wieser Mar 12, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

eo.value[:] = layout.scalar.value[:]
print("Your CGA is now somewhat PGA")

return layout.MultiVector(val_exp(B.value), copy=False)


def interpolate_TR_rotors(R_n_plus_1, R_n, interpolation_fraction):
Expand Down Expand Up @@ -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
Expand All @@ -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)