Wednesday, 8 June 2011

C vs Python Part 2

Just for record: there was this annoying issue of output matrix converging to permanent values after a very small number of iterations. Because of using one iteration's output as input on the next iteration, as with purely Lateral model, the output matrix values were very quickly becoming fixed, i.e. not changing irrelevant of how many iterations were run afterwards. This is something you don't want to happen if you have two different versions of the same simulation and want to constantly make sure they do the same thing.

To fix that we changed the model to take external inputs (as with the Afferent model) in the form of a moving dot where the value of 1.0 is moving forward each iteration down an array of zeros, and the value of 0.25 - moving backwards. For instance, on the second iteration the 3 by 3 input matrix would be:

0.0 1.0 0.0
0.0 0.0 0.0
0.0 0.25 0.0

on the fourth:

0.0 0.0 0.0
1.0 0.0 0.25
0.0 0.0 0.0

Despite using external inputs, the outputs of previous iterations are still taken into account: after the new activity matrix has been computed using dot-product response function, a hysteresis function is applied to it the following way(TC is a pre-defined time constant):

new_activity = old_activity + (new_activity - old_activity)*TC

Now the output activity values do change after each iteration, which allows us a fairly reliable test case. I also modified the python script to use the optimised (in C) version of dot-product response function, the difference in performance is quite noticeable and gets bigger for bigger matrix sizes (optimised is the first one):

[jupiter3]sXXXXXX: ./c_vs_python.sh -m 100 -i 40
Serial: 7.49869203568
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 100 -i 40
Serial: 10.6240258217


Here are the results of C vs Python runs with different densities and numbers of iterations:


