<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 After speaking to a number of academics involved with Monte Carlo simulations, I realised many of them resort to creating bash scripts to 
run a number of python instances of simulation before writing another script to piece the data back together. It turns out this is incredibly simple 
to do in IPython, and in particular in an IPython Notebook. We take a look at how to run computations in parallel as well as giving a use-case in the creation of Julia fractals.
Most modern CPUs have 2 or more cores, allowing for everyone to take full advantage of parallel computing. The IPython implementation here is probably most useful for when
you have some computationally intensive program which is not called on a regular basis. Writing a full parallel implementation can be time consuming (depending on the problem and the 
experience of the programmer), however as we'll see, writing an implementation in IPython can be quick and painless. Of course more tailored parallel implementations require more effort 
should use a more established MPI such as OpenMPI. Thankfully this can be done within the IPython framework.
To begin, we will set up the engines that we will use for our computation. 
Basics
The client is handler for all of the engines (kernels) on a cluster.
Step1:  Before we start, we first need to assign a number of engines in our cluster. This can be done through the IPython notebook interface or through the command line. The basic syntax is
$ ipcluster start -n 4
which starts a cluster of 4 engines to compute on. For more help use $ ipcluster -h
or consult the documentation to set up engines with specific settings or parameters by creating a profile.
We begin by calling the Client() function.
Step2:  If there are no engines running, or they are incorrectly set up then this will throw an error.
The way we access the engines is through client views. There are two types of view
Step3:  One can also select a subset of the engines directly using list slicing
Step4:  Load Balanced Views
A load balanced view allows the computation of a task on any one engine. The engine is chosen by the scheduler using some black box logic.
Step5:  How do we execute code?
There are many ways to execute code in a Direct view. As an example we will get each engine to return its process ID (PID).
If we were considering one engine only we would run
Step6:  As shown before, we can apply a function (with parameters if needed) to the view, which returns the result for each engine in a list.
Step7:  Alternatively we can apply a function decorator to our function and call it as normal. A function decorator is a wrapper function, or a function which has one or more functions as a parameter.
Step8:  Another method uses the cell magic %%px or line magic %px. This executes the line or cell on a direct view of the current client (specifically its current list of targets).
Step9:  Finally, and most fundamentally, one can excute a string of code using the execute function.
Step10:  Imports and Environments
Chances are if one is concerned with computation speed, you will be using numpy and scipy or any number of packages or custom scripts. As with executing code, there are multiple ways to manage imports
Step11:  The later method uses the context manager sync_imports which simultaneously sets the local and remote imports.
It is also possible for all remote engines to run a centralised script. This could be used for imports, defining application functions, or actual computation.
Step12:  Since the computational engines are all running with IPython, then %magics are supported, as well as inline plotting, for example,
Step13:  Accessing and Setting Data
Accessing and setting data can be done in multiple ways (spotting a theme here?). All variables and objects on each engine can be accessed like a dictionary using the direct view. We are not limited to accessing purely what we have assigned, anything in the global namespace on the remote engine is fair game.
Step14:  Values can be set and retrieved using the more conventional push and pull methods respectively.
Step15:  Arrays can be distributed evenly over the engines using the scatter method.
Step16:  This allows us to compute the result for each array before using the gather method to piece them back together.
Step17:  Load Balanced Views
Direct views are simple to understand and are extremely useful when we know how to split our computation up efficiently. In some cases however there may be a list of computational tasks with undetermined computational runtime which is where a load balanced view becomes useful. The load balanced view was defined earlier as
Step18:  If we apply a function to the load balanced view we receive only one result as the function has been only been executed once.
Step19:  To pass in multiple inputs (effectively creating a queue of parameters), we can use the map function. By setting bview.block=False we allow the results to stream in asynchronously.
Step20:  The map instantly returns an AsyncResult object which begins empty and collects the results as they are computed. To check whether the computation is complete we run
Step21:  which in this case is False. After waiting a sufficient amount of time, the complete result becomes available
Step22:  If we call result before the computation is complete then the system will wait until it is. Note that there is no particular ordering to the engine used - the scheduler is picking the next available engine to run the function. The load balanced view map is also iterable which allows for results to be returned as and when they become available.
Step23:  Applying Parallel Computing to Fractals (Julia sets)
Now for an interesting application. Since everyone has estimated $\pi$ a thousand times before we'll avoid that example. Instead we'll look at visualising Julia sets, closely following this well written guide to Numpy/SciPy Recipes for Image Processing.
Theory in Brief
We consider a rational function $f
Step24:  Note that Python uses the letter j rather than i for the complex unit.
We can use this function as an iterative map on the complex plane, i.e we have $z_t = f(z_{t-1}) = \dots = f^t(z_0)$. There are two possibilities which can happen, either $|z_t|$ diverges as $t \to \infty$, or $|z_t|$ remains bounded.
Points where $f^t(z)$ diverge are called the escape set of $f$, and points which do not diverge are called the prisoner set of $f$. The Julia set of $f$ is then defined to be the boundary of the prisoner set, i.e,
$$J(f) = \partial  \left{ z \in \mathbb{C} \left\vert  \frac{}{}\right . \lim_{t \to \infty} |f^t(z)|< \infty \right} $$
The interesting fact about this dynamical system is that it behaves chaotically for almost all $c$, that is to say that two nearby points $z_0$ and $z_0 + \epsilon$ will likely yield very different trajectories in the large $t$ limit.
To begin with we define the section of the complex plane we want to explore (zmin and zmax) as well as our resolution (m and n). The resolution defines how many inital conditions we are going to sample, which in this case is $mn=1024^2$.
Step25:  Using the meshgrid function we can easily create an array of complex numbers.
Step26:  Application
The Julia set itself is very difficult to be determined exactly and so we plot the escape set of the function instead.
For each initial condition we count the number of iterations until the map escapes a threshold radius, and plot this.
Step27:  The above function carries out this function in a very succinct fashion. Firstly, we initialise an array with our maximum iteration number. We then check our complex plane array to see if any values have a modulus greater than 2. This creates a boolean index array which we can use to select only the elements in our array Z which satisfy the condition.
Iterating the process tmax times we then arrive at our array J which distinguishes between escape times for each initial conditions. We can also run the algorithm in the other direction for the same result. The use of the boolean arrays means that if an iteration of the function escapes the threshold radius we do not ever apply the map again to the function - saving valuable computational time!
Below we first create a copy of Z since our julia function edits the array in-place and we will want to use the original again later.
Step28:  In parallel we first feed the functions to each engine.
Step29:  We can then scatter/distribute the array amoungst the different engines. Each engine runs the same number of loops on the same sized array so in theory the load should be evenly distributed. In reality the size of the masks may vary so one would expect some deviation in the workloads.
Step30:  Piecing together all the parts we arrive at the same fractal as before.
Step31:  So did we speed it up?
Step32:  So we see that the computational time has indeed been reduced by a factor of the number of engines we chose at the start (plus some overhead). For this small example the overhead is relatively large, however if we wanted to double or triple the resolution then this overhead would become insignificant.
Of course no speed test is complete without checking that the answers agree,
<ASSISTANT_TASK:>
Python Code:
from IPython.parallel import Client
%pylab inline
Explanation: After speaking to a number of academics involved with Monte Carlo simulations, I realised many of them resort to creating bash scripts to 
run a number of python instances of simulation before writing another script to piece the data back together. It turns out this is incredibly simple 
to do in IPython, and in particular in an IPython Notebook. We take a look at how to run computations in parallel as well as giving a use-case in the creation of Julia fractals.
Most modern CPUs have 2 or more cores, allowing for everyone to take full advantage of parallel computing. The IPython implementation here is probably most useful for when
you have some computationally intensive program which is not called on a regular basis. Writing a full parallel implementation can be time consuming (depending on the problem and the 
experience of the programmer), however as we'll see, writing an implementation in IPython can be quick and painless. Of course more tailored parallel implementations require more effort 
should use a more established MPI such as OpenMPI. Thankfully this can be done within the IPython framework.
To begin, we will set up the engines that we will use for our computation. 
Basics
The client is handler for all of the engines (kernels) on a cluster.
End of explanation
client = Client()
client.block = True # Computations run synchronously.
print client.ids
Explanation: Before we start, we first need to assign a number of engines in our cluster. This can be done through the IPython notebook interface or through the command line. The basic syntax is
$ ipcluster start -n 4
which starts a cluster of 4 engines to compute on. For more help use $ ipcluster -h
or consult the documentation to set up engines with specific settings or parameters by creating a profile.
We begin by calling the Client() function.
End of explanation
dview = client.direct_view()
def f(x):
    return x
