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

updated (and broke some) tools, add some docs. #389

Open
wants to merge 24 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
ea35eca
docfix: g3 turorial missing pi,e import
arsenovic Jun 25, 2020
da80796
Merge branch 'master' of https://github.com/pygae/clifford
arsenovic Sep 14, 2020
0ee58d2
Merge branch 'master' of https://github.com/pygae/clifford
arsenovic Dec 19, 2020
57ef74b
allow mat2Frame to take list of vectors
arsenovic Dec 19, 2020
f876b98
BREAK: changed call sig on tools.frame2Mat and mat2Frame. added test.…
arsenovic Dec 20, 2020
7d9569c
BREAK:merged erics additions
arsenovic Dec 28, 2020
97bbf64
Update docs/tutorials/g3-algebra-of-space.ipynb
arsenovic Dec 29, 2020
434dbd2
Update clifford/__init__.py
arsenovic Dec 30, 2020
e90b594
Update clifford/__init__.py
arsenovic Dec 30, 2020
78b545d
pep8d
arsenovic Dec 31, 2020
f9ff2dc
added some missing methods to docs
arsenovic Dec 31, 2020
6bd529c
Revert unwanted autopep8 changes
eric-wieser Dec 31, 2020
d283fa5
Update clifford/tools/__init__.py
arsenovic Jan 4, 2021
5d81208
docfixs
arsenovic Jan 4, 2021
fc44e78
added func2Mat
arsenovic Feb 3, 2021
eeaeb1b
added tutorial of matrices
arsenovic Feb 4, 2021
23dc2e5
added matrix tutorial to docs
arsenovic Feb 4, 2021
2595f28
fix matrix tutorial headings
arsenovic Feb 4, 2021
b36a713
updated matrix tutorial notebook
arsenovic Feb 4, 2021
0c05eeb
added randomIntMV and some docs
arsenovic Feb 5, 2021
045b006
more docs
arsenovic Feb 5, 2021
dbe768c
add randomMV and randomV to Multivector
arsenovic Feb 7, 2021
248d26e
add log import to tools
arsenovic Mar 11, 2021
bfbd85e
Merge branch 'master' into master
eric-wieser Jul 19, 2021
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
11 changes: 11 additions & 0 deletions clifford/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -354,6 +354,17 @@ def Cl(p=0, q=0, r=0, sig=None, names=None, firstIdx=1, mvClass=MultiVector):

The notation :math:`Cl_{p,q,r}` means that the algebra is :math:`p+q+r`-dimensional, with the first :math:`p` vectors with positive signature, the next :math:`q` vectors negative, and the final :math:`r` vectors with null signature.

Parameters
=============
p : int
number of positive signature basis vectors
q : int
number of negative signature basis vectors
r : int
number of zero signature basis vectors
sig, names, firstIdx
See the docs for :class:`clifford.Layout`. If ``sig`` is passed, then `p`, `q`, and `r` are ignored.

Returns
=======
layout : Layout
Expand Down
32 changes: 23 additions & 9 deletions clifford/test/test_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,17 @@
from clifford.tools import orthoFrames2Versor as of2v
from clifford._numba_utils import DISABLE_JIT

from clifford import tools


too_slow_without_jit = pytest.mark.skipif(
DISABLE_JIT, reason="test is too slow without JIT"
)


@unittest.skip("reason unknown")
# @unittest.skip("reason unknown")
eric-wieser marked this conversation as resolved.
Show resolved Hide resolved
class ToolsTests(unittest.TestCase):

@unittest.skip("reason unknown")
def checkit(self, p, q):
# p, q =4,0
N = p + q
Expand All @@ -27,9 +29,9 @@ def checkit(self, p, q):
# create frame
A = layout.randomV(n=N)
# create Rotor
R = 5.*layout.randomRotor()
R = 5. * layout.randomRotor()
# create rotated frame
B = [R*a*~R for a in A]
B = [R * a * ~R for a in A]

# find versor from both frames
R_found, rs = of2v(A, B)
Expand All @@ -38,8 +40,9 @@ def checkit(self, p, q):
self.assertTrue(R == R_found or R == -R_found)

# Determined Versor implements desired transformation
self.assertTrue([R_found*a*~R_found for a in A] == B)
self.assertTrue([R_found * a * ~R_found for a in A] == B)

@unittest.skip("reason unknown")
def testOrthoFrames2VersorEuclidean(self):
for p, q in [(2, 0), (3, 0), (4, 0)]:
self.checkit(p=p, q=q)
Expand All @@ -54,6 +57,15 @@ def testOrthoFrames2VersorBalanced(self):
for p, q in [(2, 2)]:
self.checkit(p=p, q=q)

