.. _rotating: ============================================================= Rotating particles and Python efficiency ============================================================= All the codes on this page can be found in the directory $CLASSHG/codes/particles/ (see :ref:classhg for information on obtaining the class repository). For Fortran versions of these codes, see :ref:rotating_fortran. For more about the timings reported here and the computer used, see :ref:timings. The main point of the example on this page is to show how implementing a Python code using a few basic principles can speed it up by a factor of 2000 relative to a first attempt. These principles apply equally well to Matlab codes (see :ref:numerical_python for more about Python vs. Matlab for numerical calculations). For a general introduction to Python, see :ref:python. A poorly written code --------------------- These examples solve the simple problem described in Section :ref:particle_methods. Below is a Python program that computes the rotating particles and plots them at a number of revolution angles. Note that this program is **not** written at all efficiently, in fact it was designed to be nearly as bad as possible for what it does (though it could be made even worse, see :ref:python_memory_alloc). .. literalinclude:: ../codes/particles/rotate1.py :language: python :linenos: The file$CLASSHG/codes/particles/rotate_plotter.py is not included here but can be viewed here <../../../codes/particles/rotate_plotter.py>_ and is available in the repository. (See :ref:pylab for more about plotting.) Running this code at the command line gives:: $python rotate1.py Rotating 150 particles through 40 theta values CPU time: 7.27665200 seconds It should also produce a sequence of figures similar to this animation <../../images/rotate-movie.gif>_. The output above shows that running this code took 7.2 seconds of CPU time. **You should be shocked at how long this took.** This is partly because plotting takes a lot of CPU time. This is an important issue in itself for scientific computing, see :ref:visualization. For now we do not want to be confused by this so let's make a version of the code that leaves out the plotting and just computes the locations of the particles for each value of theta. .. literalinclude:: ../codes/particles/rotate2.py :language: python :linenos: This takes much less time than the previous version::$ python rotate2.py Rotating 150 particles through 40 theta values CPU time: 0.26613700 seconds Note that this is much quicker, but is still painfully slow. All we are doing is computing 6000 particle locations on a machine that in priciple is capable of doing a few hundred million floating point operations per second (see :ref:timings). There are several major things wrong with this program. We will make some improvements by examing each in turn. Eliminating repeated computations --------------------------------- One problem with the code above is that the trigonometric functions *cos* and *sin* are being evaluated over and over again with the same arguments in the loops. Evaluating these functions takes longer than doing a single multiply or add and it is much better to evaluate the required values once and reuse them as needed. Here is an improved version of the code: .. literalinclude:: ../codes/particles/rotate3.py :language: python :linenos: This cuts the time down by a factor of 10 or so:: $python rotate3.py Rotating 150 particles through 40 theta values CPU time: 0.03809200 seconds This may seem pretty good, but note that this is for only 150 particles. Real applications often involve millions of particles. If we scaled this up to 150,000 particles we'd expect it to take at least 100 times longer, as in fact we see if we modify *N* in *rotate3.py*::$ python rotate3.py # after modifying N Rotating 150000 particles through 40 theta values CPU time: 35.73830700 seconds Vectorizing code and eliminating loops -------------------------------------- The major problem with the code above is that the inner loop on *j*, where we loop over all the particles, is written as a loop. *Loops should be avoided whenever possible* with interpreted languages like Python or Matlab. Instead, it is much better to use vectorized operations that are applied to the vectors *x* and *y* rather than to individual elements. An improved version takes the form: .. literalinclude:: ../codes/particles/rotate4.py :language: python :linenos: Note that we have kept *N = 150000* here. Even so, this code executes in under 1 second:: $python rotate4.py Rotating 150000 particles through 40 theta values CPU time: 0.27013900 seconds One might also be able to replace the loop on *theta* values by vectorized operations to compute all the locations at all times, but if one wants to add back in the plotting or other operations to be done for each value of *theta*, this may not be advantageous. Linear algebra version ---------------------- It is possible to improve this even further, however, by applying a bit of linear algebra. The location of a single particle at each theta can be computed by taking its location *(x,y)* at the previous value of *theta*, when viewed as a 2-vector, and multiplying by a *rotation matrix* of the form :math:R = \left[\begin{array}{cc}\cos(\Delta\theta) & -\sin(\Delta\theta) \\ \sin(\Delta\theta) & \cos(\Delta\theta) \end{array}\right]. where :math:\Delta\theta is the change in rotation angle, which is called *dtheta* in the code. We do not want to loop over all particles doing this, however, or we will be back to a very slow code because of the loop. Instead, notice if that we define a :math:2\times N matrix *B* whose first row consists of the *x* vector and whose second row is the *y* vector (so that the j'th column is the location vector for the j'th particle), then we can update all particle locations simultaneously simply by computing the matrix-matrix product *RB*, where *R* is again the :math:2\times 2 rotation matrix. This gives the following (note that *B* is called *xy* in the code): .. literalinclude:: ../codes/particles/rotate5.py :language: python :linenos: This gives slightly improved timing because Python takes advantage of special routines that optimize matrix-matrix multiplication::$ python rotate5.py Rotating 150000 particles through 40 theta values CPU time: 0.16931200 seconds This is about 200 times faster than the *rotate3.py* code, and probably about 2000 times faster than our original code would have been. We might now get brave and try increasing N by another factor of 100, using 15 million particles:: \$ python rotate5.py Rotating 15000000 particles through 40 theta values CPU time: 25.91272300 seconds This is more like the speed expected from this computer! Note that doing the matrix-matrix multiply for the case with *N* particles requires doing *4N* multiplies and *2N* additions, so a total of about *240N* floating point operations and hence 3.6 billion operations when *N* is 15 million. Dividing this by 25 seconds indicates that we have achieved a speed of about 144 Mflops (million floating point operations per second). Further improvements are possible. In particular, this code as written is being executed on only one core on this dual core laptop, while the other sits idle. If we can take advantage of both cores, we may be able to cut the computing time in half. See :ref:parallel for more discussion of this. We will return to this example later in the quarter.