Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    primarykey
    data
    text
    <p>I just wanted to mention that if speed is required, wrapping unutbu's code in scipy's weave.inline and passing an already existing matrix as a parameter yields a 20-fold decrease in the running time.</p> <p>The code (in rotation_matrix_test.py):</p> <pre><code>import numpy as np import timeit from math import cos, sin, sqrt import numpy.random as nr from scipy import weave def rotation_matrix_weave(axis, theta, mat = None): if mat == None: mat = np.eye(3,3) support = "#include &lt;math.h&gt;" code = """ double x = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]); double a = cos(theta / 2.0); double b = -(axis[0] / x) * sin(theta / 2.0); double c = -(axis[1] / x) * sin(theta / 2.0); double d = -(axis[2] / x) * sin(theta / 2.0); mat[0] = a*a + b*b - c*c - d*d; mat[1] = 2 * (b*c - a*d); mat[2] = 2 * (b*d + a*c); mat[3*1 + 0] = 2*(b*c+a*d); mat[3*1 + 1] = a*a+c*c-b*b-d*d; mat[3*1 + 2] = 2*(c*d-a*b); mat[3*2 + 0] = 2*(b*d-a*c); mat[3*2 + 1] = 2*(c*d+a*b); mat[3*2 + 2] = a*a+d*d-b*b-c*c; """ weave.inline(code, ['axis', 'theta', 'mat'], support_code = support, libraries = ['m']) return mat def rotation_matrix_numpy(axis, theta): mat = np.eye(3,3) axis = axis/sqrt(np.dot(axis, axis)) a = cos(theta/2.) b, c, d = -axis*sin(theta/2.) return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) </code></pre> <p>The timing:</p> <pre><code>&gt;&gt;&gt; import timeit &gt;&gt;&gt; &gt;&gt;&gt; setup = """ ... import numpy as np ... import numpy.random as nr ... ... from rotation_matrix_test import rotation_matrix_weave ... from rotation_matrix_test import rotation_matrix_numpy ... ... mat1 = np.eye(3,3) ... theta = nr.random() ... axis = nr.random(3) ... """ &gt;&gt;&gt; &gt;&gt;&gt; timeit.repeat("rotation_matrix_weave(axis, theta, mat1)", setup=setup, number=100000) [0.36641597747802734, 0.34883809089660645, 0.3459300994873047] &gt;&gt;&gt; timeit.repeat("rotation_matrix_numpy(axis, theta)", setup=setup, number=100000) [7.180983066558838, 7.172032117843628, 7.180462837219238] </code></pre>
    singulars
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. VO
      singulars
      1. This table or related slice is empty.
    2. VO
      singulars
      1. This table or related slice is empty.
    3. VO
      singulars
      1. This table or related slice is empty.
    1. This table or related slice is empty.
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload