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

[MRG] DOC: use morlet instead of multitaper #226

Merged
merged 5 commits into from Jan 4, 2021

Conversation

jasmainak
Copy link
Collaborator

@jasmainak jasmainak commented Dec 20, 2020

closes #223

I was messing around with the gamma examples and realized that HNN-GUI uses a morlet wavelet for TFR. This is an attempt to make them more similar and add a colorbar.

While I was at it, I also updated the tonic input example and made a couple of fixes there.


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
dpls[0].plot(ax=axes[0], layer='agg', show=False)
axes[0].plot(times, data)
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I was in two minds whether to add tstart and tstop to the plotting routine or just directly use matplotlib to plot. I settled with the later ... the utility of a plotting routine for Dipole object seems a bit limited to me unless we come up with more use cases

Copy link
Contributor

Choose a reason for hiding this comment

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

Are you thinking of just removing the plot() methods from Dipole? I agree it's fairly unnecessary considering most people are capable of making a line plot.

Now there's certainly a case to be made for more visualization options (multi-trial with mean dipole overlaid for example), but that would be better suited for viz.py.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

right, so perhaps we need a plot_dipole function? Still it's rather straightforward to do with 2 lines of Python code, but maybe we can add functionality later.

@jasmainak jasmainak changed the title DOC: use morlet instead of multitaper [MRG] DOC: use morlet instead of multitaper Dec 20, 2020
n_cycles=n_cycles, output='power')

im = axes[1].pcolormesh(times, freqs, power[0, 0, ...], cmap='RdBu_r',
shading='auto')
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

shading='auto' was needed to prevent a deprecation warning in matplotlib.

@codecov-io
Copy link

codecov-io commented Dec 20, 2020

Codecov Report

Merging #226 (e49b95b) into master (31b2a93) will increase coverage by 0.04%.
The diff coverage is 87.50%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #226      +/-   ##
==========================================
+ Coverage   66.84%   66.88%   +0.04%     
==========================================
  Files          21       21              
  Lines        2277     2283       +6     
==========================================
+ Hits         1522     1527       +5     
- Misses        755      756       +1     
Impacted Files Coverage Δ
hnn_core/network.py 50.36% <83.33%> (+0.36%) ⬆️
hnn_core/tests/test_feed.py 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 31b2a93...602a157. Read the comment docs.

feed_times = self.feed_times
else:
feed_times = dict()
feed_times.update({src_type: list() for src_type in src_types})
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

this is a bit of a spaghetti anti-pattern here. We should try to improve this. I think _instantiate_feeds should be called just once with n_trials rather than being overwritten by another call. It's surprising.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@cjayb hope this does not cause too much of a rebase issue for you?

Copy link
Collaborator

Choose a reason for hiding this comment

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

@cjayb hope this does not cause too much of a rebase issue for you?

don't worry about it, I'm still knee-deep in figuring out parconnect. I'll rebase when the time comes, still working on and need to simplify #221 (once I get it working, that is)

@jasmainak jasmainak force-pushed the gamma_tfr branch 2 times, most recently from 2500a63 to e34d37b Compare December 20, 2020 16:07
freqs = np.arange(20., 100., 1.)
n_cycles = freqs / 4.
n_cycles = freqs / 8.
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this to keep a consistent time resolution between morlet and multitaper? I'm struggling to remember exact impact for the two methods.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

no, not really. MNE complains about wavelet cycles being too long if n_cycles = freqs / 4. Not sure why it doesn't do the same with multitaper. Anyhow I think HNN-GUI uses n_cycles = freqs / 7 but MNE doesn't allow that, so I picked the closest possible :)

As for the difference between the two, this page gives a nice summary. Essentially, morlet is like doing multitaper with a single Gaussian taper

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Not sure why it doesn't do the same with multitaper.

maybe an opportunity to make a small PR to MNE if you are willing to dig into it? :)

@jasmainak
Copy link
Collaborator Author

@ntolley incorporated the raster plot in this example now. Take a look at the example in the CircleCI artifact tab once it is built.

It's trying to reproduce this: https://jonescompneurolab.github.io/hnn-tutorials/gamma/gamma.

The merge button is yours when you're happy with the changes/plots.

@jasmainak
Copy link
Collaborator Author

ping ping, any blockers with the reviews here?

Copy link
Collaborator

@cjayb cjayb left a comment

Choose a reason for hiding this comment

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

The new plots are beautiful! Few minor comments, otherwise do merge asap. I'll have to rebase #221 anyway, let's get this in there before.

###############################################################################
# The network requires some time to reach steady state. Hence, we omit the
# first 50 ms in our time-frequency analysis.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
trial_index = 0 # pick first trial

Comment on lines 41 to 43
mask = dpls[0].times > tstart
times = dpls[0].times[mask]
data = dpls[0].data['agg'][mask]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
mask = dpls[0].times > tstart
times = dpls[0].times[mask]
data = dpls[0].data['agg'][mask]
include_samples = dpls[trial_index].times > tstart
...

Copy link
Collaborator

Choose a reason for hiding this comment

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

Just to be explicit about what you're doing. Users may not feel comfortable with 'magic' constants like 0 before they get the hang of the objects.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

okay with adding trial_index. I actually prefer mask as opposed to include_samples though. It's easier to read the code because I know that it's a masked array as opposed to being indices or slices.

Comment on lines 86 to 91
# Notice that the Layer 5 pyramidal neurons are now firing nearly
# synchronously. They in turn synchronously activate the inhibitory basket
# neurons, which then inhibit the pyramidal neurons for ~20 ms, when the
# tonic drive outweighs the inhibition and the pyramidal neurons firing again
# creating a ~50 Hz PING rhythm. This type of synchronous rhythm is sometimes
# referred to as “strong” PING.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
# Notice that the Layer 5 pyramidal neurons are now firing nearly
# synchronously. They in turn synchronously activate the inhibitory basket
# neurons, which then inhibit the pyramidal neurons for ~20 ms, when the
# tonic drive outweighs the inhibition and the pyramidal neurons firing again
# creating a ~50 Hz PING rhythm. This type of synchronous rhythm is sometimes
# referred to as “strong” PING.
...
# neurons, which then inhibit the pyramidal neurons for ~20 ms. Once the tonic
# drive again outweighs the inhibition, the pyramidal neurons fire again,
# creating a ~50 Hz PING rhythm. This type of synchronous rhythm is sometimes
# referred to as “strong” PING.

Nitpick

@jasmainak jasmainak merged commit 93c79d9 into jonescompneurolab:master Jan 4, 2021
@jasmainak jasmainak deleted the gamma_tfr branch January 4, 2021 06:35
@jasmainak
Copy link
Collaborator Author

Gonna commit a small crime and self-merge to keep things moving with #221

Thanks for the reviews @ntolley and @cjayb

@jasmainak
Copy link
Collaborator Author

if you disagree with something, feel free to follow up with PRs :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

add example for tonic inputs
4 participants