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

Fix for multimeter origin bug. #669

Merged
merged 8 commits into from May 30, 2017
Merged

Conversation

@zbarni
Copy link

@zbarni zbarni commented Mar 3, 2017

Modified multimeter to actually consider the origin variable when sampling. Modification of the DataLoggingRequest class was also required to keep track of the origin. Main issue was that the initial sampling timepoint (relative to which the next intervals are calculated) was computed without regard to the origin.

…pling. Modification of the DataLoggingRequest class was also required to keep track of the origin.
@heplesser
Copy link
Contributor

@heplesser heplesser commented Mar 6, 2017

@zbarni As I explained in #670, multimeter works as intended, although some documentation is misleading. The main problem with your proposal here is that the purpose of origin is to allow the user to move the start-stop window along the time axis by changing origin during stimulation breaks. Your intended use of origin would conflict with this. Also, a logger will only take into account the value of origin at the first Simulate call after the multimeter is first connected. Any later changes would be difficult, because it would require resizing of logger buffers and thus risk losing or invalidating data that has been recorded already in the logger inside the neuron recorded from, but that has not yet been fetched by the multimeter.

I therefore think that you approach is not feasible.

@zbarni
Copy link
Author

@zbarni zbarni commented Mar 6, 2017

I think the discussion has moved to #670 .

@heplesser
Copy link
Contributor

@heplesser heplesser commented Mar 7, 2017

I think we should aim for a more general solution along the lines of what I sketched in #505. Therefore, I will close this PR for now.

@heplesser heplesser closed this Mar 7, 2017
@heplesser
Copy link
Contributor

@heplesser heplesser commented Mar 7, 2017

@zbarni After thinking it through a bit further, here a suggestion for how to handle subsampling with arbitrary offset:

  1. Sampling is parameterized in terms of offset and interval such that data is recorded at all times t == n h == offset + k interval, where h is the simulation resolution and n and k are non-negative integers (n>=1, k>=0).
  2. Just as interval, offset can only be set before the multimeter is connected to any neurons; we cannot change it later because that would require modifying the recording buffers, risking loss of data. Since interval already has this restriction, this should not be problematic.
  3. Windowing (anything that has to do with start, stop, origin) is orthogonal to this and handled entirely by base classes.

Could you try to implement this? It would also be excellent if you could add a test for the new offset feature, including tests that modifying offset after connecting is prohibited.

Sorry for the back-and-forth!

@heplesser heplesser reopened this Mar 7, 2017
@zbarni
Copy link
Author

@zbarni zbarni commented Mar 12, 2017

@heplesser no problem, probably the role of offset is clearer this way. Finished the changes you suggested (offset as a new variable instead of overlapping with origin), and added a simple unit test for correctness and changing the offset after connection. Let me know if anything else should be changed or extended.

Copy link
Contributor

@heplesser heplesser left a comment

@zbarni Thanks for the changes. It looks pretty good now, see the comments in the code for some detailed suggestions.

if ( updateValue< double >( d, names::offset, v ) )
{
// only throw error if given offset is > 0., otherwise it's not an error
if ( Time( Time::ms( v ) ) < Time::get_resolution() && v > 0. )

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

I think this code would be easier to read with the tests the other way around; also, we now prefer boolean operators written out, i.e.,

if ( v > 0 and Time( ... ) < ... )

It is also not necessary to write 0. here, 0 will do.

V_.new_request_ =
B_.has_targets_ && !P_.record_from_.empty(); // no targets, no request
V_.new_request_ = B_.has_targets_
&& P_.record_from_.empty() == false; // no targets, no request

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

One should not compare with true and false. The following would be better (comment moved for better line structure):

// no targets or no recordables means no request
V_.new_request_ = B_.has_targets_ and not P_.record_from_.empty();

if ( !is_active( info[ j ].timestamp ) )
if ( is_active( info[ j ].timestamp ) == false )

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

if ( not is_active( ... ) )
@@ -213,22 +256,28 @@ Multimeter::handle( DataLoggingReply& reply )

// record sender and time information; in accumulator mode only for first
// Reply in slice
if ( !device_.to_accumulator() || V_.new_request_ )
if ( device_.to_accumulator() == false || V_.new_request_ )

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

if ( not device_.to_accumulator() or V_.new_request_ )

if ( !device_.to_accumulator() )
if ( device_.to_accumulator() == false )

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

if ( not device_.to_accumulator() )
* for which (T-(origin+start)) mod interval is zero.
* data is with time stamp offset+start+1, the last recorded one
* that with time stamp offset+stop. Only such times are recorded
* for which (T-(offset+start)) mod interval is zero.

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

start, stop and origin are independent of offset and interval. Isn't the correct description now "Data is recorded at time steps T for which start < T <= stop and ( T - offset ) mod interval == 0."?

@@ -237,7 +238,8 @@ class Multimeter : public Node

struct Parameters_
{
Time interval_; //!< recording interval, in ms
Time interval_; //!< recording interval, in ms
Time offset_; //!< offset relative to which interval is calculated

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

add "in ms"

@@ -582,18 +587,22 @@ class DataLoggingRequest : public Event
inline DataLoggingRequest::DataLoggingRequest()
: Event()
, recording_interval_( Time::neg_inf() )
, recording_offset_( Time::ms( 0. ) )

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

0 would suffice

This comment has been minimized.

@zbarni

zbarni Mar 13, 2017
Author

This doesn't work.. initializing recording_offset_ with ( 0 ) leads to a compile error.

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

@zbarni But Time::ms(0) rather than (0.) should be enough.

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

Good to know -- I should open an issue on ambiguities in time.

This comment has been minimized.

@apeyser

apeyser Apr 4, 2017
Contributor

There's an explicit ms(long t). I wonder why there's an ambiguity on the constructor. The long/double constructors should probably anyway be reduced to just the double constructor.

{ ; true }
ifelse
}
Map

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

It isn't really necessary to run this test for all neuron models and all recordables for each model. Running it for a single model with two recordables would suffice and make the test run faster.

It would be useful to add another test that uses start/stop and origin (changing origin during a simulation break) to see that offset still works correctly then.

@@ -555,7 +555,7 @@ class DataLoggingRequest : public Event
DataLoggingRequest();

/** Create event for given time stamp and vector of recordables. */

This comment has been minimized.

@heplesser

heplesser Mar 13, 2017
Contributor

Could you extend this comment to explain the different parameters? That becomes more important since we have two Time arguments now.

@heplesser heplesser requested a review from apeyser Mar 13, 2017
@heplesser
Copy link
Contributor

@heplesser heplesser commented Mar 13, 2017

@zbarni Thanks for the fixes, nice work!

@zbarni
Copy link
Author

@zbarni zbarni commented Mar 13, 2017

Glad I could contribute!

"as the simulation resolution." );
"The sampling interval must be at least as long as the simulation "
"resolution." );
}

// see if we can represent interval as multiple of step

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

This could use some commenting, especially since the it's reused at line 142.
So I believe this is saying

floor(v * steps/msec) * msec/step / v = 1 +/- 10*\epsilon

which is just another way of saying that floor(v * steps/msec) - v*steps/msec = +- v * 10\epsilon
and that is just as true for 142, offset_

This comment has been minimized.

@zbarni

zbarni Apr 3, 2017
Author

This is part of the original code, I just copied the whole if block.

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

Understood -- but if possible, add some clarification to the block you copied, because otherwise the lack of commentary will just propagate through the code, leading to continuing confusion (for me).

This comment has been minimized.

@zbarni

zbarni Apr 3, 2017
Author

Can you point out which lines are you referring to exactly? 121-129? It seems to me that it's just one way to get around the numerical issues, but the main idea is present in the comment

// see if we can represent interval as multiple of step

Even if the interval is a multiple of step, std::abs( 1 - interval_.get_ms() / v ) might be very close but still != 0, hence the comparison with epsilon. Should I expand on this? not sure the details are necessary, but as you say ;)

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

The Time class provides the method is_multiple_of(), so

if ( not interval_.is_multiple_of( Time::get_resolution ) )
{
  throw ...;
}

will perform the test with all rounding logic handled in the Time class.

