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

Stochastic processes: part 4 #590

Merged
merged 6 commits into from
Dec 17, 2014
Merged

Stochastic processes: part 4 #590

merged 6 commits into from
Dec 17, 2014

Conversation

tbekolay
Copy link
Member

Here's another go at implementing stochastic processes for ensemble noise, this time using the newly refactored synapse classes. The last commit, 138ee2a, implements that new stuff (the rest is from prior iterations).

In general this works, but the Brown and white noise implementations don't pass unit tests, mostly because I don't have my eng. ring and therefore don't know how to design filters.

Brown noise

Here's plot from the Brown noise test:

test_processes test_brown_noise

The bottom subplot shows the standard deviation (along axis=1, so the standard deviation across the 5000 dimensions of Brown noise for each timestep). As can be seen, this starts off okay, but then drops off over time. I think this is because I'm using LinearFilter(num=[dt, 1], den=[1, 1]) for filtering, which, if I'm understanding things correctly, is forward Euler integration, which does a bad job.

White noise

Generating white noise works as expected, but I'm not sure how to make a filter to implement the cutoff.

test_processes test_whitenoise_high_5

My understanding is that a lowpass filter with tau = 1./ (2 * np.pi * high) where high is the cutoff frequency should do it, but it does a pretty butts job IMO.


Anyhow, these do fulfil the promise of not needing to know the signal duration beforehand, but perhaps it's also useful to also be able to generate white noise with ideal filters?

This allows for injecting noise drawn from any Distribution.
Technically this also allows for injecting arbitrary current,
as I've done for testing purposes.
@jgosmann
Copy link
Collaborator

Not sure what is going on with the Brown noise, mostly because I do not understand LinearFilter. Did you check that the plot did match in earlier iterations?

Regarding the white noise: The filter is producing an exponential shape in the power spectrum, whereas my implementation keeps the frequencies < high flat and does a hard cut-off for frequencies > high. I wouldn't be surprised if this cannot be reproduced by how we implement filters. Could have some special processes (like white noise) that do not make use of filtering?

@tbekolay
Copy link
Member Author

Brown noise looks good now, and I just pushed a fix for white noise which essentially does what utils.functions.whitenoise did, so I was able to remove it. See e723da5 for more details.

Unfortunately something's still not quite right with setting the RMS, but I'm not sure what.

It's fine for large values of high (ignore the weird font; high=500, rms=100.0 in this plot)

test_processes test_whitenoise_rms_100 0

But for small values it's messed up (high=5, rms=0.5 here):

test_processes test_whitenoise_high_5

Medium values also messed up (high=50, rms=0.5):

test_processes test_whitenoise_high_50

I'm guessing this is because I'm generating too many sine waves for low frequencies; I always generate 50 in these plots. Probably should be lower for lower cutoff frequencies? I'm also not using dt anywhere, but perhaps am missing somewhere that I should be using it.

@jgosmann
Copy link
Collaborator

Is there a reason for using the approach with the basis functions instead of the implementation I proposed in PR #523? I listed reasons why I think my implementation is better in a comment. (better performance, explicit period length of generated signal, uses as many basis function as possible given dt and duration which results in "better" noise).

Regarding the RMS: It took me quite some time to get all the normalization constants right (because they depend on the FFT normalization which can differ). By filtering the signal the relationship that the RMS in the time domain and frequency domain are the same does not hold anymore. Thus, first one has to decide which RMS one wants to measure (time domain, frequency domain over the whole domain, or frequency only over the non-zero frequencies). Then the normalization has to be figured out.

@drasmuss
Copy link
Member

I think the reason for the change was explained in the commit message (
559bcaa),
also there's some comments from @tbekolay in that comment thread you linked
to which also explain the logic for the change. To paraphrase, I think the
main motivation is to make things more simple and consistent from an API
perspective (and I agree with @tbekolay that the refactored version seems
cleaner in that regard). Is your suggestion that we could get the same API
benefits with a different implementation (also no eng ring here)?

On 16 December 2014 at 22:05, Jan Gosmann notifications@github.com wrote:

Is there a reason for using the approach with the basis functions instead
of the implementation I proposed in PR #523
#523? I listed reasons why I think
my implementation is better in a comment.
#523 (comment) (better
performance, explicit period length of generated signal, uses as many basis
function as possible given dt and duration which results in "better"
noise).

Regarding the RMS: It took me quite some time to get all the normalization
constants right (because they depend on the FFT normalization which can
differ). By filtering the signal the relationship that the RMS in the time
domain and frequency domain are the same does not hold anymore. Thus, first
one has to decide which RMS one wants to measure (time domain, frequency
domain over the whole domain, or frequency only over the non-zero
frequencies). Then the normalization has to be figured out.


Reply to this email directly or view it on GitHub
#590 (comment).

@tbekolay
Copy link
Member Author

@jgosmann honestly I don't care either way; after having worked with both I prefer your approach, so I'm adapting it to this API now (I'll be squashing my commits so I'll make sure you get the git blame). When I talked to @hunse today he suggested using the approach in utils.functions.whitenoise so I did it, but yeah, it is definitely slower so I'm in for either way. Mostly just want this done!

@jgosmann
Copy link
Collaborator

The only thing in consistency which differs in my approach is that you have to decide for a dt ahead of time. (We could interpolate to circumvent this. I was too lazy to implement that because np.interp only supports one dimension.)
Also, it might appear that utils.functions.whitenoise does not require you to know the simulation duration in advance. But that is not true as in both approaches the signal will repeat (in the utils.functions.whitenoise approach it depends on high and n_bases).

So mainly it is an implementation detail and a decision what parametrization is more natural (dt and duration vs. n_bases and a dependence on high).

