Drift update - #437
Conversation
… drops off power on the basis functions systematically. Drift is also scaled separately from the brain noise since that scale is based on a detrended version of drift
… drops off power on the basis functions systematically. Drift is also scaled separately from the brain noise since that scale is based on a detrended version of drift
|
@mihaic I didn't add any code that I didn't subsequently test. Is there any way in which this can be fixed without me having to add tests to other code? |
|
How about |
|
@CameronTEllis, I fixed the latest failure in PR #439. Update your branch from master. |
|
@mihaic Did I do something wrong? |
|
Jenkins, retest this please. |
|
That was my fault. The cleanup script was not preparing the repository properly before the test, so there were leftover files from another PR causing the error. I changed it in Jenkins. |
|
@lcnature Do you have time to review this? |
lcnature
left a comment
There was a problem hiding this comment.
@CameronTEllis Sorry for my delay. Please see my comments.
|
@mihaic Why is Travis crashing for ruby? |
| # What steps do you want to take for this basis function | ||
| random_phase = np.random.rand() * np.pi * 2 | ||
|
|
||
| timepoint_phase = (timepoints / duration * np.pi * basis_counter) |
There was a problem hiding this comment.
Here, basis_counter goes from 1 to basis_funs, which is twice the total duration of the scan, in the unit of a second. timepoints is equally spaced time points between 0 and duration. When it is divided by duration, the result is equally spaced numbers between 0 and 1. This times pi would result in equally spaced points between 0 and pi. But multiply this phase with basis_counter which is a duration seems a bit strange to me. It won't mean something as a radius anymore, but this variable is in turn added with random_phase and used as an input to cosine function. I think the unit here may be somehow wrong. Likely, the basis_funcs should be something like the ratio between the lowest frequency of the fft of the time series to the Nyquist frequency (twice of the sampling frequency of the scanner) instead, in line 1615?
There was a problem hiding this comment.
The reason for basis_funcs being equal to 2*duration is because the power drop is meant to drop off to 99% by basis_funcs at duration, with the remaining power being captured in the subsequent basis_funcs duration->2*duration.
The reason for having the basis_counter multiplier is simply to increase the number of cycles for each basis_counter. This means that for longer runs there are simulated basis functions of decreasing periodicity along with decreasing power
There was a problem hiding this comment.
If this makes sense, is there something else I could add to the comments to clarify?
There was a problem hiding this comment.
Thanks. Now I think I may have got it. You are trying to make a power spectrum that decreases exponentially from the lowest frequency (inverse of double the duration of the scan, corresponding to when basis_counter equals 1) to the highest frequency (the sampling frequency of fMRI, when basis_counter equals basis_funcs).
There was a problem hiding this comment.
But then I am confused with the description
'the power drop is meant to drop off to 99% by basis_funcs at duration, with the remaining power being captured in the subsequent basis_funcs duration->2*duration.'
Can I interpret it as, if I am integrating an exponential decaying curve from 0 up to certain number (duration) here, the integral is 99% of the integral from 0 to 2*duration? If this is the case, I would imagine the power of the numerator and denominator in the function power_drop should differ in a ratio of 2, instead of a ratio of F?
There was a problem hiding this comment.
@lcnature is right that the units aren't right there. The reasoning for having the sum be up to L isn't based on the Nyquist freq (though it is related), it is because that is the number of basis functions in the DCT. For example see the top of pg 2 here - u indexes the timepoints, m indexes the basis functions. The twos are coming from the geometric series equation, since power is proportional to r**2 (so these terms could be written, for example, (r**2)^(L/F)).
I believe the numerator is correct as is but the denominator needs to be put back into # of TRs units, e.g.
numerator = 1 - r ** (2 * L / F) # Power of this period
denominator = 1 - r ** (2 * L / tr_duration) # Power of all periods
This implies that the shortest possible period is tr_duration which also makes intuitive sense.
There was a problem hiding this comment.
Thanks for this clarification @cbaldassano , it helped my understanding. Because of it I added, these changes, some documentation, an example and a check to ensure that tr_duration < period. However, reading this did raise a concern about another line of code
brainiak/brainiak/utils/fmrisim.py
Line 1679 in e21342d
basis_weights = r ** np.arange(basis_funcs)). This line is meant to find weights for each of the basis functions before they are all added together; however, based on my new understanding shouldn't this be: basis_weights = r ** (2 * np.arange(basis_funcs)) in order to make power proportional?
There was a problem hiding this comment.
No, the power spectrum is the square of the coefficients. So the coefficients themselves are 1, r, r^2, r^3... and the power spectrum of the frequencies is 1, r^2, r^4, r^6...
There was a problem hiding this comment.
Aha, okay thanks for that clarification. The code should then be ready to go, what do you think @lcnature ?
There was a problem hiding this comment.
Looks great. I am glad to approve it!
|
|
||
| percent_retained = 0.99 # What is the percentage of drift retained | ||
| numerator = 1 - r ** (4 * L / F) # Power of this period | ||
| denominator = 1 - r ** (4 * L) # Power of all periods |
There was a problem hiding this comment.
According to my comment on Line 1658, I now guess r is the reduction rate of the weight of each frequency component from its previous one. So the power it would be raised to would be a unit-less number, not a duration. But L is the duration of the scan, which seems to be should not be a power. I am also not clear why 4 is multiplied. If the amplitude is decreased by r, then we can say the power is decreased by r^2. So I can understand if 2 is used, but not 4.
There was a problem hiding this comment.
@cbaldassano caught this issue, if you don't use 4 then the math won't work out. I may be wrong but I think it is because you need to square the increase on L, effectively 2**2
There was a problem hiding this comment.
Hmm, let me make some guess to see if I am understanding correctly. In the denominator, 4 * L is 2 * (2 * L), where 2 * L is 2 * duration. This would make sense given the codes above where you have set the highest basis_counter to be 2 * duration. And the 4 * L / F is essentially 2 * (1/F / (1/2L)), where 1/F is the frequency corresponding to the period you consider retaining 99% of the power below it, and 1/2L is the lowest frequency in your basis set. So now both of these make sense to me. My only concern is whether setting the highest frequency as 2*L makes sense. They are of different units to start with: 2L is a quantity of time, choosing its number as a frequency seems strange to me.
There was a problem hiding this comment.
Glad it makes sense, will wait on the above comment to see if this still seems workable
|
@CameronTEllis, fixed with PR #443. |
|
Thanks @lcnature for your careful review of this! |
|
@mihaic Failed because of error in the testing of FCMA, which I didn't edit. I think this was a statistical fluke, should I just retest? |
|
It's an issue we are working on, issue #448. |
…ng this wrong. Dropped hard coding but maintained spirit of the test
|
@CameronTEllis, it looks to me that the existing test was looking for some specific output, not for random output. Did I misunderstand that? |
|
@mihaic correct. The test generates random data and then runs an SVM with pseudo-random labels. The test Yida wrote is based on running this classifier once, checking what results it gets and then using those results to check that future tests get the same result. In this code there is some randomness in terms of how the CV is set, but this was controlled previously by setting the seed. However, I think the update added an additional call to the rng which means that the same code will not produce the same results. Regardless, the spirit of the test is whether the code outputs something of a certain shape and then asks whether it is at chance. I can't think of why getting the exact same result as before is necessary for testing the code: you could change the code in problematic ways and still (sometimes) get the same results or you could change the code in non-problematic ways (e.g. extra call to rng) and not get the same results. Moreover, this is not a standard we apply consistently in our tests because speaking from experience the tests in fmrisim do not require this. |
|
@CameronTEllis, the error I see is about the length of the results, not the contents. I would rather pin scikit-learn so we can resume merging and later try harder at preserving the existing comparison. What do you think? |
|
@mihaic The length discrepancy is a misdirection. The code runs a searchlight on a mask and computes the classification accuracy in each voxel within the mask mask. If the voxel is outside of the searchlight it doesn't return a |
| expected_output = [6, 6, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, | ||
| 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 1] | ||
|
|
||
| # Check that the mean accuracy with noise is close to chance |
There was a problem hiding this comment.
@CameronTEllis, you missed a line when you undid the test change. Next time, try git revert.
Update to how drift is calculated and used. The new default for drift drops off power on the basis functions systematically. Drift is also scaled separately from the brain noise since that scale is based on a detrended version of drift