Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add gamma kernel in topology #823

Merged
merged 4 commits into from Apr 18, 2018
Merged

Conversation

@babsey
Copy link
Contributor

@babsey babsey commented Sep 8, 2017

The purpose of this PR is to implement gamma kernel in topology (See #809).

Copy link
Contributor

@heplesser heplesser left a comment

@babsey This looks nice, thank you! I added some small suggestions below. Could you create a test using this parameter?

@@ -43,6 +43,9 @@
#include "topology_names.h"
#include "topologymodule.h"

// Includes for gamma kernel:
#include <gsl/gsl_sf_gamma.h>

This comment has been minimized.

@heplesser

heplesser Oct 3, 2017
Contributor

Move to very top of all includes with comment // C includes:

@@ -458,6 +461,52 @@ class Gaussian2DParameter : public TopologyParameter


/**
* Gamma parameter p(d) = d^(kappa-1)*exp(-d/theta)/(theta^kappa*gamma(kappa))

This comment has been minimized.

@heplesser

heplesser Oct 3, 2017
Contributor

Should this be Gamma(kappa)?

raw_value( double x ) const
{
return std::pow( x, kappa_ - 1 ) * std::exp( -x / theta_ )
/ ( std::pow( theta_, kappa_ ) * gsl_sf_gamma( kappa_ ) );

This comment has been minimized.

@heplesser

heplesser Oct 3, 2017
Contributor

For efficiency, 1 / ( std::pow( theta_, kappa_ ) * gsl_sf_gamma( kappa_ ) ) should be pre-computed in the constructor so that one only needs a single multiplication here. Also, precompute beta_ = 1.0 / theta_ in the constructor, so that you avoid the division in the exp() argument.

}

private:
double kappa_, theta_;

This comment has been minimized.

@heplesser

heplesser Oct 3, 2017
Contributor

Individual declarations please, with doxygen comment of parameters.

@heplesser
Copy link
Contributor

@heplesser heplesser commented Oct 30, 2017

@babsey Have you had a chance to look at the review comments and the problems occurring on Travis?

Copy link
Contributor

@hakonsbm hakonsbm left a comment

The code looks mostly good to me, once the points by @heplesser has been addressed. However, see my comments about the problems occurring on Travis.

raw_value( double x ) const
{
return std::pow( x, kappa_ - 1 ) * std::exp( -x / theta_ )
/ ( std::pow( theta_, kappa_ ) * gsl_sf_gamma( kappa_ ) );

This comment has been minimized.

@hakonsbm

hakonsbm Oct 31, 2017
Contributor

This line is underindented, causing clang to complain.

@@ -458,6 +461,52 @@ class Gaussian2DParameter : public TopologyParameter


/**
* Gamma parameter p(d) = d^(kappa-1)*exp(-d/theta)/(theta^kappa*gamma(kappa))
*/
class GammaParameter : public RadialParameter

This comment has been minimized.

@hakonsbm

hakonsbm Oct 31, 2017
Contributor

As gsl_sf_gamma() is being used in this class, it will keep NEST from compiling if it is built without GSL. Use an #ifdef HAVE_GSL directive to make it compile without the GammaParameter-class if NEST is built without GSL.

@babsey babsey force-pushed the babsey:Topology_gamma_kernel branch 2 times, most recently from b310752 to 552d132 Dec 4, 2017
@babsey
Copy link
Contributor Author

@babsey babsey commented Dec 4, 2017

Here, the updated code and test suite for Gamma kernel in topology. Note that I am not quite skilled to write in C++ but I tried to modify codes according to @heplesser points.

Parts of equation is pre-computed in the constructor for better efficiency. I also implemented factorial function but in comment in case we want to compile it without GSL. I apologize that I am not aware of Travis testing and failed to find the source of the error.

For testing suite in pynest, I implemented the python code to compare the connection distributions with Gamma PDF. Here, I do not understand why NEST connect neurons in Gamma shaped manner for kappa + 1 (See the visual comparision).

@babsey
Copy link
Contributor Author

@babsey babsey commented Dec 4, 2017

Looking at log in Travis I found that this line (3922) probably indicates fails in Travis:
../topology/libtopology.so: undefined reference to "gsl_sf_gamma"

@hakonsbm
Copy link
Contributor

@hakonsbm hakonsbm commented Dec 5, 2017

@babsey Yes, you are right. When NEST is compiled without support for GSL, functions from that library, such as gsl_sf_gamma() that's being used in the GammaParameter-class, are unavailable. As I pointed out in my review, you need to use an #ifdef HAVE_GSL directive, which will tell the compiler to use a block of code only if NEST is compiled with GSL. That way, if there is no GSL support, you can either make the whole class unavailable, or you can implement an alternative way to do the calculation.

Additionally, line 60 and line 97 in test_topology_gamma_kernel.py are too long.

@babsey babsey force-pushed the babsey:Topology_gamma_kernel branch 3 times, most recently from e8d47fa to 4c0e7a5 Dec 6, 2017
import unittest
import nest
import nest.topology as tp
from collections import Counter

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

collections isn't used in this test.

"""Tests of gamma kernel in topology"""

def test_targets(self):
"""Object targets"""

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

It should maybe have a more descriptive name?

})

gamma1 = sts.gamma.pdf(x, a=kappa, scale=theta)
gamma_normed1 = gamma1 / float(np.sum(gamma1)) / dx

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

As these are only used when plotting for debugging, they should be moved down to the debug block, or removed.

difference = np.diff([gamma_normed2[:-1], h_normed], axis=0)
error_score = np.abs(np.median(difference))
# print('Error: %s' % error_score)
assert(error_score < 10e-5)

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

I suggest using self.assertLess(error_score, 1e-4) here, as it gives better feedback if the assert fails.

ax.step(x, gamma_normed1, label='Gamma PDF: $a = \kappa$')
ax.step(x, gamma_normed2, label='Gamma PDF: $a = \kappa + 1$')
ax.legend()
pl.show()

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

As this block is only used for debugging, I suggest removing if False and either commenting the block out, or removing it completely.

x *= i;
}
return x;
}

This comment has been minimized.

@hakonsbm

hakonsbm Dec 7, 2017
Contributor

This implementation of the factorial only works for integers. However, kappa_ can be any double value, and if kappa_ is a decimal number, we get the wrong value for delta_. I suggest instead using the Spouge's approximation. You can find an example implementation here.

@babsey babsey force-pushed the babsey:Topology_gamma_kernel branch 2 times, most recently from dac0093 to 9184791 Dec 7, 2017
@babsey babsey force-pushed the babsey:Topology_gamma_kernel branch 2 times, most recently from 966132c to f3ed4e0 Mar 5, 2018
@babsey babsey force-pushed the babsey:Topology_gamma_kernel branch from f3ed4e0 to 4ddafa0 Mar 18, 2018
@heplesser heplesser added this to the NEST 2.16 milestone Mar 19, 2018
Copy link
Contributor

@heplesser heplesser left a comment

@babsey This looks quite fine now, but see my detailed comments. I do not understand your comment in the discussion from 4 Dec 2017, where you write "Here, I do not understand why NEST connect neurons in Gamma shaped manner for kappa + 1 (See the visual comparision)." Is this still relevant? If so, could you explain it?

{
updateValue< double >( d, names::kappa, kappa_ );
updateValue< double >( d, names::theta, theta_ );
if ( kappa_ < 1 )

This comment has been minimized.

@heplesser

heplesser Mar 19, 2018
Contributor

Why the limitation to >= 1? Mathematically, only kappa > 0 is required.

"topology::GammaParameter: "
"kappa >= 1 required." );
}
if ( theta_ < 1 )

This comment has been minimized.

@heplesser

heplesser Mar 19, 2018
Contributor

Why the limitation to >= 1? Mathematically, only theta > 0 is required.

"theta >= 1 required." );
}
beta_ = 1. / theta_;
delta_ = std::pow( beta_, kappa_ ) / tgamma( kappa_ );

This comment has been minimized.

@heplesser

heplesser Mar 19, 2018
Contributor

Use std::tgamma() and explicitly include <cmath> at the beginning of the file.

h = np.histogram(d, bins=x)[0]
h_normed = h / float(np.sum(h)) / dx
difference = np.diff([gamma_normed[:-1], h_normed], axis=0)
error_score = np.abs(np.median(difference))

This comment has been minimized.

@heplesser

heplesser Mar 19, 2018
Contributor

Instead of this rather ad hoc test, why not use scipy.stats.kstest() as a proper statistical test?

This comment has been minimized.

@babsey

babsey Mar 21, 2018
Author Contributor

That is my concerning. I still do not understand, why PDF of the Gamma distribution function shows not matching results with the distribution of targets as a function of distance calculated in NEST.

Using the same values of Gamma distribution function in NEST and in Python the distributions scored a low match (Fig. 1-3, blue and orange traces). Then I increased kappa by 1 for PDF in python and compared with the distribution calculated from NEST. Here, these distributions have a better matching score (Fig 1-3, blue and green traces)

In my view, the PDFs in python were calculated correctly but I do not understand why NEST Topology gave different output.

Figures

gamma_kappa_1_theta_10
Fig 1.


gamma_kappa_2_theta_5
Fig 2.


gamma_kappa_3_theta_4
Fig 3.

This comment has been minimized.

@heplesser

heplesser Mar 21, 2018
Contributor

@babsey The number of connections you will find in NEST at a given distance depends not only on the kernel value at that distance, but also on the number of potential targets at that distance; for the test, you thus need to convolve the gamma kernel with this density and then apply the KS-test. For details, see Daniel Hjertholm, Statistical tests for connection algorithms for structured neural networks , esp. Chapter 4.

heplesser and others added 2 commits Apr 17, 2018
Revisions to the gamma parameter code
@heplesser
Copy link
Contributor

@heplesser heplesser commented Apr 18, 2018

@hakonsbm Could you take another look?

@heplesser heplesser merged commit 9df6352 into nest:master Apr 18, 2018
1 check passed
1 check passed
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@babsey babsey deleted the babsey:Topology_gamma_kernel branch Apr 18, 2018
heplesser added a commit that referenced this pull request Apr 19, 2018
Added test for spatial kernels in topology (depends on #823).
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

4 participants
You can’t perform that action at this time.