SpikeSourceInhGamma

[neworderofjamie]

This is not currently implemented

[neworderofjamie]

Spike source, generating realizations of an inhomogeneous gamma process, employing the thinning method.

See: Muller et al (2007) Spike-frequency adapting neural ensembles: Beyond mean-adaptation and renewal theories. Neural Computation 19: 2958-3010.

[neworderofjamie]

If you fancy a maths-based challenge you could try and implement this! Not sure how useful it is as I've never encountered it used but it would be nice for completeness.

[neworderofjamie]

Parameters are: shape parameter a = α = k and an inverse scale parameter b = β = 1/θ, called a rate parameter.

[neworderofjamie]

scaling by b is multiplication

[tnowotny]

hum ... maybe it's not as bad as we thought to sample the gamma distribution directly: http://www.hongliangjie.com/2012/12/19/how-to-generate-gamma-random-variables/

[tnowotny]

And here is the original gsl code without the typo(s) from the above page: https://github.com/ampl/gsl/blob/master/randist/gamma.c It is distributed under GPL v3, which I think would mean we would be able to grab & modify it to fit our purposes ...?

[neworderofjamie]

Like I said, maybe add gammaDistFloat and gammaDistDouble in the style of the functions already here https://github.com/genn-team/genn/blob/master/lib/src/generateRunner.cc#L2176-L2207 then hook it to the code generator here https://github.com/genn-team/genn/blob/master/lib/include/codeGenUtils.h#L111-L125 (for consistency I'd use std::gamma_distribution for the CPU version)

[neworderofjamie]

You could maybe split it up so, for efficiency, you could provide precalculated cs and ds if you wanted to.

template<typename RNG>
__device__ float gammaDistFloatInternal(RNG *rng, float c, float d)
{
    while (true)
      {
...
     }
}

template<typename RNG>
__device__ float gammaDistFloat(RNG *rng, float a)
{
    if (a > 1)    {
        const float u = curand_uniform (rng);
        const float d = (1.0f + a) - 1.0f / 3.0f;
        const float c = (1.0f / 3.0f) / sqrtf(d);
        return gammaDistFloatInternal (rng, c, d) * powf(u, 1.0f / a);
    }
    else  {
        const float d = a - 1.0f / 3.0f;
        const float c = (1.0f / 3.0f) / sqrtf(d);
        return gammaDistFloatInternal(rng, c, d);
  }
}

[tnowotny]

ok ... have had a first go at this. I have done it in the fixed_TraubMiles branch for pynn_genn and a new python_wrapper_gamma branch forked off the python_wrapper branch. I am sure on the pynn_genn side I didn't quite get it right how to handle the three extraGlobalNeuronParameters ... but I thought I share what I have so far and maybe we can fix the rest together.

[tnowotny]

PS: I also was unsure whether to expect the same structure for a,b and tbins as for spikeTimes in the SpikeSourceArray case (individual values and bins for every neuron) but have assumed this for now.

[tnowotny]

need to also rescale b from units of ms to timesteps; I think it needs to be b/DT with the current convention of b being an ISI.

[neworderofjamie]

That looks awesome! The fixed_traubmiles branch was somewhat behind master so I cherry-picked your new changes into spike_source_gamma if you could switch over to that.

[tnowotny]

if you are at it - can you have a look at how the extraXXX parameters work? I also think the test test_SpikeSourceGamma may be missing the tbin parameters ...

[tnowotny]

Hum ... also, looking at the test, I suspect b is meant to be a rate, not a ISI ... we might need to flip it: b'= 1000.0/b ...