Reproducing the result in DOI: 10.1038/ncomms7922

Request or post your example code here
cgastald
Posts: 11
Joined: Thu Nov 12, 2015 8:54 am

Reproducing the result in DOI: 10.1038/ncomms7922

Post by cgastald »

Hello,
I am trying to reproduce the results in the paper "Diverse synaptic plasticity mechanisms orchestrated to form and retrieve memories in spiking neural networks". I think the best way to visualize whether the patterns are stored correctly is to reproduce fig.4b. In order to reproduce fig 4b, I installed Auryn-0.6 and the code in the folder 2015orchestrated (downloaded from github). I run the code, without making any modification. In order to reproduce fig 4b, I sum the 3 seconds of population activity (taken from the output file rf3.0.pact) after a specific cue is shown for every time the cue is shown and for each read-out population. Finally I divide by times the cue is shown. I did this for the 20 cues and for the 4 read-out populations and plot the result in the same fashion as in fig 4b. However the result is not quite the same, in fact I never obtain more then 35 spikes, the read-out population associated to the square doesn't "recognise" the square cues, the read-out population associated to the circle has unexpected hight activity >1s after the cue is shown and finally the activity peaks are 0.5 s large (in fig 4b they are 1 s large).
Therefore, I would like to ask you 2 questions:
1) how did you define exactly the quantity "Spikes per neuron per time [Hz]"?
2) what do the ranks (in the code) correspond to? Do they correspond to anything physical? (sorry for this last question, it is the first time I try to use MPI)
Thank you very much in advance!
Attachments
averaged_pact.png
averaged_pact.png (130.92 KiB) Viewed 13219 times
User avatar
zenke
Site Admin
Posts: 156
Joined: Tue Oct 14, 2014 11:34 am
Location: Basel, CH
Contact:

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by zenke »

Hi Chiara,

thanks for posting. You are right, the results you are showing do not seem to show delay activity. Could you post a plot of the spike raster activity that you are getting in the end of the simulation? We could use this to get a first idea whether or not there is delay activity.

Some more explanation: My best guess is that something in the analysis is going wrong. I have stored the "rf1.pat" file with the code which I used (this file holds the information about which cells in the recurrent network prefer which out of the four stimuli). This file is loaded by PatternMonitor and performs online analysis in the sense that it monitors the population firing rates of the groups specified in the file (these are the ones you are plotting). However, using my file only works if you are running exactly the same network simulation (down to the same random seed and same number of ranks), because otherwise these groups/assemblies which form spontaneously will be different. I think the Auryn version I am suggesting to use in the readme is Auryn v0.5. Since, you are using Auryn v0.6, this suggests there have been some API changes and some of them affect the random seeding mechanism (I also recall some changes to the generation of random sparse matrices). These changes will yield a different connectivity at the micro level although statistically the networks are still the same. This would explain why you would not see anything. You would effectively be averaging out the actual selective delay activity. It also means we have to regenerate the rf1.pat file for the network connectivity you are using.

However, first we need to verify if that's actually the problem. That's why we should check if there are working memory states in the network. We will be able to guess this by looking at the raw spiking data, hence the raster plots. We will know more when we see those. If I am correct, I can explain to you how to generate a rf1.pat file for your network.

To answer your specific questions:
  1. Yes the units are Hz. First column: time in s, second column: activity in Hz, third column: network activity normalized by pattern size (you can just ignore that), the two last columns repeat for each pattern, last column: again normalized activity by population size. I hope that answers it.
  2. Ranks in the code correspond to the process number you are running. MPI starts the program you write multiple times as specified with the mpirun command. Each process is identical to begin with, the only way they can figure out they are different are via their "rank" property. Auryn programs use this property to distribute the different cells on different nodes.
Best,

Friedemann
User avatar
zenke
Site Admin
Posts: 156
Joined: Tue Oct 14, 2014 11:34 am
Location: Basel, CH
Contact:

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by zenke »

Hi Chiara,

it's me again. Being curious I just ran the first part of the simulations. Specifically, the 1h intensive learning period which is executed by "1run_init.sh". I now ran this on the Auryn development branch (downstream of what you are using). Towards the end of the learning period the spike raster starts to look like this:
1run_init_spikes.png
1run_init_spikes.png (17.97 KiB) Viewed 13212 times
The high firing periods correspond to the network being stimulated with an external stimulus. The less busy timse are spontaneous activity. Although working memory states are not reliably evoked at the end of the 1run_init.sh simulation, you see that the spontaneous activity following the second stimulus has some structure which is not entire flat and that is correlates to some degree with the evoked activity. That is the signature of delay activity to look for in your data.
During the second phase (less frequent external stimulation) this type of behaviour activity firms up:
2run_learn_spikes.png
2run_learn_spikes.png (24.21 KiB) Viewed 13212 times
I haven't run the long part of the sims because I don't have my nice little cluster any more which you are presumably using now ;-)
I hope that helps!
cgastald
Posts: 11
Joined: Thu Nov 12, 2015 8:54 am

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by cgastald »

Thanks for the quick answer!
Sorry for the long time it took me to answer: I re-run the simulation and re-checked my small code for plotting the averaged read-out population activity.
First of all, I must make a clarification: I actually used v0.5 of auryn which is in the folder auryn-0.6 on github. Therefore it is supposed to be the exactly the same network simulation as yours.

Secondly, thanks for the explanation about the file "rf1.pat", I saw it is responsible for the recognition of the read-out populations, but it would be very useful for me to know how to regenerate the rf1.pat, as I will need to do it very soon. (for example: how should rf1.pat change if I what the network to learn only the first 3 pattern?)

Now let's come back on the questions.
1) I didn't formulate very clearly the first one, I actually wanted to ask how is define the quantity "Spikes per neuron per time" that appears on the y-axes of figure 4.b, in the article. In other words, I am afraid I am not plotting the same quantity. So what I do in my plot is: every time a specific cue appears, I register the population activity, starting from the moment the cue appears for 3 seconds, in time bins of 0.1 s (which is the time-resolution in the pact file).

2) thanks for the explanation, it starts to get clearer. However, when I want to analyze the output, what is the difference between the output files rf1.0.pact and rf1.1.pact?

I attach the ras plots: they are indeed very similar to yours! :-)
Attachments
rf2.0.e.png
rf2.0.e.png (115.84 KiB) Viewed 13192 times
rf1.0.e.png
rf1.0.e.png (84.49 KiB) Viewed 13192 times
User avatar
zenke
Site Admin
Posts: 156
Joined: Tue Oct 14, 2014 11:34 am
Location: Basel, CH
Contact:

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by zenke »

cgastald wrote:Thanks for the quick answer!
Sorry for the long time it took me to answer: I re-run the simulation and re-checked my small code for plotting the averaged read-out population activity.
First of all, I must make a clarification: I actually used v0.5 of auryn which is in the folder auryn-0.6 on github. Therefore it is supposed to be the exactly the same network simulation as yours.
Alright. I should check the "rf1.pat" then in my simulation. It's also possible that a different boost version or any other component in the chain is different which might be responsible for the networks not being "exactly the same". In any way the plots you are attaching below look like somethings works in principle.
cgastald wrote: Secondly, thanks for the explanation about the file "rf1.pat", I saw it is responsible for the recognition of the read-out populations, but it would be very useful for me to know how to regenerate the rf1.pat, as I will need to do it very soon. (for example: how should rf1.pat change if I what the network to learn only the first 3 pattern?)
They way I did it in the paper is to compute a PSTH for the entire excitatory population during the stimulation intervals towards the end of the second phase. I then define a cut off on the PSTH which I think was >1.0 average activity during that interval. Cells which respond with higher than mean activity are counted as active. I will check if I have some instructive code I can give you for the analysis. I think it's all Python and AWK scripts.
cgastald wrote: Now let's come back on the questions.
1) I didn't formulate very clearly the first one, I actually wanted to ask how is define the quantity "Spikes per neuron per time" that appears on the y-axes of figure 4.b, in the article. In other words, I am afraid I am not plotting the same quantity. So what I do in my plot is: every time a specific cue appears, I register the population activity, starting from the moment the cue appears for 3 seconds, in time bins of 0.1 s (which is the time-resolution in the pact file).
It is the population firing rate of the sub-population of the neurons which prefer a given stimulus. Makes sense? Different stimuli have different population sizes (ie. you need to correct for this) and divide by the time bin size.
cgastald wrote: 2) thanks for the explanation, it starts to get clearer. However, when I want to analyze the output, what is the difference between the output files rf1.0.pact and rf1.1.pact?
These files are generated on ranks 0 and 1 respectively. Any online analysis will usually only process spikes which were generated locally on that rank (ie. processes, which might be running on physically separated computers). To generate joint statistics you can do this offline at a later time. In fact you should have four of those files, because the example runs 4 processes in parallel. The same is true for spike raster files for instance. In the raster files you will see that the different files contain the spikes from different neurons. Per default they are distributed in a round robin fashion (ie. the first rank has neuron 0, the second neuron 1, the third 2, the fourth 3 -- then rank 0 neuron 4 and so forth ... )
cgastald wrote: I attach the ras plots: they are indeed very similar to yours! :-)
Jup this looks like we are on the right track!
User avatar
zenke
Site Admin
Posts: 156
Joined: Tue Oct 14, 2014 11:34 am
Location: Basel, CH
Contact:

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by zenke »

Here are the relevant marco lines from the Makefile I used for the data analyses. This should give you an idea for the syntax of the scripts I used to generate the rf1.pat file (just check the meanings of % $< $@ for gnu Makefiles so the macros make more sense to you, you can run all the steps manually from the command line - no need to completely automatize it yet).

Code: Select all

%.pat: %.sta
    ../tools/sta2pat.sh $< $@
    
%.sta: /lcncluster/zenke/reset/%.e.ras %.trig
    rassta.py -t $(subst .e.ras,.trig,$(notdir $<)) -w 0.5 -f $< -s 3000 -m 3500 -q 0.0 -o $@
    
%.trig: /lcncluster/zenke/reset/%.0.stimtimes
    ../tools/mktrigfile.awk $< > $@
Here is the logic in brief:
  1. The mktrigfile script generates a trigger file from the *.0.stimtimes file from your simulation. The trigger file contains the trigger times for your population spike triggered average (in this case a stimulus triggered average really). The file will be called rf1.trig
  2. Then rassta.py computes the population spike triggered average (STA) and stores it in a file rf1.sta
  3. Finally sta2pat.sh does the thresholding and outputs a rf1.pat.
I am attaching the three scripts. mktrigfile and sta2pat.sh I pulled from the repository of the paper. They should be in their original state and work as expected. rassta.py is part of my malleable toolkit and being actively developed. I think it did not change, but it might have. Let me know, if there are any problems with running it.

Note, in the last step that I said we take every neuron whose firing rate is higher than the average rate, but it seems in the final version I actually took 10Hz as the threshold (you can see that in the script). I guess I did both for the paper and the results should not change dramatically, because 10Hz is about where the unstable fixed point is.

So ultimately you need to run the sim. Then get the pat file rf1.pat and then run it again with that file for the online analysis. Alternatively you can also generate the PatternActivity (pact) files with your own offline analysis. I hope what I have written makes sense. Sorry, need to run now.

F
Attachments
scripts.tar.gz
(2.03 KiB) Downloaded 593 times
cgastald
Posts: 11
Joined: Thu Nov 12, 2015 8:54 am

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by cgastald »

Thanks! This clarifies a lot. :D

So in order to re-produce your rf1.pat file, I commented the line where the monitor file is passed,

Code: Select all

--monf ./data/rf1.pat \
, in 1run_init.sh and run it. Note that I don't get the same ras plot as before.

So I run the scripts (without a makefile), using the commands

Code: Select all

./mktrigfile.awk  ../output/rf1.0.stimtimes
python rassta.py -t rf1.trig -w 0.5 -f ../output/rf1.0.e.ras -s 3000 -m 3500 -q 0.0 -o rf1.sta
./sta2pat.sh  rf1.sta rf1.pat
Note that I add

Code: Select all

> "rf1.trig"
at the end of line 13 in mktrigfile.awk. The rf1.pat that I got is not the same as your: in each read-out population there are only about 50 neurons.

At this point I re-run 1run_init.sh using the rf1.pat just obtained and I get again the "right" ras plot.
Finally, re-running the scripts a second time, I obtain the rf1.pat a bit more similar to yours. In fact the number of neurons in each read-out population is about 127 (in your rf1.pat there are about 660 neurons in each population) and some values are the same as in your file (in particular the ones at the top of the file). Finally, if I run everything a third time, I still obtain the same as the second time.

I was expecting to get the "right" ras and pat files at the first shot, while I have to run everything twice (and I still don't get the same rf1.pat), I didn't give this much though yet, but it doesn't look clear to me.

Chiara
User avatar
zenke
Site Admin
Posts: 156
Joined: Tue Oct 14, 2014 11:34 am
Location: Basel, CH
Contact:

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by zenke »

Hi Chiara,

seems like we are on the right track. Yes, you should only have to run the sim twice. Since you are not changing anything to the simulation, you should get the same ras plots every time. You are only changing the readout which will affect the *.pact files. However, you probably won't get the same rf1.pat file as me, because something in your simulation is different (could be due to different library versions).
However, you should have a comparable amount of neurons in each pattern. One thing which probably contributes to not having as many neurons is that you used "../output/rf1.0.e.ras" as the input to the stimulus triggered average (STA). This file only contains the spikes of one fourth of all the neurons if you run the code as 4 parallel processes. Sorry, I should have mentioned this :roll: To avoid this problem you need to merge the multiple ras files from your sim. On the Linux command line the easiest way is:

Code: Select all

sort -g -m ../output/rf1.*.e.ras > ../output/rf1.e.ras
See also http://www.fzenke.net/auryn/doku.php?id=manual:ras

As a next step you run rassta.py on "rf1.e.ras" with the trigger file you have. This should give you a more reasonable result. Could you please post your pact plots again? At least you should see clearly different population responses during stimulation intervals in your plot when you are using the new rf1.pat as the "cells to monitor file" ("--monf").

Cheers,

F
cgastald
Posts: 11
Joined: Thu Nov 12, 2015 8:54 am

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by cgastald »

Hello! So I made my own fr1.pat file, which is very similar to yours and I attach again my pact plot (rf3.png). As you can see, there is still a problem, i.e. the network doesn't recognise the square cues. What I noticed is that, if I plot the activity of the four read-out populations at the end of 1run_init.sh (rf1.png) and 2run_learn.sh (rf2.png), the population activity associated to the square is the only one that doesn't increase. I guess this means that the consolidation failed for this pattern, doesn't it?
A interesting point is that, if I make the network learn only 3 patterns instead of 4 (I excluded the triangle) I dosen't have any more the problem with the square (see rf3_3pop.png).
I will think more about this :-)
Attachments
rf1.png
rf1.png (44.47 KiB) Viewed 13136 times
rf2.png
rf2.png (45.73 KiB) Viewed 13136 times
rf3.png
rf3.png (117.18 KiB) Viewed 13136 times
cgastald
Posts: 11
Joined: Thu Nov 12, 2015 8:54 am

Re: Reproducing the result in DOI: 10.1038/ncomms7922

Post by cgastald »

Here the last plot
Attachments
rf3_3pop.png
rf3_3pop.png (112.57 KiB) Viewed 13136 times
Post Reply