As everyone is making that disclaimer: ✋ no eng ring here (we don't get in Germany and I'm not sure whether computer science would qualify for it anyways).

@tbekolay
Copy link
Member Author

OK, pushed 7ba24e8 which reverts back to @jgosmann's original implementation. Again, don't really care which one we go with, but we now have both available so whichever the consensus likes best I'll clean up the history such that it seems like we always only wanted that one ;) But let's make it a quick consensus as I'm out of here Wednesday afternoon!

If an opinion from me is needed I'd go with @jgosmann's as it's faster and the repetition doesn't matter in most cases. When it does matter, could use StochasticProcess directly.

Oh, and it could be possible to include both, but I think that's more confusing than helpful.

@jgosmann
Copy link
Collaborator

👎 for including both

@hunse
Copy link
Collaborator

hunse commented Dec 17, 2014

I don't really care either way. I don't like the fact that you have to include simulation time and dt with the pregeneration method, since it's not like any other process. What will happen when it reaches the end of it's allotted time? Just remain constant? Jump back to the beginning (with a sudden jump to the initial value)? Play the time series backwards and then forwards again (at least that way it's continuous)? I worry that users might be surprised by what happens if they run the simulation longer than they specify for the process.

But there are cons to the sine wave method, too. As we've said, it is slower, and I guess it will repeat eventually (though it seems to me this would take a long time, since all the waves would have to get in phase again). It's also impossible to include all frequencies when the simulation can be arbitrarily long (though this is not really a practical concern IMO). I think the biggest concern is that to get a descent signal out, one should include a number of components (probably at least one each Hz in the pass-band), and for higher frequency cutoffs this could be quite a few, making it slow (though maybe it should be a set number of components no matter the cutoff).

I'm fine with either one.

@tbekolay
Copy link
Member Author

Apathy: the best consensus! I agree that the discontinuity when the signal repeats is bad, so I pushed a commit that plays it forwards then backwards again, as you suggest. This might be surprising in situations where the statistics of the noise are very important, but in that case it's hopefully not too onerous to implement a new process. But yeah, I would be surprised if there was jump in the whitenoise signal, even though reversing it isn't a perfect solution either.

Anyhow, I'll do a bit of history cleanup and make sure this is all in order this morning so we can get it merged in and move on with RC1.

jgosmann and others added 2 commits December 17, 2014 10:07
Includes implementations of a Weiner process, band-limited white noise,
and a generic SampledProcess for sampling from distributions directly.
This is necessary for implementing Brownian noise with distributions
and synapses. The matrix exponential in 'zoh' is broken, so the
'euler' method will be necessary in that case.
@tbekolay tbekolay force-pushed the ensemble_noise4 branch 2 times, most recently from 5570037 to 4eef64a Compare December 17, 2014 15:33
The overall change here is to make processes more like synapses.
That is, they don't need to know their dimensionality ahead of time;
this is instead done when sampling. Similarly, we use partials to
bind state to the function, rather than having it be part of the
object, or be explicitly passed in on each sample.

In addition to an API more consistent with synapses, the big win with
this change is that it is no longer necessary to `copy.copy` processes
when using them in multiple ensembles; `SimNoise` takes care of this
instead (or this could be relegated to the builder in other
`StochasticProcess` use cases). Using `copy` (or, even worse,
`deepcopy`) was a huge problem in earlier iterations of Nengo, so I'm
keen to never use it again. This also makes it possible to use
processes as a default with the config system.

Another behavioral change is that `WhiteNoise` now reflects instead
of loops; that is, once you reach the end of the pre-generated signal
it will replay it in reverse, until it reaches the start, then it will
start playing forward again. This is to limit the discontinuity at
the end of the generated signal; with looping, it will jump from
one value to the next. With reflection, the derivative will flip but
the signal itself will still be sort of continuous. If the derivative
happens to be small at the end of the signal, then it'll be like it's
a continuous signal.

A few more detailed changes:

- Added some docs to `Distribution`
- Naming changes: `WienerProcess` -> `BrownNoise`,
  `LimitedGuassianWhiteNoise` -> `WhiteNoise`
- Removed `GaussianWhiteNoise` to avoid confusion with `WhiteNoise`.
  If someone wants non-limited white noise, they can either pass
  `high=None` to `WhiteNoise`, or make it using `StochasticProcess`
  (as is done in `test_processes.py`)
- Added plotting to `test_processes.py` and `test_ensemble.py`
And used the `processes.WhiteNoise` instead. To make this relatively
easy from the end-user perspective, I added an `f` method to
`StochasticProcess` that returns a function that can be used as a
Node output. This basically just wraps the sample function and ignores
the incoming `t`.

Note that when used as a node function, processes maintain state even
if the simulation is reset. This is skirted around in `test_reset`
by using a white noise signal with half the duration of the simulation
so that it repeats. In general, however, we can't expect that node
functions will reset when the simulation resets.
@tbekolay
Copy link
Member Author

OK, I've cleaned up the history so I think this is ready to merge. It may not be perfect, but I think it's pretty good and definitely sufficient for 2.0!

@tbekolay tbekolay added this to the 2.0.0 release milestone Dec 17, 2014
"\n",
"model = nengo.Network(label=\"Delayed connection\")\n",
"with model:\n",
" # We'll use white noise as input\n",
" inp = nengo.Node(whitenoise(1, 5, seed=60))\n",
" inp = nengo.Node(WhiteNoise(2, high=5).f())\n",
Copy link
Collaborator

Choose a reason for hiding this comment

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

There is no WhiteNoise.f anymore, is there?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Oh my bad. It's StochasticProcess.f.

@jgosmann jgosmann merged commit 301d45c into master Dec 17, 2014
@tbekolay tbekolay deleted the ensemble_noise4 branch December 17, 2014 16:38
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants