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

New neuronal population model gif_pop_psc_exp #895

Merged
merged 28 commits into from Mar 14, 2018
Merged

Conversation

@mdeger
Copy link
Contributor

@mdeger mdeger commented Feb 28, 2018

Tilo Schwalger (@schwalger), Hesam Setareh (@hesam-setareh) and I would like to add the model gif_pop_psc_exp to NEST by means of this pull request.

The model represents the pooled activity a population of (any number of) neurons in a single NEST node, accurately treating stochastic finite size effects. The algorithm and a range of simulations of single neuron vs. population activity are published in the journal article [1]:

Schwalger T, Deger M, Gerstner W (2017) Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. PLoS Comput Biol 13(4): e1005507. https://doi.org/10.1371/journal.pcbi.1005507

The model is related to the flexible single neuron model gif_psc_exp, hence the name gif_pop_psc_exp. Using this population model, one can build models of interconnected populations, etc, as discussed in [1], which can be simulated efficiently.

@mdeger
Copy link
Contributor Author

@mdeger mdeger commented Feb 28, 2018

Somehow the travis builds are failing. As far as I can see, the problem is this:

/home/travis/build/nest/nest-simulator/models/gif_pop_psc_exp.h:315:5: error: ‘GSL_BinomialRandomDev’ in namespace ‘librandom’ does not name a type

However, I cannot reproduce this error when building on my machine. Could it be that there is a different version of GSL on the travis-CI build system? Any ideas how to proceed?

@Silmathoron
Copy link
Member

@Silmathoron Silmathoron commented Mar 1, 2018

Models that require GSL need to be registered inside the GSL guard.
The error does not occur locally because you have installed gsl; however, in the tests, we do some of them without GSL, which is where your code fails.

@terhorstd terhorstd added this to the NEST 2.16 milestone Mar 5, 2018
Copy link
Contributor

@hakonsbm hakonsbm left a comment

This looks pretty good to me, once you get the GSL guard in place. I just have a few comments, see below. Also you should PEP8ify gif_pop_psc_exp.py.

#include "compose.hpp"

#include <limits>
#include <algorithm>

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

As far as I can tell, the only includes you need here are

#include "gif_pop_psc_exp.h"
#include "universal_data_logger_impl.h"
#include "compose.hpp"
, V_m_( 0.0 )
, n_expect_( 0.0 )
, n_spikes_( 0 )
, initialized_( false )

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

theta_hat_ should also be initialized here.

B_.ex_spikes_.clear(); //!< includes resize
B_.in_spikes_.clear(); //!< includes resize
B_.currents_.clear(); //!< includes resize
B_.logger_.reset(); //!< includes resize

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

Commenting once should suffice.

double theta_tmp_;
for ( int k = 0; k < P_.len_kernel_; k++ ) // InitPopulations
{
V_.n_.push_back( 0 ); // line 3

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

There are quite a few places where you refer to a line, like here. What do they mean? And can they be made clearer or more intuitive? Or removed?

{
// must be in units voltage, as q_sfa, so no division by tau!
theta_tmp_ +=
P_.q_sfa_[ j ] * std::exp( -k * V_.h_ / P_.tau_sfa_[ j ] ); // below (89)

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

Can the comment here be made clearer?

long n_spikes_; // number of spikes
double theta_hat_; // adapting threshold for non-refractory neurons

// internal switch signalling that state vectors are initialized

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

"Signaling".

@@ -86,6 +87,8 @@
#include "gif_psc_exp_multisynapse.h"
#include "gif_cond_exp.h"
#include "gif_cond_exp_multisynapse.h"
#include "gif_pop_psc_exp.h"

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

Why the extra blank line?

@@ -201,3 +202,4 @@ M_ERROR setverbosity
assert_or_die

endusing

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

Remove the extra blank line.


This script simulate an inhibitory population with the population model. After each simulation, it calculates
the mean and variance of population rate and compare them with desired values. The simulated values should be
in a acceptable margine.

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

"margin".

/unittest using

1.0 /err_mean Set % error margine for mean rate (Hz)
6.0 /err_var Set % error margine for variance of rate (Hz^2)

This comment has been minimized.

@hakonsbm

hakonsbm Mar 7, 2018
Contributor

"margin"

@@ -86,6 +87,8 @@
#include "gif_psc_exp_multisynapse.h"
#include "gif_cond_exp.h"
#include "gif_cond_exp_multisynapse.h"
#include "gif_pop_psc_exp.h"

This comment has been minimized.

@heplesser

heplesser Mar 7, 2018
Contributor

Isn't this the same line as line 81?

This comment has been minimized.

@mdeger

mdeger Mar 9, 2018
Author Contributor

Thanks. I removed line 81 instead.

Copy link
Contributor

@heplesser heplesser left a comment

@mdeger Thank you for the fast response! I think you could improve performance by making local variables const; I also had a few minor other suggestions, see inline comments.

double h_ex_ = P_.tau_m_ * ( JNA_ex_ +
h_ex_tmpvar_ / ( P_.tau_syn_ex_ - P_.tau_m_ ) );
double h_in_ = P_.tau_m_ * ( JNA_in_ +
h_in_tmpvar_ / ( P_.tau_syn_in_ - P_.tau_m_ ) );

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

The compiler may be able to optimize better if you declare these four variables const. Also, local variable names should not end with _, that suffix is only for private members.

V_.g_[ j ] = V_.g_[ j ] * V_.Q30_[ j ]
+ ( 1. - V_.Q30_[ j ] ) * V_.n_[ V_.k0_ ]
/ ( static_cast< double >( P_.N_ ) * V_.h_ ); // line 5
double g_j_tmp_ = ( 1. - V_.Q30_[ j ] ) * V_.n_[ V_.k0_ ]

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

As above: declare local helper variables const where possible and avoid the trailing _.

S_.theta_hat_ = S_.theta_hat_ + V_.Q30K_[ j ] * V_.g_[ j ]; // line 6
}

// compute free escape rate
double lambda_tld_ = escrate( S_.V_m_ - S_.theta_hat_ ); // line 8
double lambda_tld_ = escrate( S_.V_m_ - S_.theta_hat_ ); // line 8 of [1]

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

Also this could be const

This comment has been minimized.

@mdeger

mdeger Mar 9, 2018
Author Contributor

lambda_tld_ is assigned assigned a new value in line 592, so I think it cannot be const.

double P_free_ = 1 - std::exp( -0.5 * ( V_.lambda_free_ + lambda_tld_ )
* V_.h_ / 1000. ); // line 9
* V_.h_ / 1000. ); // line 9 of [1]

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

This could probably also be const. Replacing division by 1000. with multiplication by 0.001 would improve performance. Or could you even redefine V_.h_ in a way to avoid the scaling here?

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import nest