dview.apply(f, "Hello World")
Explanation: If there are no engines running, or they are incorrectly set up then this will throw an error.
The way we access the engines is through client views. There are two types of view:
Direct Views
Direct views allows the computation of a task on all of the engines simultaneously:
End of explanation
even_dview = client[::2]
even_dview.apply(f, "Hello World")
Explanation: One can also select a subset of the engines directly using list slicing:
End of explanation
bview = client.load_balanced_view()
bview.apply(f, "Hello World")
Explanation: Load Balanced Views
A load balanced view allows the computation of a task on any one engine. The engine is chosen by the scheduler using some black box logic.
End of explanation
import os
print os.getpid()
Explanation: How do we execute code?
There are many ways to execute code in a Direct view. As an example we will get each engine to return its process ID (PID).
If we were considering one engine only we would run:
End of explanation
dview.apply(os.getpid)
Explanation: As shown before, we can apply a function (with parameters if needed) to the view, which returns the result for each engine in a list.
End of explanation
@dview.remote(block=True)
def getpid():
    import os
    return os.getpid()
getpid()
Explanation: Alternatively we can apply a function decorator to our function and call it as normal. A function decorator is a wrapper function, or a function which has one or more functions as a parameter.
End of explanation
%%px 
import os
print os.getpid()
Explanation: Another method uses the cell magic %%px or line magic %px. This executes the line or cell on a direct view of the current client (specifically its current list of targets).
End of explanation
res = dview.execute('print os.getpid()')
res.display_outputs()
Explanation: Finally, and most fundamentally, one can excute a string of code using the execute function.
End of explanation
dview.execute('import numpy')
%px import numpy
with dview.sync_imports():
    import numpy
Explanation: Imports and Environments
Chances are if one is concerned with computation speed, you will be using numpy and scipy or any number of packages or custom scripts. As with executing code, there are multiple ways to manage imports:
End of explanation
%%file script.py
def hello_world():
    print "Hello World"
hello_world()
dview.run("script.py").display_outputs()
Explanation: The later method uses the context manager sync_imports which simultaneously sets the local and remote imports.
It is also possible for all remote engines to run a centralised script. This could be used for imports, defining application functions, or actual computation.
End of explanation
%%px 
%pylab inline
z = np.random.random(10)
hist(z);
Explanation: Since the computational engines are all running with IPython, then %magics are supported, as well as inline plotting, for example,
End of explanation
dview['a'] = 5
dview['b'] = 10
print dview['a']
print dview['os.getpid']
Explanation: Accessing and Setting Data
Accessing and setting data can be done in multiple ways (spotting a theme here?). All variables and objects on each engine can be accessed like a dictionary using the direct view. We are not limited to accessing purely what we have assigned, anything in the global namespace on the remote engine is fair game.
End of explanation
dview.push(dict(a=1, b=2))
dview.pull(('a','b'))
Explanation: Values can be set and retrieved using the more conventional push and pull methods respectively.
End of explanation
a = np.arange(32)
dview.scatter('a', a)
%px print a
Explanation: Arrays can be distributed evenly over the engines using the scatter method.
End of explanation
%px a = a**2
dview.gather('a')
Explanation: This allows us to compute the result for each array before using the gather method to piece them back together.
End of explanation
def random_number(factor):
    import time
    import os
    import random
    r = random.random()
    time.sleep(5*r)
    return os.getpid(), factor*r
Explanation: Load Balanced Views
Direct views are simple to understand and are extremely useful when we know how to split our computation up efficiently. In some cases however there may be a list of computational tasks with undetermined computational runtime which is where a load balanced view becomes useful. The load balanced view was defined earlier as:
python
bview = client.load_balanced_view()
Below we have a function which will hopefully illustrate how the view behaves. The function sleeps for a random amount of time before returning the PID of the engine and a random number multiplied by an inputted factor.
End of explanation
bview.block=True
bview.apply(random_number,1)
Explanation: If we apply a function to the load balanced view we receive only one result as the function has been only been executed once.
End of explanation
bview.block = False
res = bview.map(random_number, range(10))
Explanation: To pass in multiple inputs (effectively creating a queue of parameters), we can use the map function. By setting bview.block=False we allow the results to stream in asynchronously.
End of explanation
res.ready()
Explanation: The map instantly returns an AsyncResult object which begins empty and collects the results as they are computed. To check whether the computation is complete we run
End of explanation
print res.ready()
print res.result
Explanation: which in this case is False. After waiting a sufficient amount of time, the complete result becomes available:
End of explanation
for result in bview.map(random_number, range(10)):
    print result
Explanation: If we call result before the computation is complete then the system will wait until it is. Note that there is no particular ordering to the engine used - the scheduler is picking the next available engine to run the function. The load balanced view map is also iterable which allows for results to be returned as and when they become available.
End of explanation
def f(z, c=-.065+.66j):
    return z**2 + c
Explanation: Applying Parallel Computing to Fractals (Julia sets)
Now for an interesting application. Since everyone has estimated $\pi$ a thousand times before we'll avoid that example. Instead we'll look at visualising Julia sets, closely following this well written guide to Numpy/SciPy Recipes for Image Processing.
Theory in Brief
We consider a rational function $f: \mathbb{C} \to \mathbb{C}$. In this example we'll consider the function $f(z) = z^2 + c$ where $c \in \mathbb{C}$. Written in Python this function is:
End of explanation
m=1024
n=1024
zmin = -1.3 - 1j * 1.3
zmax = 1.3 + 1j * 1.3
xs = np.linspace(zmin.real, zmax.real, n)
ys = np.linspace(zmin.imag, zmax.imag, m)
X, Y = np.meshgrid(xs, ys)
Z = X + 1j * Y
Explanation: Note that Python uses the letter j rather than i for the complex unit.
We can use this function as an iterative map on the complex plane, i.e we have $z_t = f(z_{t-1}) = \dots = f^t(z_0)$. There are two possibilities which can happen, either $|z_t|$ diverges as $t \to \infty$, or $|z_t|$ remains bounded.
Points where $f^t(z)$ diverge are called the escape set of $f$, and points which do not diverge are called the prisoner set of $f$. The Julia set of $f$ is then defined to be the boundary of the prisoner set, i.e,
$$J(f) = \partial  \left{ z \in \mathbb{C} \left\vert  \frac{}{}\right . \lim_{t \to \infty} |f^t(z)|< \infty \right} $$
The interesting fact about this dynamical system is that it behaves chaotically for almost all $c$, that is to say that two nearby points $z_0$ and $z_0 + \epsilon$ will likely yield very different trajectories in the large $t$ limit.
To begin with we define the section of the complex plane we want to explore (zmin and zmax) as well as our resolution (m and n). The resolution defines how many inital conditions we are going to sample, which in this case is $mn=1024^2$.
End of explanation
print Z[:2,:2]
Explanation: Using the meshgrid function we can easily create an array of complex numbers.
End of explanation
def julia(f, Z, tmax=256):
    J = np.ones(Z.shape) * tmax
    for t in xrange(tmax):
        mask = np.abs(Z) <= 2.
        Z[ mask] = f(Z[mask])
        J[-mask] -= 1
    return J
Explanation: Application
The Julia set itself is very difficult to be determined exactly and so we plot the escape set of the function instead.
For each initial condition we count the number of iterations until the map escapes a threshold radius, and plot this.
End of explanation
ZZ = np.copy(Z)
J = julia(f, ZZ)
figure(figsize=[10,10])
imshow(J)
axis('off');
Explanation: The above function carries out this function in a very succinct fashion. Firstly, we initialise an array with our maximum iteration number. We then check our complex plane array to see if any values have a modulus greater than 2. This creates a boolean index array which we can use to select only the elements in our array Z which satisfy the condition.
Iterating the process tmax times we then arrive at our array J which distinguishes between escape times for each initial conditions. We can also run the algorithm in the other direction for the same result. The use of the boolean arrays means that if an iteration of the function escapes the threshold radius we do not ever apply the map again to the function - saving valuable computational time!
Below we first create a copy of Z since our julia function edits the array in-place and we will want to use the original again later.
End of explanation
%%px
import numpy as np
def f(z, c=-.065+.66j):
    return z**2 + c
def julia_parallel(f,Z,tmax=256):
    J = np.ones(Z.shape) * tmax
    for t in xrange(tmax):
        mask = np.abs(Z) <= 2.
        Z[mask] = f(Z[mask])
        J[-mask] -= 1
    return J
Explanation: In parallel we first feed the functions to each engine.
End of explanation
dview.scatter('Z', Z)
%px ZZ = np.copy(Z)
%px J = julia_parallel(f,ZZ)
PJ = dview.gather('J')
Explanation: We can then scatter/distribute the array amoungst the different engines. Each engine runs the same number of loops on the same sized array so in theory the load should be evenly distributed. In reality the size of the masks may vary so one would expect some deviation in the workloads.
End of explanation
figure(figsize=[10,10])
imshow(PJ)
axis('off');
Explanation: Piecing together all the parts we arrive at the same fractal as before.
End of explanation
Z = X + 1j * Y
%%time
J = julia(f, Z)
Z = X + 1j * Y
%%time 
dview.scatter('Z', Z)
%px ZZ = np.copy(Z)
%px J = julia_parallel(f,ZZ)
PJ = dview.gather('J')
Explanation: So did we speed it up?
End of explanation
np.allclose(J,PJ)
Explanation: So we see that the computational time has indeed been reduced by a factor of the number of engines we chose at the start (plus some overhead). For this small example the overhead is relatively large, however if we wanted to double or triple the resolution then this overhead would become insignificant.
Of course no speed test is complete without checking that the answers agree,
End of explanation
<END_TASK>