Skip to content

Drift update #437

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.

Sign up for GitHub

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

Merged
CameronTEllis merged 23 commits into from
Dec 11, 2019
Merged

Drift update #437

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
a50bbef
Removed print statement
Sep 15, 2019
c82e57f
Update to how drift is calculated and used. The new default for drift…
Sep 15, 2019
13db940
Update to how drift is calculated and used. The new default for drift…
Sep 15, 2019
0dbf242
Fixed a math error and made the simulated drift include the 1% of pow…
Sep 28, 2019
9c1342b
Add command explicitly calling cos_power_drop
Sep 28, 2019
79e7457
Remove warning for short runs since this does not apply to this basis…
Oct 5, 2019
3ace12f
Remove warning for short runs since this does not apply to this basis…
Oct 5, 2019
b3140a8
Merge branch 'master' into drift_update
CameronTEllis Oct 13, 2019
ab3297f
Update for MCI's comments
Nov 5, 2019
5d229b3
Update for MCI's comments
Nov 5, 2019
c632237
Merge branch 'master' into drift_update
CameronTEllis Nov 9, 2019
192a305
Added more detail to the power drop doc string
Dec 3, 2019
d2526d9
Merge branch 'drift_update' of https://github.com/CameronTEllis/Brain…
Dec 3, 2019
5ec1c63
Change the maximum frequency to be equivalent to the TR number
Dec 5, 2019
e21342d
Make power drop unitless and add documentation. Also add test for per…
Dec 5, 2019
d80b89e
Merge branch 'master' into drift_update
lcnature Dec 6, 2019
3109a19
Fix pep8 issue
Dec 6, 2019
79f4b9b
Merge branch 'drift_update' of https://github.com/CameronTEllis/Brain…
Dec 6, 2019
debb8b2
Fix pep8 issue
Dec 6, 2019
eb82c55
Update of scikit-learn to 0.22 changed how randomizer worked, renderi…
Dec 9, 2019
88ee76e
Revert test
Dec 11, 2019
251e2bb
Merge branch 'master' into drift_update
CameronTEllis Dec 11, 2019
59346f1
Finish reverting test_mvpa_voxel_selection
mihaic Dec 11, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 124 additions & 20 deletions brainiak/utils/fmrisim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@
from scipy import signal
import scipy.ndimage as ndimage
import copy
from scipy import optimize

__all__ = [
"apply_signal",
Expand Down Expand Up @@ -1523,7 +1524,7 @@ def _generate_noise_temporal_task(stimfunction_tr,

def _generate_noise_temporal_drift(trs,
tr_duration,
basis="discrete_cos",
basis="cos_power_drop",
period=150,
):

Expand All @@ -1544,10 +1545,17 @@ def _generate_noise_temporal_drift(trs,
basis : str
What is the basis function for the drift. Could be made of discrete
cosines (for longer run durations, more basis functions are
created) or a sine wave.
created) that either have equal power ('discrete_cos') or the power
diminishes such that 99% of the power is below a specified frequency
('cos_power_drop'). Alternatively, this drift could simply be a sine
wave ('sine')

period : int
How many seconds is the period of oscillation of the drift
When the basis function is 'cos_power_drop' this is the period over
which no power of the drift exceeds (i.e. the power of the drift
asymptotes at this period). However for the other basis functions,
this is simply how many seconds is the period of oscillation of the
drift

Returns
----------
Expand Down Expand Up @@ -1597,6 +1605,86 @@ def _generate_noise_temporal_drift(trs,
phase = (timepoints / (trs - 1) * cycles * 2 * np.pi) + phaseshift
noise_drift = np.sin(phase)

elif basis == 'cos_power_drop':
Comment thread
lcnature marked this conversation as resolved.

# Make a vector counting each TR
timepoints = np.linspace(0, trs - 1, trs) * tr_duration

# Specify the other timing information
duration = trs * tr_duration

# How bases do you have? This is to adhere to Nyquist
basis_funcs = int(trs)

noise_drift = np.zeros((timepoints.shape[0], basis_funcs))
for basis_counter in list(range(1, basis_funcs + 1)):
# 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)

@lcnature lcnature Nov 11, 2019 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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?

@CameronTEllis CameronTEllis Nov 19, 2019

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

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

@CameronTEllis CameronTEllis Nov 19, 2019

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

If this makes sense, is there something else I could add to the comments to clarify?

@lcnature lcnature Dec 2, 2019

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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).

@lcnature lcnature Dec 2, 2019 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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?

@cbaldassano cbaldassano Dec 5, 2019 •
edited
Loading

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

@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.

@CameronTEllis CameronTEllis Dec 5, 2019

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

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)
(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?

@cbaldassano cbaldassano Dec 6, 2019

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

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...

@CameronTEllis CameronTEllis Dec 6, 2019

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Aha, okay thanks for that clarification. The code should then be ready to go, what do you think @lcnature ?

@lcnature lcnature Dec 6, 2019

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Looks great. I am glad to approve it!


# In radians, what is the value for each time point
timepoints_basis = timepoint_phase + random_phase

# Store the drift from this basis func
noise_drift[:, basis_counter - 1] = np.cos(timepoints_basis)

def power_drop(r, L, F, tr_duration):
# Function to return the drop rate for the power of basis functions
# In other words, how much should the weight of each basis function
# reduce in order to make the power you retain of the period's
# frequency be 99% of the total power of the highest frequency, as
# defined by the DCT.
# For an example where there are 20 time points, there will be 20
# basis functions in the DCT. If the period of the signal you wish
# to simulate is such that 99% of the power should drop off after
# the equivalent of 5 of these basis functions, then the way this
# code works is it finds the rate at which power must drop off for
# all of the 20 basis functions such that by the 5th one, there is
# only 1% of the power remaining.
# r is the power reduction rate which should be between 0 and 1
# L is the duration of the run in seconds
# F is period of the cycle in seconds It is assumed that this will
# be greater than the tr_duration, or else this will not work
# tr_duration is the duration of each TR in seconds

# Check the TR duration
if F < tr_duration:
msg = 'Period %0.0f > TR duration %0.0f' % ((F, tr_duration))
raise ValueError(msg)

percent_retained = 0.99 # What is the percentage of power retained

# Compare the power at the period frequency (in the numerator) with
# the power at the frequency of the DCT, AKA the highest possible
# frequency in the data (in the denominator)
numerator = 1 - r ** (2 * L / F) # Power of this period
denominator = 1 - r ** (2 * L / tr_duration) # Power of DCT freq.

# Calculate the retained power
power_drop = abs((numerator / denominator) - percent_retained)
return power_drop

# Solve for power reduction rate.
# This assumes that r is between 0 and 1
# Takes the duration and period as arguments
sol = optimize.minimize_scalar(power_drop,
bounds=(0, 1),
method='Bounded',
args=(duration, period, tr_duration))

# Pull out the solution
r = sol.x

# Weight the basis functions based on the power drop off
basis_weights = r ** np.arange(basis_funcs)
Comment thread
CameronTEllis marked this conversation as resolved.

# Weigh the basis functions
weighted_basis_funcs = np.multiply(noise_drift, basis_weights)

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

# Normalize so the sigma is 1
noise_drift = stats.zscore(noise_drift)

Expand Down Expand Up @@ -2012,19 +2100,6 @@ def _generate_noise_temporal(stimfunction_tr,
# Preset the volume
noise_volume = np.zeros((dimensions[0], dimensions[1], dimensions[2], trs))

# Generate the drift noise
if noise_dict['drift_sigma'] != 0:
# Calculate the drift time course
noise = _generate_noise_temporal_drift(trs,
tr_duration,
)
# Create a volume with the drift properties
volume = np.ones(dimensions)

# Combine the volume and noise
noise_volume += np.multiply.outer(volume, noise) * noise_dict[
'drift_sigma']

# Generate the physiological noise
if noise_dict['physiological_sigma'] != 0:

Expand Down Expand Up @@ -2285,6 +2360,7 @@ def _noise_dict_update(noise_dict):

def _fit_spatial(noise,
noise_temporal,
drift_noise,
mask,
template,
spatial_sd,
Expand All @@ -2306,6 +2382,9 @@ def _fit_spatial(noise,
noise_temporal : multidimensional array, float
The temporal noise that was generated by _generate_temporal_noise

drift_noise : multidimensional array, float
The drift noise generated by _generate_noise_temporal_drift

tr_duration : float
What is the duration, in seconds, of each TR?

Expand Down Expand Up @@ -2400,7 +2479,8 @@ def _fit_spatial(noise,
)

# Sum up the noise of the brain
noise = base + (noise_temporal * temporal_sd) + noise_system
noise = base + drift_noise + noise_system
noise += (noise_temporal * temporal_sd) # Add the brain specific noise

# Reject negative values (only happens outside of the brain)
noise[noise < 0] = 0
Expand All @@ -2423,6 +2503,7 @@ def _fit_temporal(noise,
spatial_sd,
temporal_proportion,
temporal_sd,
drift_noise,
noise_dict,
fit_thresh,
fit_delta,
Expand Down Expand Up @@ -2466,6 +2547,9 @@ def _fit_temporal(noise,
What is the standard deviation in time of the noise volume to be
generated

drift_noise : multidimensional array, float
The drift noise generated by _generate_noise_temporal_drift

noise_dict : dict
A dictionary specifying the types of noise in this experiment. The
noise types interact in important ways. First, all noise types
Expand Down Expand Up @@ -2581,7 +2665,8 @@ def _fit_temporal(noise,
)

# Sum up the noise of the brain
noise = base + (noise_temporal * temporal_sd) + noise_system
noise = base + drift_noise + noise_system
noise += (noise_temporal * temporal_sd) # Add the brain specific noise

# Reject negative values (only happens outside of the brain)
noise[noise < 0] = 0
Expand Down Expand Up @@ -2737,7 +2822,7 @@ def generate_noise(dimensions,
# What is the mean signal of the non masked voxels in this template?
mean_signal = (base[mask > 0]).mean()

# Generate the noise
# Generate the temporal noise
noise_temporal = _generate_noise_temporal(stimfunction_tr=stimfunction_tr,
tr_duration=tr_duration,
dimensions=dimensions,
Expand All @@ -2746,6 +2831,22 @@ def generate_noise(dimensions,
noise_dict=noise_dict,
)

# Generate the drift noise
if noise_dict['drift_sigma'] != 0:
# Calculate the drift time course
noise = _generate_noise_temporal_drift(len(stimfunction_tr),
tr_duration,
)
# Create a volume with the drift properties
volume = np.ones(dimensions[:3])

# Combine the volume and noise
drift_noise = np.multiply.outer(volume, noise) * noise_dict[
Comment thread
CameronTEllis marked this conversation as resolved.
'drift_sigma']
else:
# If there is no drift, then just make this zeros (in 4d)
drift_noise = np.zeros(dimensions_tr)

# Convert SFNR into the size of the standard deviation of temporal
# variability
temporal_sd = (mean_signal / noise_dict['sfnr'])
Expand All @@ -2765,14 +2866,16 @@ def generate_noise(dimensions,
)

# Sum up the noise of the brain
noise = base + (noise_temporal * temporal_sd) + noise_system
noise = base + drift_noise + noise_system
noise += (noise_temporal * temporal_sd) # Add the brain specific noise

# Reject negative values (only happens outside of the brain)
noise[noise < 0] = 0

# Fit the SNR
noise, spatial_sd = _fit_spatial(noise,
noise_temporal,
drift_noise,
mask,
template,
spatial_sd,
Expand All @@ -2792,6 +2895,7 @@ def generate_noise(dimensions,
spatial_sd,
temporal_proportion,
temporal_sd,
drift_noise,
noise_dict,
fit_thresh,
fit_delta,
Expand Down
14 changes: 12 additions & 2 deletions tests/utils/test_fmrisim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -534,13 +534,24 @@ def test_generate_noise():
assert power[period_freq] > power[period_freq + 1], 'Power is low'
assert power[period_freq] > power[period_freq - 1], 'Power is low'

# Check it gives a warning if the duration is too short
# Check it runs fine
drift = sim._generate_noise_temporal_drift(50,
tr_duration,
'discrete_cos',
period,
)

# Check it runs fine
drift = sim._generate_noise_temporal_drift(300,
tr_duration,
'cos_power_drop',
period,
)

# Check that when the TR is greater than the period it errors
with pytest.raises(ValueError):
sim._generate_noise_temporal_drift(30, 10, 'cos_power_drop', 5)

# Test physiological noise (using unrealistic parameters so that it's easy)
timepoints = list(np.linspace(0, (trs - 1) * tr_duration, trs))
resp_freq = 0.2
Expand Down Expand Up @@ -675,7 +686,6 @@ def test_generate_noise_spatial():

# Calculate the proportion of std relative to the mean
std_proportion = np.nanstd(fwhm3) / np.nanmean(fwhm3)
print(fwhm3)
assert std_proportion < 0.25, 'Variance is inconsistent across dim'


Expand Down