if ( updateValue< double >( d, names::offset, v ) )
{
// only throw error if given offset is > 0., otherwise it's not an error
if ( v > 0 && Time( Time::ms( v ) ) < Time::get_resolution() )

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

0 < v < t_resolution is an error? But not v > t_resolution?
Please comment.

Also, v should just be converted to time once (say v_time), rather than assuming that the compiler can optimize it out -- it's also easier to read/

This comment has been minimized.

@zbarni

zbarni Apr 3, 2017
Author

0 < v < t_resolution is an error, because if v is other than the default value (0), it has to be multiple of the resolution (this means positive and at least as large as the resolution itself.)

v > 0 was wrong, should have been v != 0. Not sure about introducing another variable v_time , or creating a new object by converting to Time(), as at that moment v is still double, and we're not comparing it against another object, just 0.

@@ -582,18 +587,22 @@ class DataLoggingRequest : public Event
inline DataLoggingRequest::DataLoggingRequest()
: Event()
, recording_interval_( Time::neg_inf() )
, recording_offset_( Time::ms( 0. ) )

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

@zbarni But Time::ms(0) rather than (0.) should be enough.

Copy link
Contributor

@apeyser apeyser left a comment

I'm sure it exists in the referenced issues, but I didn't manage to pull it out. Why is an offset for recording being added to multimeter?

@zbarni
Copy link
Author

@zbarni zbarni commented Apr 3, 2017

@apeyser as @rcfduarte put it in #670:

...we're running very long simulations where we need to sample the values of population state variables at regular intervals. For example, we would like to sample the population activity at 100.1, 200.1, 300.1, etc (where the extra 0.1 ms is added to account for a delay of 0.1 ms from an input generator to the network). As @zbarni mentioned, as it is currently implemented, the multimeter doesn't allow this (unless we set interval=0.1, which is not a very efficient solution for long simulations)...

Let me know if this answers your question.

@zbarni
Copy link
Author

@zbarni zbarni commented Apr 3, 2017

@apeyser regarding

But Time::ms(0) rather than (0.) should be enough.

Doesn't compile, yielding the error

.....event.h:591:36: error: call of overloaded ‘ms(int)’ is ambiguous
   , recording_offset_( Time::ms( 0 ) )

Same case for recording_offset_( 0 )

Copy link
Contributor

@apeyser apeyser left a comment

Thanks for the information -- if you can just add some information on the copied code, I'll accept.

"as the simulation resolution." );
"The sampling interval must be at least as long as the simulation "
"resolution." );
}

// see if we can represent interval as multiple of step

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

Understood -- but if possible, add some clarification to the block you copied, because otherwise the lack of commentary will just propagate through the code, leading to continuing confusion (for me).

@@ -582,18 +587,22 @@ class DataLoggingRequest : public Event
inline DataLoggingRequest::DataLoggingRequest()
: Event()
, recording_interval_( Time::neg_inf() )
, recording_offset_( Time::ms( 0. ) )

This comment has been minimized.

@apeyser

apeyser Apr 3, 2017
Contributor

Good to know -- I should open an issue on ambiguities in time.

@mschmidt87
Copy link

@mschmidt87 mschmidt87 commented May 29, 2017

Hi, I just stumbled upon the same issue in my simulations. I have read this PR and the related issue #670 . If I understood it correctly, this PR now aims to implement the suggestion by @heplesser , namely that the recording device should sample at times

t == n h == offset + k interval, where h is the simulation resolution and n and k are non-negative integers (n>=1, k>=0).

However, it seems to me that with the code implements k>0. Executing the following code:

import nest
nest.ResetKernel()
n = nest.Create('iaf_psc_exp', 1)
vm = nest.Create('voltmeter', 1)
nest.SetStatus(vm, {'interval': 5.,
                    'offset': 30.1})
nest.Connect(vm, n)
nest.Simulate(100.)
data = nest.GetStatus(vm, 'events')[0]['times']

gives recorded times [ 35.1 40.1 45.1], whereas I would have expected [30.1 35.1 40.1 45.1], which would be equivalent to setting k=[0,1,2,3] in the equation suggested by @heplesser (with offset=30.1, interval=5.).

Did I misunderstand something here?

"nodes." );
"The recording interval, the interval offset and the list of properties "
"to record "
"cannot be changed after the multimeter has been connected to nodes." );

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

Could you fix the line breaks to that the last line is shortest?

This comment has been minimized.

Copy link
Contributor

@heplesser heplesser left a comment

