Sparse/numba performance of radon transform

sparse_radon_performance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from clifford import g3c\n",
    "from numpy import pi, e\n",
    "import numpy as np\n",
    "import math\n",
    "import random\n",
    "\n",
    "layout = g3c.layout\n",
    "locals().update(g3c.blades)\n",
    "\n",
    "ep, en, up, down, homo, E0, ninf, no = (g3c.stuff[\"ep\"], g3c.stuff[\"en\"], \n",
    "                                        g3c.stuff[\"up\"], g3c.stuff[\"down\"], g3c.stuff[\"homo\"], \n",
    "                                        g3c.stuff[\"E0\"], g3c.stuff[\"einf\"], -g3c.stuff[\"eo\"])\n",
    "I3 = e123\n",
    "I5 = e12345\n",
    "\n",
    "def MeetSpheres(sphereA,sphereB):\n",
    "    return (sphereA*sphereB)(2)*I5\n",
    "\n",
    "def generate_sphere_from_dual(center,radius):\n",
    "    dual_sphere = center  - 0.5*radius*radius*ninf\n",
    "    sphere = I5*dual_sphere\n",
    "    return sphere\n",
    "\n",
    "def get_center_from_sphere(sphere):\n",
    "    center = sphere*ninf*sphere\n",
    "    return center\n",
    "\n",
    "def get_radius_from_sphere(sphere):\n",
    "    dual_sphere = sphere*I5\n",
    "    return math.sqrt( abs(dual_sphere*dual_sphere) )\n",
    "\n",
    "def intersect_spheres_and_plane(sphere_list, plane):\n",
    "    intersection_count = 0\n",
    "    for sphere in sphere_list:\n",
    "        intersection_circle = MeetSpheres(sphere,plane)\n",
    "        sqr_val = float( (intersection_circle*intersection_circle)[0] )\n",
    "        if sqr_val >=0:\n",
    "            intersection_count += 1\n",
    "    return intersection_count\n",
    "\n",
    "def remove_floating_point_zero(geo_obj, threshold):\n",
    "    a = 0;\n",
    "    for i in range(5):\n",
    "        if abs( geo_obj[i] ) > threshold:\n",
    "            a = a + geo_obj(i)\n",
    "    return a\n",
    "\n",
    "def generate_plane_from_dual(euclidean_normal, rho):\n",
    "    dual_of_plane = euclidean_normal + rho*ninf\n",
    "    plane = normalise_plane(-I5*dual_of_plane)\n",
    "    return plane\n",
    "\n",
    "def get_plane_origin_distance(plane):\n",
    "    # The clifford package defines no as 0.5*the no defined in \"covariant approach to geometry\"\n",
    "    return (plane*I5)|no\n",
    "\n",
    "def get_plane_normal(plane):\n",
    "    # The clifford package defines no as 0.5*the no defined in \"covariant approach to geometry\"\n",
    "    return (plane*I5 - get_plane_origin_distance(plane)*ninf)\n",
    "\n",
    "def translate_point_along_vector(point_p,vector_a):\n",
    "    return (1 + ninf*vector_a/2)*point_p*(1 + vector_a*ninf/2)\n",
    "\n",
    "def sample_sphere(npoints, radius, ndim=3):\n",
    "    vec = np.random.randn(npoints,ndim)\n",
    "    conformal_vec = []\n",
    "    for v in vec:\n",
    "        resized_v = radius*v/np.linalg.norm(v)\n",
    "        a = up(resized_v[0]*e1 + resized_v[1]*e2 + resized_v[2]*e3)\n",
    "        conformal_vec.append(a)\n",
    "    return conformal_vec\n",
    "\n",
    "def projection_point_in_plane(point,plane):\n",
    "    # https://crypto.stanford.edu/~blynn/haskell/hga.html\n",
    "    return up(down( (point|plane)*(plane.inv()) ))\n",
    "\n",
    "def rejection_point_in_plane(point,plane):\n",
    "    # https://crypto.stanford.edu/~blynn/haskell/hga.html\n",
    "    return (point^plane)*(plane.inv())\n",
    "\n",
    "def draw_points_from_plane(N_points, plane, radius):\n",
    "    sphere_points = sample_sphere(N_points,radius)\n",
    "    random_points = [ projection_point_in_plane(p,plane) for p in sphere_points ]\n",
    "    return random_points\n",
    "\n",
    "def normalise_plane(plane):\n",
    "    new_plane = plane/abs(plane)\n",
    "    return new_plane\n",
    "\n",
    "def add_sphere_to_axis(sphere, ax):\n",
    "    centre = down( get_center_from_sphere(sphere) )\n",
    "    radius = get_radius_from_sphere(sphere)\n",
    "    u, v = np.mgrid[0:2*np.pi:10j, 0:np.pi:5j]\n",
    "    x = centre[1] + radius*np.cos(u)*np.sin(v)\n",
    "    y = centre[2] + radius*np.sin(u)*np.sin(v)\n",
    "    z = centre[3] + radius*np.cos(v)\n",
    "    ax.plot_wireframe(x, y, z, color=\"g\", alpha=0.4)\n",
    "\n",
    "def add_plane_to_axis(plane, ax):\n",
    "    # The plane\n",
    "    normal = get_plane_normal(plane)\n",
    "    d = float( get_plane_origin_distance(plane) )\n",
    "\n",
    "    xx, yy = np.meshgrid(range(-15,15), range(-15,15))\n",
    "    z = (-normal[1] * xx -normal[2] * yy + d) * 1. /normal[3]\n",
    "    ax.plot_surface(xx, yy, z, alpha=0.2)\n",
    "\n",
    "def add_points_to_axis(points, ax):\n",
    "        # Now the plane_basis\n",
    "    euclidean_points = [down(p) for p in points]\n",
    "    xs =[p[1] for p in euclidean_points]\n",
    "    ys = [p[2] for p in euclidean_points] \n",
    "    zs =[p[3] for p in euclidean_points] \n",
    "    ax.scatter(xs, ys, zs, c=[(1,0,0),(1,0,0),(1,0,0)], alpha=1)\n",
    "\n",
    "def plot_points_and_plane(spheres, plane, plane_basis):\n",
    "    # Set up the axis\n",
    "    import matplotlib.pyplot as plt\n",
    "    from mpl_toolkits.mplot3d import Axes3D\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    ax.set_aspect(\"equal\")\n",
    "\n",
    "    # Draw the stuff\n",
    "    add_plane_to_axis(plane,ax)\n",
    "    for sphere in spheres:\n",
    "        add_sphere_to_axis(sphere,ax)\n",
    "    add_points_to_axis(plane_basis, ax)\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "def cga_to_euclidean_np_array(p):\n",
    "    array_point = np.float64([down(p)[1],down(p)[2],down(p)[3]])\n",
    "    return array_point\n",
    "\n",
    "def np_array_euclidean_to_cga(p):\n",
    "    cga_point = up( p[0]*e1 + p[1]*e2 + p[2]*e3 )\n",
    "    return cga_point\n",
    "\n",
    "def spherical_to_unit_vector(theta,psi):\n",
    "    x = math.sin(theta)*math.cos(psi)\n",
    "    y = math.sin(theta)*math.sin(psi)\n",
    "    z = math.cos(theta)\n",
    "    return x*e1 + y*e2 + z*e3\n",
    "\n",
    "def np_array_to_spherical(euc_point):\n",
    "    mag = np.linalg.norm(euc_point)\n",
    "    theta = math.acos(euc_point[2]/mag)\n",
    "    psi = math.atan2(euc_point[1],euc_point[0])\n",
    "    return [mag,theta,psi]\n",
    "\n",
    "\n",
    "def planar_radon_transform(rho,theta,psi,point_list,beam_radius):\n",
    "\n",
    "    euclidean_normal = spherical_to_unit_vector(theta,psi)\n",
    "\n",
    "    # Our query plane!\n",
    "    plane = generate_plane_from_dual(euclidean_normal, rho)\n",
    "\n",
    "    # Generate a sphere around each point \n",
    "    sphere_list = [generate_sphere_from_dual(center, beam_radius) for center in point_list]\n",
    "\n",
    "    # Show them on a figure\n",
    "    #plot_points_and_plane(sphere_list,plane,plane_base_points)\n",
    "\n",
    "    # Do the intersection\n",
    "    intersection_count = intersect_spheres_and_plane(sphere_list, plane)\n",
    "    #print('Intersection count ', intersection_count)\n",
    "    return intersection_count\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " "
     ]
    }
   ],
   "source": [
    "%%prun -s cumulative\n",
    "\n",
    "# Seed\n",
    "random.seed(1)\n",
    "\n",
    "# Define a plane to generate spheres at\n",
    "plane_base_points = [up(5.0*e1),up(5.0*e2),up(5.0*e3)]\n",
    "plane =  normalise_plane(plane_base_points[0]^plane_base_points[1]^plane_base_points[2]^ninf)\n",
    "\n",
    "# Where we expect to find the plane in our radon space\n",
    "rho = get_plane_origin_distance(plane)\n",
    "euclidean_normal = get_plane_normal(plane)\n",
    "\n",
    "spherical_coords = np_array_to_spherical( np.float64( [euclidean_normal[1],euclidean_normal[2],euclidean_normal[3] ]) )\n",
    "target_spherical = spherical_coords[:]\n",
    "\n",
    "# Draw points from a plane\n",
    "N_points = 255\n",
    "point_list = draw_points_from_plane(N_points, plane, 10)\n",
    "\n",
    "# Check this position of plane\n",
    "\n",
    "#print(spherical_to_unit_vector(phi,psi))\n",
    "intersection_rows = []\n",
    "for theta in np.arange(0,np.pi,np.pi/20.0):\n",
    "    intersection_column = []\n",
    "    for psi in np.arange(0,2*np.pi,np.pi/20.0):\n",
    "        intersection_column.append( planar_radon_transform(rho,theta,psi,point_list,1) )\n",
    "    intersection_rows.append(intersection_column[:])\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "24093363 function calls (24039567 primitive calls) in 19.056 seconds\n",
    "\n",
    "   Ordered by: cumulative time\n",
    "\n",
    "   ncalls  tottime  percall  cumtime  percall filename:lineno(function)\n",
    "        1    0.000    0.000   19.056   19.056 {built-in method builtins.exec}\n",
    "        1    0.036    0.036   19.056   19.056 <string>:4(<module>)\n",
    "      800    0.002    0.000   18.086    0.023 <ipython-input-1-0bd74762afc1>:153(planar_radon_transform)\n",
    "   821216    2.156    0.000   11.805    0.000 __init__.py:740(__mul__)\n",
    "      800    0.431    0.001   11.670    0.015 <ipython-input-1-0bd74762afc1>:32(intersect_spheres_and_plane)\n",
    "   204000    0.458    0.000    8.091    0.000 <ipython-input-1-0bd74762afc1>:16(MeetSpheres)\n",
    "  1237895    0.894    0.000    7.935    0.000 __init__.py:706(_checkOther)\n",
    "      800    0.089    0.000    6.296    0.008 <ipython-input-1-0bd74762afc1>:161(<listcomp>)\n",
    "   204000    0.431    0.000    6.207    0.000 <ipython-input-1-0bd74762afc1>:19(generate_sphere_from_dual)\n",
    "  1029158    0.346    0.000    5.035    0.000 __init__.py:364(__ne__)\n",
    "  1029158    1.077    0.000    4.688    0.000 __init__.py:361(__eq__)\n",
    "  1029158    0.822    0.000    3.611    0.000 fromnumeric.py:1973(all)\n",
    "  1029158    0.352    0.000    2.258    0.000 {method 'all' of 'numpy.ndarray' objects}\n",
    "   821221    2.198    0.000    2.198    0.000 __init__.py:206(gmt_mult)\n",
    "  2792159    0.586    0.000    2.068    0.000 {built-in method builtins.isinstance}\n",
    "   204001    0.253    0.000    2.037    0.000 __init__.py:836(__sub__)\n",
    "  1444774    0.615    0.000    2.019    0.000 __init__.py:729(_newMV)\n",
    "  1029158    0.247    0.000    1.906    0.000 _methods.py:40(_all)\n",
    "   204000    1.377    0.000    1.703    0.000 __init__.py:1158(__call__)\n",
    "  1029302    1.660    0.000    1.660    0.000 {method 'reduce' of 'numpy.ufunc' objects}\n",
    "  1242899    0.896    0.000    1.482    0.000 abc.py:178(__instancecheck__)\n",
    "  1444774    0.807    0.000    1.404    0.000 __init__.py:666(__init__)\n",
    "   206403    0.467    0.000    1.034    0.000 __init__.py:756(__rmul__)\n",
    "  2477128    0.692    0.000    0.692    0.000 {built-in method numpy.core.multiarray.array}\n",
    "        6    0.000    0.000    0.622    0.104 dispatcher.py:294(_compile_for_args)\n",
    "        6    0.000    0.000    0.621    0.104 dispatcher.py:554(compile)\n",
    "        6    0.000    0.000    0.621    0.103 dispatcher.py:70(compile)\n",
    "        6    0.000    0.000    0.620    0.103 compiler.py:753(compile_extra)\n",
    "        6    0.000    0.000    0.617    0.103 compiler.py:345(compile_extra)\n",
    "        6    0.000    0.000    0.616    0.103 compiler.py:717(_compile_bytecode)\n",
    "        6    0.000    0.000    0.616    0.103 compiler.py:668(_compile_core)\n",
    "        6    0.000    0.000    0.615    0.103 compiler.py:229(run)\n",
    "  2276586    0.587    0.000    0.587    0.000 _weakrefset.py:70(__contains__)\n",
    "  1029206    0.440    0.000    0.531    0.000 numeric.py:534(asanyarray)\n",
    "        6    0.000    0.000    0.510    0.085 compiler.py:639(stage_nopython_backend)\n",
    "        6    0.000    0.000    0.510    0.085 compiler.py:588(_backend)"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

24093363 function calls (24039567 primitive calls) in 19.056 seconds

Ordered by: cumulative time

ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 19.056 19.056 {built-in method builtins.exec}
1 0.036 0.036 19.056 19.056 :4()
800 0.002 0.000 18.086 0.023 :153(planar_radon_transform)
821216 2.156 0.000 11.805 0.000 init.py:740(mul)
800 0.431 0.001 11.670 0.015 :32(intersect_spheres_and_plane)
204000 0.458 0.000 8.091 0.000 :16(MeetSpheres)
1237895 0.894 0.000 7.935 0.000 init.py:706(_checkOther)
800 0.089 0.000 6.296 0.008 :161()
204000 0.431 0.000 6.207 0.000 :19(generate_sphere_from_dual)
1029158 0.346 0.000 5.035 0.000 init.py:364(ne)
1029158 1.077 0.000 4.688 0.000 init.py:361(eq)
1029158 0.822 0.000 3.611 0.000 fromnumeric.py:1973(all)
1029158 0.352 0.000 2.258 0.000 {method 'all' of 'numpy.ndarray' objects}
821221 2.198 0.000 2.198 0.000 init.py:206(gmt_mult)
2792159 0.586 0.000 2.068 0.000 {built-in method builtins.isinstance}
204001 0.253 0.000 2.037 0.000 init.py:836(sub)
1444774 0.615 0.000 2.019 0.000 init.py:729(_newMV)
1029158 0.247 0.000 1.906 0.000 _methods.py:40(_all)
204000 1.377 0.000 1.703 0.000 init.py:1158(call)
1029302 1.660 0.000 1.660 0.000 {method 'reduce' of 'numpy.ufunc' objects}
1242899 0.896 0.000 1.482 0.000 abc.py:178(instancecheck)
1444774 0.807 0.000 1.404 0.000 init.py:666(init)
206403 0.467 0.000 1.034 0.000 init.py:756(rmul)
2477128 0.692 0.000 0.692 0.000 {built-in method numpy.core.multiarray.array}
6 0.000 0.000 0.622 0.104 dispatcher.py:294(_compile_for_args)
6 0.000 0.000 0.621 0.104 dispatcher.py:554(compile)
6 0.000 0.000 0.621 0.103 dispatcher.py:70(compile)
6 0.000 0.000 0.620 0.103 compiler.py:753(compile_extra)
6 0.000 0.000 0.617 0.103 compiler.py:345(compile_extra)
6 0.000 0.000 0.616 0.103 compiler.py:717(_compile_bytecode)
6 0.000 0.000 0.616 0.103 compiler.py:668(_compile_core)
6 0.000 0.000 0.615 0.103 compiler.py:229(run)
2276586 0.587 0.000 0.587 0.000 _weakrefset.py:70(contains)
1029206 0.440 0.000 0.531 0.000 numeric.py:534(asanyarray)
6 0.000 0.000 0.510 0.085 compiler.py:639(stage_nopython_backend)
6 0.000 0.000 0.510 0.085 compiler.py:588(_backend)