def testframe2Mat(self):
for N in [2, 3, 4]:
l, b = Cl(N)
X = np.random.rand((N**2)).reshape(N, N)
I = l.pseudoScalar
B, I = tools.mat2Frame(X, I=I)
X_, I = tools.frame2Mat(B=B, I=I)
testing.assert_almost_equal(X, X_)


class G3ToolsTests(unittest.TestCase):

Expand Down Expand Up @@ -104,9 +116,10 @@ def test_generate_rotation_rotor_and_angle(self):
euc_vector_n = random_unit_vector()
theta = angle_between_vectors(euc_vector_m, euc_vector_n)

rot_rotor = generate_rotation_rotor(theta, euc_vector_m, euc_vector_n)
rot_rotor = generate_rotation_rotor(
theta, euc_vector_m, euc_vector_n)
v1 = euc_vector_m
v2 = rot_rotor*euc_vector_m*~rot_rotor
v2 = rot_rotor * euc_vector_m * ~rot_rotor
theta_return = angle_between_vectors(v1, v2)

testing.assert_almost_equal(theta_return, theta)
Expand All @@ -124,7 +137,7 @@ def test_find_rotor_aligning_vectors(self):
u_list = [random_euc_mv() for i in range(50)]
for i in range(100):
r = random_rotation_rotor()
v_list = [r*u*~r for u in u_list]
v_list = [r * u * ~r for u in u_list]
r_2 = rotor_align_vecs(u_list, v_list)
print(r_2)
print(r)
Expand Down Expand Up @@ -194,7 +207,8 @@ def test_GADelaunay_facets(self):
from clifford.tools.g3c import random_conformal_point, project_points_to_plane
from clifford.tools.point_processing import GADelaunay
point_list = [random_conformal_point() for i in range(100)]
point_list_flat = project_points_to_plane(point_list, (up(0)^up(e1)^up(e2)^einf).normal())
point_list_flat = project_points_to_plane(
point_list, (up(0) ^ up(e1) ^ up(e2) ^ einf).normal())
hull = GADelaunay(point_list_flat, hull_dims=2)
facets = hull.conformal_facets()

Expand Down
74 changes: 53 additions & 21 deletions clifford/tools/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,11 @@
-----------------------------------------------------------

Given two frames that are related by a orthogonal transform, we seek a rotor
which enacts the transform. Details of the mathematics and pseudo-code used the
which enacts the transform. Details of the mathematics and pseudo-code used to
create the algorithms below can be found at Allan Cortzen's website,
:cite:`ctz-frames`.


There are also some helper functions which can be used to translate matrices
into GA frames, so an orthogonal (or complex unitary) matrix can be directly
translated into a Versor.
Expand All @@ -45,7 +46,9 @@
orthoFrames2Versor
orthoMat2Versor
mat2Frame
frame2Mat
omoh
rotor_decomp