@zbarni Sorry for the long delay. I added a suggestion of how to address the issues raised by @apeyser. Could you also check the issue raise by @mschmidt87? His expectation is correct, recording in his case should start with 30.1, not 35.1. Could you add his case to the unittest?

"as the simulation resolution." );
"The sampling interval must be at least as long as the simulation "
"resolution." );
}

// see if we can represent interval as multiple of step

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

The Time class provides the method is_multiple_of(), so

if ( not interval_.is_multiple_of( Time::get_resolution ) )
{
  throw ...;
}

will perform the test with all rounding logic handled in the Time class.

offset_ = Time::step( Time( Time::ms( v ) ).get_steps() );
if ( std::abs( 1 - offset_.get_ms() / v ) > 10
* std::numeric_limits< double >::epsilon() )
{

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

Same as above here with not offset_.is_multiple_of( Time::get_resolution() ).

This comment has been minimized.

@zbarni

zbarni May 29, 2017
Author

applied both changes along with the fix for @mschmidt87 's example and edited the unit test to account for this. Let me know if anything else needs to be done.

@heplesser
Copy link
Contributor

@heplesser heplesser commented May 29, 2017

@zbarni The conflict reported by Github seem to be just line break/empty line/whitespace changes and should be easy to fix. The easiest is probably to pull the current master into your branch.

Copy link
Contributor

@heplesser heplesser left a comment

@zbarni Looks good, just some tiny details left (see individual comments). And could you add one more test case, just to check that DataLogger initiatialization works correctly: Create network with multimeter so that if records between 20 and 30 ms, e.g.. Then simulate 10 ms, and add another multimeter that also should record between 20 and 30 ms. Both should return the same data.

@@ -103,8 +103,8 @@ nest::Multimeter::Parameters_::set( const DictionaryDatum& d,
{
throw BadProperty(
"The recording interval, the interval offset and the list of properties "
"to record "
"cannot be changed after the multimeter has been connected to nodes." );
"to record cannot be changed after the multimeter has been connected "

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

This needs a space at the end of the string, since the message output mechanism may change line breaks when the message is printed.

( kernel().simulation_manager.get_time().get_steps() / rec_int_steps_
+ 1 ) * rec_int_steps_
- 1;
}

This comment has been minimized.

@heplesser

heplesser May 29, 2017
Contributor

Could you add a comment on why you need to add 1 and move the -1 up one line?

This comment has been minimized.

@zbarni

zbarni May 30, 2017
Author

clang requires the -1 to be on a separate line...


MSGBLD0175: [DIFF] <     ( kernel().simulation_manager.get_time().get_steps() / rec_int_steps_
MSGBLD0175: [DIFF] <       + 1 ) * rec_int_steps_ - 1;
MSGBLD0175: [DIFF] ---
MSGBLD0175: [DIFF] >     ( kernel().simulation_manager.get_time().get_steps() / rec_int_steps_ + 1 )
MSGBLD0175: [DIFF] >       * rec_int_steps_
MSGBLD0175: [DIFF] >     - 1;
@zbarni
Copy link
Author

@zbarni zbarni commented May 29, 2017

@heplesser I'm actually unsure whether the offset is supposed to be relative to 0. or relative to the elapsed simulation time since 0.? For instance,

nest.Simulate(10.)
nest.Create("multimeter", 1, {"interval": 3., "offset": 5.})
(...)
nest.Simulate(10.)

should give [15. 18.] (relative to elapsed time, 10. ms) or [11. 14. 17. 20.] (relative to 0.)? Currently, the intervals are calculated relative to 0. when a new multimeter is created and not the elapsed time..

@heplesser
Copy link
Contributor

@heplesser heplesser commented May 30, 2017

@zbarni We should stick to the current approach (relative to 0.), so that the recording times are given by t_k = offset + k * interval. This makes offset/interval semantics consistent with the origin/start/stop semantics which also use 0 as a reference point. Furthermore, it ensures that recording times remain fixed even if you change simulation intervals or move code lines around.

I just suggested to add the test to ensure that recording times are correctly based on 0 even if a multimeter is created after the simulation has run for some time.

@heplesser
Copy link
Contributor

@heplesser heplesser commented May 30, 2017

Add possiblity to set offset for multimeter recording. Thanks, @zbarni!

@heplesser heplesser merged commit 08e830b into nest:master May 30, 2017
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.