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

Implement recording of current generated by stimulating devices #663

Merged
merged 14 commits into from
Apr 18, 2017
Merged

Implement recording of current generated by stimulating devices #663

merged 14 commits into from
Apr 18, 2017

Conversation

appukuttan-shailesh
Copy link
Contributor

@appukuttan-shailesh appukuttan-shailesh commented Feb 20, 2017

@heplesser : I was running the tests (make installcheck) and got the following output:

ERROR: test_synapse_deletion_one_to_one_no_sp (nest.tests.test_sp.test_disconnect.TestDisconnectSingle)
TypeError: allgather() takes exactly 1 positional argument (2 given)
FAILED (SKIP=11, errors=1)

PyNEST tests: 548
Passed: 547
Failed: 1

NEST Testsuite Summary

NEST Executable: /home/shailesh/opt/nest/bin/nest
SLI Executable : /home/shailesh/opt/nest/bin/sli
Total number of tests: 772
Passed: 771
Failed: 1 (1 PyNEST)

The link to the reports archive.
Do let me know how to proceed.

@heplesser heplesser added ZC: Model DO NOT USE THIS LABEL I: No breaking change Previously written code will work as before, no one should note anything changing (aside the fix) ZP: PR Created DO NOT USE THIS LABEL S: Normal Handle this with default priority T: Enhancement New functionality, model or documentation labels Feb 20, 2017
@heplesser
Copy link
Contributor

@appukuttan-shailesh I will provide a more detailed review later. But not that you cannot record from the noise generator as you do in your PR. The noise generator provides a different current for each of its targets. Your code records the current provided to the "last" target, which in the current implementation would be the one to which the generator is connected last on a given thread. The only consistent thing possible for the noise generator would be to return the mean current.

@heplesser
Copy link
Contributor

@appukuttan-shailesh The test failure seems strange, since allgather() according to all documentation I found should take two arguments. The test also works on Travis. Could you check your Python implementation?

The travis checks failed because of formatting errors. To see the detailed report, click on "Details" above and then look at the end of the report for, e.g., the last test case (....7).

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : I was wondering if we could do either of the following for the noise generator:

  1. If only one target (B_.amps_.size()==1), then have the actual current, else the mean current (as you suggested)
  2. The instantaneous average of the current evaluated for all the targets

Both the above approaches would have the benefit of giving the actual current when only one target is present, though I am not totally sure if the second approach is technically sound (with regards to being true to gaussian white noise). Let me know, and I shall proceed accordingly.

The advantage of having a true current profile (atleast for just a single target) is with regards to comparisons across simulators (e.g. using PyNN; my main focus currently) and debugging.

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : I was checking out the formatting errors, and I believe most of it was applicable to the original code. For example:

models/noise_generator.cpp:

  • vera++:
    • negation operator used in its short form : 4

All the 4 instances correspond to statements in the original code and not modified here. I am not sure whether I should (try) to fix all of those. Do let me know. Also, about the approach to follow for the noise generator (described in the previous post).

@heplesser
Copy link
Contributor

@appukuttan-shailesh NEST's static code analysis checks all files that differ from master in your pull request. Thus, if you even just correct a spelling error in a comment in abc.cpp, the entire file will be checked and you will need to fix all issues reported. In this way, we will slowly bring the entire code base into shape. We don't have the resources to do this centralized in one go. We have just improved the format of the report provided (#546), so if you pull recent changes from master into your branch, you can enjoy the new format.

By running extras/check_code_style.sh in your source directory, you can run the static code checks on your own computer before pushing to Github. You should commit before running the check, since check_code_style.sh by default checks files differing between master and head.

@heplesser
Copy link
Contributor

@appukuttan-shailesh Now to your questions about the code itself.

  1. If only one target (B_.amps_.size()==1), then have the actual current, else the mean current (as you suggested)
  2. The instantaneous average of the current evaluated for all the targets

Doing different things depending on the number of targets would not be a good thing, the risk of user errors (or rather, misinterpretations) is too big and I would discourage this. The instantaneous average would be safer (and for just one target it would be just the current for the one target).

Both the above approaches would have the benefit of giving the actual current when only one target is present, though I am not totally sure if the second approach is technically sound (with regards to being true to gaussian white noise).

The noise generator does not generate truly gaussian white noise, since it only switches amplitude once per time step; see doc/model_details/noise_generator.ipynb for more).

The advantage of having a true current profile (atleast for just a single target) is with regards to comparisons across simulators (e.g. using PyNN; my main focus currently) and debugging.

One way to make currents fully recordable would be to extend parrot_neuron to support CurrentEvent as well as SpikeEvent. One could then pass the current through one parrot neuron per actual target and record from the parrot neuron. We do this in some tests when we need complete information about poisson spike trains injected into neurons. But this would hit performance, so it would make sense mostly for testing and debugging.

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : Thanks for those inputs. I shall go ahead with the instantaneous average approach for the noise generator. I am away this coming week, but shall complete this fix along with the formatting issues in the following week.

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : I believe this is done. The formatting has been taken care of for the files involved.

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : I was wondering if you had chance to check this out. I wasn't able to join in for the NEST developers conference as we have another project call at the same time; so not sure if there has been any update.

Copy link
Contributor

@heplesser heplesser left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@appukuttan-shailesh Well done and sorry for the delay. I added a number of comments in the code, mostly small stuff that repeats for most models.

for ( long lag = from; lag < to; ++lag )
{
B_.logger_.record_data( origin.get_steps() + lag );
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

record_data should be called at the end of the time step, because NEST works on left-open, right-close intervals. So you need to move this line after the if-block below.

if ( device_.is_active( Time::step( start + lag ) ) )
{
const double y_0 = S_.y_0_;
S_.y_0_ = V_.A_00_ * y_0 + V_.A_01_ * S_.y_1_;
S_.y_1_ = V_.A_10_ * y_0 + V_.A_11_ * S_.y_1_;
S_.I_ = S_.y_1_ + P_.offset_;

CurrentEvent ce;
ce.set_current( S_.y_1_ + P_.offset_ );
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both for clarity and efficiency, you could use S_.I_ as argument here.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But since S_.I_ can be removed entirely, you can just drop line 225.

@@ -108,6 +125,7 @@ class ac_generator : public Node
{
double y_0_;
double y_1_;
double I_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a doxygen comment on the parameter (either /** ... */ above or //!< ... at end of line. It is important to make clear that you need I_ to be able to record the current. But see for another idea below, I think you do not need I_ at all.

double
get_I_() const
{
return S_.I_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you replace this by

return S_.y1_ + P_.offset_;

you don't need S_.I_. Avoiding an additional state variable would be good.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser isn't it contradictory to the remark on line 228?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gtrensch No. If you remove S_.I_ entirely, you would also delete line 225 above.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser this is true, if you leave line 228 as is, correct?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, 228 then stays as it is.

Copy link
Contributor Author

@appukuttan-shailesh appukuttan-shailesh Apr 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser : I am not sure this would work properly (without additional changes). The problem would be before t_start and after t_stop, where S_.y_1_ is not being set to zero.

The option might be to explicitly set S_.y_1_ to zero outside the required time period. But then we have the problem of P_.offset_, which will always be reflected in the recorded current. Unless we put in a check for time in get_I_().

If so, it might be better to go ahead with an extra state variable (S_.I_), as the implementation would also be more uniform across the various electrode classes (as others have this variable).

Copy link
Contributor Author

@appukuttan-shailesh appukuttan-shailesh Apr 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried incorporating the changes for AC and DC (removing S_.I_). This is what I got:
output_all_bad
There are issues outside of the current injection period.
(ignore the noisy input - that was due to a mistake)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@appukuttan-shailesh I am sorry, I had not thought about the fact that the device might be inactive. To account for that, we need S_.I_. It might be a good idea to explain this in the doxygen comment for that variable.

long start = origin.get_steps();

CurrentEvent ce;
ce.set_current( P_.amp_ );

for ( long offs = from; offs < to; ++offs )
{
B_.logger_.record_data( origin.get_steps() + offs );
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should be called after if block.

Copy link
Contributor Author

@appukuttan-shailesh appukuttan-shailesh Apr 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was wondering if the same should also be true for 'sinusodial_poisson_generator.cpp', line 268:. That's what I had based these changes on.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In sinusoidal_poisson_generator, the recorded variable S_.rate_ is updated in lines 256-265, so once we get to 268, it has its value at the end of the time step; but the logging in that model ignores the possibility that the rate is zero while the device is inactive. Could you correct that there (and in the corresponding gamma generator?).

The line needs to move down to after the if so that it records the correct value (0 or amp) for the end of the time step.

@@ -125,16 +126,29 @@ class noise_generator : public Node
* @see Technical Issues / Virtual Functions: Overriding, Overloading, and
* Hiding
*/
using Node::handle;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the user documentation above (I cant comment on those lines), you need to explain the recordable I, since it is an average.

@@ -186,6 +201,7 @@ class noise_generator : public Node
{
double y_0_;
double y_1_;
double I_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe use I_avg_ as variable name and add a Doxygen comment.


for ( long offs = from; offs < to; ++offs )
{
const long curr_time = t0 + offs;

B_.logger_.record_data( origin.get_steps() + offs );
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Move this after the if ( device_.is_active( ... ) ) block.

struct State_
{
double y_0_;
double y_1_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need for y_0_ and y_1_ here.

{
double y_0_;
double y_1_;
double I_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add doxygen comment explaining why we need this and what it is.

Copy link
Contributor

@gtrensch gtrensch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did some functional tests by running the 'nest_test_records_all.py' script which is mentioned in the related issue. Everything works well. I added some minor comments. If those are addressed, I approve the PR.

@@ -97,6 +112,8 @@ class ac_generator : public Node
double phi_deg_; //!< Phase of sine current (0-360 deg)

Parameters_(); //!< Sets default parameter values
Parameters_( const Parameters_& );
Parameters_& operator=( const Parameters_& p ); // Copy constructor EN
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does "EN" have a special meaning?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The initials of the person who wrote the comment in 2009 which got here by copy-paste from another model (really looking forward to NESTML ...). The comment is also incorrect---this is the assignment operator, not the copy constructor. @appukuttan-shailesh You can just delete the comment, the code speaks for itself. Same applies in other models.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removing this comment from sinusoidal_poisson_generator and sinusoidal_gamma_generator as well.

@@ -97,6 +112,8 @@ class ac_generator : public Node
double phi_deg_; //!< Phase of sine current (0-360 deg)

Parameters_(); //!< Sets default parameter values
Parameters_( const Parameters_& );
Parameters_& operator=( const Parameters_& p ); // Copy constructor EN

void get( DictionaryDatum& ) const; //!< Store current values in dictionary
void set( const DictionaryDatum& ); //!< Set values from dicitonary
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo: dicitonary -> dictionary

@@ -116,6 +134,24 @@ class ac_generator : public Node

// ------------------------------------------------------------

// These friend declarations must be precisely here.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also this comment arrived from another model, where I put it in 2009. I am not sure why I wrote it then, I believe, looking at other models, that the declarations could also be placed above the class Parameter {}; declaration. iaf_psc_alpha.h provides a more sensible comment:

// The next two classes need to be friends to access the State_ class/member

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am updating the original comments to match the one from iaf_psc_alpha.h

double
get_I_() const
{
return S_.I_;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@heplesser isn't it contradictory to the remark on line 228?

@@ -106,15 +121,57 @@ class dc_generator : public Node
double amp_; //!< stimulation amplitude, in pA

Parameters_(); //!< Sets default parameter values
Parameters_( const Parameters_& );
Parameters_& operator=( const Parameters_& p ); // Copy constructor EN
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does "EN" have a special meaning?

@@ -174,6 +188,7 @@ class noise_generator : public Node

Parameters_(); //!< Sets default parameter values
Parameters_( const Parameters_& );
Parameters_& operator=( const Parameters_& p ); // Copy constructor EN
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does "EN" have a special meaning?

@@ -106,15 +121,57 @@ class dc_generator : public Node
double amp_; //!< stimulation amplitude, in pA

Parameters_(); //!< Sets default parameter values
Parameters_( const Parameters_& );
Parameters_& operator=( const Parameters_& p ); // Copy constructor EN

void get( DictionaryDatum& ) const; //!< Store current values in dictionary
void set( const DictionaryDatum& ); //!< Set values from dicitonary
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo: dicitonary -> dictionary


// ------------------------------------------------------------

// These friend declarations must be precisely here.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why?

if ( device_.is_active( Time::step( start + lag ) ) )
{
const double y_0 = S_.y_0_;
S_.y_0_ = V_.A_00_ * y_0 + V_.A_01_ * S_.y_1_;
S_.y_1_ = V_.A_10_ * y_0 + V_.A_11_ * S_.y_1_;
S_.I_ = S_.y_1_ + P_.offset_;

CurrentEvent ce;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line was originally intentionally outside the for-loop so that the event object needs to be constructed only once. I believe that the compiler cannot automatically pull this definition out of the loop since ce is manipulate inside the loop. Therefore, the line should move back before the loop.

Copy link
Contributor Author

@appukuttan-shailesh appukuttan-shailesh Apr 5, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had moved it inside, thinking that the object be constructed only if it was going to be used (device active). This, I think, is also how it was coded in sinusoidal_poisson_generator.cpp.
But there I suppose a choice had to be made between DSSpikeEvent and SpikeEvent, and thus was probably a preferable approach.

I shall revert this change in ac_generator.cpp

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented Apr 5, 2017

Have incorporated all the changes discussed above. Do let me know if I have missed anything (despite double-checking). I have re-tested for the four stimulating devices described in the initial issue, but haven't run any specific tests for the minor edits for sinusoidal_gamma_generator and sinusoidal_poisson_generator.

@appukuttan-shailesh
Copy link
Contributor Author

Have fixed the formatting errors.

@gtrensch
Copy link
Contributor

gtrensch commented Apr 5, 2017

Looks good to me!

Copy link
Contributor

@heplesser heplesser left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@appukuttan-shailesh This looks very fine now, just one typo in a comment. But could I ask you to also create a test checking that recording works for all these generators? This would in particular ensure that we do not break things by future changes. We prefer test written in SLI (in testsuite/unittests), but since SLI is a bit exotic tests written in Python (in pynest/nest/tests) are fine as well.

@@ -96,6 +96,10 @@ frequency double - Frequency of sine modulation in Hz

To obtain comparable results for different values of dt, you must
adapt std.
- As The noise generator provides a different current for each of its targets,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"the" with lower-case "t"

Copy link
Contributor Author

@appukuttan-shailesh appukuttan-shailesh Apr 5, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah... I missed that typo! ... I can try my hand at quickly writing a test case. I shall commit the typo change and the test case in a couple of hours. Keen to have this done and dusted today itself.

@heplesser : For the test, I presume a simple check to ensure that the current vectors are being populated should suffice? With the length of the current vector being equal to the length of the membrane potential vector. Does that seem sufficient?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

naming the test file: test_current_recording_generators.py

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

p.s. and maybe just a check to ensure that current = 0 before t_start and after t_end

@appukuttan-shailesh
Copy link
Contributor Author

In addition to the earlier comments, I had a few other concerns:

  1. In the test, I have set time: origin = 5.0 ms, start = 2.5 ms. I would have expected the current to be injected at 7.5 ms. But on the plot, it initiates at 7.4 ms. What could be the reason?
  2. Also, the corresponding change in Vm initiates at 8.5 ms, i.e. after a delay of 1 ms. This delay corresponds to the min_delay for the simulation. I suppose this is fine?

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : Travis has finally given the all clear. If the comments in the earlier posts don't demand additional changes, this is good to go.

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented Apr 7, 2017

@heplesser Not sure why it shows these conflicts today... (was all ok yesterday). Is this related to the "Code style systematized" email that was sent out today? Do I need to do anything about this? Most of the issues seem to be related to spacing and location of comments.

@heplesser
Copy link
Contributor

@appukuttan-shailesh This is due to the merge of #691, which changed code style in many places (see my mail to NEST User). The best way to fix this is to pull master into your branch. This will probably lead to some (minor) conflicts in the files listed above, which you need to fix. Let me know if you run into any trouble!

@heplesser
Copy link
Contributor

@appukuttan-shailesh Concerning the time at which you first observe output current from the dc generator, it is indeed a bit peculiar, but with good reason. Consider the following example (I drop the voltmeters here to keep the code short):

n_dc = nest.Create('iaf_psc_alpha')
n_ie = nest.Create('iaf_psc_alpha')

d = nest.Create('dc_generator', params={'amplitude': 100., 'start': 5., 'stop': 10.})
nest.Connect(d, n_dc)

nest.Simulate(5.+1.)
nest.SetStatus(n_ie, {'I_e': 100.})
nest.Simulate(5.)
nest.SetStatus(n_ie, {'I_e': 0.})
nest.Simulate(5.)
  • As reference, we consider neuron n_ie, which receives no external input.
    • This neuron is simulated for 6 ms, then an internal current is turned on at time 6.0 ms, and turned off again at time 11.0 ms.
    • Since that current "kicks in" at 6.0 ms, we expect to see the first effect on the membrane potential at 6.1 ms when recording at 0.1 ms sampling intervals.
    • Similarly, we expect the first effect of the current being turned off at 11.1 ms.
    • This is indeed the case.
  • A neuron simulated by an external generator should show exactly the same response.
    • n_dc receives input from a dc_generator, starting at 5 ms and stopping at 10 ms.
    • Since the generator is connected with default delay 1 ms, the effect on the neuron should last from 6 ms to 11 ms, just as for the reference case
    • This is the case
  • We actually have a unittest for this

Since NEST models always handle current input at the end of a time step (the precise reason escapes me right now, it is related to working on left-open, right-closed intervals), this requires that any current injecting devices must be pro-active, i.e., turning on/off their output one time step earlier than one would expect. Therefore you see the current in the recording already at 7.4, I believe (I haven't tested with your branch).

@appukuttan-shailesh
Copy link
Contributor Author

@heplesser Thanks for the detailed explanation. Makes sense now.

I resolved the conflicts and committed the changes. Travis gave me an all clear for six of the jobs, but the last one has failed with:
"The command "sudo apt-get install -y build-essential cmake libltdl7-dev libreadline6-dev libncurses5-dev libgsl0-dev python-all-dev python-numpy python-scipy python-matplotlib ipython" failed and exited with 100 during .
Your build has been stopped."

Not sure how to proceed. I don't imagine it has anything to do with the submitted changes.

@heplesser
Copy link
Contributor

@appukuttan-shailesh It looks like Travis had some trouble downloading software packages from some servers. I have restarted the 7th job now. We are currently looking at more robust solutions.

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented Apr 11, 2017

@heplesser ... finally an all clear from Travis.

@heplesser
Copy link
Contributor

This PR adds the possibility to record current emitted by current generators; see model documentation for details. This is mainly added for PyNN support.

@heplesser heplesser merged commit c0cb633 into nest:master Apr 18, 2017
@appukuttan-shailesh
Copy link
Contributor Author

@heplesser : A small correction in what I posted earlier regarding the timing. When I have set time: origin = 5.0 ms, start = 2.5 ms, the electrode current does infact initiate at 7.5 ms (and not 7.4 ms as I reported earlier.... there was a minor mistake in my code), with the effect on the target neuron at 8.5 ms (after min_delay).

Further to this, I had a few minor queries. On recording the currents, it is seen that we get no data for t = 0. The recorded data begins from t = 0.1 ms (timestep = 0.1 ms).

  1. If in a model, I wish to inject current (using a stimulating device) from t = 0... would this be possible in NEST? Or would the current injection into the cell be actually initiated only at t = timestep + min_delay? (... the query is not for practical purposes, but just to be sure how it works).

  2. Would it be fine to define that at t = 0, all injected currents are zero?

  3. For my final query, it might be easier to present an example.
    Consider the following:
    ac = nest.Create('ac_generator', 1, params = {'amplitude': 550.0, 'offset': 100.0, 'frequency': 100.0, 'phase' : 0.0, 'origin' : 2.5, 'start' : 2.5, 'stop' : 40.0})

The values that are computed are:
t (ms) <-> I (pA)
4.9 <-> 0.0
5.0 <-> 134.53
5.1 <-> 168.93
5.2 <-> 236.77

Shouldn't the value at t = 5.0, be I = 100.00.... i.e. offset + sin(t=0). This is what I might have expected intuitively, and also seems to be the case in Brian and NEURON. The absence of this value (100.0), I believe, causes a more abrupt change in current at the start of current injection as seen below, which is absent in the next cycle (starting at 15 ms).
ac_jump

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented May 3, 2017

I had one other concern regarding my present implementation of recording for noisy current. For a simulation with timestep = 0.1 ms and noise_generator's dt = 1.0 ms, we get:
noisy_before
<< Code for testing >>
i.e. values are recorded only at the instants when the current changes (here every 1.0 ms), and the plot basically joins these points via straight lines. This gives a false impression of the current value in between these changes (due to linear interpolation). The correct representation, I believe, would have been to keep the current value constant for the duration of noise_generator's dt (here 1.0 ms).

Also, the present implementation causes the length of the recorded current and membrane potential vectors to be different when simulation timestep is not equal to the electrode's dt.

Changing https://github.com/nest/nest-simulator/blob/master/models/noise_generator.cpp#L328 to https://github.com/nest/nest-simulator/blob/master/models/noise_generator.cpp#L348 as follows will achieve this:

// >= in case we woke from inactivity
if ( now >= B_.next_step_ )
{
  // compute new currents
  for ( AmpVec_::iterator it = B_.amps_.begin(); it != B_.amps_.end();
        ++it )
  {
    *it = P_.mean_
      + std::sqrt( P_.std_ * P_.std_ + S_.y_1_ * P_.std_mod_ * P_.std_mod_ )
        * V_.normal_dev_( kernel().rng_manager.get_rng( get_thread() ) );
  }
  // use now as reference, in case we woke up from inactive period
  B_.next_step_ = now + V_.dt_steps_;
}

// record values
for ( AmpVec_::iterator it = B_.amps_.begin(); it != B_.amps_.end();
      ++it )
{
  S_.I_avg_ += *it;
}
S_.I_avg_ /= B_.amps_.size();
B_.logger_.record_data( origin.get_steps() + offs );

DSCurrentEvent ce;
kernel().event_delivery_manager.send( *this, ce, offs );

This gives:
noisy_after

@heplesser : If you agree with this, I shall make the changes.

@gtrensch
Copy link
Contributor

gtrensch commented May 4, 2017

@appukuttan-shailesh @heplesser

Please can you have a look at the following branch/fix ? https://github.com/gtrensch/nest-simulator/tree/bugfix_get_recordables

I went into an issue with getting the "recordables" from the generators.
e.g. print(nest.GetDefaults('step_current_generator')['recordables'][0])
I fixed this in the branch mentioned above. I also enhanced your unit test test_current_recording_generators.pyaccordingly.
My changes caused the SLI test test_multimeter_support.sli to fail. The test checks all models providing "recordables" and verifies the size of the recorded data vector. This failed for the noise_generator due to an invalid time regime.
'times': array([ 1.1, 2.1, 3.1, 4.1, 5.1, 6.1, 7.1, 8.1]) instead of
'times': array([ 1., 2., 3., 4., 5., 6., 7., 8., 9.])
The problem is exact the issue you describe. So I added our change as well.

You may pull the changes from my branch and re-open this PR or we can also create a new PR.

I discovered another issue which also occurs with the noise generator only. The data recorded from the noise generator (no params set up) is:
{'times': array([ 1., 2., 3., 4., 5., 6., 7., 8., 9.]), 'I': array([ nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'senders': array([1, 1, 1, 1, 1, 1, 1, 1, 1])}
If the noise generator is additionally connected to a neuron it becomes:
{'senders': array([1, 1, 1, 1, 1, 1, 1, 1, 1]), 'I': array([ 0., 0., 0., 0., 0., 0., 0., 0., 0.]), 'times': array([ 1., 2., 3., 4., 5., 6., 7., 8., 9.])}

I would expect zeros in all cases for the current.

@appukuttan-shailesh
Copy link
Contributor Author

@gtrensch : Thanks for looking into this. I shall try to fix that issue and re-open the PR.

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented May 5, 2017

@gtrensch : I have made a small change. See here: https://github.com/appukuttan-shailesh/nest-simulator/blob/record-stim/models/noise_generator.cpp [Lines 343 to 356]. I believe this fixes the above mentioned issue. Could you check and confirm.

@heplesser
Copy link
Contributor

@appukuttan-shailesh It is a bit unfortunate to continue discussion in a merged (and thus closed PR), it would be better to open a new issue or PR for discussion. That also allows for line-by-line comments on the code. Would you open a new PR?

You are right that the B_.logger_.record_data() call must come after and outside the if ( now >= ... ) block. But you do not need to recalculate S_.I_avg_ for every time step, if you treat the variable as a proper state member; you just need to be careful then not to set it to 0 at the beginning of each pass through the loop and to set it to 0 when the device becomes inactive.

You can also avoid the division by zero problem more easily:

if ( not B_.amps_.empty() )
{
  S_.I_avg_ /= B_.amps_.size();
}

or

S_.I_avg_ /= std::max( 1, B_.amps_.size() ); 

which is correct because S_.I_avg_ == 0 if B_.amps_ is empty.

Since all transmission in NEST has a delay (except gap junctions, but that is another topic), the earliest time at which a current from a current generator can begin to affect a neuron is after min_delay.

Concerning the start of the sine, this may again be related to the need of devices to be proactive, see my comment from April 7.

@gtrensch
Copy link
Contributor

gtrensch commented May 8, 2017

@appukuttan-shailesh, @heplesser:
Even though I suggested to re-open this PR, I agree with @heplesser. We should follow up this in a new issue/pull request. This is not only because one meanwhile told me that merging a second time might not be possible (I never tried this before), it also turned out that there is another issue with this implementation.

While looking at the latest changes (initializing B_.amps_ did not work well in my tests) and after a discussion with @jougs, we came to the conclusion that the recording is not correct. To get the array initialized properly in init_buffers(), num_targets needs to be increased. This explains why it works in my tests when neurons are connected (B_.amps_.size() will then be > 0). A multimeter connection does not increase num_targets. This could be easily changed but lead us to the question: Which connection do we actually record ? The current implementation does not take into accout that for the noise_generator a buffer per connection exists. We need to specify the connection we want to record or record everything.

I would suggest the following:
My branch already contains the bugfix for the recordables and @appukuttan-shailesh fix for the time regime. Pls. see https://github.com/gtrensch/nest-simulator/tree/bugfix_get_recordables
I can create a PR for that. We open an issue for the problem described above to make a design descission and solve this seperatly.

@appukuttan-shailesh
Copy link
Contributor Author

appukuttan-shailesh commented May 8, 2017

@heplesser : Sorry about that. With my initial post, I wasn't certain if this required a fix, and so was intended more as a query. But I see that this now requires a new PR. I suppose @gtrensch plans to create it. I like the fix that you have suggested... its way neater than what I had in place.

Regarding the AC current, I believe this requires a fix. The current implementation causes an initial jump in Vm as well. If I adjust the phase to compensate for this, the resultant Vm (and I_inj) is as expected (and also matching outputs from other simulators). I don't have my work pc currently, so cannot post examples. Let me know if you think I should open a new issue for this.

@gtrensch : Yes, regarding what should be recorded, might require further discussion. Currently, for the noise generator, its simply the average of all instantaneous currents it generates; without any option to specify the current generated for a specific target (we don't have the problem for other generators, I believe, because its a single common value for each target). We would need to add the fix suggested above by @heplesser in your branch. Let me know if you want me to work on it.

@gtrensch
Copy link
Contributor

gtrensch commented May 9, 2017

@appukuttan-shailesh : It would be wonderful if you could proceed with this. I needed the step_current_generator recording functionality for a project. Unfortunately, I went into some difficulties and started resolving them. However, for the sake of simplicity I would suggest to pull this from my repository: master...gtrensch:bugfix_get_recordables
As you already mentioned L348 needs to be replaced by S_.I_avg_ /= std::max( 1, B_.amps_.size() ); as @heplesser suggested. I think this line additionally requires a cast.

This will add the following:

  • Reading the recordables: GetDefaults('..._generator')['recordables']
  • Enhances your Python unit test by testing reading the recordables.
  • SLI test test_multimeter_support.sli will recognize the generator recordables and run the test against the generators as well and the test will not fail anymore.

The behavior would be:

  • The average value is recorded.
  • If there is no neuron connected, the values will be zero.

This is probaly a point for discussion. Hard to say for me in how far this behavior is useful.

@jougs jougs changed the title Fixes #658: Implements recording of current generated by stimulating devices Implement recording of current generated by stimulating devices Jun 6, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
I: No breaking change Previously written code will work as before, no one should note anything changing (aside the fix) S: Normal Handle this with default priority T: Enhancement New functionality, model or documentation ZC: Model DO NOT USE THIS LABEL ZP: PR Created DO NOT USE THIS LABEL
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants