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

Adding inhomogeneous poisson generator model #671

Merged
merged 15 commits into from
Mar 10, 2018

Conversation

zbarni
Copy link

@zbarni zbarni commented Mar 4, 2017

Model for generating Poisson spikes with varying rate. Authorship shared with Renato Duarte (@rcfduarte).

@terhorstd terhorstd 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 Mar 6, 2017
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.

@zbarni This looks quite fine, but I have some comments in the code, including one which is serious, I believe (in the update() method). You should also include a test, including corner cases.

const bool ut =
updateValue< std::vector< double_t > >( d, names::rate_times, rate_times_ );
const bool uv = updateValue< std::vector< double_t > >(
d, names::rate_values, rate_values_ );
Copy link
Contributor

Choose a reason for hiding this comment

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

ut and uv are rather too short and uninformative.

Copy link
Author

Choose a reason for hiding this comment

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

done

}

// ensure amp times are strictly monotonically increasing
if ( rate_times_.empty() == false )
Copy link
Contributor

Choose a reason for hiding this comment

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

if ( not rate_times_.empty() )

Copy link
Author

Choose a reason for hiding this comment

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

done

// random number generator
librandom::RngPtr rng = kernel().rng_manager.get_rng( get_thread() );


Copy link
Contributor

Choose a reason for hiding this comment

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

Remove empty line.

Copy link
Author

Choose a reason for hiding this comment

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

done

== Time( Time::ms( P_.rate_times_[ B_.idx_ ] ) ).get_steps() )
{
B_.rate_ = P_.rate_values_[ B_.idx_ ] / 1000.0; // scale the rate to ms^-1
B_.idx_++;
Copy link
Contributor

Choose a reason for hiding this comment

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

++B_.idx_;, we generally prefer prefix operators since they are faster for complex data types.

Copy link
Author

Choose a reason for hiding this comment

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

done

// on class SimulatingDevice.
if ( B_.idx_ < P_.rate_times_.size()
&& curr_time + 1
== Time( Time::ms( P_.rate_times_[ B_.idx_ ] ) ).get_steps() )
Copy link
Contributor

Choose a reason for hiding this comment

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

What will happen here if the rate changes twice within a time step? Since you can only increase B_.idx_ by one for each value of offs, you would probably get out of step for the remainder of the simulation.

Note that this is a deep problem: which rate should you use if the rate changes multiple times within an interval. Indeed, the question is also how to handle rate changes inside an interval correctly. This is especially critical to ensure consistent results when simulating the same network with different resolutions. The sinusoidal_poisson_generator always uses the rate at the end of the time step, so a consistent solution would be to always use the rate valid at then end of the step for the entire step. This needs to be documented, and you need code here that skips far enough.

Copy link
Contributor

Choose a reason for hiding this comment

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

I agree to @heplesser, you will get out of step if the rate_times are not on the time-grid of the simulation. For example you will increase B_.idx_ only once within the dmin loop, if you specify the rate_times with a time-interval equal to half of the simulation time-step. I anyway suggest that the user can only specify rate_times which lie on the simulation time grid. However, the simulation time grid can change after the device has been specified and parameters have already been given. Then a warning should be raised. @heplesser, do you have an idea how to handle this?

Copy link
Contributor

Choose a reason for hiding this comment

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

@jschuecker NEST prohibits changes in resolution after the first network node has been created, so resolution changes are not an immediate problem. Simulating the same network at different resolutions may lead to different rate profiles, though, unless all times at which the rates change are multiples of all resolutions used.

Forcing the times of rate changes onto the resolution time grid is a sensible idea. We do this in spike_generator, see its Parameters::set() method for an implementation. That model also shows how to optionally disable the grid-forcing.

Choose a reason for hiding this comment

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

So far, when using this model, we have always made sure that all the rate_times are multiples of the simulation resolution. I guess the simplest way to address this issue would be to raise an error if this is not the case, or, as @heplesser suggests, force the rate_times onto the time grid (issuing a warning, so that the user knows the chosen rate_times will be changed)...

Copy link
Contributor

Choose a reason for hiding this comment

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

When writing code for a widely used simulator, one always needs to keep an eye on protecting users from doing things with untoward consequences. I would suggest to implement the grid-forcing as in spike_generator. That would also provide consistency between these two generators that can be "fed" arrays of times.

Copy link
Author

Choose a reason for hiding this comment

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

Sorry for the long delay, will try to implement all requests / suggestions as soon as possible.. @heplesser should grid enforcement be optional as in spike_generator, controlled via the precise_times and allow_offgrid_times flags or not?

if ( B_.rate_ > 0 && device_.is_active( Time::step( t0 + offs ) ) )
{
DSSpikeEvent se;
kernel().event_delivery_manager.send( *this, se, offs );
Copy link
Contributor

Choose a reason for hiding this comment

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

This code sends different spike trains to each target. To send the same train to all targets, you need additional code, see, e.g., sinusoidal_poisson_generator.

Choose a reason for hiding this comment

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

Should this be included? It is generally desirable that each target receives an independent realization..

Copy link
Contributor

Choose a reason for hiding this comment

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

@rcfduarte Sending the same train to all targets is not strictly necessary. One can always achieve the same effect by connecting the generator to a single parrot_neuron and then connecting that to all neurons that should receive identical spike trains. It was just a suggestion, since we have this feature in sinusoidal_{poisson,gamma}_generator. The reason for implementing it there was that we use those generators often to represent retinal ganglion cells, i.e., individual neurons, not background populations projecting noise. But if that use case is not relevant for you, there is not need to add support for it.



/*BeginDocumentation
Name: inh_poisson_generator - provides Poisson spike trains at a piecewise
Copy link
Contributor

Choose a reason for hiding this comment

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

I am not happy with the name, inh is too little informative. I would suggest inhomogeneous_poisson_generator, even though it is rather long, as the very least inhom_... to avoid any confusion with "inhibitory". The class name should change accordingly.

Copy link
Author

Choose a reason for hiding this comment

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

done


Description:
The inhomogeneous Poisson generator provides Poisson spike trains at a
piecewise constant rate
Copy link
Contributor

Choose a reason for hiding this comment

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

Beautify line breaks.

Copy link
Author

Choose a reason for hiding this comment

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

done

CopyModel
before a generator node is created, the generator will send the same spike
train
to all 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.

Beautify line breaks.

Copy link
Author

Choose a reason for hiding this comment

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

done

Buffers_( const Buffers_&, inh_poisson_generator& );
UniversalDataLogger< inh_poisson_generator > logger_;
};
*/
Copy link
Contributor

Choose a reason for hiding this comment

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

Delete unused code.

Copy link
Author

Choose a reason for hiding this comment

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

done


// Keep the amplitude up-to-date at all times.
// We need to change the amplitude one step ahead of time, see comment
// on class SimulatingDevice.
Copy link
Contributor

Choose a reason for hiding this comment

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

Here is a typo, should read StimulatingDevice right?

Copy link
Author

Choose a reason for hiding this comment

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

I think so, changed.


// create spikes
if ( B_.rate_ > 0 && device_.is_active( Time::step( t0 + offs ) ) )
{
Copy link
Contributor

@jschuecker jschuecker Mar 15, 2017

Choose a reason for hiding this comment

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

instead of t0 + offs you could use curr_time from line 179.

Copy link
Author

Choose a reason for hiding this comment

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

done

Copy link
Contributor

@jschuecker jschuecker left a comment

Choose a reason for hiding this comment

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

@zbarni, nice work! I added comments in the code and also commented on the critical issue raised by @heplesser.

@@ -70,7 +70,7 @@ set( models_sources
music_message_in_proxy.h music_message_in_proxy.cpp
noise_generator.h noise_generator.cpp
parrot_neuron.h parrot_neuron.cpp
inh_poisson_generator.h inh_poisson_generator.cpp
inhomogeneous_poisson_generator.h inhomogeneous_poisson_generator.cpp
Copy link
Contributor

Choose a reason for hiding this comment

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

This looks like a spurious indent.

Copy link
Author

Choose a reason for hiding this comment

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

done

@@ -93,7 +93,7 @@ nest::inh_poisson_generator::Parameters_::set( const DictionaryDatum& d,
}
}

if ( ut && uv )
if ( times && rates )
Copy link
Contributor

Choose a reason for hiding this comment

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

&& -> and

Copy link
Author

Choose a reason for hiding this comment

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

done

@heplesser
Copy link
Contributor

@zbarni This concerning the times at which rates change has now also popped up in #740. We must move all times to the grid, using the same mechanism as in spike_generator. precise_times are irrelevant here, since this generator is strictly interval-based (poisson_generator_ps is for poisson processes with precise times, but that is another story). It makes sense to handle allow_offgrid_times as in spike_generator. It is important to check that the spike times after conversion to Time objects on the grid are strictly increasing.

@heplesser
Copy link
Contributor

@zbarni Please run extras/check_code_style.sh on your code and fix all formatting issues!

@heplesser
Copy link
Contributor

@zbarni Please address the problems in this PR by 12 November or we will close this PR due to inactivity.

@zbarni
Copy link
Author

zbarni commented Nov 8, 2017

@heplesser sorry for the extreme delays with this PR, but I have been working on my master's thesis in the last months and didn't really get a chance to look at the revisions. Would it be possible the keep the PR open for a few more weeks, maybe early December? I would implement the requested changes and some tests after my deadline, which is end of this month. It would be a pity not to have this PR go through as most of the work is already done.

@heplesser
Copy link
Contributor

@zbarni Thank you for the explanation. We will keep the PR open for now, but at most until the end of the year.

@zbarni
Copy link
Author

zbarni commented Dec 27, 2017

I'm not sure why the current commit(s) fails with cppcheck.. it complains about not finding all the include files. Can somebody take a quick look at the build? It might be related to #871.

@jougs
Copy link
Contributor

jougs commented Jan 12, 2018

The problem is not caused by cppcheck, but rather by an error in the new test test_inhomogeneous_poisson_generator.sli:

Dec 24 14:01:51 dup [Error]: StackUnderflow

Probably the error can even be reproduced locally using make installcheck. Could you please check? Thanks!

@zbarni
Copy link
Author

zbarni commented Jan 15, 2018

Thanks for taking a look! Yes, the error is indeed there, but I've no idea why it happens. It's a very simple dummy test case, and when I run the same code (create a generator, set the parameters, and simulate for 10 ms) in a regular python script it runs without errors.

@stinebuu
Copy link
Contributor

@zbarni I don't think the problem with the test is that there is an error per se, it is because with assert_or_die you have to have something to actually test, usually with the eq function, which checks for equality. A simple example can be found here, which asserts that when you apply round to 1.4 it equals 1. A bit more complex example is here in test_spike_detector.sli, which test that precision is what we expect it to be.

Your test runs, but then it comes to assert_or_die and gives the StackUnderflow because it is missing something to assert.

Hope this helps!

@zbarni
Copy link
Author

zbarni commented Feb 5, 2018

@stinebuu thanks for the hint, it works now!
@heplesser I might have overlooked something, but I guess all the changes that were requested are now implemented. Please let me if anything has to be corrected / modified. Thank you!

EDIT: maybe one question: what should be the expected behavior when trying to set a rate at time 0.?
Currently this silently ignored as the rate will be set/stored but there will be no spikes generated..

@heplesser
Copy link
Contributor

@zbarni The code looks almost fine now, except for the problem with setting a rate at 0. Since that rate is silently ignored, this will lead to errors: Users will expect the generator to run with a given rate from time 0, and will likely not notice that the generator isn't emitting any spikes. That is dangerous and we need to prevent it.

The simplest approach would be to allow only rate times that are in the future. We could also allow setting the rate "now"; to support this, maybe calibrate() can check what rate should apply in the beginning.

@zbarni
Copy link
Author

zbarni commented Feb 26, 2018

@heplesser I added a check to make sure that rate times are strictly in the future, including a test case. This avoids the problem at 0., and I think it's an acceptable limitation for the users. However, I can also implement support for 'now' if you think it's necessary.

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

@jschuecker Could you re-check this PR and approve it if all is in order? Then it could go into 2.16.

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.

None yet

8 participants