Today we'll look at performing parallel computations with IPython. Much of this material will draw from resources available online, including
Parallel computing is generally very difficult, and IPython provides primitives for a number of scenarios, including:
In addition, IPython parallel includes tools for executing parallel jobs interactively (this is the "I" in "IPython").
One common complaint levied against IPython is that it's too bloated. Why does a project which started as an enhanced Python interpreter include Parallel tools? You can hear the answer from IPython's creators in This Scipy 2012 Keynote. In it, Fernando Perez says IPython is
"A tool for managing the entire lifecycle of a scientific idea"
This life-cycle includes the following:
Within this context, you can see how the IPython interpreter, parallel tools, and the notebook all fit the vision.
Additionally, it is a natural fit because the underlying message-passing architecture necessary for the IPython notebook is also a natural platform for parallel computing.
The IPython engine is a Python kernel which takes commands over a network connection. There is one engine per core.
An important property of an IPython engine is that it blocks while user code is being executed, which allows asynchronous execution.
The IPython controller provides an interface for working with a set of engines.
At a general level, the controller is a collection of
processes to which IPython engines and clients can connect. The
controller is composed of a Hub and a collection of
Schedulers, which may be in processes or threads.
The controller provides a single point of contact for users who
wish to utilize the engines in the cluster. There is a variety of
different ways of working with a controller, but all of these
models are implemented via the View.apply method, after
constructing View objects to represent different collections engines.
The two primary models for interacting with engines are:
Advanced users can readily extend the View models to enable other styles of parallelism.
The center of an IPython cluster is the Hub. The Hub can be viewed as an über-logger, which keeps track of engine connections, schedulers, clients, as well as persist all task requests and results in a database for later use.
All actions that can be performed on the engine go through a Scheduler. While the engines themselves block when user code is run, the schedulers hide that from the user to provide a fully asynchronous interface to a set of engines. Each Scheduler is a small GIL-less function in C provided by pyzmq (the Python load-balanced scheduler being an exception).
The first thing required for using IPython Parallel is to start the cluster. By default, this will start a local cluster on your own computer, utilizing your CPU's multiple cores if available.
There are two ways to do this.
On the command-line, you can run the following command to start an IPython cluster:
[~]$ ipcluster start -n 4
will start a cluster with 4 cores (probably what you have available on a modern laptop). To use as many cores as are available, leave out the n:
[~]$ ipcluster start
This spins-up a Python kernel on each core, along with a hub that connects to them.
Within the IPython notebook, you can use the Clusters tab of the dashboard, and press Start with the desired number of cores, under the desired profile (more on profiles later). This will automatically run the correct commands to start your IPython cluster.
In [1]:
from IPython import parallel
clients = parallel.Client()
clients.block = True # use synchronous computations
print clients.ids
Now we can use this object to run some code. Let's define a simple function that we will execute on the cores:
In [2]:
def mul(a, b):
return a * b
mul(5, 6)
Out[2]:
Now let's execute mul on the first engine:
In [3]:
clients[0].apply(mul, 5, 6)
Out[3]:
Notice that the function becomes the first argument of apply, and any
additional arguments or keyword arguments are passed after it.
We can similarly execute mul in parallel on all the engines at once:
In [4]:
clients[:].apply(mul, 5, 6)
Out[4]:
But say we want to pass different arguments to each instance of mul. In normal Python, we could do this with the map() function:
In [5]:
map(mul, [5, 6, 7, 8], [8, 9, 10, 11])
Out[5]:
The client views also have a map function which does the execution in parallel. Let's use a load-balanced view so that the IPython scheduler takes care of farming the tasks for us:
In [6]:
view = clients.load_balanced_view()
view.map(mul, [5, 6, 7, 8], [8, 9, 10, 11])
Out[6]:
This is the basic interface to parallel computation in IPython!
In [7]:
clients = parallel.Client(profile='default')
print len(clients)
In [8]:
clients.block = True
dview = clients.direct_view()
dview.block = True
print dview
In [9]:
dview.apply(sum, [1, 2, 3])
Out[9]:
You can also obtain a Direct view by slicing the Clients object:
In [10]:
clients[::2]
Out[10]:
In [11]:
clients[::2].apply(sum, [1, 2, 3])
Out[11]:
In [12]:
lview = clients.load_balanced_view()
print lview
In [13]:
lview.apply(sum, [1, 2, 3])
Out[13]:
One way to introspect what's going on is to get the Process ID of the engines. This gives a unique label of any running process on your computer. If you have a cluster running on multiple computers, you can also print the hostname to see which host each cluster is running on:
In [14]:
import os
import socket
print os.getpid()
print socket.gethostname()
In [15]:
dview.apply(os.getpid)
Out[15]:
In [16]:
for i in range(5):
print lview.apply(os.getpid)
Notice that the load-balanced view chooses engines somewhat randomly (it's not actually random, but is based on a process which is opaque to the user).
The direct view is an interface in which each engine is directly exposed to the user.
The direct view is a very flexible wrapper of the client. Let's look at a few ways of managing remote execution and remote data with direct views. We'll look at the following:
In [17]:
clients = parallel.Client()
dview = clients[:]
In a DirectView interface, we can either use Blocking (synchronous) execution, in which all results must finish computing before any results are recorded, or non-blocking (asynchronous) execution, where we receive results as they finish.
Here's an example of the two:
In [18]:
def get_pid_slowly():
# imports in the function, otherwise these imports
# don't happen on the engines
import os
import time
import random
# sleep up to 10 seconds
time.sleep(10 * random.random())
return os.getpid()
In [19]:
dview.block = True
dview.apply(get_pid_slowly)
Out[19]:
Notice that all three results come back at the same time, after all are finished computing.
Let's now see what happens if we use non-blocking execution:
In [20]:
dview.block = False
dview.apply(get_pid_slowly)
Out[20]:
What is returned is an AsyncResult object.
In [21]:
res = dview.apply(get_pid_slowly)
In [22]:
res.ready()
Out[22]:
We can ask if the result is ready using the ready() command. Once it's ready, we can get the result with the get() command
In [25]:
res.ready()
Out[25]:
In [26]:
res.get()
Out[26]:
For convenience, you can also use the apply_sync and apply_async commands
In [27]:
dview.apply_sync(get_pid_slowly)
Out[27]:
In [28]:
dview.apply_async(get_pid_slowly)
Out[28]:
In [36]:
import numpy
def normalize_vector(v):
return numpy.asarray(v) / numpy.linalg.norm(v)
In [37]:
x = numpy.random.random(5)
normalize_vector(x)
Out[37]:
In [38]:
dview.block = True
dview.apply(normalize_vector, x)
Though we've imported numpy on our local engine, we haven't imported numpy on any of our other engines, so the code won't work.
There are a few ways to execute imports on other engines. We can do this directly with the execute method:
In [39]:
dview.execute('import numpy')
Out[39]:
Alternatively, we can use the sync_imports context manager to make this happen:
In [40]:
with dview.sync_imports():
import numpy
Or, we can use the px (parallel execute) magic function to execute the contents of the cell on every engine:
In [41]:
%px import numpy
After doing any of these, we're able to run our function and see the results:
In [42]:
dview.apply(normalize_vector, x)
Out[42]:
In [43]:
%%file myscript.py
import numpy
a = 5
def square(x):
return x * x
Now we can run this script on each of our clients:
In [44]:
dview.run("myscript.py")
Out[44]:
To look at the data, let's use the dictionary interface:
In [45]:
print dview['a']
In [46]:
print dview['square']
We can also use the dictionary interface to assign data:
In [47]:
dview['b'] = 4
print dview['b']
Now the variables a and b are defined globally on all our engines.
In [48]:
print dview['a']
print dview['b']
Now we'll use execute() to add a and b and save the result to c:
In [49]:
dview.execute('c=a+b')
Out[49]:
In [50]:
dview['c']
Out[50]:
We can also execute code on any subset of the engines:
In [51]:
clients[:2].execute('c=a*b')
Out[51]:
In [52]:
dview['c']
Out[52]:
This is the basic way to push and pull data to and from the client and engines.
In [53]:
%px c = numpy.random.random(3)
%px print c
If you want to execute multiple lines at once, you can use the %%px cell magic to do this:
In [54]:
%%px
c = numpy.random.random(5)
c *= a
print c
Again, we can access the variables directly using the dictionary interface:
In [55]:
dview['c']
Out[55]:
In [56]:
dview.scatter('a', np.arange(24))
dview['a']
Out[56]:
Notice that the contents of a are partitioned and scattered among the instances. To bring them back, we can use gather():
In [57]:
a = dview.gather('a')
a
Out[57]:
So, for example, you might do an operation something like this:
In [58]:
x = np.random.random(400)
dview.scatter('x', x)
dview.execute('y = numpy.sqrt(x)')
y = dview.gather('y')
print np.allclose(y, np.sqrt(x))
This is a somewhat trivial computation, but if you had a larger computation that could be done on partitions of the data, it would allow you to speed your calculation according to the number of cores available.
In [59]:
A = np.random.random((8, 3))
dview.scatter('A', A)
print A
print
for chunk in dview['A']:
print chunk
So scattering a matrix breaks it up into blocks... we can perform a matrix multiplication on each of these blocks, and gather the individual results to get the final result!
In [60]:
def parallel_dot(dview, A, B):
dview.scatter('A', A)
dview['B'] = B
dview.execute('C = numpy.dot(A, B)')
return dview.gather('C')
A = np.random.random((10, 3))
B = np.random.random((3, 5))
np.allclose(parallel_dot(dview, A, B),
np.dot(A, B))
Out[60]:
Direct views are a nice low-level interface to the clients, but for most tasks it can be more helpful to use Load-balanced views.
A load-balanced view is a fault-tolerant way to distribute tasks among the engines. While there is no direct access to individual engines, it provides a simple and powerful interface.
In [61]:
lview = clients.load_balanced_view()
In [62]:
def sleep_and_return(factor=10):
import time
import random
r = factor * random.random()
time.sleep(r)
return r
In [63]:
lview.block = True
lview.apply(sleep_and_return, 1)
Out[63]:
Using apply on a load balanced view simply executes the function once.
We can use the map function to execute multiple inputs. We'll use block=False to make sure the results come asynchronously:
In [67]:
lview.block = False
res = lview.map(sleep_and_return, range(10))
In [68]:
res.ready()
Out[68]:
The progress attribute tells us how many of the jobs have been completed up to this point:
In [69]:
import time
while not res.ready():
time.sleep(1)
print res.progress
Once all the results are finished, we can use the result attribute to get at them:
In [70]:
res.result
Out[70]:
The load-balanced view also includes an iterator which will return results as they come in:
In [71]:
for result in lview.map(sleep_and_return, range(10)):
print result
In [72]:
%pylab inline
x, y = 2 * np.random.random((2, 500)) - 1
ax = plt.axes(aspect='equal', xticks=[], yticks=[], frameon=False)
ax.scatter(x, y, c=(x ** 2 + y ** 2 < 1), cmap=plt.cm.binary)
ax.add_patch(mpl.patches.Circle((0, 0), 1, zorder=0, fc='#FFAAAA'))
Out[72]:
The probability that those darts will fall within 1 unit of the center is
$$ p = \frac{\rm Area~of~circle}{Area~of~square} = \frac{\pi\cdot 1^2}{2^2} = \frac{\pi}{4} $$So if you throw $N$ darts, you can estimate $\pi$ to be equal to $$ \pi = 4\frac{N_{in}}{N} $$
Here we'll write a function to estimate $\pi$ using a random number generator:
In [73]:
def estimate_pi(N):
r = np.random.random((N, 2))
inside = np.sum(r ** 2, 1) <= 1
return 4 * inside.sum() * 1. / N
In [74]:
estimate_pi(1E6)
Out[74]:
Write a function which uses a Load-balanced View to parallelize this estimation of pi:
In [ ]:
def estimate_pi(N, lview, N_per_trial=1E6):
# your code here:
# distribute the N trials across the load-balanced view,
# using lview.map() and return the average of the results
Note that this is certainly not the best way to estimate $\pi$, but it shows how you can use IPython parallel to do distributed Monte Carlo simulations, which can be very useful in practice.
For the homework today, we'll look at two Astronomy applications. You may choose one (or both, if you want more practice) to turn in. In both of them, we'll return to a problem that we've seen before, but extend our computation using IPython parallel.
As usual, you can turn-in your results via GitHub: use either an IPython notebook, or a stand-alone .py file.
In the Scikit-learn homework a few weeks ago, we looked at different regression schemes for computing photometric redshifts (see Exercise 2). Return to this, and determine the minimum rms you can achieve given combinations of max_depth and n_estimators in the random forest regressor
Here's how we'll download and prepare the data:
In [ ]:
from astroML.datasets import fetch_sdss_specgals
import numpy as np
data = fetch_sdss_specgals()
# put magnitudes in a matrix
X = np.vstack([data['modelMag_%s' % f] for f in 'ugriz']).T
y = data['z']
# down-sample the data for speed
X = X[::10]
y = y[::10]
Here's an example of using cross-validation to compute the r2 score for a given C:
In [ ]:
from sklearn.ensemble import RandomForestRegressor
from sklearn.cross_validation import cross_val_score, StratifiedKFold
# define the model with the given value of C
model = RandomForestRegressor(n_estimators=10, max_depth=3)
# compute the scores via cross-validation
scores = cross_val_score(model, X, y,
scoring='r2', cv=3)
# print the mean of the cross-validation scores
score = np.mean(scores)
print score
Write a script that uses a parallel computation to explore what the best value of max_depth and n_estimators is. What's the maximal score you can obtain?
In the lecture on optimization, we looked at computing the periods of LINEAR stars using the lomb-scargle algorithm.
Write a script which uses IPython parallel to compute the periods of ~1000 stars from the LINEAR dataset. Plot the log of the period you found vs the g-i color of the star. You can compare your results to this paper, though the Lomb-Scargle period-finder might lead to slightly different results.
Here's how to get the data:
In [ ]:
from astroML.datasets import fetch_LINEAR_sample
lcs = fetch_LINEAR_sample()
id_ = lcs.ids[0]
time, flux, dflux = lcs[id_].T
And here's how to compute the lomb-scargle periodogram:
In [ ]:
from astroML.time_series import lomb_scargle
periods = np.logspace(-2, 0, 10000)
periodogram = lomb_scargle(time, flux, dflux, omega=2 * np.pi / periods, generalized=True)
Now we find the maximum period:
In [ ]:
idx = np.argsort(periodogram)
print periods[idx[-1]]
And we get the $g-r$ color:
In [ ]:
g_r = lcs.get_target_parameter(id_, 'gr')
print g_r
You'll have to automate the task of computing the best period, and use a LoadBalancedView to distribute this across as many cores as you have available.