diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 8283a17b8..94c51d31f 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -90,6 +90,7 @@ from scipy import signal import scipy.ndimage as ndimage import copy +from scipy import optimize __all__ = [ "apply_signal", @@ -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, ): @@ -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 ---------- @@ -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': + + # 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) + + # 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) + + # 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) @@ -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: @@ -2285,6 +2360,7 @@ def _noise_dict_update(noise_dict): def _fit_spatial(noise, noise_temporal, + drift_noise, mask, template, spatial_sd, @@ -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? @@ -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 @@ -2423,6 +2503,7 @@ def _fit_temporal(noise, spatial_sd, temporal_proportion, temporal_sd, + drift_noise, noise_dict, fit_thresh, fit_delta, @@ -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 @@ -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 @@ -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, @@ -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[ + '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']) @@ -2765,7 +2866,8 @@ 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 @@ -2773,6 +2875,7 @@ def generate_noise(dimensions, # Fit the SNR noise, spatial_sd = _fit_spatial(noise, noise_temporal, + drift_noise, mask, template, spatial_sd, @@ -2792,6 +2895,7 @@ def generate_noise(dimensions, spatial_sd, temporal_proportion, temporal_sd, + drift_noise, noise_dict, fit_thresh, fit_delta, diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 2cfb3b8e3..6b52986f6 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -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 @@ -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'