Skip to content

Instantly share code, notes, and snippets.

@hugohadfield

hugohadfield/SparseMultiplicationTest.ipynb

Created May 29, 2018 15:47
Show Gist options
  • Save hugohadfield/f847a19aaf654f39b75f94da9b4e6fb2 to your computer and use it in GitHub Desktop.
Save hugohadfield/f847a19aaf654f39b75f94da9b4e6fb2 to your computer and use it in GitHub Desktop.
Download ZIP
SparseMultiplicationTest.ipynb
Display the source blob
Display the rendered blob
Raw
SparseMultiplicationTest.ipynb
{
"cells": [
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import numpy as np\nimport numba\nimport clifford as cf\nfrom clifford.g3c import *",
"execution_count": 1,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def get_mult_function_dense(mult_table,n_dims):\n ''' \n Returns a function that implements the mult_table on two input multivectors\n '''\n non_zero_indices = mult_table.nonzero()\n k_list = non_zero_indices[0]\n l_list = non_zero_indices[1]\n m_list = non_zero_indices[2]\n mult_table_vals = np.array([mult_table[k,l,m] for k,l,m in np.transpose(non_zero_indices)],dtype=int)\n\n @numba.njit\n def mv_mult(value,other_value):\n output = np.zeros(n_dims)\n for ind,k in enumerate(k_list):\n l = l_list[ind]\n m = m_list[ind]\n output[l] += value[k]*mult_table_vals[ind]*other_value[m]\n return output\n return mv_mult\n",
"execution_count": 2,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def get_mult_function_graded(mult_table,n_dims,gradeList,grades_a=None,grades_b=None,filter_mask=None):\n ''' \n Returns a function that implements the mult_table on two input multivectors\n '''\n non_zero_indices = mult_table.nonzero()\n k_list = non_zero_indices[0]\n l_list = non_zero_indices[1]\n m_list = non_zero_indices[2]\n mult_table_vals = np.array([mult_table[k,l,m] for k,l,m in np.transpose(non_zero_indices)],dtype=int)\n \n if filter_mask is not None:\n # We can pass the sparse filter mask directly\n k_list = k_list[filter_mask]\n l_list = l_list[filter_mask]\n m_list = m_list[filter_mask]\n mult_table_vals = mult_table_vals[filter_mask]\n \n elif ((grades_a is not None) and (grades_b is not None)):\n # We can also specify sparseness by grade\n filter_mask = np.zeros(len(k_list), dtype=bool)\n for i in range(len(filter_mask)):\n if gradeList[k_list[i]] in grades_a:\n if gradeList[m_list[i]] in grades_b:\n filter_mask[i] = 1\n\n k_list = k_list[filter_mask]\n l_list = l_list[filter_mask]\n m_list = m_list[filter_mask]\n mult_table_vals = mult_table_vals[filter_mask]\n \n @numba.njit\n def mv_mult(value,other_value):\n output = np.zeros(n_dims)\n for ind,k in enumerate(k_list):\n m = m_list[ind]\n l = l_list[ind]\n output[l] += value[k]*mult_table_vals[ind]*other_value[m]\n return output\n return mv_mult",
"execution_count": 3,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "def get_mult_function_full(mult_table,n_dims,gradeList,grades_a=None,grades_b=None,filter_mask=None):\n ''' \n Returns a function that implements the mult_table on two input multivectors\n '''\n non_zero_indices = mult_table.nonzero()\n k_list = non_zero_indices[0]\n l_list = non_zero_indices[1]\n m_list = non_zero_indices[2]\n mult_table_vals = np.array([mult_table[k,l,m] for k,l,m in np.transpose(non_zero_indices)],dtype=int)\n \n if filter_mask is not None:\n # We can pass the sparse filter mask directly\n k_list = k_list[filter_mask]\n l_list = l_list[filter_mask]\n m_list = m_list[filter_mask]\n mult_table_vals = mult_table_vals[filter_mask]\n @numba.njit\n def mv_mult(value,other_value):\n output = np.zeros(n_dims)\n for ind,k in enumerate(k_list):\n m = m_list[ind]\n l = l_list[ind]\n output[l] += value[k]*mult_table_vals[ind]*other_value[m]\n return output\n return mv_mult\n \n elif ((grades_a is not None) and (grades_b is not None)):\n # We can also specify sparseness by grade\n filter_mask = np.zeros(len(k_list), dtype=bool)\n for i in range(len(filter_mask)):\n if gradeList[k_list[i]] in grades_a:\n if gradeList[m_list[i]] in grades_b:\n filter_mask[i] = 1\n\n k_list = k_list[filter_mask]\n l_list = l_list[filter_mask]\n m_list = m_list[filter_mask]\n mult_table_vals = mult_table_vals[filter_mask]\n \n @numba.njit\n def mv_mult(value,other_value):\n output = np.zeros(n_dims)\n for ind,k in enumerate(k_list):\n m = m_list[ind]\n l = l_list[ind]\n output[l] += value[k]*mult_table_vals[ind]*other_value[m]\n return output\n return mv_mult\n\n # This case we specify no sparseness in advance, the algorithm checks for zeros\n @numba.njit\n def mv_mult(value,other_value):\n output = np.zeros(n_dims)\n for ind,k in enumerate(k_list):\n v_val = value[k]\n if v_val!=0.0:\n m = m_list[ind]\n ov_val = other_value[m]\n if ov_val!=0.0:\n l = l_list[ind]\n output[l] += v_val*mult_table_vals[ind]*ov_val\n return output\n return mv_mult",
"execution_count": 4,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "# Lets make a sphere\na = (up(-e1)^up(e1)^up(e2)^up(e3)).normal()\n# And a line\nb = (up(e1)^up(-e1)^einf).normal()",
"execution_count": 5,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "d_gp = get_mult_function_dense(layout.gmt,layout.gaDims)\n\ngraded_dual_gp = get_mult_function_graded(layout.gmt,layout.gaDims,layout.gradeList,grades_a=[5],grades_b=[0,1,2,3,4,5])\ngraded_gp = get_mult_function_graded(layout.gmt,layout.gaDims,layout.gradeList) # This will be the same as d_gp\ngraded_gp_34 = get_mult_function_graded(layout.gmt,layout.gaDims,layout.gradeList,grades_a=[4],grades_b=[3])\n\nfull_gp = get_mult_function_full(layout.gmt,layout.gaDims,layout.gradeList)",
"execution_count": 6,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "# Now lets see how they perform"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nc = a*b",
"execution_count": 7,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "4.67 µs ± 11.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nc_val = d_gp(a.value,b.value)",
"execution_count": 8,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "2.56 µs ± 9.06 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nc_val = graded_gp(a.value,b.value)",
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "2.58 µs ± 10.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nc_val = graded_gp_34(a.value,b.value)",
"execution_count": 10,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "715 ns ± 1.64 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nc_val = full_gp(a.value,b.value)",
"execution_count": 11,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "1.5 µs ± 2.48 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "As expected when we know the sparsity of our objects in advance we can generate and compile functions that are far more efficient for multiplying them than the standard dense operation. This however requires us to know the grade of all our objects in advance. Our full multiplication operation is a trade off. We can still exploit the sparsity of the input multivectors by checking for zero values but at the expense of iterating over the values anyway and checking them against 0."
},
{
"metadata": {},
"cell_type": "markdown",
"source": "# How about on something that is a bit more realistic"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "I5 = e12345\nI5_val = e12345.value\n\ndef no_optimisation_func(A, B):\n return I5*((I5*A)*(I5*B))\n\n@numba.njit\ndef dense_func_val(A, B):\n return d_gp(I5_val,d_gp(d_gp(I5_val,A) , d_gp(I5_val,B)))\n\ndef dense_func(A, B):\n return cf.MultiVector(layout, dense_func_val(A.value, B.value))\n\n@numba.njit\ndef graded_func_val(A, B):\n return graded_dual_gp(I5_val,graded_gp(graded_dual_gp(I5_val,A), graded_dual_gp(I5_val,B)))\n\ndef graded_func(A, B):\n return cf.MultiVector(layout, graded_func_val(A.value, B.value))\n\n@numba.njit\ndef full_func_val(A, B):\n return full_gp(I5_val,full_gp(full_gp(I5_val,A) , full_gp(I5_val,B)))\n\ndef full_func(A, B):\n return cf.MultiVector(layout, full_func_val(A.value, B.value))",
"execution_count": 12,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "# First check that they are all equivalent"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "np.testing.assert_almost_equal(dense_func(a,b).value, graded_func(a,b).value)\nnp.testing.assert_almost_equal(dense_func(a,b).value, full_func(a,b).value)",
"execution_count": 13,
"outputs": []
},
{
"metadata": {},
"cell_type": "markdown",
"source": "# Now look at their speed"
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nno_optimisation_func(a,b)",
"execution_count": 14,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "18.6 µs ± 104 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\ndense_func(a,b)",
"execution_count": 15,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "9.74 µs ± 206 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\ngraded_func(a,b)",
"execution_count": 16,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "4.19 µs ± 53 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "%%timeit\nfull_func(a,b)",
"execution_count": 17,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": "5.11 µs ± 37.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
}
]
},
{
"metadata": {},
"cell_type": "markdown",
"source": "In this case the we are still much faster than the dense function but without having to know the grades of our objects in advance."
}
],
"metadata": {
"language_info": {
"nbconvert_exporter": "python",
"mimetype": "text/x-python",
"pygments_lexer": "ipython3",
"file_extension": ".py",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"name": "python",
"version": "3.5.2"
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3",
"language": "python"
},
"gist": {
"id": "",
"data": {
"description": "SparseMultiplicationTest.ipynb",
"public": true
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment