title: Writing a fast 3d matrix and vector library in cython
body: <p>As part of rewriting my game engine in cython, I am trying to improve the performance of my python+numpy classes for matrix and vector math as this was one of the major bottlenecks I had run into previously. This group of modules had defined classes for types such as <code>Vector2/3/4</code>, <code>Matrix2/3/4</code>, and <code>Quaternion</code>. </p>
  
  <p>Taking a page from the <a href="http://glmatrix.net/" rel="nofollow noreferrer">glMatrix javascript library</a>, I thought that one thing I could do this time around is switch away from a class-based system to a module with just a bunch of math functions to reduce more overhead. That way, instead of returning a new object every time I, say, added two vectors together, I would not have to construct a custom object.</p>
  
  <p>To test this out, I wrote a benchmark demo for creating two <code>Vec2</code> objects <code>a</code> and <code>b</code> summing them component-wise to get a <code>Vec2</code> object <code>out</code>. The code for this is broken down into a <code>main.py</code> that does the timing, a <code>vec2.pyx</code> for the cython code, and a <code>pyvec2.py</code> for the python code. Here is the code for each of those components:</p>
  
  <p><strong>main.py</strong></p>
  
  <pre><code>import time
  import array
  import math3d.vec2 as vec2
  import math3d.pyvec2 as pyvec2
  
  def test(n, func, param_list):
      start = time.time()
      for i in range(n):
          func(*param_list)
      end = time.time()
      print func, end-start
  
  test(1000000, pyvec2.pyadd, [[1, 2], [3, 4]])
  test(1000000, pyvec2.pyadd2, [[0, 0], [1, 2], [3, 4]])
  test(1000000, vec2.add, [[1, 2], [3, 4]])
  test(1000000, vec2.add2, [array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add3, [array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add4, [array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add5, [[0, 0], [1, 2], [3, 4]])
  test(1000000, vec2.add6, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add7, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add8, [array.array("f", [0, 0]), array.array("f", [1, 2]), array.array("f", [3, 4])])
  test(1000000, vec2.add9, [[0, 0], [1, 2], [3, 4]])
  </code></pre>
  
  <p><strong>vec2.pyx</strong></p>
  
  <pre><code>from libc.stdlib cimport malloc, free
  from cpython cimport array
  import array
  
  def add(list a, list b):
      cdef float[2] out = [0, 0]
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add2(float[:] a, float[:] b):
      cdef float[2] out = [0, 0]
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add3(array.array a, array.array b):
      cdef float[2] out = [0, 0]
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add4(array.array a, array.array b):
      cdef array.array out = array.array("f", [0, 0])
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add5(list out, list a, list b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add6(float[:] out, float[:] a, float[:] b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add7(array.array out, array.array a, array.array b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add8(array.array out, array.array a, array.array b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  
  def add9(out, a, b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  </code></pre>
  
  <p><strong>pyvec2.py</strong></p>
  
  <pre><code>def pyadd(a, b):
      out = [a[0] + b[0], a[1] + b[1]]
  
  def pyadd2(out, a, b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  </code></pre>
  
  <p>Here the results after running <code>main.py</code>:</p>
  
  <pre><code>&lt;function pyadd at 0x0000000003354828&gt; 0.380000114441
  &lt;function pyadd2 at 0x0000000003354908&gt; 0.31299996376
  &lt;built-in function add&gt; 0.261000156403
  &lt;built-in function add2&gt; 0.680999994278
  &lt;built-in function add3&gt; 0.268000125885
  &lt;built-in function add4&gt; 0.601000070572
  &lt;built-in function add5&gt; 0.144999980927
  &lt;built-in function add6&gt; 1.06299996376
  &lt;built-in function add7&gt; 0.241000175476
  &lt;built-in function add8&gt; 0.237999916077
  &lt;built-in function add9&gt; 0.141000032425
  </code></pre>
  
  <p>From this, it looks like just using python lists are faster than typed arrays! Furthermore, just blindly compiling my python function without typing into cython yielded the best results! It looks like a significant portion of the time spent executing the code is used to convert to and from python types.</p>
  
  <p>I am therefore curious to know if there are even faster ways to perform the math on the cython side while minimizing python overhead passing in arguments. I am not really interested in necessarily exposing lists or array.array objects for direct access to the contents of my vectors or matrices. It would be ideal to just pass a pointer to and from python into my cython math module, but this does not seem possible as pointers are not python objects. Any suggestions would be appreciated.</p>
  
  <hr>
  
  <p><strong>UPDATE:</strong></p>
  
  <p>Below is the code for my <code>Vec2</code> class. It is composed of two files: the first is a base <code>_Vec</code> class, and the second inherits from it as a specific <code>Vec2</code> class. </p>
  
  <p><strong>vec.py</strong></p>
  
  <pre><code>import math
  import numpy as np
  import random
  
  class _Vec(object):
  
      def __init__(self, *args):
          try:
              data, = args
              data = np.array(data, dtype=np.float32)
              cls_name = self.__class__.__name__
              vec2_check = cls_name == "Vec2" and len(data) != 2
              vec3_check = cls_name == "Vec3" and len(data) != 3
              vec4_check = cls_name == "Vec4" and len(data) != 4
              if any([vec2_check, vec3_check, vec4_check]) == True:
                  raise TypeError("{0} is not a valid {1}".format(data, cls_name))
          except ValueError:
              data = np.array(args, dtype=np.float32)
          self._data = data
  
      def __add__(self, other):
          if isinstance(other, self.__class__):
              return self.__class__(self._data + other._data)
          return self.__class__(self._data + other)
  
      def __radd__(self, other):
          return self.__class__(other + self._data)
  
      def __sub__(self, other):
          if isinstance(other, self.__class__):
              return self.__class__(self._data - other._data)
          return self.__class__(self._data - other)
  
      def __rsub__(self, other):
          return self.__class__(other - self._data)
  
      def __mul__(self, other):
          if isinstance(other, self.__class__):
              return self.__class__(self._data * other._data)
          return self.__class__(self._data * other)
  
      def __rmul__(self, other):
          return self.__class__(other * self._data)
  
      def __div__(self, other):
          if isinstance(other, self.__class__):
              return self.__class__(self._data / other._data)
          return self.__class__(self._data / other)
  
      def __rdiv__(self, other):
          return self.__class__(other / self._data)
  
      def __neg__(self):
          return self.__class__(-self._data)
  
      def __pos__(self):
          return self.__class__(+self._data)
  
      def __eq__(self, other):
          return np.array_equal(self._data, other._data)
  
      def __ne__(self, other):
          return not self.__eq__(other)
  
      def __lt__(self, other):
          return self.square_length() &lt; other.square_length()
  
      def __le__(self, other):
          return self.square_length() &lt;= other.square_length()
  
      def __gt__(self, other):
          return self.square_length() &gt; other.square_length()
  
      def __ge__(self, other):
          return self.square_length() &gt;= other.square_length()
  
      def __repr__(self):
          return "{0}(data={1})".format(self.__class__.__name__, self.get_data())
  
      def __str__(self):
          return np.array_str(self._data)
  
      def ceil(self):
          return self.__class__(np.ceil(self._data))
  
      def floor(self):
          return self.__class__(np.floor(self._data))
  
      def get_data(self):
          return self._data.flatten().tolist()
  
      def inverse(self):
          return self.__class__(1.0/self._data)
  
      def length(self):
          return float(np.linalg.norm(self._data))
  
      def negate(self):
          return self.__class__(-self._data)
  
      def normalize(self):
          length = self.length()
          if length == 0.0:
              return self.__class__(np.zeros(self._data.shape()))
          return self.__class__(self._data/length)
  
      def round(self, decimal=0):
          return self.__class__(np.round(self._data, decimal))
  
      def square_length(self):
          return float(np.sum(np.square(self._data)))
  
      @classmethod
      def distance(cls, a, b):
          c = b - a
          return c.length()
  
      @classmethod
      def dot(cls, a, b):
          return float(np.dot(a._data, b._data))
  
      @classmethod
      def equals(cls, a, b, tolerance=0.0):
          diffs = np.fabs((a - b)._data)
          pairs = zip(list(np.fabs(a._data)), list(np.fabs(b._data)))
          tolerance_calcs = [tolerance * max(1, a_val, b_val) for (a_val, b_val) in pairs]
          tests = [d &lt;= t for (d, t) in zip(diffs, tolerance_calcs)]
          return all(tests)
  
      @classmethod
      def lerp(cls, a, b, t):
          return a*(1-t) + b*t
  
      @classmethod
      def max_components(cls, a, b):
          return cls(np.maximum(a._data, b._data))
  
      @classmethod
      def min_components(cls, a, b):
          return cls(np.minimum(a._data, b._data))
  
      @classmethod
      def random(cls, n):
          return cls(np.random.rand((n)))
  
      @classmethod
      def square_distance(cls, a, b):
          c = b - a
          return c.square_length()
  </code></pre>
  
  <p><code>vec2.py</code></p>
  
  <pre><code>from vec import _Vec
  from vec3 import Vec3
  from vec4 import Vec4
  import math
  import numpy as np
  import random
  
  class Vec2(_Vec):
  
      @property
      def x(self):
          return float(self._data[0])
  
      @x.setter
      def x(self, value):
          self._data[0] = float(value)
  
      @property
      def y(self):
          return float(self._data[1])
  
      @y.setter
      def y(self, value):
          self._data[1] = float(value)
  
      def __repr__(self):
          return "Vec2(x={0}, y={1})".format(self.x, self.y)
  
      def transform_mat2(self, a):
          prod = np.dot(a._data.T, self._data.T).T
          return Vec2(prod)
  
      def transform_mat3(self, a):
          v3 = Vec3(self.get_data() + [1])
          prod = np.dot(a._data.T, v3._data.T).T
          return Vec2(prod[0:2])
  
      def transform_mat4(self, a):
          v4 = Vec4(self.get_data() + [0, 1])
          prod = np.dot(a._data.T, v4._data.T).T
          return Vec2(prod[0:2])
  
      @classmethod
      def random(cls):
          return super(Vec2, cls).random(2)
  </code></pre>
  
  <hr>
  
  <p><strong>UPDATE 2:</strong></p>
  
  <p>As kazemakase pointed out, some of the values in <code>main.py</code> were integers. Defining everything as floating point by appending a .0, I got the following times instead:</p>
  
  <pre><code>&lt;function pyadd at 0x0000000002FF4828&gt; 0.384000062943
  &lt;function pyadd2 at 0x0000000002FF4908&gt; 0.332000017166
  &lt;built-in function add&gt; 0.227999925613
  &lt;built-in function add2&gt; 0.640000104904
  &lt;built-in function add3&gt; 0.258999824524
  &lt;built-in function add4&gt; 0.556999921799
  &lt;built-in function add5&gt; 0.145999908447
  &lt;built-in function add6&gt; 0.983999967575
  &lt;built-in function add7&gt; 0.217000007629
  &lt;built-in function add8&gt; 0.236000061035
  &lt;built-in function add9&gt; 0.131000041962
  </code></pre>
  
  <p>These appear to be similar to the original times in that case 5 and 9 are faster.</p>
  
  <hr>
  
  <p><strong>UPDATE 3:</strong></p>
  
  <p>As BrenBarn pointed out, it would probably be useful to explain how I used my original python+numpy class and why I am struggling with its performance. Originally, my entire 3d game library project was in pure python, using PyOpenGL for rendering the graphics. In order to position meshes/models in my 3d world, every frame, I need to calculate a Mat4 transformation matrix that defines that object's position, rotation, and scale in the world that is then uploaded to the GPU. When I was positioning many objects (> 1000), the frame rate in my app would slow down. While temporarily disabling 3d did improve performance somewhat, my app was still lagging below 60 fps. That is when I realized that simply calculating the <code>Mat4</code> matrices in python was to blame. When I removed these calculations and just drew all of my 3d objects at the origin, performance went back up.</p>
  
  <p>Wondering then if ALL of my math3d library was slow, I thought I would start benchmarking all of these classes. I did the following benchmark for my <code>Vec2</code> class's add function as it was the simplest and computationally cheapest function to test:</p>
  
  <pre><code>import time
  from vec2 import Vec2
  
  def add(a, b):
      out = a + b
      return out
  
  def test(n, func, param_list):
      start = time.time()
      for i in range(n):
          func(*param_list)
      end = time.time()
      print func, end-start
  
  test(1000000, add, [Vec2([1.0, 2.0]), Vec2([3.0, 4.0])])
  #&lt;function add at 0x00000000041B1CF8&gt; 2.81699991226 (using the real def add(a, b) function)
  #&lt;function add at 0x000000000362CCF8&gt; 0.168999910355 (just passing in values to def add(a, b): pass)
  </code></pre>
  
  <p>Those comments showed the times I had got just for adding two <code>Vec2</code> together. From that, I concluded that the math was slow, that there was significant python function call + class overhead, and that this was a performance bottleneck to consider in my 3d library. I hope this provides some rationale for this question.</p>
  
  <hr>
  
  <p><strong>Update 4:</strong></p>
  
  <p>What Paul Cornelius said about incurring the penalty of translation from python to cython only once got me thinking: why not just fetch a "pointer" to a <code>float *</code> object whenever I create a "Vec2 object"? Then, I can just pass in pointers for future math operations, cython can dereference those pointers to get the actual data, and the math operation can be performed. The result of that is the following code:</p>
  
  <p><strong>main.py</strong></p>
  
  <pre><code>import time
  import array
  import math3d.vec2 as vec2
  
  def make_list(a, b):
      out = [a, b]
  
  def test(n, func, param_list):
      start = time.time()
      for i in range(n):
          func(*param_list)
      end = time.time()
      print func, end-start
  
  test(1000000, vec2.create, [1, 2])
  test(1000000, make_list, [1, 2])
  a = vec2.create(1, 2)
  #b = vec2.get_data(a)
  b = vec2.create(3, 4)
  c = vec2.create(0, 0)
  test(1000000, vec2.add, [c, a, b])
  test(1000000, vec2.add2, [[0, 0], [1, 2], [3, 4]])
  print vec2.get_data(c)
  </code></pre>
  
  <p><strong>vec2.pyx</strong></p>
  
  <pre><code>def add(uintptr_t out, uintptr_t a, uintptr_t b):
      cdef float *a_data = &lt;float *&gt;a
      cdef float *b_data = &lt;float *&gt;b
      cdef float *out_data = &lt;float *&gt;out
      out_data[0] = a_data[0] + b_data[0]
      out_data[1] = a_data[1] + b_data[1]
  
  def create(float x, float y):
      cdef float* a = &lt;float *&gt;malloc(sizeof(float))
      a[:] = [x, y]
      cdef uintptr_t u_ptr = &lt;uintptr_t&gt; a
      return u_ptr
  
  def get_data(uintptr_t u_ptr):
      cdef float *b = &lt;float *&gt;u_ptr
      return b[0], b[1]
  
  def add2(out, a, b):
      out[0] = a[0] + b[0]
      out[1] = a[1] + b[1]
      return out
  </code></pre>
  
  <p><strong>Timing Results</strong></p>
  
  <pre><code>&lt;built-in function create&gt; 0.19000005722
  &lt;function make_list at 0x0000000002994828&gt; 0.269999980927
  &lt;built-in function add&gt; 0.111000061035
  &lt;built-in function add2&gt; 0.141999959946
  (4.0, 6.0)
  </code></pre>
  
  <p>Of course, handling <code>uintptr_t</code> integers in python is pretty unsafe, as they could mistakenly be added together on the python side with the <code>+</code> operator. Also, it is unclear if the slight performance benefit of this (0.03 seconds for 1 million operations) is really worth it.</p>