Add monitor to estimate rates of each neuron in group

[awillats]

It's often useful to estimate a smoothed estimate of firing rate for each neuron in a group.

However, the current options are:

1. PopulationRateMonitor which has a smooth_rate() method, but averages across the population
   • in theory you could connect one of these to each neuron in a group individually but this is somewhat clunky
2. SpikeMonitor which can return an array of spike_trains() one for each neuron in a group, but doesn't have a built-in smooth_rate() method

it would be great to extend SpikeMonitor to be able to call smooth_rate() for each neuron.

as a bonus, it would be nice for an option of returning a vector of spike counts binned at a certain interval (this could be thought of as smooth_rate(window='flat', width=dt, stride=dt) but may also be clearer as it's own function)

[mstimberg]

This would indeed be a nice feature. If anyone wants to work on this, please let us know – I think this is a great issue for external contributors, since it does not touch on anything related to code-generation (which is code that is quite hard to read).
Here are some thoughts on how I think this could be implemented:

1. We introduce a new class RateMonitor that provides two functions: smooth_rate and binned (I agree that it would be better to keep them separate). Both PopulationMonitor and EventMonitor (SpikeMonitor is just a special case) inherit from this class.
2. The binned function takes a single argument bin_size, stating the desired bin size in seconds (which should be a multiple of dt). I think this should probably be an abstract function without implementation in RateMonitor, since the implementation is quite different between EventMonitor and PopulationRateMonitor. In particular, for bin_size=dt, the PopulationRateMonitor will simply return the recorded rates without doing any binning. I think the function should return a pair (bins, binned_values), where the first part are the midpoints of the bins, and the second part are the actual values as a 2D array for EventMonitor, and as a 1D array for PopulationRateMonitor. I guess for consistency with StateMonitor, the 2D array should have shape neurons × bins. This has the disadvantage that for plotting you need to do transpose it (plt.plot(bins, binned_values.T)), but on the other hand it is more straightforward to plot the binned values for a single neuron (plt.plot(bins, binned_values[idx])).
3. The smooth_rate function can probably have only one implementation (which is useful if we want to add other window functions, etc.). Instead of referring to something like self.rate, as it is currently done in PopulationRateMonitor, we'd refer to self.binned(bin_size=self.clock.dt) – this would create the binned array for EventMonitor, but would simply return self.rate for a PopulationRateMonitor.

[lysea-haggie]

Hello! I am interested in working on this since I submitted the comment on the forum about it. I'm a graduate student not a software developer but am keen to try and learn. Let me know how I can get involved. My email is [email protected].

[mstimberg]

Hi @MunozatABI Apologies for the very late reply! It would be great if you could work on it. The first step would be to fork the repository and start working on it in the way I outlined earlier. As soon as you have something that is worth discussing (can be quite an early state), you can open a pull request. I'd prefer to have discussions here (or in the new pull request) instead over email. Please let us know if things about the general development process, the code itself, or anything related to the issue are unclear, or if you run into any unexpected issues.
https://docs.github.com/en/get-started/quickstart/contributing-to-projects

[lysea-haggie]

Thanks @mstimberg. I forked the repository and found the SpikeMonitor and PopulationRateMonitor modules. I was just wondering what is a good way of testing it? I saw the test folder and test_monitor.py file but it does not run. Apologies if I am missing something obvious in this process.

[mstimberg]

Hi. You should be able to run test_monitors.py directly, it will call all the test functions. What kind of error do you get when you try that? This way of running it is only a convenience for running tests during development, though. It will not run the tests for all the code generation targets, etc. The "official" way to run the test suite is to call the brian2.test() function, or the run_tests.py script in dev/tools (which does the same thing). This will run all the tests though, which will take some time. You can pass additional arguments to brian2.test in the additional_args argument, which will be forwarded to pytest. E.g. with additional_args=['-k', 'monitor'] it should only run tests with the word monitor in their name.

[lysea-haggie]

When I run the test_monitor.py file I get the error below (with no changes or added files).

The brian2.test(additional_args=['-k', 'monitor'] completes successfully. But what is this running?

[mstimberg]

Hi @MunozatABI , sorry for the late reply, I somehow overlooked your comment. From the stack trace, it seems that you are running tests in Brian's checked out repo (C:/.../Documents/Repos/Brian2), but it imports brian2 from an existing installation in your Anaconda root environment. The test and the code it tests therefore do not come from the same version, leading to the error. If you run brian2.test(...), it runs both from the same location (the installed package). Also see https://brian2.readthedocs.io/en/latest/introduction/install.html#development-install

[mstimberg]

Hi @MunozatABI. Just as a follow-up: would you still be interested in working on this?

[lysea-haggie]

Hello, yes I would be interested in giving it a go but also don’t mind if there is someone else who can do it faster/better as I’m also currently writing my thesis. Was there a date you wanted it done by?

[mstimberg]

There's no rush and there is no specific date. I was just curious since this is one of the few open issues that could be very useful for users which can be implemented by someone outside the "core" developer team. Either way, if you or anyone else starts to work on this, please let us know here in this issue.

[lysea-haggie]

Hi, I can run the test_monitor.py file and have been going through the code. I've been looking at the test_rate_monitor_1 function to get an idea of how the PopulationRateMonitor is built. Just a few questions:

• When I step into the rate_mon = PopulationRateMonitor(G) line it takes me to tracking.py and doesn't seem to go through ratemonitor.py, so how is it calling the PopulationRateMonitor?
• In the PopulationRateMonitor class in ratemonitor.py I can see 'rate' as a variable added as a dynamic array, but where is it actually calculating the rate?
• As mentioned above, would the way to do this to be to add a new class 'RateMonitor' to the ratemonitor.py file? Would I use a similar structure with the magic network lines

 invalidates_magic_network = False
    add_to_magic_network = True

What do these do?

• Also would I use a similar __init__ structure? I'm currently unsure how to access the spiking data from the source group.

    def __init__(self, source, name='ratemonitor*', codeobj_class=None,
                 dtype=np.float64):

        #: The group we are recording from
        self.source = source

        self.codeobj_class = codeobj_class
        CodeRunner.__init__(self, group=self, code='', template='ratemonitor',
                            clock=source.clock, when='end', order=0, name=name)

        self.add_dependency(source)

        self.variables = Variables(self)

Anyways, that's where I'm at! I may be out of my depth here.

[mstimberg]

Hi @MunozatABI, here a few comments/suggestions:

• When I step into the rate_mon = PopulationRateMonitor(G) line it takes me to tracking.py and doesn't seem to go through ratemonitor.py, so how is it calling the PopulationRateMonitor?

This is because PopulationRateMonitor inherits from Trackable (via Group→BrianObject→Nameable) and this class defines a __new__ method which is called before __init__. This system is used to keep track of all created Brian objects so that we can give them unique names (neurongroup_1, etc.). You do not have to worry about this here.

• In the PopulationRateMonitor class in ratemonitor.py I can see 'rate' as a variable added as a dynamic array, but where is it actually calculating the rate?

This is happening in generated code (Python, Cython, or C++ code). The code generation machinery is quite complex, but you do not need to worry about it here.

• As mentioned above, would the way to do this to be to add a new class 'RateMonitor' to the ratemonitor.py file? Would I use a similar structure with the magic network lines

Yes, we'd like to have a new RateMonitor class in ratemonitor.py, but this should be an abstract class. Its only role is to provide the smooth_rate and binned classes to PopulationRateMonitor and EventMonitor, it never gets instantiated. The class attributes related to magic networks are therefore not needed. They are related to the fact that you can run a "magic network" (i.e. a network where you don't manually add all the objects) for a certain time and then add a monitor and continue to run it. This is ok, but adding a new group like a NeuronGroup would not be ok, since it is unclear whether the user actually wants to continue a simulation or start a new one. See https://brian2.readthedocs.io/en/stable/user/running.html#multiple-magic-runs for more details. But again, feel free to skip it, it is not needed for this issue :-)

Also would I use a similar __init__ structure? I'm currently unsure how to access the spiking data from the source group.

See my comment above, we don't need an __init__.

My initially suggested approach was mostly a "note to self" (or a note to someone that knows Brian internals well). Let me suggest a more step-wise approach:
Do not worry about a new RateMonitor class for now, instead:

1. Add a new function def binned(bin_size) to PopulationRateMonitor, which takes a desired bin_size (as a Quantity, i.e. something like 0.5*ms) and returns the corresponding binned spike rate. This function should:
   • raise an error for a bin_size < self.clock.dt
   • return self.rate if bin_size == self.clock.dt
   • return an appropriately binned rate for bin_size that are multiples of self.clock.dt , e.g. if the PopulationRateMonitor record the rates [0, 1, 1, 1, 0, 2]*Hz with a dt of 0.1ms, then asking for rate_mon.binned(0.2*ms) should return [0.5, 1, 1]*Hz.
2. If that seems to work, add def binned(bin_size) to EventMonitor. The implementation here will be different, since the EventMonitor does not record any rates. Instead, it has to calculate them from the individual event times. You can get the event times for each neuron by calling self.event_trains(). This returns a dictionary mapping each neuron index to its spike times. E.g. if this dictionary is:

{0: [0]*ms,
 1: [0.5]*ms,
 2: [1.5, 2.5]*ms}

then binned(1*ms) should return

[[1000, 0, 0],
 [1000, 0, 0],
 [0, 1000, 1000]]*Hz

and binned(2*ms) should return

[[500, 0]
 [500, 0]
 [0, 1000]]*Hz

i.e. the firing rate for each neuron in each of the time bins.

If you get that far (or at any earlier time when you need specific feedback on your code), please open a draft pull request with the changes.

[lysea-haggie]

Thanks for the step by step.

How do I test my binned function? I added a function to test_monitor.py, basically something like:

def test_rate_monitor_binned():
    G = NeuronGroup(1, 'v : 1', threshold='v>1') 
    G.v = 1.1 # All neurons spike every time step
    r_mon = PopulationRateMonitor(G)
    run(1*ms)
    binned_rate = r_mon.binned(0.5*ms)

but I get this error.

The test function does work if I do the smooth_rate function and I think I added the binned function to the PopulationRateMonitor class in the same way that smooth_rate is defined. So what am I missing?

Also, in regards to your 3rd step:

return an appropriately binned rate for bin_size that are multiples of self.clock.dt , e.g. if the PopulationRateMonitor record the rates [0, 1, 1, 1, 0, 2]Hz with a dt of 0.1ms, then asking for rate_mon.binned(0.2ms) should return [0.5, 1, 1]*Hz.

Does the bin_size have to be a multiple of self.clock.dt? So would I find the width of the window in terms of steps of dt (ie. bin_size/self.clock.dt) then find the average rate in the window?

[mstimberg]

Hi @MunozatABI

How do I test my binned function? I added a function to test_monitor.py, basically something like:

the error indicates that it is running things in an installed version of Brian, which does not have the changes you applied. If you installed brian2 with conda, uninstall it with conda remove brian2; if you installed it with pip, uninstall it with pip uninstall brian2. You can then install the version you are working on in "development mode" by using pip install -e . in the C:\Users\lmun373\Documents\Repos\Brian2 folder. The "installation" will then just install a link to your development repository instead of copying all the files. Changes will therefore automatically been taken into account when you run the tests for example.

Does the bin_size have to be a multiple of self.clock.dt? So would I find the width of the window in terms of steps of dt (ie. bin_size/self.clock.dt) then find the average rate in the window?

Yes and yes :) Checking whether a certain window is a multiple of dt is a bit tricky in practice due to floating point issues (e.g. 0.3 != 3*0.1), but don't worry about it for now and use something like bins = int(round(bin_size/self.clock.dt)).

[mstimberg]

Hi @MunozatABI just a ping to see whether you are still interested in working on this and to see whether you encountered any issues you need help with?

[lysea-haggie]

So sorry! Been a bit busy with writing my thesis and getting COVID... Think I'm a bit over committed right now so this isn't a priority at the moment but I'll have more time next year after I submit? Totally happy for someone else to take this on if I am being more of a hinderance than a help.

[mstimberg]

No worries, @MunozatABI. If someone else wants to pick this up, please let me know. If not, happy to work with you on this task when you have more available time.

[mstimberg]

@lysea-haggie It's a bit of a shot in the dark, but please let me know if you are still interested in working on this and have the bandwidth for it.

[lysea-haggie]

Hello! I do have some bandwidth these days now that I've submitted my thesis. I'll have a look over the next few days =)

…

On Thu, Nov 30, 2023 at 2:12 AM Marcel Stimberg ***@***.***> wrote: @lysea-haggie <https://github.com/lysea-haggie> It's a bit of a shot in the dark, but please let me know if you are still interested in working on this and have the bandwidth for it. — Reply to this email directly, view it on GitHub <#1368 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMK6HOGC3XLYFO6T67SF5ALYG4YCLAVCNFSM5HMNFMYKU5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCOBTGE4DOMJUGE2Q> . You are receiving this because you were mentioned.Message ID: ***@***.***>