does this sound like an overkill?

100 combined LoC for a python script that builds graphs...

Note: if you're planning to build something like this with matplot don't use pyplot.contour or .contourf unless you want to get some serious headache. It's easier to use pyplot.fill (fillx, etc.). I should've discovered it earlier


[update] 250-something LoC

Current Plan

Quoting Jim:


1. Measure effect of Bare-C timing code by comparing real times for
   versions with and without timing code.  As long as this difference
   is small (under 1%) for sufficiently large problem sizes, then we
   can ignore the timing code.  Otherwise we'll have to account for it
   in any plots with total run time.


2. Finish reporting and plotting, making it very clear (a) how the
   components add up to the total time (and highlighting any
   differences between the sum of the components and the total time),
   and (b) which numbers and plots are comparable between the
   Topographica (hybrid Python+C) and the bare C versions.


3. Prepare good, clear plots of the scaling behavior for us to
   discuss: density versus time (on a single processor), and then
   number-of-CPUs versus time for various (relatively large) density
   values.
   
4. Analyze the plots to understand the scaling behaviour of the bare C
   code with number of CPUs (each with different memory subsystems).
   Is the communication overhead (and any other overheads) low enough
   that it scales well, up to a hundred processors or so?  


5. If the bare-C scaling with CPU is bad, we probably need to focus on the
   communication pattern, reducing time CPUs spend waiting on data.


6. If the bare-C scaling with CPU is good, rejoice and start focusing
   on the Topographica implementation, first making it have the same
   pattern of communication and C computation. Once it appears to be
   doing the same thing as the bare-C version, analyse the scaling
   with CPU for Topographica -- does the Python overhead limit the
   number of processors that are useful?  Is that limit high enough
   for our practical purposes?  If not, we'll need to focus on
   optimizing other parts of Topographica until the scaling behavior
   with CPU is good.


7. Once we have good scaling with CPU for Topographica, for
   sufficiently large problem sizes, we'll be ready to release that to
   everyone to start using.
   
8. Then we can focus on scaling with cores -- can we do anything to
   exploit other cores that share memory hierarchies?  Right now we do
   get some speedup using OpenMP on some shared-memory machines, but
   it's far from ideal. We'll at least want to exploit what speedup we
   already have from using all the cores, and possibly investigate how
   to make it better.  But we're unlikely to get quite that far in
   this project.


Suggested time output:

Total time (CPU usage 96%): 14.301s

Components: init=8.387s + sim.run=5.48s = 13.867s

Unaccounted-for: 0.434s (3.0%)

I'm hoping to finish my plotting script in about 2 days. As soon as that happens I will start using ECDF instead of Jupiter.

Hand-crafted profiling

Just finished placing timers around all the ins and outs in my C code and formatting the outputs. Now the result of a single C vs Python run looks like

./c_vs_python.sh -m 100 -n 100 -i 30 -s examples/matchmodel.ty
TOPOGRAPHICA RUNS

### Python+C ###
Total script init time: 8.38710808754
Serial: 5.48633885384

real    0m14.301s
user    0m13.097s
sys     0m1.021s

### BARE C ###

====================== main() ======================
Total:                                  10.7808
----------------------------------------------------
Utils:                                  0.0557
Create and normalise CFs:               1.4893
Serial activate:                        4.6100
Barrier between serial and parallel:    0.0000
Parallel activate:                      4.6257
=================== activateSeq() ==================
Total:                                  4.6100
----------------------------------------------------
Utils:                                  0.0008
Dot product:                            4.6083
Output function:                        0.0008
=================== activateMPI() ==================
Total:                                  4.6257
----------------------------------------------------
Utils:                                  0.0065
Distributing weights:                   0.0001
Broadcasting inputs:                    0.0000
Dot product:                            4.6174
Gathering activity:                     0.0008
Output function:                        0.0008
====================================================

real    0m10.840s
user    0m10.484s
sys     0m0.303s

### Comparing results ###
Results matched

Timings are printed to screen only for the first node, although in this run I only used just one. For more detailed per-processor results (only for bare C code yet) the timings are printed into separate files (appending if the files exist) called timings.X.out where X corresponds to the rank of the node. A sample file might look like this:

#size 4
#m 100
#n 100
#i 20
main_tot 7.87351489
main_utils 1.07544136
main_create_cfs 0.00000095
main_seq_act 0.00000095
main_mpi_act 2.45419502
main_barrier 4.34387183
act_seq_tot 0.00000000
act_seq_utils 0.00000000
act_seq_dp 0.00000000
act_seq_of 0.00000000
act_mpi_tot 2.45419288
act_mpi_utils 0.00215340
act_mpi_distr_weights 1.10045886
act_mpi_bc 0.00546145
act_mpi_dp 1.34518266
act_mpi_gather 0.00084949
act_mpi_of 0.00001097
#end

I planned this to be more machine-readable for later use in two plotting scripts that I'm currently working on: one would take filename1 filename2 ... "param1 param2 ... >paramN" as arguments to produce N bars where params 1 through N-1 make up the total of paramN (e.g. timings.0.out timings.1.out "act_seq_utils act_seq_dp act_seq_of > act_seq_tot" should show two bars with utils, dot-product and output function displayed as percentages of the total serial activation time). The second script would work with scaling of one parameter against another. For example the total parallel activation time against the number of runs, the number of processors, or the size of the input sheet (m*n). That will be displayed in form of line graphs. Sounds like a lot of work, but in the end I will have a hand-crafted and very customisable set of profiling tools that do exactly what I need, should be worth the effort.

Un-squaring simulations

Got a bit distracted by other work in the previous few days, but thankfully it's all sorted now. Today I fixed a very annoying bug that wasn't allowing me to use rectangular sheets (i.e. with sides of different lengths, as opposed to square sheets): forgot to multiply one value by -1 and there we go, hours of debugging.

Other than that, work in progress on detailed timings and pretty graphics...

Week 3

Another weekly meeting. The aim for this week is to improve my testing script to do graph plotting and measure CPU load. Apart from that I will perform series of parallel scaling experiments on Python and bare C topographicas. If everything goes well, there will be a bunch of colourful graphs here by the end of this week.

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.

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

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.