LinearGeometricFunctions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5efae962",
   "metadata": {},
   "source": [
    "## Linear Geometric Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "id": "71106e12",
   "metadata": {
    "code_folding": [
     11,
     28,
     43,
     66,
     71
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "from functools import reduce,wraps\n",
    "import warnings\n",
    "from clifford import op \n",
    "from clifford.tools import mat2Frame,frame2Mat,orthoMat2Versor,func2Mat,log_rotor\n",
    "from clifford.invariant_decomposition import bivector_split,log\n",
    "import numpy as np\n",
    "from scipy import e, pi \n",
    "\n",
    "cayley =  lambda x:(1-x)/(1+x)\n",
    "\n",
    "def outermorphic(f):\n",
    "    '''\n",
    "    convert a geometric function `f` which operates on vectors into a \n",
    "    geometric function which operates on multivectors, using outermorphism.\n",
    "    \n",
    "    useful as decorator\n",
    "    '''\n",
    "    def F(M): # M is a multivector\n",
    "        N = []\n",
    "        for blade in M.blades_list:            # decompose MV into blades\n",
    "            factors,scale  = blade.factorise() # decompose each blade into factors \n",
    "            y = [f(x) for x in factors]        # f(x) each factor \n",
    "            Y = reduce(op,y)*scale             # compile factors back into blades\n",
    "            N.append(Y)                        # compile blades back into MV\n",
    "        return sum(N)\n",
    "    return F \n",
    "\n",
    "class MatrixOperator(object):\n",
    "    def __init__(self, I,M):\n",
    "        '''\n",
    "        An abstract base class inhereted by the various matrix operators\n",
    "        '''\n",
    "        self.I = I \n",
    "        self.M = M\n",
    "        self.validate_M()\n",
    "    \n",
    "    @property\n",
    "    def rank(self):\n",
    "        return self.M.shape[0]\n",
    "    def validate_M(self):\n",
    "        pass\n",
    "      \n",
    "class Linear(MatrixOperator):\n",
    "    @property\n",
    "    def Me(self):\n",
    "        return (self.M + self.M.T)/2.\n",
    "    \n",
    "    @property\n",
    "    def Mo(self):\n",
    "        return (self.M - self.M.T)/2.\n",
    "    \n",
    "    @property\n",
    "    def symmetric(self):\n",
    "        return Symmetric(I= self.I ,M = self.Me)\n",
    "    \n",
    "    @property\n",
    "    def skewmetric(self):\n",
    "        return Skewmetric(I= self.I ,M = self.Mo)\n",
    "    \n",
    "    def as_f(self):\n",
    "        raise NotImplementedError('works if self.normal?') \n",
    "        g = self.symmetric.as_f()\n",
    "        h = self.skewmetric.as_f()\n",
    "        return lambda x: g(x) + h(x)\n",
    "    \n",
    "class Skewmetric(MatrixOperator):\n",
    "    def validate_M(self):\n",
    "        if not np.allclose(self.M,-self.M.T):\n",
    "            warnings.warn('M is not skewmetric')\n",
    "        \n",
    "    def as_bivector(self):\n",
    "        '''\n",
    "        convert a skewmetric matrix to its curling bivector\n",
    "        '''\n",
    "        B = mat2Frame(self.M,I=self.I)[0]\n",
    "        A = mat2Frame(np.eye(self.rank),I=self.I)[0]\n",
    "        F     = sum([a*b for a,b in zip(A,B)])/2\n",
    "        return F\n",
    "\n",
    "    def as_f(self):\n",
    "        '''\n",
    "        convert a skewmetric (antisymmetric) matrix into a\n",
    "        geometric function `f`\n",
    "        '''\n",
    "        F   = self.as_bivector()\n",
    "        f_x = lambda x: x|F\n",
    "        f   = outermorphic(f_x)\n",
    "        return f\n",
    "    \n",
    "class Symmetric(MatrixOperator):\n",
    "    def validate_M(self):\n",
    "        if not np.allclose(self.M,self.M.T):\n",
    "            warnings.warn('M is not symmetric')\n",
    "    def as_eigenframe(self):\n",
    "        w,v = np.linalg.eig(self.M)\n",
    "        return mat2Frame(np.diag(w)@v,I=self.I)[0]\n",
    "    \n",
    "    def as_vector_and_versor(self):\n",
    "        w,v = np.linalg.eig(self.M)\n",
    "        d   = sum([a*k for a,k in zip(w,self.I.basis())])\n",
    "        R,rs  = orthoMat2Versor(v,I= self.I)\n",
    "        return d,R\n",
    "        \n",
    "    def as_f(self):\n",
    "        '''\n",
    "        convert a symmetric matrix into a geometric function `f`\n",
    "        '''\n",
    "        #V = self.symmetric_mat_2_eigenframe(M)\n",
    "        #f_x = lambda x: sum([(x|k)/k*abs(k) for k in V])\n",
    "        w,v = np.linalg.eig(self.M)\n",
    "        V   = mat2Frame(v,I=self.I)[0]\n",
    "        f_x = lambda x: sum([(x|k)/k*a for k,a in zip(V,w)])\n",
    "        f   = outermorphic(f_x)\n",
    "        return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 334,
   "id": "9f4036eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from clifford import Cl\n",
    "\n",
    "n     = 3\n",
    "M     = np.random.randint(0,100,n**2).reshape(n,n)\n",
    "\n",
    "layout,blades = Cl(n)\n",
    "I      = layout.pseudoScalar\n",
    "linear = Linear(I=I, M=M)\n",
    "\n",
    "\n",
    "assert(np.allclose(func2Mat(linear.symmetric.as_f(),  I=I)[0], linear.Me))\n",
    "assert(np.allclose(func2Mat(linear.skewmetric.as_f(), I=I)[0], linear.Mo))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "id": "f18a2dac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1628.25]"
      ]
     },
     "execution_count": 335,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f   = linear.skewmetric.as_f()\n",
    "F   = linear.skewmetric.as_bivector()\n",
    "Fks = bivector_split(F)\n",
    "[f(Fk)/Fk for Fk in Fks] # should all be scalar. def of skewsymmetry and invariant eigenbivector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "id": "fd846922",
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1.65271^e1) + (0.32157^e2) + (0.10082^e3),\n",
       " (0.83076^e1) + (0.39892^e2) + (0.31929^e3) - (0.22079^e123))"
      ]
     },
     "execution_count": 362,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def symmetric_2_skewup1(M):\n",
    "    '''convert a symmetric matrix of rank `n` to a skewmetric matrix of rank `n+1`'''\n",
    "    rank = M.shape[0]\n",
    "    out= np.zeros([rank+1,rank+1])\n",
    "    out[0,1:]= np.diagonal(M)\n",
    "    out[1:,1:]=M\n",
    "    out = np.triu(out,1)\n",
    "    return out-out.T\n",
    "\n",
    "from clifford import Cl\n",
    "n=3\n",
    "\n",
    "lay,b = Cl(sig=[1]+list(ones(n)),firstIdx=0)\n",
    "locals().update(b)\n",
    "M       = np.random.rand(n,n)\n",
    "M       = .5*(M+M.T)\n",
    "lo      = Symmetric(I=lay.pseudoScalar*e0, M=M)\n",
    "hi      = Skewmetric(I=lay.pseudoScalar, M=symmetric_2_skewup1(lo.M))\n",
    "\n",
    "d,G = lo.as_vector_and_versor()\n",
    "#(d*e0),log_rotor(G)\n",
    "d,G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 375,
   "id": "5b62c348",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((array([[-0.19422567,  0.47846007,  0.85635994],\n",
       "         [ 0.28645756,  0.8625971 , -0.4169752 ],\n",
       "         [ 0.93819958, -0.16432349,  0.304597  ]]),\n",
       "  array([[0.87981623, 0.15912039, 0.19389111],\n",
       "         [0.15912039, 1.04978986, 0.43166266],\n",
       "         [0.19389111, 0.43166266, 0.88848037]])),\n",
       " array([[0.34719408, 0.47173146, 0.1823311 ],\n",
       "        [0.47173146, 0.06083838, 0.8554263 ],\n",
       "        [0.1823311 , 0.8554263 , 0.087329  ]]))"
      ]
     },
     "execution_count": 375,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy \n",
    "M = np.random.rand(n,n)\n",
    "l = Linear(M=M,I=I)\n",
    "scipy.linalg.polar(np.random.rand(n,n)), l.symmetric.M\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "id": "dbc076a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.21675 + (0.13881^e01) + (0.41163^e02) + (0.61551^e03) + (0.12185^e12) + (0.43234^e13) + (0.34528^e23) + (0.25391^e0123)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-1.43558 - (0.46538^e0123)"
      ]
     },
     "execution_count": 347,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    ",F = hi.as_bivector()\n",
    "G = log_rotor(cayley(F))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "id": "56bdc67e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.98771 - (0.66897^e01) + (0.1071^e02) + (0.12932^e12),\n",
       " 0.88617 - (0.6405^e01) - (0.0395^e02) + (0.32402^e12)]"
      ]
     },
     "execution_count": 236,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = hi.as_f() \n",
    "Ks = [k*e0 for k in lo.as_eigenframe()]\n",
    "\n",
    "[f(K)/K for K in Ks]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "id": "bc403f2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Layout([-1, 1, 1],\n",
       "        ids=BasisVectorIds.ordered_integers(3, first_index=0),\n",
       "        order=BasisBladeOrder.shortlex(3),\n",
       "        names=['', 'e0', 'e1', 'e2', 'e01', 'e02', 'e12', 'e012']),\n",
       " {'': 1,\n",
       "  'e0': (1^e0),\n",
       "  'e1': (1^e1),\n",
       "  'e2': (1^e2),\n",
       "  'e01': (1^e01),\n",
       "  'e02': (1^e02),\n",
       "  'e12': (1^e12),\n",
       "  'e012': (1^e012)})"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Cl(sig=[-1]+list(ones(2)),firstIdx=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "537891d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## commutator patterns\n",
    "def symmetric_M(n):\n",
    "    M     = np.random.randint(0,20,n**2).reshape(n,n)\n",
    "    M = (M+M.T)/2.\n",
    "    return M\n",
    "def skewmetric_M(n):\n",
    "    M     = np.random.randint(0,20,n**2).reshape(n,n)\n",
    "    M = (M-M.T)/2.\n",
    "    return M\n",
    "\n",
    "M1 = symmetric_M(3)    \n",
    "M2 = symmetric_M(3)    \n",
    "N1 = skewmetric_M(3)\n",
    "N2 = skewmetric_M(3)\n",
    "\n",
    "com = lambda x,y: x@y-y@x\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}