- James Bednar, professor at Edinburgh University, for a cool project proposal - Chris, PhD student at Edunburgh University CSNG group for all the help he offered from the very beginning of this project till, literally, last minutes. - mpi4py developers for mpi4py and letting keep the address name for this blog - Olaf Lenz for PMI - this module saved me a lot of time! - Abe for thesis proof-reading and commenting (at 6am on Wednesday morning!) - That guy in Chris's office who had a thesis cover page template - As always, family and friends for support
P.S. One funny (-ish) fact is that actually, until september 2009 I was going to work on a completely different (all projects are allocated in May), and then suddenly changed my mind, came to Jim and asked if he had anything I could work on, and he offered me to optimise Topographica with MPI.
Well, as I have finally submitted my thesis (Have to admit, this was the first paper of such kind that I have ever written in my life and, bloody hell, that was a nightmare!) the project is officially over now, and it's time to summarise the results. Perhaps, it would be useful to reiterate the objectives first (I'll just copy bits from my thesis):
Purpose of the project:
Optimising the Open-Source Neural Map Simulator Topographica by making it capable of distributing its heaviest computations over CPU cores of one machine or over a network of independent machines would allow users of Topographica to run experiments on a bigger scale than it was previously possible. This is expected to help neuroscientists have a better understanding of brain activity at the cortical level. Such understanding, apart from expanding the pool of human knowledge, could potentially assist in the creation of medicines to combat and, hopefully, cure mental diseases, and, possibly benefit the human race in other, as yet undiscovered ways.
Primary Objective Improve the performance of Topographica simulations by using the High Performance Computing approach, allowing Topographica to take advantage of the computational power of both Multiprocessing and Cluster architectures.
Other Objectives - Research and summarise how distributing computations with the Message Passing Interface can benefit Topographica - Research what MPI software packages are currently available and can be utilised for optimising Topographica simulations - Create a strong base for future development and optimisation
Was it successful with regards to that? Can't tell, it's really not for me to decide, let the markers do their job=) Am I satisfied? Yes. Anyway, a quick summary of achievements:
- We proved that it is possible to reach high performance with MPI in Topographica simulations - x3.2 maximum speed-up factor, relative to the original single-processed implementation, has been recorded at realistic simulations with Lissom model at 142 cortex density, 24 LGN and retina densities (values used in most of examples for the CMVC book) -We analysed the current range of available MPI tools and found the most suitable up to this date modules and packages - Pinpointed the slow-down factors of both shared-memory and distributed models, on which further development efforts should be concentrated - Created a strong base for further optimisations and development work
I remember once Jim (my supervisor) has mentioned that he had posted a proposal for the next year's Honours Projects - "High Performance Cortical Simulator", basically taking this project to the next level. He also said that someone was very interested and there was a fair chance this work would be picked from up from the point where I finished. Well, all I can say about this as someone who knows the current state of things with Topographica and MPI: it's not going to be easy in any way. However, having said this, I hope that the person won't change his or her mind and will continue my work. I wish that person all the best of luck and I would be very excited to hear about their progress. Maybe this blog will even get a second author, who knows!=)
What about myself? Well, in the short run, I still have some data from the experiments that I want to post here in the next few days, so stay tuned. In the long-run, I'm applying to EPCC (Edinburgh Parallel Computing Centre), for the MSC in High Performance Computing and to some other universities that do similar courses. If I get a place (fingers crossed), I will have more stuff to write about.
Day 114 (actually, it's still Day 113 for me now, biologically, since I haven't went to bed since the last post): Well, it's time to wrap it up
Just a quick update to the last note: hardware architectures of all nodes have to be the same. I.e. making a 64-bit machine and a 32 one run one parallel task together won't work (at least I don't see how they could do it). How can I explain this? Well, before running a parallel task, you need to have your programme and all libraries that it uses (including OpenMPI) compiled. It is obvious that an executable compiled on a 64-bit machine won't run on a 32-bit machine, but MPI requires each node to be able to run the same executables, so there you go - you cannot combine machines with different architectures into one cluster to run parallel jobs, it will not work.
[update added MUCH later] While theoretically MPI can work on a heterogeneous cluster, I still have no clue how to do that and if it would make any sense with Topographica
--------------
Anyway, I found out that running topographica on 3 Jupiter nodes is faster than on one Jupiter node with the same number of processes launched, which was a nice surprise right at the end of this project. I got the maximum performance optimisation factor of 3 (i.e. MPI versoin can get up to 3 times faster than the original sequential model) and I believe that at this stage I will not make it any faster since this is almost the end of the project and I haven't even started writing my report yet.
So, this is it...
P.S. I will do some benchmarking in the following 1.5 weeks and post the results here, so stay tuned. It's not like I'm abandoning this blog now anyway ;)
Oh God, I did it! I figured out how to run topographica on a cluster, and it wasn't that easy. A very good guide that helped me was: Setting up a Beowulf Cluster Using Open MPI on Linux. Running MPI jobs points 1-4 are also useful. Just to summarise, in order to run an MPI job on a cluster, you need:
- Hostfile with hostnames of machines in the cluster - OpenMPI installed on each machine in the cluster, in the same directory - system PATH needs to know where OpenMPI's bin/ directory is - system LD_LIBRARY_PATH needs to know where OpenMPI's lib/ directory is - All common libraries that are going to be used during the runtime need to be installed in the same directory (i.e. path to them should be the same on all nodes). In my case, all topographica libraries are local and therefore I had to have it installed on each node in the same directory, e.g. ~/CODE/topographica
One note: if using an mpirun command with an absolute path to it (e.g. /usr/bin/openmpi/mpirun), then OpenMPI bin/ and lib/ directories will be added automatically to PATH and LD_LIBRARY_PATH respectively on all nodes.
A small (more or less) example of programming with MPI in C. This is a copy of my solution to the second practical for Parallel Programming Languages and Systems, which our course organiser Dr Murray Cole has kindly allowed me to post in this blog.
The algorithm calculates adaptive quadrature using MPI and the "bag of tasks" approach.
"Adaptive Quadrature is a recursive algorithm that computes an approximation of the integral of a function F(x), using static quadrature rules on adaptively reļ¬ned sub-intervals of the integration domain."
To compile and run it you need some kind of MPI libraries installed on your system (I used OpenMPI) and the following commands:
mpicc -o
to compile and
mpirun -c 5
to run. Here's the source code (a bit lengthy, need to find a way to minimise code listings...):
#define EPSILON 1e-3 #define F(arg) cosh(arg)*cosh(arg)*cosh(arg)*cosh(arg) #define A 0.0 #define B 5.0
#define SLEEPTIME 1
/*************************************************************************** Notes on implementation: Tested on my own machine with 2 CPU cores. MPI interface: OpenMPI 1.3.2-3ubuntu1.1
Description: The implementation is based on a standard "bag of tasks" technique. The farmer (aka controller) and worker functions initiate two loops:
Farmer:
After initial declarations, farmer loop is initiated with exit conditions of: - All workers have finished computing area (maintained by the int idle variable which increments every time a worker returns results and decrements each time some worker is passed task to perform) - Stack is empety (additional method isEmpty was implemented in stack.c to check that condition) The loop starts with wild-card synchronous MPI receive function. When data is received, controller checks whether it is a new task or computed area (indicated by a tag). If it is a new task, it is pushed into the stack, otherwise, adds the received value to the total area. Then idle process counter is incremented and the value of the proc_waiting array slot, corresponding to the worker's ID, changes to 1, indicating that this worker is idle. Then, if stack is not empty, the controller iterates over proc_waiting, starting from the point where it last finished (which ensures that farmer would not be hijacked by any worker), and if the current value in the array is 1, sends new task to the worker and continues the main loop. When stack is empty and all workers are idle it means that the area has been computed, farmer breaks the main loop and sends exit signal to all workers.
Worker:
Worker loop starts with synchronous MPI receive comand waiting for input from controller. The received input is processed according to the algorithm provided and the results are sent back to the controller. When exit signal (message with specific tag) is received, worker breaks the loop and terminates.
MPI primitives:
Synchronous blocking wild-card MPI receive is used on the controller in order to avoid useless iterating over all workers with asynchronous receive scanning for results. Blocking receive ensures synchronisation with all workers. MPI gather is not used since workers can finish their tasks with different speed, and waiting for input from all workers together on every iteration of the controller loop would have decreased the overall performance, leaving alone the fact that implementation would have not been so straight forward in case of using MPI gather.
Semi-synchronous (but blocking) send is used since neither the farmer nor the workers need to wait till data is received by other side as long as they know it is going to be received eventually, which is ensured by the blocking property of the mechanism.
double farmer(int numprocs) { int n_workers = numprocs - 1; //total number of idle workers int idle=0; //iterator over the list of workers int iter=0; // list of workers. values: 1 if waiting for input, 0 if computing int* proc_waiting = (int*) malloc(sizeof(int)*(n_workers)); double result = 0; MPI_Status status; double* temp = (double*) malloc(sizeof(double)*2);
stack* bag; bag = new_stack();
temp[0] = A; temp[1] = B; push(temp,bag);
int i=0; for (i;i<n_workers;i++){ proc_waiting[i] = 0; }
//Controller loop do{ // Receiving data from workers. MPI_Recv(temp, 2, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status); idle++; proc_waiting[status.MPI_SOURCE - 1] = 1; if (status.MPI_TAG == 1){ result += temp[0]; }else{ push(temp,bag); int tag = status.MPI_TAG; int source = status.MPI_SOURCE; MPI_Recv(temp, 2, MPI_DOUBLE, source, tag, MPI_COMM_WORLD, &status); push(temp,bag); } //Iterating over proc_waiting list to find and idle worker and send a task to it while(!isEmpty(bag) && idle>0){ if (proc_waiting[iter]){ MPI_Send(pop(bag),2,MPI_DOUBLE,iter+1,0,MPI_COMM_WORLD); proc_waiting[iter]=0; idle--; tasks_per_process[iter+1]++; } iter = (iter+1)%n_workers; } }while(!(isEmpty(bag) && idle==n_workers));
i = 0; // Sending exit signal to all workers for (i;i<n_workers;i++){ temp[0] = 0; temp[1] = 0; MPI_Send(temp,2,MPI_DOUBLE,i+1,1,MPI_COMM_WORLD); }
return result; }
// mypid argument is not used in this implementation void worker(int mypid) { MPI_Status status; double buf[]={0,0}; // Sending 0 to the farmer to indicate that the worker has initialized and // ready to receive a task. MPI_Send(buf,2,MPI_DOUBLE,0,1,MPI_COMM_WORLD); // Worker loop while (1){ MPI_Recv(buf,2,MPI_DOUBLE,0,MPI_ANY_TAG,MPI_COMM_WORLD,&status); int tag = status.MPI_TAG; // exit signal if (tag == 1) break; double left = buf[0]; double right = buf[1]; double lrarea = (F(left) + F(right)) * (right - left) / 2; double fleft = F(left); double fright = F(right);
int isEmpty(stack * s){ if (s==NULL || s->top == NULL){ return 1; }else{ return 0; } }
// Simple type for stack of doubles
// creating a new stack stack * new_stack() { stack *n;
n = (stack *) malloc (sizeof(stack));
n->top = NULL;
return n; }
// cleaning up after use void free_stack(stack *s) { free(s); }
// Push data to stack s, data has to be an array of 2 doubles void push (double *data, stack *s) { stack_node *n; n = (stack_node *) malloc (sizeof(stack_node)); n->data[0] = data[0]; n->data[1] = data[1];
Ok, I got lissom working and giving correct results. It's time now to see how fast it is. Damn, I think over the course of working on this project I probably learned as much as I did over the past 4 years at edinburgh uni...
Finally, I made Lissom work. Can't say this is it yet, since I haven't tested for correctness of results, but if you imagine that this project is a long journey by train to a city far away, then the view from the window has now changed from seemingly infinite countryside landscapes to urban scenery. Ok, metaphors aside, a quick summary of what's left to do:
- Modify lissom script in order to be able to check the correctness of results produced and fix whatever needs to be fixed (1 day) - Modify plotting script so that it draws graphs of multi-CPU relative to single-CPU performance improvements (a few hours at most) - Run series of tests on shared- distributed-memory machines (might take a few days, but I don't need to interact with the tests in any way, just need to start them and collect results when finished) - Write report (... 5-8 pages per day? ...)
def __set_min_matrix_radius(self,min_matrix_radius): pmi.call(self.pmiobj,'set_min_matrix_radius',min_matrix_radius) def __get_min_matrix_radius(self): ## !KKUFOALERT: replacing the following invoke with call, localcall or simply removing the whole property ## from MPI_CFProjection (as it does not need to be here really) and trying to simulate tiny_mpi with cortex_density ## between 25 and 40 sometimes (!) freezes the simulation on distributing connection fields with set_flatcfs (WTF?!) return pmi.invoke(self.pmiobj,'get_min_matrix_radius')[0] def __del_min_matrix_radius(self): pmi.call(self.pmiobj,'set_min_matrix_radius',None) min_matrix_radius = property(__get_min_matrix_radius,__set_min_matrix_radius, __del_min_matrix_radius)
Hey! I tweaked my code a little bit to improve the speed of communications by replacing PMI invokes in some places with MPI scatter. The thing is, PMI invokes distribute data using MPI bcast, which is perfectly fine if all workers need to work with the same sets of data. However, if each worker has to process only a chunk of this data that corresponds to the node, it would be more logical to use MPI scatter that sends chunks of data to nodes (instead of the full data set, with each node taking its chunk from it) - that's what I thought. PMI does not have mechanisms for scattering, so I had to reimplement my invocation as a call, that breaks data down into chunks on node 0 and scatters it around the rest of nodes. Wondering if my tweak was at all useful, I constructed a simple test that checks how long it would take to send the same chunk of data using MPI bcast and scatter (I was only interested in communication time):
from mpi4py import MPI from numpy.random import beta from time import time
comm = MPI.COMM_WORLD size = comm.Get_size() rank = comm.Get_rank()
big_array = [beta(a=1,b=1000,size=(500,500)) for x in range(size)]
if rank == 0: data = big_array else: data = None
scatter_time = 0 bcast_time = 0
tries = 3
for j in range(tries): t0 = time() for i in range(100): x = comm.scatter(data, root=0) scatter_time += time() - t0
t0 = time() for i in range(100): x = comm.bcast(data, root=0) bcast_time += time() - t0
I picked this course by mistake. At the beginning of semester 2 I decided to take every module that was somehow concerned with high performance computing and thus signed up for the two that were on offer: Parallel Programming Languages and Systems (PPLS) and Parallel Architectures (PA). While PPLS was a hit on target for obvious reasons, with PA I failed to realise that this was a hardware course (probably haven't processed what was written on the module description page well enough). I am a software person. By saying this I mean that I have a more or less adequate low to high - level programming experience, familiarity with various software engineering techniques and a bit of working experience as a software developer. However, my general idea of what's happening inside the laptop I'm using used to be not very far from a mythical wonderland where fairies and elves are flying around in a silicon forest and carry data (probably in some sort of boxes or even bags) in and out of the Castle of Princess Ursula where it gets sorted out somehow. Well, I got 2 CPU's actually, so there must be two Ursula's there, dunno how they get along together... A bit exaggerating, but you get the point.
Anyway, while finishing the last assignment for the course I realised that I actually quite liked it. I mean, I managed to get 84% for the first assignment where I had to write about some weird multi-processing architecture designed in 1999 as a pure concept and never implemented since then. I had a very rough idea of what I was writing about (Wikipedia did neither!) - and there you go, 6% or 7% higher than class average! Then I managed to complete the second assignment of implementing a simulator of a multi-cpu cache in just one night (I had a really tough week and simply forgot about the deadline and almost missed it. Ah, shit happens...), give that prior to this assignment I had no idea what cache was all about! And got 65% for it - maybe not that impressive, but given the circumstances... However, the point of this is not how cool I am or something, the point is: I'm glad I made that mistake and took the course. The fairies finally flew away for good. With elves and Ursula's.
Dr. Randal Smith (big name in Sun Microsystems, the inventor of SunSPOTs and probably not only), while giving a presentation at our uni last year said: "... I'm a software person, and all I can say about this piece of technology as a software person is, electrons are involved...". It's still one of my favourite quotes, but I can't refer this to myself anymore. I mean, I know how cache works, yeah!=))
