Tuesday, June 12, 2012

Parallel computing in Python for the masses


Case scenario: you wrote a python routine that does some kind of time-consuming computation. Then you think, wow, my computer has N cores but my program is using only one of them at a time. What a waste of computing resources. Is there a reasonably easy way of modifying my code to make it exploit all the cores of my multicore machine?

The answer is yes and there are different ways of doing it. It depends on how complex your code is and which method you choose to parallelize your computation.

I will talk here about one relatively easy way of speeding up your code using the multiprocessing python package. I should mention that there are many other options out there but the multiprocessing package comes pre-installed with any python distribution by default.

I am assuming that you really need to make your code parallel. You will have to stop and spend time thinking about how to break your computation in smaller parts that will be sent to the different cores. And I should mention that debugging is harder for parallel code compared to serial code, obviously.

Parallelization is one way of optimizing your code. Other ideas for optimizing your code is using Cython or f2py. Both these approaches may imply >10x speedup and are worth exploring depending on your situation. But both will involve using the C or Fortran languages along with your python code.

The ideal case is when your problem "embarassingly parallel". What I mean by this is: your problem  can be made parallel in a reasonably easy way since the computations which correspond to the bottleneck of the code can be carried out independently and do not need to communicate between each other. Examples:

  • You have a "grid" of parameters that you need to pass to a time-consuming model (e.g., a 1000x1000 matrix with the values of two parameters). Your model needs to evaluate those parameters and provide some output.
  • Your code performs a Monte Carlo simulation with 100000 trials which are carried out in a loop. You can then easily "dismember" this loop and send it to be computed independently by the cores in your machine.

Instead of giving code examples myself, I will point out the material I used to learn parallelization. I learned the basics of parallel programming by reading the excellent tutorial "introduction to parallel programming" written by Blaise Barney.

The next step was learning how to use the multiprocessing package. I learned this with the examples posted in the AstroBetter blog. I began by reading the example implemented with the pprocess package. The caveat here is that 'pprocess' is a non-standard package. The multiprocessing package which comes with python should be used instead. Somebody posted the original example discussed in the blog ported to the multiprocessing package.

As the posts above explain, the basic idea behind using 'multiprocessing' is to use the parallel map method to evaluate your time-consuming function using the many cores in your machine. Once you figure out a way of expressing your calculation in terms of the 'map' method, the rest is easy.

In my experience doing parallel programming in python using 'multiprocessing' I learned a few things which I want to share:

  1. Do not forget to close the parallel engine with the close() method after your computation is done! If you do not do this, you will end up leaving a lot of orphan processes which can quickly consume the available memory in your machine.
  2. Avoid using lambda functions when passing arguments to the parallel 'map' at all costs! Trust me,  multiprocessing does not play well with lambda constructs.
  3. Finally, as I mentioned before, parallelizing a code increases development time and the complexity of debugging your code. Only resort to parallelization if you really need it, i.e. if you think you will get a big speedup in your code execution. For example, if you code takes 24 hours to execute and you think you can get a 6x speedup by resorting to 'multiprocessing', then the execution time can be reduced to 4 hours which is not bad.

Thursday, May 24, 2012

Hidden features of Python

I learned about this link with many useful hidden features of Python via Eduardo.

I particularly like:

  • the use of enumerator in loops: for i,x in enumerate(array)
  • decorators as a simple way of enhancing methods: @method
  • one-line swapping of variables: a,b=b,a

Wednesday, May 9, 2012

Parallel programming with Python: coming soon

Coming soon here: parallel programming in python for the masses!

I will write a tutorial as soon as I have time from my work duties.

Stay tuned.

Thursday, April 26, 2012

Using git to manage source code and more

Recently I learned how to use Git to manage source code (thanks to this guy). Let me tell you, it is such a fantastic tool! Especially when you have thousands of line of source code constantly evolving and you need to keep track of what changes.

In my case, I have been using it to manage the source code I wrote for my different scientific projects. And I will soon begin using it even to manage the writing of one paper.

Let me list the tutorials that I read and have been very useful in getting me started quickly:

  • Git Magic: I began learning git with this one. It goes straight to the point and illustrates the most important commands.
  • Pro Git: need more detailed information and have more time to spend learning git? Have a look at this one.

Quick reference for gitds:

I use SourceTree, a GUI on Mac, to check the evolution of the source code.



Changelog
May 24th 2012: replaced suggestion of GUI GitX -> SourceTree.

Wednesday, April 11, 2012

How to easily do error propagation with Python

You can easily do error propagation using the uncertainties package in Python, without having to estimate analytically the propagated error or doing Monte Carlo simulations.

Example: Suppose you have two arrays x and y which have associated uncertainties errx and erry. Using x and y, you calculate some function z = f(x, y) which can be arbitrarily complicated (but can be expressed in an analytical form) and you want to estimate the resulting uncertainty errz in z from errx and erry.

The script below

  • defines the arrays x, y, errx and erry using numpy
  • defines a function z=log10(x+y^2) for illustration purposes
  • demonstrates how to invoke 'uncertainties' in order to estimate the uncertainty in z from errx and erry (i.e., errz = f(errx, erry) )


Requirements:



 import uncertainties as unc  
 import uncertainties.unumpy as unumpy  
 import numpy  
 import nemmen  
   
 # Defines x and y  
 x=numpy.linspace(0,10,50)  
 y=numpy.linspace(15,20,50)  
   
 # Defines the error arrays, values follow a normal distribution  
 # (method random_normal defined in http://astropython.blogspot.com/2012/04/how-to-generate-array-of-random-numbers.html)  
 errx=nemmen.random_normal(0.1,0.2,50);     errx=numpy.abs(errx)  
 erry=nemmen.random_normal(0.3,0.2,50);     erry=numpy.abs(erry)  
   
 # Defines special arrays holding the values *and* errors  
 x=unumpy.uarray(( x, errx ))  
 y=unumpy.uarray(( y, erry ))  
   
 """  
 Now any operation that you carry on xerr and yerr will   
 automatically propagate the associated errors, as long  
 as you use the methods provided with uncertainties.unumpy  
 instead of using the numpy methods.  
   
 Let's for instance define z as   
 z = log10(x+y**2)  
 and estimate errz.  
 """  
 z=unumpy.log10(x+y**2)  
   
 # Print the propagated error errz  
 errz=unumpy.std_devs(z)  
 print errz  

Update Oct. 23rd 2014: code snippet available on github.

How to generate an array of random numbers following a normal distribution with arbitrary mean and standard deviation

Here is a recipe to generate an array of random numbers following a normal distribution with the supplied mean and standard deviation in Python:

 def random_normal(mean,std,n):  
      """  
 Returns an array of n elements of random variables, following a normal   
 distribution with the supplied mean and standard deviation.  
      """  
      import scipy  
      return std*scipy.random.standard_normal(n)+mean  

Update Oct. 23rd 2014: code snippet available on github. 

Friday, March 23, 2012

BCES code ported to Python!

Just ported the good old BCES linear regression algorithm (the four methods described in the Akritas & Bershady 1996 paper), originally written in Fortran 77, to Python!!

I am testing the python code and making sure it produces exactly the same answers as the old code.

More news soon...

Update (March 26th 2012): implemented the bootstrapping BCES!

Update (May 19th 2014): the BCES python code is now available as a Github project. Feel free to download it and even contribute improvements.