Compare right inverse jit

CompareRightInverseJit.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "from clifford.g3c import *\nimport numpy as np\nfrom numpy import linalg\nimport numba\nfrom clifford import _myDot\nimport clifford as cf\n_eps = 1e-12",
      "execution_count": 9,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "def get_right_inverse_function_solve(mult_table, n_dims, gradeList):\n    '''\n    Returns a function that implements the right_inverse\n    '''\n\n    identity = np.zeros((n_dims,))\n    identity[gradeList.index(0)] = 1\n\n    tempAxes = (1,0,2)\n    newB = np.transpose(mult_table, tempAxes)\n\n    non_zero_indices = newB.nonzero()\n    k_list = non_zero_indices[0]\n    l_list = non_zero_indices[1]\n    m_list = non_zero_indices[2]\n    newB_vals = np.array([newB[k, l, m] for k, l, m in np.transpose(non_zero_indices)], dtype=int)\n\n    @numba.njit\n    def rightLaInv_func(value):\n        intermed = np.zeros((n_dims,n_dims))\n        for ind,l in enumerate(l_list):\n            m = m_list[ind]\n            k = k_list[ind]\n            intermed[k, m] += value[l]*newB_vals[ind]\n        sol = linalg.solve(intermed, identity)\n        return sol\n    return rightLaInv_func\n\n\ndef get_right_inverse_function_lstsq(mult_table, n_dims, gradeList):\n    '''\n    Returns a function that implements the right_inverse\n    '''\n\n    identity = np.zeros((n_dims,))\n    identity[gradeList.index(0)] = 1\n\n    tempAxes = (1,0,2)\n    newB = np.transpose(mult_table, tempAxes)\n\n    non_zero_indices = newB.nonzero()\n    k_list = non_zero_indices[0]\n    l_list = non_zero_indices[1]\n    m_list = non_zero_indices[2]\n    newB_vals = np.array([newB[k, l, m] for k, l, m in np.transpose(non_zero_indices)], dtype=int)\n\n    @numba.njit\n    def rightLaInv_func(value):\n        intermed = np.zeros((n_dims,n_dims))\n        for ind,l in enumerate(l_list):\n            m = m_list[ind]\n            k = k_list[ind]\n            intermed[k, m] += value[l]*newB_vals[ind]\n        sol = linalg.lstsq(intermed, identity)[0]\n        return sol\n    return rightLaInv_func\n\n\ndef rightLaInv(value):\n    \"\"\"Return right-inverse using a computational linear algebra method\n    proposed by Christian Perwass.\n     -1              -1\n    M    where M * M  == 1\n    rightLaInv() --> MultiVector\n    \"\"\"\n\n    identity = np.zeros((layout.gaDims,))\n    identity[layout.gradeList.index(0)] = 1\n\n    intermed = _myDot(value, layout.gmt)\n\n    if abs(linalg.det(intermed)) < _eps:\n        warn('Inverse might be inaccurate')\n        #raise ValueError(\"multivector has no right-inverse\")\n\n    sol = linalg.solve(intermed, identity)\n\n    return cf.MultiVector(layout,sol)\n",
      "execution_count": 10,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "rinv_solve = get_right_inverse_function_solve(layout.gmt,layout.gaDims,layout.gradeList)\nrinv_lstsq = get_right_inverse_function_lstsq(layout.gmt,layout.gaDims,layout.gradeList)",
      "execution_count": 11,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "value = layout.randomMV().value",
      "execution_count": 12,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "%%timeit\nrightLaInv(value)",
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": "208 µs ± 8.34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "%%timeit\nrinv_solve(value)",
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": "33 µs ± 1.15 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "%%timeit\nrinv_lstsq(value)",
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": "165 µs ± 909 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "a = rinv_solve(value)\nb = rightLaInv(value).value\nc = rinv_lstsq(value)\nnp.testing.assert_almost_equal(a,b)\nnp.testing.assert_almost_equal(a,c)",
      "execution_count": 19,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "mimetype": "text/x-python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Compare right inverse jit",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}