Had a few free hours this morning, and decided to write a bit more code for my simulation program. Very well, now the first parallel version that uses the same algorithms that's currently implemented in topographica, is finally working! So, two thirds of the task are done=))) the only thing left to work on is new parallel algorithm, hopefully faster than the current.
At the moment stuck with creating a derived datatype for the connection field matrix. The problem is somewhere between three lines of code and I'm out of ideas. However, I'm having a Message Passing Programming practical today - think I'll ask there for help.
Today I'm going to try out my ideas for hiding communications on a piece of C code that's doing something similar to what's happening inside the Topographica's DotProduct_opt. Now, because this piece of code obviously does not exist yet, I'll have to write it, so the plan is:
Write the sequential version
Then modify this to run in parallel, the way it does in the current parallel topographica
Final stage would be trying to modify it further so that communication happens in parallel with computations and see if it gives any speed-ups
...The sequential version is now working. Ew, salt and vinegar crisps are awful... God I hate vinegar! Anyway, on to the next part.
FUUUUUUUUUUUUUUUUUUUUUUUUUUUCK!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! accidentally destroyed half the code I've written.... so stupid=( rewriting, what else...
ok, finally recovered, so now back to where I've been. Probably I'll call it a night now, it's a bit late. Tomorrow - working on the MPI bit.
More or less done with optimising the activation function in python, just realised that I should probably roll back to the original code now. While I could potentially get some speed-ups, the whole approach seems to be a bit pointless since I will ANYWAY have to take the optimisation down to the C level and undo all the changes I will have done in Python. It sort of made sense when I though it was going to be a quick and easy hack, but turned out not to be so and any further work in that direction would be a waste of time.
Next step will be writing a program in C that works roughly the same way as DotProduch_opt, and then trying out several approaches at optimising it to pick the best one.
Working on implementing the idea I had earlier today. Doesn't seem to be hard so far. I got topographica with all the MPI stuff running on my laptop as well, no clue why it's been so difficult the first time... Hope to get some results today, though I don't know how long it's going to take me... Anyway, I've still got 3/4 of a Bud crate.
Ok, braking the connection field array into parts and computing each part works now. Time to work on sending/receiving. 1/2 the crate...
well, it's not as easy as it seems, maybe I could finish it later. Figured what my problem was tho.... almost finished the crate
Got the activation to work with blocking send/receive. Time to improve it to use the non-blocking versions.
Tried testing performance with the single-run test_script-short test file, everything ok. Now need to find out where the communication happens.
wondering if I should move all the MPI code to a new directory topo/mpi. Good idea, but I'll save this for later.
Yes, it is a very good idea, because I already spent about 15 mins just trying to find bits of my parallel code, especially the PMI part. Where the hell did I put it?!!!......
ok, found everything now.
Will have to add new response functions, capable of exchanging data between nodes iteratively. This might mean that I'll have to do some MPI programming in C for the inline functions.
So, to summarise what I found out today: parallelising communication with dot-product computations at the C-level (i.e. modifying the Dot_Product_opt) is probably the way to go, however, this implies a lot of new changes, so I'll put that away for now. I'll need to speak to Chris or Jim about it first. What I could actually try doing for now is modify the activation function so that it breaks it's own chunk of the activation job into smaller pieces and applies the response_fn to each one, on every iteration sending the computed activity to other nodes in a non-blocking way. This is the first step in hiding away the communications, and it will hopefully help reducing the overhead of gathering the activity matrix on master node after each activation. Then we'll see how it goes. Hope it all makes at least some minimal amount of sense...
As far as I remember, communications were the biggest slowdown factor with the distributed model. Perhaps if I could hide them away, I could get a substantial increase in performance. At the moment I'm thinking of the following trick: instead of gathering the entire projection matrix after each activation (which is the main cause of the communication bottle-neck), why don't I make each node send smaller chunks while processing the rest of the matrix. What I mean by this is: say each time process A computes X number of neuron activities, where X is some small portion of A's chunk, it will broadcast (in a non-blocking way) these X entries to all other nodes, until it reaches the end of its chunk. Hopefully, this will help to hide away the communication overheads. My main concerns, however, are:
- How do we figure out what this X should be? Pick it too small, and broadcast mechanism, whatever it is, won't be able to control the flood. Pick it too big, and the whole point of this optimisation tweak is lost.
- Synchronising this will be a massive pain in the neck: there's no such thing as non-blocking broadcast in MPI-2. As for now, I'll go with bSend or iSend, which, tbh, looks a bit more of a temporary hack to me rather than a long-term solution. I'll come back to it later, probably will have to research this a bit further.
- Lastly, we'll most definitely have to use OpenMP for doing computations on each shared-memory system, since this will help to reduce the number of MPI nodes, and as a result decrease the amount of communication: simply speaking, there will be less packages flying around, because I can really see how the whole thing could get very messy otherwise.
At the moment I got Topographica up and running again, MPI code seems to be working fine and I aslo created a new repository: http://gitorious.org/parallel-topographica/parallel-topographica .
Starting to get an urge to do some coding. Hm... what should I start with, interesting=) Probably do some work on Topographica again, especially since I got quite a few new ideas.
- 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 refined 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!=))
There are a few ways for the Controller to communicate with the Workers. In most cases I'm using pmi.invoke. Invoke call takes function name to be executed on workers (as a string value), a list of arguments to be passed to workers and returns a list of results returned from workers (if any). For example, let's say we have a function "get_data" of "SomeClass" on workers that takes "data" as an argument and returns "modified_data". First we need to initialise SomeClass on all workers and get reference to it. pmi.create does the job. Let's say the reference is stored in "pmiobj" object. Now, calling
pmi.invoke(pmiobj,"get_data",data)
Would return (assume we have 4 nodes, and all run the same code) a list:
How that works: pmi.invoke simply passes data to the worker on node 0, which saves it (by reference!) and broadcasts to all other nodes, which involves pickling it on node 0 and unpickling on receiving by other nodes. If this data object does not implement Single Segment Buffer Interface (like numpy arrays, for instance) it sent, generally, 10 times slower than if it does. Ok, nothing new, I guess
However, the communication can be organised in different, slightly lower-level, fashion. Let's say, now instead of passing data and processing it on the nodes, we want workers to store the reference to "some_object" and then use data stored in that object later in some computations. This could be done by issuing pmi.invoke again. However, what if data in "some_object" gets modified on controller? All workers apart from the one on node 0 would have out-dated (which in many situations would mean useless) copy of "some_object", and would need to somehow get a new one, e.g. from the worker on node 0. In order to avoid redundant pickling/unpickling and network overheads we can use pmi localcall (or was it local_call?), that just passes data to the Worker on node 0. Node 0 can now store the reference and when needed by all nodes, broadcast data from the object using MPI bcast. This would ensure that all nodes will work with up to date copy of the data. Simple trick, but very useful, especially since sometimes pickling objects can have unpredictable results,(like for example, in my case, bloody dest..... spent hours debugging this!).
A side note on PMI terminology I'm using: "Serial" can be also referred to as "the Controller" as well a "Parallel" or "Nodes" can be referred to as "Workers". Tonight while debugging my code I had to, for the first time, have a look at PMI source code. That actually helped... Ok, sleep now. I have to take care of other uni work, so probably I'll resume the Project not sooner than Friday. I really wish I will finish it in one go then...
In this example I will demonstrate how changing numpy printoptions threshold affects the communication speed between Serial and Parallel. Honestly, I have no idea what printoptions threshold has to do with MPI or PMI communication speed and I am really tempted to investigate this issue further, but sadly don't have time for this at the moment. Anyway, output first (using mpirun with 2 nodes), then source code:
numpy_printoptions threshold set to 0 sending 200x200 numpy matrix with random elements Sending x 10 times took 1.40654301643 numpy_printoptions threshold set to 201*201 (greater than the size of communicated matrix) sending 200x200 numpy matrix with random elements Sending x 10 times took 12.252920866
The numbers speak for themselves. The source code:
if __name__ != 'pmi': ################################################## ## Serial code ################################################## import pmi from time import time from numpy.random import beta from numpy import array, set_printoptions
pmi.setup()
pmi.execfile_(__file__)
x = beta(a=1,b=1000,size=(200,200))
# create a frontend class class Hello(object): def __init__(self, name): self.pmiobj = pmi.create('HelloLocal', name) def __call__(self): return pmi.invoke(self.pmiobj, '__call__') def time_send(self): t0 = time() for i in range(0,10): pmi.invoke(self.pmiobj,'recv',x) print "Sending x 10 times took",(time() - t0) def time_recv(self): t0 = time() for i in range(0,10): pmi.invoke(self.pmiobj,'send') print "Receiving x 10 times took",(time() - t0)
# use the class hello = Hello('Konstantin')
print "numpy_printoptions threshold set to 0" set_printoptions(threshold=0) print "sending 200x200 numpy matrix with random elements" hello.time_send() set_printoptions(threshold=201*201) print "numpy_printoptions threshold set to 201*201 (greater than the size of communicated matrix)" print "sending 200x200 numpy matrix with random elements" hello.time_send()
else: ################################################## ## Parallel code ################################################## from mpi4py import MPI from numpy.random import beta
y = beta(a=1,b=1000,size=(100,100))
class HelloLocal(object): def __init__(self, name): self.name = name def __call__(self): return 'Hello %s, this is MPI task %d!' % (self.name, MPI.COMM_WORLD.rank) def recv(self,x): "do nothing"
In this example I will be testing if PMI passes objects from Serial into Parallel by reference or by value. First I will give the output of my little test script, and then the source code.
############ TEST 1: PASSING FROM SERIAL INTO PARALLEL ############ Testing PMI pass by reference: an array is created in Serial and then passed into Parallel array before mpi.invoke: [1 2 3] If the array was passed by _value_ then after modifying it in Parallel it should remain unchanged in Serial array after mpi.invoke: [1000 2 3] As you can see, the array has been changed in serial as well, which proves that objects are passed by reference in pmi.
############ TEST 2: PASSING FROM PARALLEL INTO SERIAL ############ Now we'll check if modifying the same array in Serial will make it change in Parallel The array: [1000 2 3] Changed the first element to be 999: [999 2 3] Now calling pmi.invoke that will make Parallel print out the array (it has been stored as a class atribute during the previous experiment) Node 0 array: [999 2 3] Node 1 array: [1000 2 3] Oops! It seems like objects are passed by reference only between Serial and node 0 of parallel. Let's verify that: the next pmi.invoke will call a method in Parallel that will attempt to set the second array element to be 1000 on node 0 (and only on node 0) and the third element to be 1000 on node 1 (and node 1 only)
############ TEST 3: CHANGING DIFFERENT ELEMENTS ON DIFFERENT NODES ############ Once again, array before pmi.invoke: [999 2 3] Array after pmi.invoke: [ 999 1000 3]
Ok, that proves the case that objects are indeed passed by reference in pmi between Serial and Parallel,but only for the first node. For all other nodes they are passed by reference. This is a bit tricky, but might be useful.
and the source code:
if __name__ != 'pmi': ################################################## ## Serial code ################################################## import pmi from time import time from numpy import array
pmi.setup()
pmi.execfile_(__file__)
# create a frontend class class Hello(object): def __init__(self, name): self.pmiobj = pmi.create('HelloLocal', name) def __call__(self): return pmi.invoke(self.pmiobj, '__call__') def check_pass_by_ref1(self,arr): pmi.invoke(self.pmiobj,"modify_array",arr) def check_pass_by_ref2(self): pmi.invoke(self.pmiobj,"print_array") def check_pass_by_ref3(self): pmi.invoke(self.pmiobj,"modify_array2")
# use the class hello = Hello('Konstantin')
a = array([1,2,3])
print "############ TEST 1: PASSING FROM SERIAL INTO PARALLEL ############" print "Testing PMI pass by reference: an array is created in Serial and then passed into Parallel" print "array before mpi.invoke:",a
hello.check_pass_by_ref1(a)
print "If the array was passed by _value_ then after modifying it in Parallel it should remain" print "unchanged in Serial" print "array after mpi.invoke:",a print "As you can see, the array has been changed in serial as well, which proves that objects" print "are passed by reference in pmi."
# test 2
print "" print "############ TEST 2: PASSING FROM PARALLEL INTO SERIAL ############" print "Now we'll check if modifying the same array in Serial will make it change in Parallel" print "The array:", a
a[0] = 999
print "Changed the first element to be 999:", a print "Now calling pmi.invoke that will make Parallel print out the array (it has been stored" print "as a class atribute during the previous experiment)"
hello.check_pass_by_ref2()
print "Oops! It seems like objects are passed by reference only between Serial and node 0 of parallel." print "Let's verify that: the next pmi.invoke will call a method in Parallel that will attempt to set" print "the second array element to be 1000 on node 0 (and only on node 0) and the third element to be" print " 1000 on node 1 (and node 1 only)"
print "" print "############ TEST 3: CHANGING DIFFERENT ELEMENTS ON DIFFERENT NODES ############" print "Once again, array before pmi.invoke:",a hello.check_pass_by_ref3() print "Array after pmi.invoke:",a print "" print "Ok, that proves the case that objects are indeed passed by reference in pmi between Serial" print "and Parallel,but only for the first node. For all other nodes they are passed by reference." print " This is a bit tricky, but might be useful." else: ################################################## ## Parallel code ################################################## from mpi4py import MPI from numpy.random import beta
y = beta(a=1,b=1000,size=(210,210))
class HelloLocal(object): def __init__(self, name): self.name = name self.rank = MPI.COMM_WORLD.Get_rank() def __call__(self): return 'Hello %s, this is MPI task %d!' % (self.name, MPI.COMM_WORLD.rank) def modify_array(self,arr): self.arr = arr self.arr[0] = 1000 def modify_array2(self): if self.rank== 0: self.arr[1] = 1000 if self.rank== 1: self.arr[2] = 1000 def print_array(self): print "Node",self.rank,"array:",self.arr
I think, the reason for this is very simple: PMI seems to be broadcasting the data to all nodes, when pmi.invoke is issued, (perhaps, using MPI bcast()) from node 0 and keep the broadcast data on node 0 as it is, i.e. reference to the object passed into pmi.invoke, and not the copy of this object.
Today I fixed learning. This means that tiny_mpi.ty simulation can now work in parallel mode with all of its original features turned back on. At this point, the only thing left for me to do, apart from writing the final report (40-60 pages btw.), is running lissom_oo_or.ty simulation in parallel. If When I manage to do that, it will be the successful end of this project. So far, we managed to achieve speed-ups on any cortex/retina densities, on any number of CPU's (although the performance improvement scales down with each new CPU added into the system) and with learning and optimisation turned on, - not so bad. We might have expected greater speed-ups and the system to be scalable (i.e. adding one more CPU increases performance accordingly), but all in all, topographica is now twice (at least!) faster then it used to be and there's a whole lot of potential for improvement. Also, we have identified the new problems, bottlenecks and fixed a few old bugs along the way. Personally, I have developed a lot of understanding what parallel programming's all about, so it's all good.
Anyway, I haven't finished the project yet and there's still hell of a lot of work to do, so let's move on.
At the moment I am applying to EPCC (Edinburgh Parallel Computing Centre) for Masters of High Performance Computing degree (God, the application form is taking me ages to complete!). I just finished writing personal statement and decided that it would be nice to post it in my blog. I apoligise in advance for the manner in which it is written and for the things like "would be an honour for me", but this is a personal statement, after all, it has to be like that. Also, I do believe in every word I have written in this letter:
My interest in High Performance Computing comes mainly from the Honours Project that I have chosen to work on during my fourth year of studies at Edinburgh University. The goal of the Project is to optimise a neural map simulator Topographica by distributing some of its heaviest computations using Message Passing Interface (MPI). Topographica is an open-source software that has been continually developed over the past ten years and is used by many scientists around the world. Distributing its computations would allow experimenting with neural networks of sizes and densities much greater than it has been possible to simulate ever before. Hopefully, this would give scientists a better understanding of how human brain functions.
Due to the complexity of Topographica, this Project was probably the most challenging programming exercise that I have ever undertaken. More importantly, it gave me an opportunity to realise the significance of Parallel Computing for both industry and science. Nowadays software requirements for processing power continue to increase while processors are not getting faster even close to as much as they used to in the past century and early 2000's. Multi-core and multi-processor machines already dominate the hardware market. However, such systems require a new approach in software development, in many aspects very different from the paradigms we are accustomed to in modern-world IT. This makes me believe that Parallel Programming is the future of Computer Science and IT Consultancy and this is why it is so important for me to acquire the skills and experience in that area. I hope that having such skills and experience will allow me to work on world-class scientific projects in areas ranging from Neuroscience to Aerospace and in my own way contribute to the bank of human knowledge.
In shared-memory models the most interesting thing for me is algorithms. I find the complexity of synchronisation problems to be very challenging and exciting. Apart from that, I find the power of modern supercomputers with hundreds or even thousands of processors to be nothing less than fascinating and having access to such machines would be an honour for me. However, not all scientific institutions around the world can afford to have such machines and this is where, I believe, the importance of distributed-memory models and Cluster Computing comes from. Being able to make a network of machines function as a “DIY” supercomputer would allow computer scientists around the world to carry out research that otherwise would have been impossible due to budget limitations. This is why it is so important for me to have knowledge and skills allowing to develop software that can take full advantage of Parallel Computing and this is what makes me believe that specialising in the area of HPC is a perfectly right choice for me.
Good Lord, so much has happened over the last week and I have so much to write about, but unfortunately can't afford to spend my time on this at the moment... I will update this blog with detailed results of experiments and my thoughts and ideas in regards to the project as soon as I get a chance, promise. Anyway, a short summary of what has happened:
- MPI is faster than serial when learning and optimisation are turned off
- In general, MPI seems to be slower than Serial when optimised mode is on, however, at some very high densities of cortex and retina MPI seems to behave quite differently and seems to outperform Serial. Also, in some cases processing bigger sets of data takes less time than processing smaller sets of data. I find this to be very bizarre, and I really want to find out what's causing this.
- I've built a python script that plots performance graphs and will upload pictures from it soon
- Turning learning on breaks simulations and that's my next thing to fix. Working on it at the moment.
- We also had some ideas in regards to distributing initialisation stage, i.e. generation of Connection Fields, so that all CFs are generated on nodes instead of serially and then distributed. This will give a great performance boost, however this is a bit tricky.
- MPI performance seems to be greatly affected by necessity to communicate the Activity matrix into Serial code on each iteration, and we still have no idea what to do with this.
- At the last skype meeting Jim, my supervisor, has pointed out that successful run of parallel topographica using lissom_oo_or.ty (standard simulation script) and not simplified and trimmed down tiny.ty I'm suing at the moment will signify the end of the project. This seems to be the light in the end of the tunnel. However, at the moment I can't estimate how far I am from it...
Benchmarking. The time has come. This evening I will test the performance of my implementation and update this post as I get results with each test. The machine I'm going to use is a shared-memory multicore with 8 Intel(R) Xeon(R) X5450 CPU's, 3GHz, 6MB cache, 4 cores each. This little monster is called Jupiter III and together with Jupiters I and II forms a 3-piece cluster.
So, at the moment I am testing the test-bench and soon will be ready to launch. Don't switch the channel=)
[Update]
Phew! Preparation for testing took me way more time than I expected... However, look at this beautiful script! (./topographica line's were created by Chris, everything else is mine)
#!/bin/bash # Argument = -c cortex density -n number of cpu's -i number of iterations
usage() { cat << EOF usage: $0 options
This script tests the performance of parallel Topographica implementation
OPTIONS: -c cortex density -n number of cpu's -i number of iterations EOF }
C_DENSITY= N_CPU= ITERATIONS= while getopts “c:n:i:” OPTION do case $OPTION in c) C_DENSITY=$OPTARG ;; n) N_CPU=$OPTARG ;; i) ITERATIONS=$OPTARG ;; ?) usage exit ;; esac done
if [[ -z $C_DENSITY ]] || [[ -z $N_CPU ]] || [[ -z $ITERATIONS ]] then usage exit 1 fi
# Mock run to make sure that all C code has compiled before proceeding to actual testing ./topographica -p cortex_density=1 examples/tiny.ty -c "topo.sim.run(1)"
# Serial run ./topographica -p cortex_density=$C_DENSITY examples/tiny.ty -c "import timeit; print 'Serial:', timeit.Timer('topo.sim.run($ITERATIONS)','import topo').timeit(number=1)" -c 'import pickle; pickle.dump(topo.sim["V1"].activity,open("results.pickle","wb"))'
Now to do testing you just need to run this script with 3 parameters: -c for cortex density, i for the number of iterations and -n for the number of processes to use. I already tried this out on my machine and so far the results are not that inspiring:
On the bright side, the results computed in MPI are still correct. Now I'm going to try this out on Jupiter.
[update 2] Ok, found a bug in the script: it ran tiny.ty instead of tiny_mpi.ty. How the hell did it get correct results? Anyway, new output (my laptop): ./test_script -c 25 -n 2 -i 100 Not using MPI Not using MPI Serial: 1.02229189873 Using MPI NODE 0 initialized NODE 1 initialized MPI: 3.37540221214
[update 3] Ok, I need to make a break now. For some reason Jupiter won't let me run the tests, so I'll have to deal with that somehow. But later.
Hi, Mr Internet. Yesterday I was finally able to run my first experiment with PMI implementation! This means that now I can start benchmarking and see where I really am. At the moment I am working on setting up my first cluster network - something I started over Christmas Holidays, spent a couple of days working on, hadn't managed to make it and decided to shift my attention to more important things. However, now is the time. I want to see for myself that it is possible to make a scalable network of arbitrary PC's (not specialized hardware) with different architectures to work as a single unit doing computations in parallel and then gathering results on a single node. Currently setting up ssh-agent to allow password-less access from my computer to DICE (university) network, which is absolutely necessary for running parallel jobs between these machines, is giving me the most headache. Anyway, should be manageable.
Also, since it's finally woring now, here's the (part of) actual implementation of MPI_CFProjection:
def __set_activity(self,activity): #flattening activity matrix self.activity_shape = activity.shape #reshaping into one-dimensional matrix (2d matrix that has only one row) activity = activity.reshape(1,self.activity_shape[0] * self.activity_shape[1]) pmi.invoke(self.pmiobj,'_set_activity',list(activity[0])) def __get_activity(self): activity_list = pmi.invoke(self.pmiobj,'_get_activity') activity = [] for activity_row in activity_list: activity.extend(activity_row) activity = numpy.array(activity) return activity.reshape(self.activity_shape[0],self.activity_shape[1]) def __del_activity(self): pmi.invoke(self.pmiobj,'_set_activity', None) activity = property(__get_activity,__set_activity,__del_activity)
This is the serial part of implementation. This means that this code is run as if mpirun command wasn't issued. Members of the super-class (CFProjection) that have to be distributed in order to do computations are implemented as properties. For each property, get and set methods call pmi.invoke which, in turn, calls the MPI method specified as the parameter to mpi.invoke call. If something else is passed as parameter to pmi.invoke, it goes through to parallel method as a parameter, one copy per node. Thus, if you do pmi.invoke(self.pmiobj, "some_method", "abc"), then some_method will be called in the parallel mode and an instance of "abc" string will be passed to it. Obviously, some_method(self, some_string) has to exist somewhere.
Here's the (part of) parallel code that is being run every time pmi.invoke call occurs:
Hi! It's been a really busy week, had to submit 3 assignments. Now some virus has struck me down, I asked for extension for this week's assignment and have at least a couple of days to work on the project. So, the status is: I am really close to my first (real) distributed simulation with pmi. There is one thing that stands between me and this target: tight coupling of one of the components I need to distribute. Hopefully, as soon as I deal with it, I'll be able to topo.sim.run(1) *hwang! rock'n'roll*... ah, whatever.
A quick update. PMI has proven to be very useful, the code now looks a lot better. I've rewritten the MPI_CFProjection so that now it is entirely serial and all the content of this class that has to be parallelized now has properties that link to the MPI_CFProjection_node parallel implementation of CFProjection. This is what I'm doing my best to debug at the moment. I'll post some code as soon as it starts working.
# use the class hello = Hello('Olaf') print('\n'.join(hello()))
a = hello.func()
print "Wow, PMI rules!"
else: ################################################## ## Parallel code ################################################## from mpi4py import MPI
class HelloLocal(object): def __init__(self, name): self.name = name def __call__(self): return 'Hello %s, this is MPI task %d!' % (self.name, MPI.COMM_WORLD.rank) def func(self): if MPI.COMM_WORLD.rank == 0: return raw_input() else: return "Slave node"
Things to pay attention to:
1) First, name object is passed as a parameter from the serial code into the parallel class, and than a separate call to parallel code uses this object. This means that all objects persist between separate calls to parallel code.
2) While calling the func function, only the master node (rank 0) is waiting for user input. However, as opposed to using just mpi, all other nodes are waiting for the master node to complete operation before all nodes join again in the serial code.
3) Notice that "Wow, PMI rules!" is printed out only once, no matter how many nodes you create with mpirun. Now comment out everything else, leave just this print statement, and mpirun this file once again, creating say 4 processes. What will you get?
I just thought I haven't posted for a while, so this is a quick update on what's currently going on.
First, this Monday my laptop has refused to switch on and now seems to be pretty much dead, but hopefully this is just coma and anyways my warranty hasn't expired yet, so this Wednesday it's going to be picked up and hopefully returned in a week, fully functional. Not a big issue, just a bit annoying.
Second, some time ago Chris suggested I have a look at PMI - "a pure python module that allows libraries to provide functions that are parallelized using MPI but that can nonetheless be called from serial Python scripts." I'll write about it if I find the way to utilize this thing, but at the moment it seems to be quite a handy module. However, here's what I want to get rid of in my implementation by using PMI:
''' Created on 20 Jan 2010
@author: megatelevizor ''' from mpi4py import MPI import time
from MPI_client import client
""" Vocabulary
a = assign flatcfs c = compute activity i = Initialize i_a = Get a copy of input_activity matrix k = Kill p = Print r = Return activity s = Set strength
""" Set a copy of input activity matrix """ if command_prefix == "i_a": input_activity = self.comm.recv(source=0,tag=1) self.client_node.set_input_activity(input_activity)
""" Compute activity """ if command_prefix == "c": self.client_node.compute_activity()
""" Set strength""" if command_prefix == "s": strength = self.comm.recv(source=0,tag=1) self.client_node.set_strength(strength)
""" TAGS: 0 - command prefix 1 - input data 2 - output data """
This class allows master node to communicate with slave nodes. Potentially, PMI will eliminate all need to have such class.
[update] No, it will not. What PMI allows you to do is hide the parallel implementation, so that your programme looks like it is a serial (single-processed) code. However, it does not enhance the underlying MPI implementation in any way, so it still needs to be implicit. What does that mean in my case? Well, if I was to introduce PMI into my code as it is now, then there would have been a class - something like CFProjection_proxy, linked to the real (i.e. parallel) implementation of CFProjection_slave using PMI mechanism. CFPRojection_slave's would still have to communicate with the master node CFProjection using MPI_io (the code listed above) because PMI does not influence or change the way different nodes of an MPI application communicate. In simple words, it just makes your code look like it's not using MPI, when it actually does.
Well, so what's so bad about it? Why not use PMI anyway, since from all that it looks like it's actually doing a good job? Actually, there's nothing wrong with it and it might be a good idea to use it. No, it won't add any new functionality, but it will definitelly make things look a lot nicer. However, I prefer dealing with things one at a time, and also finishing what I started, so I think I'll concentrate on looking for the solution to my initial problem.
To summarise what I'm looking for, have a look at this: Java RMI . Without going into too much detail, RMI is a mechanism that allows Java objects calling each other's methods accross several different JVM's (which can be run on different hosts!) and it stands for Remote Method Invocation. Combined with serialization (something similar to Python's pickling), it allows to hide all explicit communication-level code away from the user, which is great! Thus, my assumption is that if something exists in Java, it has a chance to exist in Python. As simple as that.
Hi! I remember myself promising to post some info on simulations here and perhaps give a better idea of what I'm working on. I think, I'll break it into several parts: the first one will be just a brief intro, a very top-level description. In part II we'll go through the simplest simulation script: tiny.ty. Maybe there will be even part III and IV, but let's not go that far yet.
Also, it's kinda funny that I do that. I don't believe there is a single person in the world who reads this apart from my project supervisor Jim (hopefully) and his Ph.D. student Chris, who's offered me an incredible amount of help on this project, and neither of these two people need to be explained how their code works (especially since I'm not entirely sure my explanation is going to be correct). Anyway, there you go:
This is what a sample Topographica vision system model might look like:
However, the simplest model would only include an image from the retina of an eye connected to V1, the primary visual cortex. In Topographica's implementation an image from the retina is represented as the Input Activity matrix where each element corresponds to the brightness value of each pixel of the image. This matrix could be generated from an image file or randomly, using various distributions.
V1, or the primary visual cortex is a matrix of Connection Fields, which are sets of weights on one input Sheet. Each Connection Field contributes to the activity of one unit on the output sheet, and is normally used as part of a Projection including many other Connection Fields.
The simplest simulation would involve generating an output activity matrix using an input activity matrix and Projection, which can be described as a connection between two sheets of activities: input and output. More specifically, a Projection is a Connection that can produce an Activity matrix when given an input Activity matrix, which will typically be used by the destination Sheet when it computes its activation. [2]
Computing the activation requires applying a response function to the input matrix and weights on the Connection Fields. In our example this function would simply dot-product slices of the input activity matrix with the corresponding weight matrices.
I'm in the process of writing interim report at the moment. God, this is so boring and useless... Well, actually there is one good thing about it: after creating a more or less detailed timeline I realised that for the next couple of months my personal life has seized to exist. Reminds last year's System Design Project a little bit, except for there's no team this time, just me. I loved SDP! There's a saying that wise people learn from their own mistakes and smart people learn from others'. Well, looking from that perspective I'm definitely a wise person. I think that all the f*** ups we had made as a team and each one of us personally during this project were very well worth getting the relatively low mark (~60%). I mean, I don't like getting low marks, but this one brought me a hell of a lot of experience. And it counts for much. Also, I think I've done a good job overall, made new friends and learned a lot, not only from mistakes.
Well, anyway I guess this is where the fun part begins... except for this report isn't much fun, to be honest. Allright, time to get back to it.
Hi. Here's my to-do list for the first implementation. Everything in green is successfully implemented (accomplished), yellow - not required any more or not my responsibility and red - failures.
Modify the topographica script so that it runs with mpiexec
Make a new subclass of CFProjection, e.g. MultiProcessorCFProjection , that will be a plug-in replacement for CFProjection but with parallel capabilities. At first, MultiProcessorCFProjection does exactly the same thing as CFProjection. Make your test .ty file(s) use this class instead of CFProjection.
To make the rest of the steps easier, turn off optimization (remove _opt from the various components in your copy of tiny.ty or other testing file), so that you don't have to worry about the C code's implicit assumptions. Also turn off learning, for now. Also turn off randomness in the initial weights -- use Gaussian patterns for now.
Find out how to instantiate CFProjection
Once that works, change the constructor or other initialization routines, if any, for MultiProcessorCFProjection to divide the work between compute nodes from the start, so that each node only has a copy of the weights relevant to the processing that it has to do. At first, it's probably easiest if the master node (rank 0) has a copy of all the weights, while the others have only a subset. While in this stage, we can debug everything using the usual plotting and analysis mechanisms that use the data from the master node, without having to reimplement all of that just yet. Implementation plan:
[Chris] Refactor CFProjection to replace the list of lists of cfs (self._cfs) with a flattened list (self.flatcfs?), to remove the assumption that all CFs are stored in 2D arrays. This will require changing all of the C code, which is a pain, but it should actually significantly simplify it.
[Chris] Make the MaskedCFiter use self.flatcfs instead of the self.cfs numpy object array. Again, this should simplify the implementation of MaskedCFIter , while making it simpler to keep everything coherent. Actually, we should probably just make the mask update the flatcfs list whenever it changes, so that only unmasked cfs are in the list; otherwise we'll have a problem relating the 2D mask to the 1D list. Up to this point, none of the changes should involve multiprocessing at all.
When multiprocessing (or on a single processor just for testing), split up the list and distribute it to multiple Python processes, so that any operation invoked on MultiprocessorCFProjection is actually distributed across all processors and run on MultiprocessorCFProjection running on that processor, using its local list of cfs. It's probably best just to assert numpy.alltrue(mask) (or whatever is syntactically correct), and not worry about masks at all.
Find out what exactly CFSheet consists of, and how to distribute it
Figure out what topo.sim['V1'].activity is
Read through the rest of the class
Read through the ProjectionSheet class
Figure out how simulation.connect and simulation.run work
Find out how exactly CFProjection is dotproducting cfs
Figure out what input_sheet_slice.submatrix() is
Figure out how CFIter works and how to iterate from and to a certain point in matrix
Create the s/r mechanism
Using the s/r mechanism, send cfs to slave nodes and store. Check if sending and receiving work
Create new activate() method in MPI_CFProjection
At each call broadcast the input_activity matrix
Create according method in MPI_io
Create according method in MPI_client
Add broadcasting to MPI_CFProjection.activate()
Send an appropriate portion of CF's to each node
Create according method in MPI_io
Create according method in MPI_client
Add broadcasting to MPI_CFProjection.activate()
Gather resulting activities into one matrix
Compute activity on all nodes
Create methods (io,client?) for sending back resulting activities
Gather activities into new matrix at MPI_CFProjection
Transform the non-zero part into a list
Send
Gather at master node and put into the activity array, calculating appropriate position for gathered activity values
Ok, now I can officially announce that a few days ago I've reached a major milestone: the slowest ever and probably buggiest as well, however distributed implementation of topographica has been created. I've reached the bottom, and now it's time to start climbing to the top, creating the fastest (and coolest!) version of topographica the World has seen. It's a bit too early to dring champaigne or whatever, but now I have no doubts: this project will be delivered.
Thank you. Questions?
P.S. I have a mountain of work to do on this weekend, so I really don't know whether I will be able to write about simulations, but I haven't forgotten my promise, so stay tuned!
[update] Either I don't get something, or there is indeed something wrong with mpi4py's non-blocking communication. If non-blocking communication of arbitrary objects was introduced in version 1.2 then why does Irecv() still not work the same way as recv(). True, Irecv does not require it's buffer to be tupled with MPI basic types, which allows pickling. But why not get rid of buffer parameter completely? Anyway, I've switched back to blocking communication already, so it's fine.
Wrong results again... Still trying to bring the chaos I've created in my new implementation back to order. I would love to share my ideas and code with you, but at the moment that would be absolutely pointless, and, besides, embarrassing. Going to bed now, maybe tomorrow will be more productive. Think my next informative post will be about looking closer at Topographica's simulations - .ty files, simulation.py etc. In my opinion, it's quite a sophisticated mechanism, that I still don't have a full understanding of, but regard as a very fine piece of Python code. If everything goes well, I'll write about it on the weekend.
Hi! I'm back in Edinburgh and finally (almost) done dealing with all the small to medium issues that stacked up while I was away, so I think I'll resume writing pretty soon. In the meanwhile, want to say that I'm glad to be back. I've got about 2.5 month left till the deadline for this project and I'm not quite where I want to be with it, so it's about time to accelerate. Watch me ;)
A few years ago, when I was a first year student at Edinburgh University (I'm 4th now) I had an Opportunity. Two of my friends were dealing Georgian Wine in Edinburgh (God, that was some wonderful wine!). They'd been doing it for a few years by then already, and were only selling wine to local shops in bulk. Once I asked one of them why they were not dealing wine on the Web, and he said they didn't know how. Then I offered my help and they agreed to let me into their business. However, I never managed to proceed any further than that. Given the amount of free time I had in my first year at uni, I definitely could've made it work out, if not for profit (It was a really small business and both guys left UK a year ago anyway), then at least for the extremely valuable entrepreneurial experience.
What was my excuse then? Very simple - I thought I had too much time. Have you ever had anything similar?
On the positive side, back then I've learned a lesson: grab it as soon as you see it, which, interestingly enough, in real life works for me more like start caring about it a moment before it's too late. Anyway, whatever gets the job done.
How much of the pre-paid "too much" time do you think you still have in your life? What is your allowance for f*** ups like the one I made 3 years ago?
