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

Update to deal with issues for convolution and drift (ready for review) #309

Merged
merged 19 commits into from Jan 13, 2018

Conversation

Projects
None yet
3 participants
@CameronTEllis
Copy link
Contributor

CameronTEllis commented Dec 11, 2017

Two open issues reflect problems in fmrisim. One problem was with the convolution incorrectly downsampling — it now takes the mean and also increases the speed of analysis. Another problem was with the drift: it was a sine wave and that didn't make much sense. It now defaults to a combination of cosines.

CameronTEllis added some commits Dec 10, 2017

Changed the order of downsampling in convolve_hrf to accelerate compu…
…tation. Changed drift to default to cosine functions
@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Dec 11, 2017

@CameronTEllis, are you in touch with anyone for a review?

@CameronTEllis

This comment has been minimized.

Copy link
Contributor

CameronTEllis commented Dec 11, 2017

Given that they are small changes, are you able to review them? Or should I ask someone more familiar with the code?

@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Dec 11, 2017

Please ask someone more familiar.

@CameronTEllis

This comment has been minimized.

Copy link
Contributor

CameronTEllis commented Dec 11, 2017

@lcnature Would it be possible for you to find time to review this please?

@CameronTEllis CameronTEllis changed the title Update to deal with issues for convolution and drift Update to deal with issues for convolution and drift (ready for review) Dec 12, 2017

@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Dec 18, 2017

@lcnature, we would really appreciate your help reviewing this PR.

@lcnature
Copy link
Contributor

lcnature left a comment

Hi @CameronTEllis @mihaic
Sorry for my delayed review! Was in conference and travelling. Hope this review is helpful.

@@ -469,6 +469,16 @@ def generate_stimfunction(onsets,
# Pull out the onsets, weights and durations, set as a float
for line in text:
onset, duration, weight = line.strip().split()

# Check if the onset is more precise than the temporal resolution
dp = len(onset[onset.find('.') + 1:])

This comment has been minimized.

@lcnature

lcnature Dec 19, 2017

Contributor

There may be issue if some digits are written in scientific format, e.g., 1.0234e+03. We may want to comment in the docstring that such formatting can cause problem. Or, maybe we can solve it by converting the string to float and find the resolution in floating number.

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

Wow astute observation! This was a great comment and I have hopefully fixed it now

1000). If temporal_resolution is 1 then the output will be the same
length as stimfunction.
Be aware that if scaling is on and events are very short
(> temporal_resolution * tr_duration) then the hrf may or may not come

This comment has been minimized.

@lcnature

lcnature Dec 19, 2017

Contributor

What does it mean to be (> temporal_resolution * tr_duration) and "very short"? It seems that for tr_duration=2 and temporal_resolution=1000, this means something is larger than 2000. I suppose this is not the event duration, is it?

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

You are right, this wasn't clear. Hopefully it is more obvious now. The point is that if the events last less than a TR you may get weird results

absolute response after convolution but if there are only short events
and you scale then this will look identical to a convolution with longer
events. In general scaling is useful, which is why it is the default,
but be aware of this edge case

This comment has been minimized.

@lcnature

lcnature Dec 19, 2017

Contributor

Sorry I forget about this. Can user avoid this issue by setting scaling to false?

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

Yes they can, I added a comment to say as much


# Generate the hrf to use in the convolution
if hrf_type == 'double_gamma':
hrf = _double_gamma_hrf(temporal_resolution=temporal_resolution)
hrf = _double_gamma_hrf(temporal_resolution=1 / tr_duration)

This comment has been minimized.

@lcnature

lcnature Dec 19, 2017

Contributor

Although down-sampling stim_function helps speed up the convolution, I worry that this defeats the original purpose of building a smooth double_gamma_hrf. Consider a case of TR=2s, a stimulus of condition A starts at 1.5s and lasts for 0.5s, and a stimulus of condition B starts at 2s and lasts for 0.5s. I think the original implementation would make the design matrix of the two conditions appear as down-sampled version of two HRF with onsets 0.5 sec apart (they would have slightly different shapes, but very close to each other). But the new implementation might make them exactly the same shape, but 2 sec apart. Is this true? If so, maybe the old implementation is more desirable. You may set the default temporal resolution to some lower number such as 100 or 50 to achieve some speed-up with the old implementation. (most often the displays in scanners have refresh rate of just 60 Hz)

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

Fantastic point. I had been testing this in cases where the event duration was greater than the TR (and so this isn't a problem). But as you point out when the events are shorter than a TR then you get this exact problem. I have shifted the downsampling step back to after convolution. I also changed the default temporal_resolution: it was 1000 so that it was measuring milliseconds, but you are right >100Hz is excessive

@@ -1420,6 +1448,11 @@ def _generate_noise_temporal_drift(trs,
tr_duration : float
How long in seconds is each volume acqusition
basis : str
What is the basis function for the drift. Could be made of discrete
cosines (number of bases scale with duration) or a sine wave with

This comment has been minimized.

@lcnature

lcnature Dec 19, 2017

Contributor

You probably meant tr_duration by duration? And I did not quite understand with the temporal order

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

Rewrote this to clarify. The point is that the drift gets more basis functions (more kinks) when the run is longer

This comment has been minimized.

@CameronTEllis

CameronTEllis Dec 22, 2017

Contributor

Thank you so much for the review again, these were fantastic comments and very insightful

CameronTEllis added some commits Dec 22, 2017

In response to Mingbo's review: now downsample after convolution (as …
…before) in convolve_hrf, changed some defaults and updated the test. Also added some more plots to the example
@lcnature
Copy link
Contributor

lcnature left a comment

Looks pretty good! I added a few comments for your consideration.
Maybe in the final commit, use nbstripout to strip out the output in the example notebook (Is this still required @mihaic to pass the test?)

noise_drift[:, basis_counter - 1] = np.cos(timepoints_basis)

# Average the drift
noise_drift = np.mean(noise_drift, 1)

This comment has been minimized.

@lcnature

lcnature Jan 3, 2018

Contributor

Suggestion for future improvement: It is possible that the power of low frequency drift is higher than higher frequency drift. There may be some literature quantifying this property of fMRI data. It may be possible to weight fluctuations at different frequencies according to such pattern if it is reported in the literature, instead of putting equal weight. I am not expert on this. But I remember learning that fMRI power spectrum looks like 1/f.

This comment has been minimized.

@CameronTEllis

CameronTEllis Jan 4, 2018

Contributor

Good suggestion, I will keep it in mind

@@ -97,6 +99,9 @@ def test_generate_stimfunction():
signal_function = sim.convolve_hrf(stimfunction=stimfunction,
tr_duration=tr_duration,
)

max_response = np.where(signal_function != 0)[0].max()
assert 25 < max_response <= 30, "HRF is incorrect length"

This comment has been minimized.

@lcnature

lcnature Jan 3, 2018

Contributor

is -> has?

This comment has been minimized.

@CameronTEllis

CameronTEllis Jan 4, 2018

Contributor

Thanks!

# Down sample the stim function so that it only has one element per
# TR. This accelerates the convolution greatly
signal_vox = np.zeros((duration,))
for sample in list(range(duration)):

This comment has been minimized.

@lcnature

lcnature Jan 3, 2018

Contributor

I think the suggestion in this post might do the same as the for loop but may be more efficient.
https://stackoverflow.com/questions/15956309/averaging-over-every-n-elements-of-a-numpy-array
Just be careful to make sure that the reshape works, i.e., the length of signal_vox is dividable by stride

@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Jan 3, 2018

The recommendation to use nbstripout is not automatically enforced. While we accepted non-stripped notebooks, we very much prefer not to pollute Git with non-code data, such as images. The long-term solution to using notebooks is something like sphinx-gallery (see PR #270).

if upsampled_onset - np.round(upsampled_onset) != 0:
# Because of float precision, there can be issues. E.g.
# float('1.001') * 1000 = 1000.99
if np.allclose(upsampled_onset, np.round(upsampled_onset)):

This comment has been minimized.

@lcnature

lcnature Jan 8, 2018

Contributor

This is great. You can even set an additional tolerance for the allclose to make it more permissive.
Actually I think it may be also helpful to add an optional input flag to decide whether such check is necessary. Maybe the user is fine to ignore a few millisecond difference? Such flag can save them trouble of regenerating timing files that round exactly to the temporal resolution.

This comment has been minimized.

@CameronTEllis

CameronTEllis Jan 8, 2018

Contributor

Good suggestion, I have made it a warning now so the user can decide. Thanks for all the help!

@lcnature

This comment has been minimized.

Copy link
Contributor

lcnature commented Jan 9, 2018

@CameronTEllis @mihaic I think the code looks pretty good to me. I checked the error message. It is from my test code of the gen_design function in utils. I think it arises because in the past, the down-sampling after convolution in fmrisim was a direct down-sampling, while it is an averaging in a time window in the new code. I included that original test code based on the old fmrisim implementation. I was not to sure exactly which approach match the MRI physics better so I asked mrihelp yesterday. Mark told me that the acquisition time of a slice (ignoring multi-band) is the TR divided by the number of slices. So for any single slice, it seems that the old approach of direct down-sampling makes sense. But for slice-timing corrected image, interpolation is involved to infer the BOLD signal at one particular time for all slices. To me, this would be a bit close to the average signal within that TR, but not exactly. So I am still not sure whether direct down-sampling or averaging within a TR is better for your code -- both make some senses. Maybe it is up to you to choose one of them. If you stick with averaging, then I think deleting Lines 124-128 (

design5 = gen_design(stimtime_files=[files['FSL2']],
) in test/utils/test_utils.py would resolve the error message. If you opt for direct down-sampling as you did before, then I think the error message will go away automatically.

@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Jan 9, 2018

@CameronTEllis, when implementing @lcnature's suggestion, please document your choice in a docstring.

@CameronTEllis

This comment has been minimized.

Copy link
Contributor

CameronTEllis commented Jan 9, 2018

@lcnature,

This is an interesting nuance. As I see it, the contrast is between saying the signal within a TR is either represented by the activity at single time point or it reflects the average of the activity within a TR. Practically, these are small differences because the HRF is slow relative to the length of a TR. However, as you point out there are substantive theoretical differences.

In an EPI we are indeed measuring voxel intensity in a slice by slice manner and this acquisition time is very brief. So in some sense the signal that a voxel represents should come from a single timepoint of the HRF. I believe the first approach is correct, but inaccurate if applied by downsampling as suggested.

If we downsampled by taking every nth timepoint, we would (re-)introduce a substantial bias (remember this is what we previously had as the default). For instance, if you have a signal function (estimate of the HRF convolved with the event sequence) that is 100,000 elements, and you want to down sample to 100 TRs, each of 1s duration, then you might take the 0th element, 100th element, 200th ... etc to represent your expected activity for each TR. This would effectively mean you are treating this and every voxel like it as something like the first acquired slice. You could pick a different nth element (the 50th or the 99th) but either way you are making an assumption that all signal is being acquired at the same time. This is also a problem because it can effectively shift the timing of events: sampling the 0th point is a TR different that sampling the 99th time point, so your peak response would come later in the 0th version than the 99th version.

As mentioned, we are in fact sampling data at multiple times throughout a TR. This means that taking the average of the signal function over that time window is approximately the same as if the voxels containing signal come from a distribution of voxels uniformly distributed across the slices. However, this is also often wrong, especially when we are doing whole brain acquisitions, since there is often non brain in the first and last slices. If we are simulating interleaved or multiband slice acquisition sequences this would be less of a problem but still present.

The most optimal way of designing the simulator would be to sample at a different time within a TR depending on the slice of the voxel being modelled. In other words, if the voxel is in the slice acquired first then the signal function should be sampled early in the TR, if the slice is last it should be sampled late in the TR. I could add this option to the code (provide a slice time of the voxel for each signal that is being simulated). This would be good additional functionality that would provide a lot of flexibility; although as mentioned before it would make little practical difference compared to averaging (I would default to take the middle time point).

If you agree with this logic I will go ahead and make this change.

Thanks,
Cameron

@mihaic mihaic added this to the v0.7 milestone Jan 9, 2018

@lcnature

This comment has been minimized.

Copy link
Contributor

lcnature commented Jan 9, 2018

Hi @CameronTEllis
Yes, I agree that for most of the usage, the two would not make much difference. The only case that they may differ a lot is when people use super long TR in some auditory task in order to make the auditory stimuli and the scanner noise far apart in time.

And yes, if the simulator capture the acquisition time difference between slices, it would be much closer to the reality. But of course this can be left for future extension if you want to do it.

But I think the worry about choosing to sample at 0th or 99th interval is less of concern. This is because if the simulator assume the same acquisition time for all slices, the user should make sure the simulated data are used to represent slice-timing corrected data. And in practice when doing slice-timing correction, you always specify which slice all other slices are aligned to. As long as the same timing choice is made for both the simulator and the real data, choosing the 0th or 99th slice are both correct.
The down-sampling here is different from the issue with down-sampling in the previous PR. Because last time the down-sampling was done to the event-timing series, in which sparse binary events may be entirely missed by the down-sampling. But here we are down-sampling at the simulated fMRI signal, which is already a smooth signal. (and that is why down-sampling and averaging are probably very close for most of the usage cases)

Maybe another (minor) argument for using direct down-sampling here is that respiratory or heart beat noise are sampled at discrete time instead of averaged within a TR in your code, if I understand correctly. Maybe it makes slightly more sense to use the same approach for both signal and noise.

Other people behind mrihelp just stopped by my office and some suggestion they gave include checking what other packages like FSL choose to do... They are happy to answer any questions you have, btw.

@CameronTEllis

This comment has been minimized.

Copy link
Contributor

CameronTEllis commented Jan 10, 2018

@lcnature,

Thank you for the follow up. I think your suggestion to simulate the data 'as if' it is slice time corrected is a good one. FSL does slice time correction by aligning to the middle time point in a TR so if I downsample at the middle time point then this will resolve this. This is a much simpler solution than providing a slice time option.

You are right about the difference in downsampling the event time course rather than the convolved response, my mistake.

Fortunately this will be a very minor change since the mean and middle time point are equal when the trend within a TR is linear (which is not always the case but is often close).

I will implement this now and re-commit the code. I will make sure to note this in the docstring

Thanks,
Cameron

CameronTEllis added some commits Jan 10, 2018

@@ -123,7 +123,7 @@ def test_gen_design():
'gen_design does not treat missing values correctly')
design5 = gen_design(stimtime_files=[files['FSL2']],
scan_duration=[48, 20], TR=1)
assert np.all(np.isclose(design4, design5[::2])), (
assert (design4 - design5[::2]).mean() < 0.1, (

This comment has been minimized.

@lcnature

lcnature Jan 11, 2018

Contributor

Is this change to threshold of 0.1 still necessary given you have used direct down-sampling?
If it is still necessary. Maybe it should be the absolute value of difference that is smaller than some threshold, and 0.1 on the average difference may be too loose.

This comment has been minimized.

@CameronTEllis

CameronTEllis Jan 11, 2018

Contributor

Up to you to decide whether it is still necessary.

I tested this value and it seemed in an appropriate range for the average difference. Do you have an alternative suggestion?

This comment has been minimized.

@lcnature

lcnature Jan 11, 2018

Contributor

Ah I see, I thought they will become exactly the same so np.all(np.isclose...) can actually work. But maybe I am wrong.
Oh I think I know why, you are using the time of the middle slice. So then design4 and design5 are guaranteed to be different (they are off by 0.5 s). Maybe we should just remove this test between design4 and design5, given your current implementation.

This comment has been minimized.

@CameronTEllis

CameronTEllis Jan 12, 2018

Contributor

Yes that is exactly right, the middle slice makes it different than if it was the first or last slice. I actually like this test though since they should be similar, just not the exact same. I vote to keep it in so that we know future changes to these functions don't have a big effect

This comment has been minimized.

@lcnature

lcnature Jan 12, 2018

Contributor

OK! Let's keep it. Does assert np.abs(design4 - design5[::2]).mean() < 0.1 work? I just think the mean should be applied to some metric of distance, rather than just the subtraction. If not, I think decreasing TR would make it work.

@lcnature
Copy link
Contributor

lcnature left a comment

Pretty good! Just not sure about the testing code.

CameronTEllis and others added some commits Jan 12, 2018

@lcnature

This comment has been minimized.

Copy link
Contributor

lcnature commented Jan 12, 2018

@mihaic The PR looks fine for me now! Thanks for your effort @CameronTEllis !

@mihaic

This comment has been minimized.

Copy link
Contributor

mihaic commented Jan 12, 2018

@lcnature, you have the right to squash and merge the commits to master yourself. First you have to update your review to formally approve the PR.

@lcnature lcnature merged commit eefb1de into brainiak:master Jan 13, 2018

2 of 3 checks passed

continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
linux Build finished.
Details
macos Build finished.
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment