Add Poisson Spike Train neuron function

[GuillaumeLam]

Commit that adds the necessary feature is 753a486.

[SBuercklin]

Could you provide an example use case for poissonST, like an example piece of code? It would help me contextualize this addition within the larger scope of WaspNet.jl. I'm curious if this function used for driving a Network, intended to slot into a Layer, or some other use

[GuillaumeLam]

Thanks for leaving such helpful comments! Most comments are very clear to me, I just had a couple questions/remarks that I'll put into context with a piece of code for poissonST:

simulate!(lsm, poissonST([1,1,1]), 0.001, 0.1)

Thanks to the design, I was able to implement what I needed in this fashion. So here the vector [1,1,1] denotes the input to the network for this input period.

In the tutorials, the simulation is driven by a function where the input is t. When I would simulate a Network with a function that takes no input, I get an error at the line 61 in simulate.jl [ input_matrix = hcat(input.(tsteps[1:(end-1)])...) ]. :/ This didn't seem like an easy fix so I added a throwaway variable. Could I add a `be a valid alternative? Also, I noticed that the function for simulating aNeuronvsLayervsNetworkwas different, so I decided to only write the functionality forNetworkas I would rarely be simulating a layer only. Would there be an easy fix for this issue or should it remain only utilizable byNetwork`?

After re-reading your comments, I think I have in mind what you meant with closure over l and Bernoulli distribution. So I will push the code and please lmk if it fits with what you had in mind.

[GuillaumeLam]

So while implementing the method with closure, I wanted to see the speed up. It so happens that the original method is pretty quick. Furthermore, since there is closure over the Bernoulli variable but it needs to be generated for every val of the inputs (converted to prob) and for each simulation time, I can't seem to find the type to solve the code instability.

using Distributions

l = [1,1,1]
d = Bernoulli

@btime getPoissonST(l,d)
@btime getPoissonST_old(l)

@btime getPoissonST(l,d)(1)
@btime getPoissonST_old(l)(1)

@code_warntype getPoissonST(l,d)(1)
@code_warntype getPoissonST_old(l)(1)

Here's some code that will time and show the code instability. I changed the original implementation for which the timing is the same but it's a bit cleaner and spits out a BitVector. Lmk if I implemented the way you thought of it and which we should go with!

[SBuercklin]

Thanks for giving this a shot, and I like your addition of a d field to allow for more general input distributions. It might help lay the foundation for more interesting signalling down the line :)

Regarding the implementation and timing, there are a couple things we can look at.

On the implementation, we're allocating arrays on every time step with the broadcasting. This gets us the answer, but it's not quite the fastest way. Every time you need to make a new array instead of re-using an old one, and when those arrays get big (as they would for large input layers), the allocations add up for longer calculations. An example of this is in WaspNet.jl here. We have the input signals in a vector input, and then we do the matrix multiplication in-place to a preallocated l.input vector with the mul! call, then this l.input vector drives every AbstractNeuron in the Layer.

Below, we can preallocate the vector we store the signal output in, and it'd be good if we could sample of the distribution without needing to broadcast either since that would also allocate an array. It turns out Distributions.jl has a distribution type which can do exactly that, so we don't need to broadcast the rand calls anymore. My implementation is below.

using BenchmarkTools
using Distributions
using LinearAlgebra

function getPoissonST_old(l::AbstractVector)
    return (_) -> rand.(Bernoulli.(normalize(l)))
end

function getPoissonST(l::AbstractVector, d)
    norm = normalize(l)
    return (_)->rand.(d.(norm))
end

function getPoissonST_propose(l::AbstractVector, d)
    norm = normalize(l)                       # normalize the array
    ds = product_distribution(d.(norm))       # make a distribution we can sample with rand
    sample = rand(ds)                         # pre-allocate the output array
    return (_) -> rand!(ds, sample)
end

Now, let's compare how each of them performs. I did two different tests:

1. We construct the function and then generate a random output one time
2. We have the function already generated, we just generate the output multiple times

If we need a different distribution at every time step, then the first number is the important one. If the distribution is unchanging over time, then we just need to generate the output quickly, and this corresponds to the second case below.

arr = rand(256) # Let's test with a bigger input layer

julia> @btime getPoissonST_old($arr)(1);
  2.111 μs (4 allocations: 6.45 KiB)

julia> @btime getPoissonST($arr, $Bernoulli)(1);
  2.589 μs (9 allocations: 6.59 KiB)

julia> @btime getPoissonST_propose($arr, $Bernoulli)(1);
  2.978 μs (7 allocations: 4.70 KiB)

old = getPoissonST_old(arr)
current = getPoissonST(arr, Bernoulli)
proposed = getPoissonST_propose(arr, Bernoulli)

julia> @btime old(1);
  2.100 μs (4 allocations: 6.45 KiB)

julia> @btime current(1);
  2.133 μs (7 allocations: 4.44 KiB)

julia> @btime proposed(1);
  943.333 ns (0 allocations: 0 bytes)

So if we pre-allocate everything as in the proposed solution, we can sample an input distribution of size 256 in half the time, and we don't have to allocate any space. But, if you need the distribution to change at every time step, then my implementation is the worst performing.

The benefits of pre-allocating will vary with the size of that layer, and the gains here aren't going to be noticeable unless you're running a long simulation. This was probably a case of premature/unnecessary optimization on my part now that I tried it out.

[GuillaumeLam]

Thank you so much for taking the time the explain this to me! Your proposed method makes alot of sense so I'll commit it and push it to the pr!

I did have a couple small questions: Would an allocation be saved by not having to assign the normalized vector? (ie ds = product_distribution(d.(normalize(l))) ) Is this idea to pre-allocate a feature readily available through Julia? (rand! seems like a part of the key)

[SBuercklin]

Thank you so much for taking the time the explain this to me! Your proposed method makes alot of sense so I'll commit it and push it to the pr!

Is this idea to pre-allocate a feature readily available through Julia? (rand! seems like a part of the key)

Pre-allocation is generally available, but it means that you have to work in-place, meaning that the same memory (array) is used repeatedly. The alternative is every time we call that function, we need to block out a piece of memory the same size as the array each time, which takes time and eats up memory that we don't need to necessarily use.

Here's a section from the Julia manual about pre-allocating outputs: https://docs.julialang.org/en/v1/manual/performance-tips/#Pre-allocating-outputs

The "bang" or ! at the end of a function name implies to the user that the function is mutating or operates in-place. Many functions have an out-of-place and in-place version, like rand vs rand!, where the former will allocate space but the latter will do it in-place, in an array that you provide. Matrix multiplication is another operation which we usually like to do in-place:

LinearAlgebra.mul!(Y, B, C)

carries out A*B and stores the result in Y, so if you have Y allocated already in the right shape, you can save some time not needing to allocate.

I did have a couple small questions: Would an allocation be saved by not having to assign the normalized vector? (ie ds = product_distribution(d.(normalize(l))) )

We might be able to reduce allocations further by some trick with how we normalize the vector of probabilities and construct the distribution. However, this is an operation that happens once, so if it's already pretty fast then we don't save much time. If we needed to do this allocation on every step, and it was significantly slowing us down, then we'd want to look into it. I just checked and it looks like normalize! is a function that exists, so we could conceivably normalize l in-place and save on an allocation. If we were really chasing performance, we could do that, but I don't think it's necessary in this instance.

[SBuercklin]

We also have some rebasing to do on this PR before it can merge with master as well. Let's get the inhib-neuron PR merged, then we can rebase this branch onto master and get this available as well.

[SBuercklin]

I've merged the proper commits into master. Can you check that the behavior on the master branch of this repo matches what you expect? Once that's confirmed, I'll close this PR and bump the version available on the registry

[GuillaumeLam]

Hello Sam! Thank you for merging in the Inhib PR! I had an issue as I needed to control the rng of the spike trains. As I was trying to rewrite the method, I had an idea for a more general struct. (mininal code idea on my repo: link). This struct could act (and be named) as encoder to the network. This could allow for more general and different spike train generation based on an input. Additionally, a decoder struct could also be made to dictate the behaviour of how to interpret output of the spiking nn! If you think this fits in the scope of your repo, I would gladly make the required changes and we can keep the pr open. If not, I don't mind keeping the code separate. Lmk what you think!

[SBuercklin]

Hi Guillaume, it's been a while since I looked at this, but I think it'd be better to leave this as it is. Part of the point of WaspNet.jl was to be very extensible. If you have a lot of utilities that you're finding useful, I think a great idea would be to make another package (e.g. WaspNetUtils.jl or something similar) which extends the base library and can be loaded on top of this. You could add this to the Julia General Registry and have it widely available.

Another reason to do this is that when I find the time, I'm interested in giving the DiffEq.jl backend another shot as mentioned on this Discourse thread: https://discourse.julialang.org/t/spiking-neural-networks/63270/12

I'm not sure if we'll be able to maintain compatibility with the existing interface in all cases, so there will likely be some changes to make to any utilities. The DiffEq.jl backend would also constitute a version change that can be used to set compatibility with these utilities; if your utilities can't be updated to the new interface, setting a [compat] bound in the Project.toml for your utilities will ensure that you remain on the last compatible version of WaspNet.jl.

If you have any questions about getting a project on the registry, feel free to ping me or ask on the community on Discourse. I'm also on the Julia Slack if you want to ask questions directly, I suspect finding me on there by name shouldn't be too hard. I'm very happy to help answer questions about Julia development or point you to resources that I found useful.

All of this being said, I'm going to close this PR and bump the WaspNet.jl version on the registry with all of the work that you've helped integrate.

[GuillaumeLam]

Hey Samuel!

So lately i was trying to add more stuff and realized that having it in a separate package would allow for faster iterations! All this to say that you absolutely right! I've still been adapting to the functional and reusable nature of functional programming of Julia. Just a few days ago I created a package named LiquidStateMachine.jl.

Thank you for all your help and for the great piece of software you've written! Thanks for offering your help with package registry! I did have a question about what is the idea behind adding in DiffEq.jl? The neuron model currently are solving diff eq through the Euler method, would the new backend allow for more seamless solving of diff eqs?

Thanks again and hope to ttys,

Guillaume Lam

---

From: Samuel Buercklin @.***> Sent: Saturday, October 30, 2021, 11:04 AM To: leaflabs/WaspNet.jl Cc: Guillaume Lam; Author Subject: Re: [leaflabs/WaspNet.jl] Add Poisson Spike Train neuron function (#17)

Hi Guillaume, it's been a while since I looked at this, but I think it'd be better to leave this as it is. Part of the point of WaspNet.jl was to be very extensible. If you have a lot of utilities that you're finding useful, I think a great idea would be to make another package (e.g. WaspNetUtils.jl or something similar) which extends the base library and can be loaded on top of this. You could add this to the Julia General Registry and have it widely available.

Another reason to do this is that when I find the time, I'm interested in giving the DiffEq.jl backend another shot as mentioned on this Discourse thread: https://discourse.julialang.org/t/spiking-neural-networks/63270/12

I'm not sure if we'll be able to maintain compatibility with the existing interface in all cases, so there will likely be some changes to make to any utilities. The DiffEq.jl backend would also constitute a version change that can be used to set compatibility with these utilities; if your utilities can't be updated to the new interface, setting a [compat] bound in the Project.toml for your utilities will ensure that you remain on the last compatible version of WaspNet.jl.

If you have any questions about getting a project on the registry, feel free to ping me or ask on the community on Discourse. I'm also on the Julia Slack if you want to ask questions directly, I suspect finding me on there by name shouldn't be too hard. I'm very happy to help answer questions about Julia development or point you to resources that I found useful.

All of this being said, I'm going to close this PR and bump the WaspNet.jl version on the registry with all of the work that you've helped integrate.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/leaflabs/WaspNet.jl/pull/17#issuecomment-955284650, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AFQJINOOXNMZM7QFMWNRZETUJQCN5ANCNFSM5AKCCPAA. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.