[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 50 -i 50
Serial: 0.688660860062
=== Sequential Time ===
Dot Product:    0.471745
-----------------------
Total:          0.472502
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 50 -i 200
Serial: 2.88342404366
=== Sequential Time ===
Dot Product:    1.92014
-----------------------
Total:          1.92326
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 100 -i 50
Serial: 9.23764681816
=== Sequential Time ===
Dot Product:    7.99027
-----------------------
Total:          7.99298
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 100 -i 100
Serial: 18.9039599895
=== Sequential Time ===
Dot Product:    16.5311
-----------------------
Total:          16.5367
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 150 -i 20
Serial: 17.1115911007
=== Sequential Time ===
Dot Product:    15.8554
-----------------------
Total:          15.858
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 150 -i 40
Serial: 34.4383950233
=== Sequential Time ===
Dot Product:    32.8897
-----------------------
Total:          32.8947
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 170 -i 20
Serial: 28.8636169434
=== Sequential Time ===
Dot Product:    26.6929
-----------------------
Total:          26.6962
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 170 -i 40
Serial: 59.0316541195
=== Sequential Time ===
Dot Product:    53.9727
-----------------------
Total:          53.9792

Notice that the timings for both versions of the code are very close with some amount of relatively static overhead to the Python code which decreases for bigger problem sizes: ~50% for m=50, ~10% for m=170

This looks like a very good starting base for my comparative performance investigations.

Friday, 3 June 2011

C vs Python

First of all, I built a script to automate topographica and c runs and compares the results. Called it c_vs_python.sh. Because of the output stabilising issue (after a certain number of iterations the output matrix ) I can't say I'm 100% sure we've built matching models, but I'd say my confidence in their correctness is about 90%. I ran it with different matrix sizes and numbers of iterations, and the results were matching in each case up to the 8th digit of precision (please correct me if I'm wrong, but isn't that the most you can get out of 32 bit float?).

The script also measures the simulation time (not the total run time!) for both runs. Both serial runs, no MPI involved yet. I don't know, either I'm timing things in a wrong way, or C is indeed that much faster than Python. I truly hope it's the former, because otherwise it will be really hard to compare the parallel performance...

-m - matrix height/width (it's square)
-i - the number of iterations

./c_vs_python.sh -m 120 -i 20
Serial: 9.2740278244
=== Sequential Time ===
Dot Product:    6.20771
Normalise:      0.00157809
Copy Activity:  0.000202894
-----------------------
Total:          6.20958

Toporaphica's time is at the top (9.27s), C - at the bottom (6.2s). 3:2. This ratio doesn't depend on the number of iterations, but seems to decrease slightly for bigger matrices:

[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 170 -i 20
Serial: 31.0484178066
=== Sequential Time ===
Dot Product:    25.015
Normalise:      0.00315642
Copy Activity:  0.000400066
-----------------------
Total:          25.0187
[jupiter3]sXXXXXXX: ./c_vs_python.sh -m 170 -i 40
Serial: 62.2372758389
=== Sequential Time ===
Dot Product:    50.2376
Normalise:      0.00620961
Copy Activity:  0.000801563
-----------------------
Total:          50.2448

Thursday, 2 June 2011

vampirtrace

Oh, and I finally got Vampirtrace to work, only to realise that it isn't that great for tracing Python code: it only seems to display communication routine calls, while it would be very useful for me to track other events, e.g. function calls. Ideally I would want to use some tool that combines profiling based on PC sampling together with MPI tracing and visualises that in some human-friendly way. However, even being able to manually set events or labels and track when and how many times the application triggers them would do. MPE seems to be capable of that, I'll try it out later on.

Week 2

Had my second project meeting yesterday, the plan for this week is to modify my C code and Topographica (hopefully only model scripts) to match in terms of the results produced and internal operations. Starting off with lateral-only square all-to-all projections with no learning. This means that the output matrix is produced by dot-producting an input matrix with a 4-dimensional matrix of weights: first two dimensions contain an M*M grid of neurons (at this stage the size and shape of the grid are the same as that of the input and output matrices).  At this stage, for simplicity each neuron is connected to all other neurons in the sheet, and therefore is represented by an M*M matrix of connection weights.Thus, the grid of neurons together with their weights can be expressed as M*M*M*M matrix, input and output - as M*M. Output is generated by dot-producting input with each neuron's weights, the result of each operation is a single number. After M*M such operations full output matrix is produced. "Lateral" term in our case means that the system uses the output matrix from previous iteration as input to produce new output, as opposed to using some external feed. No learning means that the neuron sheet is static at all points of the simulation.

At some point we want the neuron sheet to be M*K*L*N, i.e. use non-square patterns as well as for inputs and outputs; have limited connection fields (i.e. non-all-to-all), some form of learning and, perhaps, afferent projections (taking inputs from the outside). All that will be used to isolate any performance problems and see the potential for parallel scalability.

Tuesday, 31 May 2011

URL for the experimental MPI code

http://gitorious.org/parallel-topographica/parallel-topographica/trees/master/mpi_experiments

Another update

Spent the last two days trying to do some tracing with vampirtrace so that I don't have to hardcode timers anymore but haven't managed to get anywhere yet. Hoping to get some help on the mpi4py google group, because at the moment I'm in the dead end and out of ideas...

Done some benchmarking with the new python code I built a few days ago and the main conclusion I have to make is: there is indeed a late reduce problem causing communication times to soar, but it's definitely not due to load imbalance. In fact, in my app the load is perfectly balanced. It's just that when several SMP nodes are computing data in parallel, and the computation speed is ultimately limited by memory bandwidth (which is exactly the case with topographica and both my testing apps), the nodes will slow each other down. An interesting observation is that for some reason this behaviour is asymmetric in a way that some nodes will have advantage all the time during one run (next run some other process could dominate). I'm not quite sure why this is happening, and could dig into that if needed, but both my experimental apps demonstrate this behaviour. Here's an output of my C app, after I ran it on the same SMP machine subscribing 6 CPU cores:

=== Sequential Time ===
Dot Product:    14.4671
Normalise:      0.00240159
Copy Activity:  0.000317574
-----------------------
Total:          14.4698

=== Parallel v.1 Time ===
Parallel
dot product:    7.45656
        Broadcasting:   0.00737047
        Dot Product:    4.6
        Gathering:      2.84918
Normalise:      0.00228691
Copy Activity:  0.000320673
-------------------------
Total:          7.45918
######### P0 #########
dot product:    7.45656
        Broadcasting:   0.00737047
        Dot Product:    4.6
        Gathering:      2.84918
######### P1 #########
dot product:    7.45768
        Broadcasting:   0.0125265
        Dot Product:    7.44423
        Gathering:      0.000930071
######### P2 #########
dot product:    7.45741
        Broadcasting:   0.00942779
        Dot Product:    7.44338
        Gathering:      0.00460601
######### P3 #########
dot product:    7.45694
        Broadcasting:   0.00952554
        Dot Product:    7.43021
        Gathering:      0.017205
######### P4 #########
dot product:    7.45681
        Broadcasting:   0.00902605
        Dot Product:    7.05847
        Gathering:      0.389309
######### P5 #########
dot product:    7.45687
        Broadcasting:   0.00766182
        Dot Product:    4.95379
        Gathering:      2.49543


Processes 0 and 5 were dominating this run in that the computations on these have taken 30%-40% less time than on the other nodes. And on the same two processes we can see Gathering taking as much as about 50% of total computation time, a lot more than on the other nodes (notice on process 1 the total Gathering time is less than a millisecond). This is a very clear early reduce pattern that can be observed here: some nodes finish their computations quicker than others and have to wait inside the collective routine. Why isn't Broadcasting taking that much? Well, this collective operation is called very shortly after Gathering, which synchronises all processes, so they are entering the routine roughly at the same time.

My Python app was giving similar results, so there's no point in discussing them. Oh, except one thing: pickling Gather turned out to be not much slower than it's non-pickling version on SMP, and sometimes even faster on MPP!

Now let's see what happens if we run the app with the same problem size on 3 boxes, subscribing only 1 CPU on each machine:


=== Sequential Time ===
Dot Product:    15.3458
Normalise:      0.00247908
Copy Activity:  0.000368357
-----------------------
Total:          15.3486

=== Parallel v.1 Time ===
Parallel
dot product:    5.17208
        Broadcasting:   0.00858498
        Dot Product:    5.15076
        Gathering:      0.0127301
Normalise:      0.00232387
Copy Activity:  0.000324488
-------------------------
Total:          5.17474
######### P0 #########
dot product:    5.17208
        Broadcasting:   0.00858498
        Dot Product:    5.15076
        Gathering:      0.0127301
######### P1 #########
dot product:    5.16623
        Broadcasting:   0.158481
        Dot Product:    5.00647
        Gathering:      0.00128365
######### P2 #########
dot product:    5.16298
        Broadcasting:   0.192611
        Dot Product:    4.96899
        Gathering:      0.0013752


Almost perfectly linear speed-up and communication time kept to minimum. QED.

Now, I'm not saying there's no load imbalance problem in Topographica, to be honest I'd be quite happy if there was because that would give me an interesting and challenging problem to solve and write about in my dissertation. I only want to admit that massive communication overheads that I observed were to a big extent not due to them but rather to SMP memory access bottleneck.

It all seems very obvious now, but I think I just had to demonstrate and prove that. Tomorrow I'm profiling Topographica to see what else is slowing it down and trying to get the effing vampirtrace to work. Also, I now have my very own account on Eddie and will be able to use it as soon as I get get my miserable 1GB disk quota extended to something a bit less depressing. I think I might still have my HECToR account, need to check, but the service is down at the moment anyway.

4 am ffs....

Thursday, 26 May 2011

Load Balance

While testing my code, I had a very interesting observation. I added computations into the loops with collective routines: now on each iteration nodes were not only communicating data, but also computing dot-products of local sets of distributed data with rows of a local array loc_arr (each row of loc_arr is dot-product-ed with the local set). I defined it as

loc_arr = beta(5,20,size=(ceil(N/100),N))

Initially it was a square NxN matrix but I divided the number of elements along the first dimension by a hundred to reduce the overall runtime (I just needed some computations there, not necessarily too heavy): fewer elements in first dimension => fewer dot-products => less time spent doing them. While changing the value of the divisor I noticed something interesting: computation time was changing, as expected, but communication time was changing as well! Furthermore, the communication time differed between application runs. I couldn't get my head around this for some time until I realised: of course, Early Reduce (well, scatter/... in my case) problem! Since I'm measuring time only on the master node, it could arrive to the communication routine while some other node hasn't finished computing! To check this I put barriers before each communication call, and - here we are, comm times went down, and became stable.

Now this made me think of the following: we were always assuming that data distribution in Topographica is load-balanced. Frankly, I wasn't even thinking about it until now. And suddenly, it seems I found the source of the problem with Topographica's massive communication overheads: poor load balance! I really have to look into the distribution patterns, and test this assumption, but I have a strong reason to believe it.

Apart from that, why the hell did my test code un-balance?! All nodes have exactly the same amount of data to compute, plus the computation time taken by all nodes was nearly the same! Some OS-related issues? I don't think so: in my benchmarking code in C there are no barriers before communication routines, and the timings are always stable... I should probably start using ECDF/HECToR to exclude any OS-related overheads anyway.

hacking on...

Quick summary for the day

Upgrade to new versions of mpi4py and OpenMPI went successfully. Then I decided to check whether my assumption that I could get some performance out of replacing the currently used in Topographica collective routines with different versions was correct. In particular, I wanted to see whether Scatter would work faster than scatter and Gather - than gather. The assumption was based on this conversation on mpi4py googlegroup: http://groups.google.com/group/mpi4py/browse_thread/thread/a1741896517c7d12/5c5e604047365348?lnk=gst&q=konstantin#5c5e604047365348. This way I though I would re-familiarise myself with topographica code and also hope to get some performance improvement.

To summarise, scatter/gather involve pickling while Scatter/Gather just make use of the underlying C routines MPI_Gather() and MPI_Scatter(). Thus, the advantage of the first way is that it can communicate arbitrary Python objects while the latter is meant to be faster.

I decided to write a simple Python application that gathers and scatters (later on added broadcasting) large arrays of floats in a loop, 100 times each, and measures the execution  time for each routine. N controlled the size of individual chunks, so the master node (rank 0) would gather/scatter N*P elements, where P is the number of processes used. After having finished the coding I played around with different values of N and P, also running my code in SMP mode and MPP (i.e. using a cluster of 3 machines). Here are some results:


N=10,000
SMP, 3 nodes:


scatter time: 0.0482320785522
Scatter time: 0.0166721343994
bcast time: 0.0553240776062
Bcast time: 0.041867017746
gather time: 0.0940809249878
Gather time: 0.0110039710999


MPP, 3 nodes:

scatter time: 0.177076816559
Scatter time: 0.176827907562
bcast time: 0.211786031723
Bcast time: 0.213672876358
gather time: 0.22430896759
Gather time: 0.168394088745

N=50,000
SMP, 3 nodes:
scatter time: 0.220118045807
Scatter time: 0.0812599658966
bcast time: 0.321494102478
Bcast time: 0.180558919907
gather time: 0.190404176712
Gather time: 0.0511260032654

MPP, 3 nodes:
scatter time: 3.98046803474
Scatter time: 0.681133985519
bcast time: 1.60409212112
Bcast time: 1.39482522011
gather time: 0.849080085754
Gather time: 0.734710931778

N=200,000
SMP, 3 nodes
scatter time: 1.79250884056
Scatter time: 0.383674144745
bcast time: 2.98570203781
Bcast time: 0.743601083755
gather time: 1.25536489487
Gather time: 0.288822889328


MPP, 3 nodes
scatter time: 3.97047996521
Scatter time: 2.73933315277
bcast time: 6.40952301025
Bcast time: 5.25424790382
gather time: 3.54560399055
Gather time: 3.84726691246


As we can see, SMP benefits from non-pickling routines the most, and the difference in performance seems to increase with the amount of data in each packet sent. On the contrary, MPP hasn't shown any significant difference between the performance of two routines. Except in the Scatter/scatter case. I think there must be a hint there: what if the performance is limited by the interconnect bandwidth? Thing is, in my experiment, data is sent in loops, without any interruption or work done between the communication routine calls. This is quite far from what a real application would look like, so I decided to put some work there. Now each node on each gathering/scattering iteration will do some computations. To make it a bit similar to topographica, each node will compute a dot-product of the distributed array and a local set of data. I'll measure both the computation and gathering/scattering time.

Wednesday, 25 May 2011

Parallel Topographica part 2

This is the first post for my MSc project. After the first meeting today my task for now is to review my old TODO list (which meant quick skim through, shift-select and delete), re-familiarise with the code and report on the current state.

Topographica install went ok, parallel stuff works as well. I checked if there were new versions of OpenMPI and mpi4py released, and yes they were: mpi4py 1.2.2 and OpenMPI 1.4.3. I successfully upgraded both, everything still seems to run. Think I'll go ahead and run a quick simulation in parallel just to make sure nothing has been broken.