"""

Expand All @@ -54,14 +57,16 @@
from typing import Union, Optional, List, Tuple
from math import sqrt
from numpy import eye, array, sign, zeros, sin, arccos

import numpy as np
from .. import Cl, gp, Frame, MultiVector, Layout
from .. import Cl, gp, Frame, MultiVector
from .. import eps as global_eps

from warnings import warn


def omoh(A: Union[Frame, List[MultiVector]], B: Union[Frame, List[MultiVector]]) -> np.ndarray:
def omoh(A: Union[Frame, List[MultiVector]],
B: Union[Frame, List[MultiVector]]) -> np.ndarray:
r'''
Determines homogenization scaling for two :class:`~clifford.Frame`\ s related by a Rotor

Expand Down Expand Up @@ -112,8 +117,8 @@ def omoh(A: Union[Frame, List[MultiVector]], B: Union[Frame, List[MultiVector]])


def mat2Frame(A: np.ndarray,
layout: Optional[Layout] = None,
is_complex: bool = None) -> Tuple[List[MultiVector], Layout]:
I: Optional[MultiVector] = None,
is_complex: bool = None) -> Tuple[List[MultiVector], MultiVector]:
'''
Translates a (possibly complex) matrix into a real vector frame

Expand All @@ -130,12 +135,19 @@ def mat2Frame(A: np.ndarray,
A : ndarray
MxN matrix representing vectors

I: None, pseudoscalar of the frame
arsenovic marked this conversation as resolved.
Show resolved Hide resolved
if none we generate an algebra of Gn, if layout we take the
vector basis from that, and if its a list we will assume its
a vector basis.


Returns
-------
a : list of clifford.MultiVector
The resulting vectors
layout : clifford.Layout
The layout of the vectors in ``a``.
I : clifford.MultiVector
The blade holding the vectors in ``a``.

'''

# TODO: could simplify this by just implementing the real case and then
Expand All @@ -155,10 +167,10 @@ def mat2Frame(A: np.ndarray,
N = N * 2
M = M * 2

if layout is None:
if I is None:
layout, blades = Cl(M)

e_ = layout.basis_vectors_lst[:M]
I = layout.pseudoScalar
e_ = I.basis()

a = [0 ^ e_[0]] * N

Expand All @@ -179,19 +191,37 @@ def mat2Frame(A: np.ndarray,
a[n_ + 1] = (a[n_ + 1]) \
+ ((-A[m, n].imag) ^ e_[m_]) \
+ ((A[m, n].real) ^ e_[m_ + 1])
return a, layout
return a, I


def frame2Mat(B, A=None, is_complex=None):
def frame2Mat(B, A=None, I=None, is_complex=None):
'''
convert a list of vectors to a matrix

B : list
eric-wieser marked this conversation as resolved.
Show resolved Hide resolved
a list of vectors that have been transformed
A : None, list of vectors
a list of vectors in their initial state. if none we assume
orthonormal basis given by B.pseudoScalar, or I
I : Multivector, None
pseudoscalar of the space. if None, we use B.pseudoScalar
is_complex: Bool
do you want a complex matrix?

'''
if is_complex is not None:
raise NotImplementedError()

if I is None:
I = B[0].pseudoScalar
if A is None:
# assume we have orthonormal initial frame
A = B[0].layout.basis_vectors_lst
A = I.basis()

# you need float() due to bug in clifford
M = [float(b | a) for b in B for a in A]
M = [float(b | a) for a in A for b in B]
M = array(M).reshape(len(B), len(B))
return M, I


def orthoFrames2Versor_dist(A, B, eps=None):
Expand All @@ -207,7 +237,8 @@ def orthoFrames2Versor_dist(A, B, eps=None):

'''
# TODO: should we test to see if A and B are related by rotation?
# TODO: implement reflect/rotate based on distance (as in:cite:`ctz-frames`)
# TODO: implement reflect/rotate based on distance (as
# in:cite:`ctz-frames`)

# keep copy of original frames
A = A[:]
Expand Down Expand Up @@ -401,7 +432,7 @@ def orthoFrames2Versor(B, A=None, delta: float = 1e-3,
return R, r_list


def orthoMat2Versor(A, eps=None, layout=None, is_complex=None):
def orthoMat2Versor(A, eps=None, I=None, is_complex=None):
'''
Translates an orthogonal (or unitary) matrix to a Versor

Expand All @@ -413,16 +444,17 @@ def orthoMat2Versor(A, eps=None, layout=None, is_complex=None):
------------

'''
B, layout = mat2Frame(A, layout=layout, is_complex=is_complex)
B, layout = mat2Frame(A, I=I, is_complex=is_complex)
N = len(B)

# if (A.dot(A.conj().T) -eye(N/2)).max()>eps:
# warn('A doesnt appear to be a rotation. ')
A, layout = mat2Frame(eye(N), layout=layout, is_complex=False)
A, dum = mat2Frame(eye(N), I=I, is_complex=False)
return orthoFrames2Versor(A=A, B=B, eps=eps)


def rotor_decomp(V: MultiVector, x: MultiVector) -> Tuple[MultiVector, MultiVector]:
def rotor_decomp(
V: MultiVector, x: MultiVector) -> Tuple[MultiVector, MultiVector]:
'''
Rotor decomposition of rotor V

Expand Down Expand Up @@ -456,7 +488,7 @@ def rotor_decomp(V: MultiVector, x: MultiVector) -> Tuple[MultiVector, MultiVect


def sinc(x):
return sin(x)/x
return sin(x) / x


def log_rotor(V):
Expand All @@ -470,4 +502,4 @@ def log_rotor(V):
return log(float(V(0)))
# numpy's trig correctly chooses hyperbolic or not with Complex args
theta = arccos(complex(V(0)))
return V(2)/sinc(theta).real
return V(2) / sinc(theta).real
3 changes: 1 addition & 2 deletions docs/tutorials/euler-angles.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -50,8 +50,7 @@
"2. rotate about the rotated $e_1$-axis, call it $e_1^{'}$\n",
"3. rotate about the twice rotated axis of $e_3$-axis, call it $e_3^{''}$\n",
"\n",
"So the elemental rotations are about $e_3, e_{1}^{'}, e_3^{''}$-axes, respectively.",
"\n",
"So the elemental rotations are about $e_3, e_{1}^{'}, e_3^{''}$-axes, respectively.\n",
"![](../_static/Euler2a.gif)\n",
"\n",
"(taken from wikipedia)"
Expand Down
12 changes: 11 additions & 1 deletion docs/tutorials/g3-algebra-of-space.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -861,7 +861,17 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.0"
"version": "3.7.3"
eric-wieser marked this conversation as resolved.
Show resolved Hide resolved
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"toc_cell": false,
"toc_position": {},
"toc_section_display": "block",
"toc_window_display": false
}
},
"nbformat": 4,
Expand Down