diff --git a/clifford/__init__.py b/clifford/__init__.py index ffcba8b9..5372948a 100644 --- a/clifford/__init__.py +++ b/clifford/__init__.py @@ -70,6 +70,8 @@ grade_obj randomMV + randomIntMV + randomV """ @@ -393,6 +395,9 @@ def randomMV( -------- >>> randomMV(layout, min=-2.0, max=2.0, grades=None, uniform=None, n=2) # doctest: +SKIP + See Also + ---------- + randomIntMV """ if n > 1: @@ -421,6 +426,21 @@ def randomMV( return mv +def randomIntMV(layout, min=-10, max=10, uniform=np.random.randint, + *args,**kw): + ''' + Random MultiVectors with given layout. + + Coefficients are between min and max, and if grades is a list of integers, + only those grades will be non-zero. + + This is just a wrapper of `randomMV`, which sets `uniform` and `min`/`max` + + See Also + ---------- + randomMV + ''' + return randomMV(layout=layout, min=min,max=max,uniform=uniform, *args, **kw) def conformalize(layout: Layout, added_sig=[1, -1], *, mvClass=MultiVector, **kwargs): ''' diff --git a/clifford/_layout.py b/clifford/_layout.py index ebb5a43a..428534cb 100644 --- a/clifford/_layout.py +++ b/clifford/_layout.py @@ -764,6 +764,15 @@ def randomMV(self, n=1, **kwargs) -> MultiVector: ''' return cf.randomMV(layout=self, n=n, **kwargs) + def randomIntMV(self, n=1, **kwargs) -> MultiVector: + ''' + Convenience method to create a random multivector wiht integer + values for components. Nice for reading, + + see `clifford.randomIntMV` for details + ''' + return cf.randomIntMV(layout=self, n=n, **kwargs) + def randomV(self, n=1, **kwargs) -> MultiVector: ''' generate n random 1-vector s diff --git a/clifford/_multivector.py b/clifford/_multivector.py index e8ecb1ac..93da0aaa 100644 --- a/clifford/_multivector.py +++ b/clifford/_multivector.py @@ -897,6 +897,38 @@ def conjugate(self) -> 'MultiVector': """ return (~self).gradeInvol() + def randomMV(self,*args, **kw): + ''' + create a random MV within space of self + + This is a wrapper for `randomMV` + + See Also + --------- + Layout.randomMV + randomMV + ''' + if not self.isBlade(): + raise ValueError('I must be a blade') + return cf.randomMV(layout=self.layout, *args, **kw)(self) + + def randomV(self,*args, **kw): + ''' + create a random MV within space of self + + This is a wrapper for `Layout.randomMV` + + See Also + --------- + Layout.randomMV + ''' + if not self.isBlade(): + raise ValueError('I must be a blade') + kw.update(grades=[1]) + return cf.randomMV(layout=self.layout, *args, **kw)(self) + + + # Subspace operations def project(self, other) -> 'MultiVector': r"""Projects the multivector onto the subspace represented by this blade. diff --git a/clifford/test/test_tools.py b/clifford/test/test_tools.py index 77a6e7f3..09503cae 100644 --- a/clifford/test/test_tools.py +++ b/clifford/test/test_tools.py @@ -7,14 +7,15 @@ from clifford.tools import orthoFrames2Versor as of2v from clifford._numba_utils import DISABLE_JIT +from clifford import tools from . import rng # noqa: F401 + too_slow_without_jit = pytest.mark.skipif( DISABLE_JIT, reason="test is too slow without JIT" ) -@pytest.mark.skip(reason="unknown") class TestTools: def checkit(self, p, q, rng): # noqa: F811 @@ -39,6 +40,7 @@ def checkit(self, p, q, rng): # noqa: F811 # Determined Versor implements desired transformation 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) @@ -53,6 +55,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 TestG3Tools: diff --git a/clifford/tools/__init__.py b/clifford/tools/__init__.py index c25ee0e4..2427352e 100644 --- a/clifford/tools/__init__.py +++ b/clifford/tools/__init__.py @@ -45,7 +45,10 @@ orthoFrames2Versor orthoMat2Versor mat2Frame + frame2Mat + func2Mat omoh + rotor_decomp """ @@ -53,15 +56,16 @@ from typing import Union, Optional, List, Tuple from math import sqrt -from numpy import eye, array, sign, zeros, sin, arccos +from numpy import eye, array, sign, zeros, sin, arccos,log 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 @@ -112,8 +116,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 @@ -130,12 +134,17 @@ def mat2Frame(A: np.ndarray, A : ndarray MxN matrix representing vectors + I : MultiVector + 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 @@ -155,10 +164,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 @@ -179,20 +188,59 @@ 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, I=None, is_complex=None): + ''' + convert a list of vectors to a matrix -def frame2Mat(B, A=None, is_complex=None): + Parameters + ------------ + B : list + 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 func2Mat(f,I): + ''' + Convert a function to a matrix by acting on standard basis + + Parameters + --------------- + f : function + function that maps vectors to vectors + I : MultiVector + psuedoscalar of basis + + See Also + --------- + frame2Mat + ''' + A = I.basis() + B = [f(a) for a in A] + return frame2Mat(B=B, A=A,I=I) def orthoFrames2Versor_dist(A, B, eps=None): ''' @@ -401,7 +449,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 @@ -411,14 +459,22 @@ def orthoMat2Versor(A, eps=None, layout=None, is_complex=None): Parameters ------------ + A : matrix + matrix to be transformed + eps : number + tolerance + I : MultiVector + GA of A + is_complex : boolean + is A complex? ''' - 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) diff --git a/docs/index.rst b/docs/index.rst index fc7942fa..78827186 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -52,6 +52,7 @@ Scalars, vectors, and higher-grade entities can be mixed freely and consistently tutorials/PerformanceCliffordTutorial tutorials/cga/index tutorials/linear-transformations + tutorials/MatrixRepresentationsOfGeometricFunctions tutorials/apollonius-cga-augmented .. toctree:: diff --git a/docs/tutorials/MatrixRepresentationsOfGeometricFunctions.ipynb b/docs/tutorials/MatrixRepresentationsOfGeometricFunctions.ipynb new file mode 100644 index 00000000..cc35ac68 --- /dev/null +++ b/docs/tutorials/MatrixRepresentationsOfGeometricFunctions.ipynb @@ -0,0 +1,480 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "valuable-satin", + "metadata": {}, + "source": [ + "# Matrix Representations of Geometric Functions\n", + "More info can be found in Chapter 5 of \"New foundations for classical mechanics\", by David Hestenes.\n", + "\n", + "This notebook shows how some matrix groups can be represented in geometric algebra. Not as spinors in CGA, or PGA, just as functions in plain old Geometric Algebra. \n", + "This is done by: \n", + "\n", + " 1. Creating a geometric function\n", + " 2. Apply it to an orthonormal frame\n", + " 3. Convert the resultant frame into a matrix \n", + "\n", + "The matrix is defined as the inner product of each basis element of original and transformed frame. \n", + "\n", + "$$M_{ij} = a_i\\cdot b_j = a_i\\cdot f(a)_j $$ \n", + "\n", + "(or vice-versa with the i,j, you get the point). Since we are going to do this repeatedly, define a `func2Mat()`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "exterior-cutting", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np \n", + "\n", + "def func2Mat(f,I):\n", + " '''\n", + " Convert a function acting on a vector into a matrix, given \n", + " the space defined by psuedoscalar I\n", + " '''\n", + " A = I.basis()\n", + " B = [f(a) for a in A]\n", + " M = [float(b | a) for a in A for b in B]\n", + " return np.array(M).reshape(len(B), len(B)) " + ] + }, + { + "cell_type": "markdown", + "id": "green-medicare", + "metadata": {}, + "source": [ + "Start with initializing a euclidean N-dimensional algebra and assign our pseudoscalar to $I$, pretty standard. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "absent-executive", + "metadata": {}, + "outputs": [], + "source": [ + "from clifford import Cl\n", + "from math import * \n", + "\n", + "l,b = Cl(3) # returns (layout,blades). you can change dimesion here\n", + "I = l.pseudoScalar " + ] + }, + { + "cell_type": "markdown", + "id": "opposite-footage", + "metadata": {}, + "source": [ + "## Anti-symmetric\n", + "This is so easy,\n", + "\n", + " $$x \\rightarrow x\\cdot B$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "seeing-miniature", + "metadata": {}, + "outputs": [], + "source": [ + "B = l.randomIntMV()(2) # we use randIntMV because its easier to read\n", + "f = lambda x:x|B\n", + "func2Mat(f,I=I)" + ] + }, + { + "cell_type": "markdown", + "id": "bridal-coordinator", + "metadata": {}, + "source": [ + "Whats the B? you can read its values straight off the matrix. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "healthy-sacramento", + "metadata": {}, + "outputs": [], + "source": [ + "B" + ] + }, + { + "cell_type": "markdown", + "id": "animal-conviction", + "metadata": {}, + "source": [ + "## Diagonal ( Directional Scaling)\n", + "\n", + "A bit awkward this one, but its made by projection onto each basis vector, then scaling the component by some amount. \n", + "\n", + "$$ x \\rightarrow \\sum{\\lambda_i (x\\cdot e_i) e_i} $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "internal-desperate", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "ls = range(1,len(I.basis())+1) # some dilation values (eigenvalues) \n", + "A = I.basis()\n", + "\n", + "\n", + "d = lambda x: sum([(x|a)/a*l for a,l in zip(A,ls)])\n", + "func2Mat(d,I=I)" + ] + }, + { + "cell_type": "markdown", + "id": "electrical-community", + "metadata": {}, + "source": [ + "## Orthgonal, Rotation\n", + "\n", + "$$ x\\rightarrow Rx\\tilde{R}$$\n", + "\n", + "where \n", + "\n", + "$$R=e^{B/2}$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "quarterly-arnold", + "metadata": {}, + "outputs": [], + "source": [ + "B = l.randomMV()(2)\n", + "R = e**(B/2)\n", + "r = lambda x: R*x*~R\n", + "func2Mat(r,I=I)" + ] + }, + { + "cell_type": "markdown", + "id": "authorized-tuesday", + "metadata": {}, + "source": [ + "The inverse of this is , \n", + "$$ x\\rightarrow \\tilde{R}xR $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "single-lodging", + "metadata": {}, + "outputs": [], + "source": [ + "rinv = lambda x: ~R*x*R # the inverse rotation \n", + "func2Mat(rinv,I=I)" + ] + }, + { + "cell_type": "markdown", + "id": "danish-focus", + "metadata": {}, + "source": [ + "## Orthogonal, Reflection \n", + "\n", + "$$ x \\rightarrow -axa^{-1} $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "atmospheric-lender", + "metadata": {}, + "outputs": [], + "source": [ + "a = l.randomIntMV()(1)\n", + "n = lambda x: -a*x/a\n", + "func2Mat(n,I=I)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "focal-velvet", + "metadata": {}, + "outputs": [], + "source": [ + "a" + ] + }, + { + "cell_type": "markdown", + "id": "premier-cooking", + "metadata": {}, + "source": [ + "Notice the determinant for reflection is -1, and for rotation is +1. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "proper-startup", + "metadata": {}, + "outputs": [], + "source": [ + "from numpy.linalg import det \n", + "det(func2Mat(n,I=I)), det(func2Mat(r,I=I))" + ] + }, + { + "cell_type": "markdown", + "id": "powered-intelligence", + "metadata": {}, + "source": [ + "## Symmetric \n", + "This can be built up from the functions we just defined ie Rotation\\*Dilation/Rotation\n", + "\n", + "$$ x \\rightarrow r(d(r^{-1}(x)) $$\n", + "\n", + "which if you write it out, looks kind of dumb\n", + "\n", + "$$ x \\rightarrow R[\\sum{\\lambda_i (\\tilde{R}x R\\cdot e_i) e_i]R} $$\n", + "\n", + "So, the antisymmetric matrix is interpreted as a set dilations about a some orthogonal frame rotated from the basis (what basis,eh? exactly what basis!). \n", + "\n", + "\n", + "More generally we could include reflection in the $R$ too." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "postal-crazy", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "g = lambda x: r(d(rinv(x)))\n", + "func2Mat(g,I=I)" + ] + }, + { + "cell_type": "markdown", + "id": "friendly-rebel", + "metadata": {}, + "source": [ + "## Eigen stuffs\n", + "By definition the eigen-stuff is the invariants of the transformation, sometimes this is a vector, and other times it is a plane. " + ] + }, + { + "cell_type": "markdown", + "id": "integrated-front", + "metadata": {}, + "source": [ + "### Rotation" + ] + }, + { + "cell_type": "markdown", + "id": "hearing-principle", + "metadata": {}, + "source": [ + "The eigen blades of a rotation are really the axis and plane of rotation. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "entitled-clearing", + "metadata": {}, + "outputs": [], + "source": [ + "from numpy.linalg import eig\n", + "\n", + "vals, vecs = eig(func2Mat(r,I=I))\n", + "np.round(vecs,3)" + ] + }, + { + "cell_type": "markdown", + "id": "independent-terminology", + "metadata": {}, + "source": [ + "If you checkout the real column, and compare this to the bivector which generated this rotation (aka the generator), after its been normalized " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "insured-karma", + "metadata": {}, + "outputs": [], + "source": [ + "B/(abs(B))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "median-planet", + "metadata": {}, + "outputs": [], + "source": [ + "B" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "resident-detroit", + "metadata": {}, + "outputs": [], + "source": [ + "vals" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "frank-alexandria", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "cos(abs(B)), sin(abs(B))" + ] + }, + { + "cell_type": "markdown", + "id": "spare-fleece", + "metadata": {}, + "source": [ + "### Symmetric\n", + "For the symmetric matrix, the invariant thing the orthonormal frame along which the dilations take place" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "turkish-brooklyn", + "metadata": {}, + "outputs": [], + "source": [ + "vals, vecs = eig(func2Mat(g,I=I))\n", + "np.round(vecs,5).T" + ] + }, + { + "cell_type": "markdown", + "id": "rolled-hebrew", + "metadata": {}, + "source": [ + "This is easily found by using the rotation part of the symmetric operator, " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "atomic-controversy", + "metadata": {}, + "outputs": [], + "source": [ + "[R*a*~R for a in I.basis()]" + ] + }, + { + "cell_type": "markdown", + "id": "variable-honor", + "metadata": {}, + "source": [ + "## Primitive Visualization in 2D" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fatal-french", + "metadata": {}, + "outputs": [], + "source": [ + "from pylab import linspace, plot,axis,legend\n", + "\n", + "def plot_ps(ps,**kw):\n", + " x = [p[e1] for p in ps]\n", + " y = [p[e2] for p in ps]\n", + " plot(x,y, marker='o',ls='',**kw)\n", + "\n", + " \n", + "l,b = Cl(2)\n", + "locals().update(b)\n", + "I = l.pseudoScalar\n", + "\n", + "## define function of interest \n", + "B = l.randomMV()(2)\n", + "R = e**(B/2)\n", + "f = lambda x: R*x*~R\n", + "\n", + "## loop though cartesian grid and apply f, \n", + "ps,qs=[],[]\n", + "for x in linspace(-1,1,11):\n", + " for y in linspace(-1,1,11):\n", + " p = x*e1 + y*e2\n", + " q = f(p)\n", + " ps.append(p)\n", + " qs.append(q)\n", + "\n", + "plot_ps(ps,label='before')\n", + "plot_ps(qs,label='after')\n", + "axis('equal')\n", + "legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "specific-myrtle", + "metadata": {}, + "outputs": [], + "source": [ + "func2Mat(f,I =I )" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.7" + }, + "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, + "nbformat_minor": 5 +} diff --git a/docs/tutorials/euler-angles.ipynb b/docs/tutorials/euler-angles.ipynb index 4ebe0e25..cd5caa82 100644 --- a/docs/tutorials/euler-angles.ipynb +++ b/docs/tutorials/euler-angles.ipynb @@ -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)"