[Numpy-discussion] Faster way to generate a rotation matrix?

[Hoyt Koepke]

Hello,

> def rotation(theta, R = np.zeros((3,3))):
>    cx,cy,cz = np.cos(theta)
>    sx,sy,sz = np.sin(theta)
>    R.flat = (cx*cz - sx*cy*sz, cx*sz + sx*cy*cz, sx*sy,
>        -sx*cz - cx*cy*sz, -sx*sz + cx*cy*cz,
>        cx*sy, sy*sz, -sy*cz, cy)
>    return R
>
> Pretty evil looking ;) but still wouldn't mind somehow getting it faster

I would definitely encourage you to check out cython.  I have to write
lots of numerically intensive stuff in my python code, and I tend to
cythonize it a lot.  In cython, the above would be (something like):

from numpy cimport ndarray

cdef extern from "math.h":
   double cos(double)
   double sin(double)

def rotation(ndarry[double] theta, ndarray[double, ndim=2] R = np.zeros((3,3))):
     cdef double cx = cos(theta[0]), cy = cos(theta[1]), cz = cos(theta[2])
     cdef double sx = sin(theta[0]), sy = sin(theta[1]), sz = sin(theta[2])

    R[0,0] = cx*cz - sx*cy*sz
    R[0,1] = cx*sz + sx*cy*cz
    R[0,2] = sx*sy
    ...
    R[2,2] = cy

    return R

And that will be probably be orders of magnitude faster than what you
currently have, as everything but the function call and the return
statement would become C code.  Compilers these days are very good at
optimizing that kind of thing too.

--Hoyt

++++++++++++++++++++++++++++++++++++++++++++++++
+ Hoyt Koepke
+ University of Washington Department of Statistics
+ http://www.stat.washington.edu/~hoytak/
+ hoytak at gmail.com
++++++++++++++++++++++++++++++++++++++++++