'''
Next, we set the parameters of the microscopic model

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

"Next" doen't seem to fit well for the very first step in the example.

@@ -23,20 +23,17 @@
/* Point process population model with exponential postsynaptic currents and
* adaptation */

#include "exceptions.h"
/* [1]: line numbers in comments refer to the algorithm pseudocode in
Figures 11 and 12 of the paper

This comment has been minimized.

@heplesser

heplesser Mar 9, 2018
Contributor

Please indent all lines so that they line up under "line".

Copy link
Contributor

@hakonsbm hakonsbm left a comment

Thanks for the changes. You are still missing the GSL guards, so it's still failing on Travis. And I also have a few additional remarks.

for ( long lag = from; lag < to; ++lag )
{

// main update routine, Fig. 10

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

Should be Fig. 11. Also I suggest referring to the paper.

// compute free escape rate
double lambda_tld_ = escrate( S_.V_m_ - S_.theta_hat_ ); // line 8 of [1]
double P_free_ = 1 - std::exp( -0.5 * ( V_.lambda_free_ + lambda_tld_ )
* V_.h_ / 1000. ); // line 9 of [1]

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

This exponential can be simplified to std::exp( -500. * ( V_.lambda_free_ + lambda_tld_ )* V_.h_ ).

double g_j_tmp_ = ( 1. - V_.Q30_[ j ] ) * V_.n_[ V_.k0_ ]
/ ( static_cast< double >( P_.N_ ) * V_.h_ );
V_.g_[ j ] = V_.g_[ j ] * V_.Q30_[ j ] + g_j_tmp_; // line 5 of [1]
S_.theta_hat_ = S_.theta_hat_ + V_.Q30K_[ j ] * V_.g_[ j ]; // line 6

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

You can use += here.

V_.lambda_free_ = lambda_tld_; // line 10
S_.theta_hat_ -= V_.n_[ 0 ] * V_.theta_tld_[ 0 ]; // line 11

for ( int l = 0; l < P_.len_kernel_; ++l )

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

Don't use l as a variable.

// which is a recordable
double theta_hat_ = S_.theta_hat_;
double theta_;
for ( int l = 0; l < P_.len_kernel_ - V_.k_ref_; ++l ) // line 13 of [1]

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

Don't use l as a variable. I suggest using something like k_marked here to also be closer to the pseudocode from the paper.

stochastic population rate dynamics derived in the paper
[Schwalger et al. PLoS Comput Biol. 2017]. The stochastic population dynamics
is implemented in the NEST model gif_pop_psc_exp. We demonstrate this model
using the example of a Brunel network of two coupled populations, one exhitory

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

"excitatory"


nest_pops = nest.Create('gif_pop_psc_exp', M)

C_m = 250. # irrelavant value for membrane capacity, cancels out in simulation

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

"irrelevant"

# record all spikes from population to compute population activity
nest.Connect(nest_pops[i], nest_sd[i], 'all_to_all')

# for each population i the first Nrecord[i] neurons are recorded

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

"in"

ResetKernel
0 <<
/resolution res
>> SetStatus

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

Shift >> SetStatus one space to the right, to align it with its matching <<.

err_mean leq assert_or_die

var_rate expected_var sub abs
err_var leq assert_or_die

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

Remove trailing whitespace.

, I_syn_ex_( 0.0 )
, I_syn_in_( 0.0 )
, V_m_( 0.0 )
, theta_hat_( 0.0 ) // initialization value has no effect for theta_hat_

This comment has been minimized.

@hakonsbm

hakonsbm Mar 9, 2018
Contributor

theta_hat_ should be initialized after n_spikes_ to avoid warnings during compilation.

@mdeger
Copy link
Contributor Author

@mdeger mdeger commented Mar 9, 2018

Thank you very much for your comments. I have just committed some changes in order to address them.

Copy link
Contributor

@heplesser heplesser left a comment

@mdeger Thanks a lot! I just have one detail left, namely skipping GSL-dependent test using our new skip_if_without_gsl SLI function.

@@ -38,6 +38,9 @@ SeeAlso: gif_pop_psc_exp, testsuite::test_pp_pop_psc_delta

*/

% This test should only run if we have GSL
statusdict/have_gsl :: not {statusdict/exitcodes/success :: quit_i} if

This comment has been minimized.

@heplesser

heplesser Mar 10, 2018
Contributor

Please move this check to after (unittest) run and replace it with skip_if_without_gsl. This is not only simpler, but it will also mark the test as skipped due to lack of GSL instead of implicitly counting it as passed.

@hakonsbm
Copy link
Contributor

@hakonsbm hakonsbm commented Mar 12, 2018

@mdeger Thank you for the changes! This looks good now.

Sorry for the erroneous comment. I have created a PR in your fork reinstating the triple quotes. As soon as you merge that, I will approve this PR.

Reinstated triple quotes for documentation comments
@heplesser heplesser removed the request for review from Silmathoron Mar 14, 2018
@heplesser heplesser merged commit 9230ddd into nest:master Mar 14, 2018
1 check passed
1 check passed
continuous-integration/travis-ci/pr The Travis CI build passed
Details
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

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