From c002be12dd29e4b937c3a0a0ea2de7159a925cfd Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 29 Mar 2017 12:58:06 -0400 Subject: [PATCH 1/8] Substantial edits to how noise is calculated. Some edits to how noise is created. Altered test and example to suit update --- brainiak/utils/fmrisim.py | 337 ++++++++++++++++++++---------- examples/utils/fmrisim_example.py | 15 +- tests/utils/test_fmrisim.py | 32 ++- 3 files changed, 250 insertions(+), 134 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index e387518ce..df6204556 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -50,17 +50,21 @@ Show the brain as it unfolds over time with a given opacity. - Authors: Cameron Ellis (Princeton) 2016 + Authors: + Cameron Ellis (Princeton) 2016-2017 + Chris Boldassano (Princeton) 2016-2017 """ import logging from itertools import product +import os import nitime.algorithms.autoregressive as ar import math import numpy as np from pkg_resources import resource_stream from scipy import stats import scipy.ndimage as ndimage +import nibabel __all__ = [ "generate_signal", @@ -630,8 +634,9 @@ def double_gamma_hrf(stimfunction, signal_function = signal_function[0::decimate_interval] # Cut off the HRF - signal_function = signal_function[0:int(len(stimfunction) / tr_duration - / temporal_resolution)] + signal_function = signal_function[0:int((len(stimfunction) / + tr_duration) / + temporal_resolution)] # Scale the function so that the peak response is 1 if scale_function == 1: @@ -678,9 +683,10 @@ def apply_signal(signal_function, def _calc_fwhm(volume, mask, + voxel_size=[1.0, 1.0, 1.0], ): """ Calculate the FWHM of a volume - Takes in a volume and mask and outputs the FWHM of each TR of this + Takes in a volume and mask and outputs the FWHM (mm) of each TR of this volume for the non-masked voxels Parameters @@ -691,11 +697,14 @@ def _calc_fwhm(volume, mask : 4 dimensional array A mask of the voxels to have the FWHM measured from + voxel_size : length 3 list, float + Millimeters per voxel for x, y and z. + Returns ------- float, list - Returns the FWHM of each TR """ + Returns the FWHM of each TR in mm""" # What are the dimensions of the volume dimensions = volume.shape @@ -716,8 +725,13 @@ def _calc_fwhm(volume, for tr_counter in list(range(0, trs)): # Preset - sum_sq_der = [0.0, 0.0, 0.0] - countSqDer = [0, 0, 0] + v_count = 0 + v_sum = 0 + v_sq = 0 + + d_sum = [0.0, 0.0, 0.0] + d_sq = [0.0, 0.0, 0.0] + d_count = [0, 0, 0] # If there is more than one TR, just pull out that volume if trs > 1: @@ -737,31 +751,55 @@ def _calc_fwhm(volume, # Is this within the mask? if mask[x, y, z, tr_counter] > 0: + # Find the the volume sum and squared values + v_count += 1 + v_sum += np.abs(volume[x, y, z]) + v_sq += volume[x, y, z] ** 2 + + # Get the volume variance + v_var = (v_sq - ((v_sum ** 2) / v_count)) / (v_count - 1) + + for i in list(range(len(coordinates))): + + # Pull out this coordinate + x, y, z = coordinates[i] + + # Is this within the mask? + if mask[x, y, z, tr_counter] > 0: # For each xyz dimension calculate the squared # difference of this voxel and the next - if x < dimensions[0] - 1 and mask[x + 1, y, z, - tr_counter] > 0: - temp = (volume[x, y, z] - volume[x + 1, y, z]) ** 2 - sum_sq_der[0] += temp - countSqDer[0] += + 1 - - if y < dimensions[1] - 1 and mask[x, y + 1, z, - tr_counter] > 0: - temp = (volume[x, y, z] - volume[x, y + 1, z]) ** 2 - sum_sq_der[1] += temp - countSqDer[1] += + 1 - - if z < dimensions[2] - 1 and mask[x, y, z + 1, - tr_counter] > 0: - temp = (volume[x, y, z] - volume[x, y, z + 1]) ** 2 - sum_sq_der[2] += temp - countSqDer[2] += 1 - - # What is the FWHM for each dimension - fwhm3 = np.sqrt(np.divide(4 * np.log(2), np.divide(sum_sq_der, - countSqDer))) - fwhm[tr_counter] = np.prod(fwhm3) ** (1 / 3) # Calculate the average + inrange = (x < dimensions[0] - 1) + inmask = inrange and (mask[x + 1, y, z, tr_counter] > 0) + included = inmask and (~np.isnan(volume[x + 1, y, z])) + if included: + d_sum[0] += volume[x, y, z] - volume[x + 1, y, z] + d_sq[0] += (volume[x, y, z] - volume[x + 1, y, z]) ** 2 + d_count[0] += 1 + + inrange = (y < dimensions[1] - 1) + inmask = inrange and (mask[x, y + 1, z, tr_counter] > 0) + included = inmask and (~np.isnan(volume[x, y + 1, z])) + if included: + d_sum[1] += volume[x, y, z] - volume[x, y + 1, z] + d_sq[1] += (volume[x, y, z] - volume[x, y + 1, z]) ** 2 + d_count[1] += 1 + + inrange = (z < dimensions[2] - 1) + inmask = inrange and (mask[x, y, z + 1, tr_counter] > 0) + included = inmask and (~np.isnan(volume[x, y, z + 1])) + if included: + d_sum[2] += volume[x, y, z] - volume[x, y, z + 1] + d_sq[2] += (volume[x, y, z] - volume[x, y, z + 1]) ** 2 + d_count[2] += 1 + + # Find the variance + d_var = np.divide((d_sq - np.divide(np.power(d_sum, 2), + d_count)), (np.add(d_count, -1))) + + o_var = np.divide(-1, (4 * np.log(1 - (0.5 * d_var / v_var)))) + fwhm3 = np.sqrt(o_var) * 2 * np.sqrt(2 * np.log(2)) + fwhm[tr_counter] = np.prod(np.multiply(fwhm3, voxel_size)) ** (1 / 3) return fwhm @@ -800,24 +838,35 @@ def _calc_snr(volume, # Divide by the mask (if the middle slice was low in grey matter mass # then you would have a lower mean signal by default) - mean_signal = (slice_volume[slice_mask > 0] / slice_mask[slice_mask > - 0]).mean() + mean_signal = (slice_volume[slice_mask > 0]).mean() - # What is the overall noise - overall_noise = float(slice_volume[slice_mask == 0].std()) + # What is the noise in the non brain + background_noise = float(slice_volume[slice_mask == 0].std()) + + # What is the max activation of this volume + max_activity = volume.max() + + # What is the brain voxel by time activity. + brain_voxels = volume[mask[:, :, :, 0] > 0] + + # Find the noise to brain voxels + temporal_noise = np.std(brain_voxels, 1).mean() + + # Convert temporal noise into percent signal change + temporal_noise = temporal_noise / mean_signal * 100 # Calculate snr - snr = mean_signal / overall_noise + snr = mean_signal / background_noise # Convert from memmap snr = float(snr) - return overall_noise, snr + return temporal_noise, snr, max_activity def _calc_autoregression(volume, mask, - ar_order=2, + auto_reg_order=1, ): """ Calculate the autoregressive sigma of the data. @@ -829,7 +878,7 @@ def _calc_autoregression(volume, mask : 4d array, float What voxels of the input are within the brain - ar_order : int + auto_reg_order : int What order of the autoregression do you want to pull out Returns @@ -837,13 +886,17 @@ def _calc_autoregression(volume, float A sigma of the autoregression in the data + float + Sigma of the drift in the data + """ # Calculate the time course timecourse = np.mean(volume[mask[:, :, :, 0] > 0], 0) # Pull out the AR values (depends on order) - auto_reg_sigma = ar.AR_est_YW(timecourse, ar_order)[0][1] + auto_reg_sigma = ar.AR_est_YW(timecourse, auto_reg_order)[1] + auto_reg_sigma = np.sqrt(auto_reg_sigma) # What is the size of the change in the time course drift_sigma = timecourse.std().tolist() @@ -853,6 +906,7 @@ def _calc_autoregression(volume, def calc_noise(volume, mask=None, + noise_dict={'voxel_size': [1.0, 1.0, 1.0]}, ): """ Calculates the noise properties of the volume supplied. This estimates what noise properties the volume has. For instance it @@ -879,7 +933,6 @@ def calc_noise(volume, """ # Preset - noise_dict = {} # Create the mask if mask is None: @@ -888,21 +941,37 @@ def calc_noise(volume, # Since you are deriving the 'true' values then you want your noise to # be set to that level - # Calculate the overall noise and SNR of the volume - noise_dict['overall'], noise_dict['snr'] = _calc_snr(volume, mask) + # Calculate the temporal_noise noise and SNR of the volume + noise_dict['temporal_noise'], noise_dict['snr'], noise_dict[ + 'max_activity'] = _calc_snr(volume, mask) - # Calculate the fwhm of the data and use the average as your spatial sigma - noise_dict['spatial_sigma'] = _calc_fwhm(volume, mask).mean() + # Calculate the fwhm on a subset of volumes + + if volume.shape[3]>100: + # Take only 100 shuffled TRs + trs = np.arange(volume.shape[3]) + np.random.shuffle(trs) + trs = trs[0:100] + else: + trs = list(range(0, volume.shape[3])) + + noise_dict['fwhm'] = _calc_fwhm(volume[:, :, :, trs], + mask[:, :, :, trs], + noise_dict['voxel_size'], + ).mean() # Calculate the autoregressive and drift noise auto_reg_sigma, drift_sigma = _calc_autoregression(volume, mask) - noise_dict['auto_reg_sigma'] = auto_reg_sigma / noise_dict['overall'] - noise_dict['drift_sigma'] = drift_sigma / noise_dict['overall'] + # Calibrate for how sigma is originally calculated + auto_reg_sigma = auto_reg_sigma / noise_dict['temporal_noise'] - # Calculate the white noise - noise_dict['system_sigma'] = float(volume[mask == 0].std() / noise_dict[ - 'overall']) + # Total temporal noise, since these values only make sense relatively + total_temporal_noise = auto_reg_sigma + drift_sigma + + # What proportion of noise is accounted for by these variables? + noise_dict['auto_reg_sigma'] = auto_reg_sigma / total_temporal_noise + noise_dict['drift_sigma'] = drift_sigma / total_temporal_noise # Return the noise dictionary return noise_dict @@ -930,18 +999,18 @@ def _generate_noise_system(dimensions_tr, """ # Generate the Rician noise - noise_rician = stats.rice.rvs(1, 1, size=dimensions_tr) - - # Apply the gaussian noise - noise_gaussian = np.random.normal(0, 1, size=dimensions_tr) + noise_rician = stats.rice.rvs(b=0, loc=0, scale=1.527, size=dimensions_tr) + # + # # Apply the gaussian noise + # noise_gaussian = np.random.normal(0, 1, size=dimensions_tr) + # + # # Combine these two noise types + # noise_system = noise_rician + noise_gaussian + # + # # Normalize + # noise_system = stats.zscore(noise_system) - # Combine these two noise types - noise_system = noise_rician + noise_gaussian - - # Normalize - noise_system = stats.zscore(noise_system) - - return noise_system + return noise_rician def _generate_noise_temporal_task(stimfunction_tr, @@ -1130,7 +1199,7 @@ def _generate_noise_temporal_phys(timepoints, def _generate_noise_spatial(dimensions, mask=None, - spatial_sigma=-4.0, + fwhm=4.0, ): """Generate code for Gaussian Random Fields. @@ -1150,16 +1219,22 @@ def _generate_noise_spatial(dimensions, mask : 3d array The mask describing the boundaries of the brain - spatial_sigma : float - Approximately, what is the FWHM of the gaussian fields being - created. Use calc_fwhm to check this. - This is approximate for two reasons: firstly, although the + fwhm : float + What is the full width half max of the gaussian fields being created. + This is converted into a sigma which is used in this function. + However, this conversion was found empirically by testing values of + sigma and how it relates to fwhm values. The relationship that would be + found in such a test depends on the size of the brain (bigger brains + can have bigger fwhm). + However, small errors shouldn't matter too much since the fwhm + generated here can only be approximate anyway: firstly, although the distribution that is being drawn from is set to this value, this will manifest differently on every draw. Secondly, because of the masking and dimensions of the generated volume, this does not - behave simply. For instance, values over 10 do not work well because - the size of the volume. Moreover, wrapping effects matter (the - outputs are closer to this value if you have no mask). + behave simply- wrapping effects matter (the outputs are + closer to this value if you have no mask). + Use _calc_fwhm on this volume alone if you have concerns about the + accuracy of the fwhm. Returns ---------- @@ -1168,6 +1243,17 @@ def _generate_noise_spatial(dimensions, Generates the spatial noise volume for these parameters """ + if len(dimensions) == 4: + return + + # Convert from fwhm to sigma (relationship discovered empirical, only an + # approximation up to sigma=0-5 which corresponds to fwhm=0-8, relies + # on an assumption of brain size). + def logcurve(x, a, b, c): + return (np.log(x + a) / np.log(b)) + c + + spatial_sigma = logcurve(fwhm, -0.36778719, 2.10601011, 2.15439247) + # Set up the input to the fast fourier transform def fftIndgen(n): a = list(range(0, int(n / 2 + 1))) @@ -1214,16 +1300,17 @@ def _generate_noise_temporal(stimfunction_tr, tr_duration, dimensions, mask, - spatial_sigma, + fwhm, motion_sigma, drift_sigma, auto_reg_sigma, physiological_sigma, ): - """Generate the signal dependent noise + """Generate the temporal noise - This noise depends on things like the signal or the timing of the - experiment. + To increase or decrease the amount of total noise change the + temporal_noise noise_dict entry. To change the relative mixing of the + noise components, change the sigma's specified below. Parameters ---------- @@ -1282,18 +1369,19 @@ def _generate_noise_temporal(stimfunction_tr, # Generate the volumes that will differ depending on the type of noise # that it will be used for - volume_drift = _generate_noise_spatial(dimensions=dimensions, - spatial_sigma=spatial_sigma, - ) + volume_drift = np.ones(dimensions) + # volume_drift = _generate_noise_spatial(dimensions=dimensions, + # fwhm=fwhm, + # ) volume_phys = _generate_noise_spatial(dimensions=dimensions, mask=mask, - spatial_sigma=spatial_sigma, + fwhm=fwhm, ) volume_autoreg = _generate_noise_spatial(dimensions=dimensions, mask=mask, - spatial_sigma=spatial_sigma, + fwhm=fwhm, ) # Multiply the noise by the spatial volume @@ -1316,11 +1404,14 @@ def _generate_noise_temporal(stimfunction_tr, ) volume_task = _generate_noise_spatial(dimensions=dimensions, mask=mask, - spatial_sigma=spatial_sigma, + fwhm=fwhm, ) noise_task_volume = np.multiply.outer(volume_task, noise_task) noise_temporal = noise_temporal + (noise_task_volume * motion_sigma) + # Finally, z score each voxel so things mix nicely + noise_temporal = stats.zscore(noise_temporal, 3) + return noise_temporal @@ -1350,7 +1441,7 @@ def mask_brain(volume, What is the threshold (0 -> 1) for including a voxel in the mask? mask_self : bool - If set to 1 then it makes a mask from the volume supplied (by + If set to true then it makes a mask from the volume supplied (by averaging across time points and changing the range). Returns @@ -1371,7 +1462,7 @@ def mask_brain(volume, mask_raw = np.load(mask_name) # Is the mask based on the volume - if mask_self is 1: + if mask_self is True: mask_raw = np.zeros([volume.shape[0], volume.shape[1], volume.shape[ 2], 1]) @@ -1380,8 +1471,13 @@ def mask_brain(volume, else: mask_raw[:, :, :, 0] = np.array(volume) - # Make sure the mask values range from 0 to 1 - mask_raw = mask_raw / mask_raw.max() + mask_max = volume.max() + else: + mask_max = 1 + + # Make sure the mask values range from 0 to 1 (make out of max of volume + # so that this is invertible later) + mask_raw = mask_raw / mask_max # If there is only one brain volume then make this a forth dimension if len(volume.shape) == 3: @@ -1420,7 +1516,12 @@ def _noise_dict_update(noise_dict): ---------- noise_dict : dict - A dictionary specifying the types of noise in this experiment + A dictionary specifying the types of noise in this experiment. The noise + types interact in important ways. First, all noise types ending with + sigma (e.g. motion sigma) are mixed together in + _generate_temporal_noise. These values describe the proportion of mixing + of these elements. However critically, temporal_noise is the parameter + that describes how much noise these components contribute to the brain. Returns ------- @@ -1431,32 +1532,24 @@ def _noise_dict_update(noise_dict): # Check what noise is in the dictionary and add if necessary. Numbers # determine relative proportion of noise - if 'overall' not in noise_dict: - noise_dict['overall'] = 0.1 + if 'temporal_noise' not in noise_dict: + noise_dict['temporal_noise'] = 5 if 'motion_sigma' not in noise_dict: - noise_dict['motion_sigma'] = 0.1 + noise_dict['motion_sigma'] = 0 if 'drift_sigma' not in noise_dict: - noise_dict['drift_sigma'] = 0.5 + noise_dict['drift_sigma'] = 0.45 if 'auto_reg_sigma' not in noise_dict: - noise_dict['auto_reg_sigma'] = 1 + noise_dict['auto_reg_sigma'] = 0.45 if 'physiological_sigma' not in noise_dict: noise_dict['physiological_sigma'] = 0.1 - if 'system_sigma' not in noise_dict: - noise_dict['system_sigma'] = 1 if 'snr' not in noise_dict: noise_dict['snr'] = 30 - if 'spatial_sigma' not in noise_dict: - if noise_dict['overall'] == 0: - noise_dict['spatial_sigma'] = 0.015 - else: - noise_dict['spatial_sigma'] = 0.015 / noise_dict['overall'] - - # Print the mixing of the noise - total_var = noise_dict['motion_sigma'] ** 2 - total_var += noise_dict['drift_sigma'] ** 2 - total_var += noise_dict['auto_reg_sigma'] ** 2 - total_var += noise_dict['physiological_sigma'] ** 2 - total_var += noise_dict['system_sigma'] ** 2 + if 'max_activity' not in noise_dict: + noise_dict['max_activity'] = 1000 + if 'voxel_size' not in noise_dict: + noise_dict['voxel_size'] = [1.0, 1.0, 1.0] + if 'fwhm' not in noise_dict: + noise_dict['fwhm'] = 4 return noise_dict @@ -1488,8 +1581,9 @@ def generate_noise(dimensions, The mask of the brain volume, using noise_dict : dictionary, float - This is a dictionary that must contain the key: overall. If there - are no other variables provided then it will default values + This is a dictionary which describes the noise parameters of the + data. If there are no other variables provided then it will default + values Returns ---------- @@ -1521,8 +1615,8 @@ def generate_noise(dimensions, tr_duration=tr_duration, dimensions=dimensions, mask=mask[:, :, :, 0], - spatial_sigma=noise_dict[ - 'spatial_sigma'], + fwhm=noise_dict[ + 'fwhm'], motion_sigma=noise_dict[ 'motion_sigma'], drift_sigma=noise_dict[ @@ -1533,17 +1627,34 @@ def generate_noise(dimensions, 'physiological_sigma'], ) - noise_system = (_generate_noise_system(dimensions_tr=dimensions_tr) - * noise_dict['system_sigma']) + # Create the base (this inverts the process to make the mask) + base = mask * noise_dict['max_activity'] + + # Set the amount of background based on the snr value + + # What are the midpoints to be extracted + mid_x_idx = int(np.ceil(base.shape[0] / 2)) + mid_tr_idx = int(np.ceil(base.shape[3] / 2)) + + # Pull out the slices + slice_volume = base[mid_x_idx, :, :, mid_tr_idx] + slice_mask = mask[mid_x_idx, :, :, mid_tr_idx] + + # What is the mean signal of the non masked voxels in this slice? + mean_signal = (slice_volume[slice_mask > 0]).mean() + + # What is the standard deviation of the background activity + system_sigma = mean_signal / noise_dict['snr'] - # Find the outer product for the brain and nonbrain and add white noise - noise = noise_temporal + noise_system + # Set up the machine noise + noise_system = (_generate_noise_system(dimensions_tr=dimensions_tr) * + system_sigma) - # Z score the volume and multiply it by the overall noise - noise = stats.zscore(noise) * noise_dict['overall'] + # Convert temporal noise (in percent) to real numbers + abs_change = noise_dict['temporal_noise'] * mean_signal / 100 - # Set the SNR of the brain - noise = noise + (mask * noise_dict['snr'] * noise_dict['overall']) + # Sum up the noise of the brain + noise = base + (noise_temporal * abs_change) + noise_system return noise diff --git a/examples/utils/fmrisim_example.py b/examples/utils/fmrisim_example.py index cd0ee4036..f436228c9 100644 --- a/examples/utils/fmrisim_example.py +++ b/examples/utils/fmrisim_example.py @@ -36,8 +36,7 @@ [[32, 32, 18], [26, 32, 18], [32, 26, 18], [32, 32, 12]]) coordinates_B = np.array( [[32, 32, 18], [38, 32, 18], [32, 38, 18], [32, 32, 24]]) -signal_magnitude = [1, 0.5, 0.25, -1] - +signal_magnitude = [1, 0.5, 0.25, -1] # In percent signal change # Inputs for generate_stimfunction onsets_A = [10, 30, 50, 70, 90] @@ -65,7 +64,6 @@ signal_magnitude=signal_magnitude, ) - # Visualize the signal that was generated for condition A fig = plt.figure() sim.plot_brain(fig, @@ -119,13 +117,21 @@ # matter likelihood) signal *= mask +# Generate original noise dict for comparison later +orig_noise_dict = sim._noise_dict_update({}) + # Create the noise volumes (using the default parameters noise = sim.generate_noise(dimensions=dimensions, stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, mask=mask, + noise_dict=orig_noise_dict, ) +# Standardize the signal activity to make it percent signal change +mean_act = (mask * orig_noise_dict['max_activity']).sum() / (mask > 0).sum() +signal = signal * mean_act / 100 + # Combine the signal and the noise brain = signal + noise @@ -178,5 +184,6 @@ noise_dict=noise_dict, ) -brain_noise = nibabel.Nifti1Image(noise, affine_matrix) # Create a nifti brain +# Create a nifti brain +brain_noise = nibabel.Nifti1Image(noise, affine_matrix) nibabel.save(brain_noise, 'examples/utils/example2.nii') # Save diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index f853dfc9b..81a9ffab0 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -209,11 +209,12 @@ def test_generate_noise(): stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, mask=mask, - noise_dict={'overall': 0}, + noise_dict={'temporal_noise': 0, 'snr': 10000}, ) - assert np.sum(noise) == 0, "Noise strength could not be manipulated" - assert np.std(noise) == 0, "Noise strength could not be manipulated" + temporal_noise = np.std(noise[mask[:,:,:,0]>0],1).mean() + + assert temporal_noise <= 0.1, "Noise strength could not be manipulated" def test_mask_brain(): @@ -235,7 +236,7 @@ def test_mask_brain(): ) # Mask the volume to be the same shape as a brain - mask = sim.mask_brain(volume)[:,:,:,0] + mask = sim.mask_brain(volume)[:, :, :, 0] brain = volume * (mask > 0) assert np.sum(brain != 0) == np.sum(volume != 0), "Masking did not work" @@ -270,12 +271,12 @@ def test_calc_noise(): dimensions_tr = np.array([10, 10, 10, tr_number]) # Preset the noise dict - nd_orig = {'auto_reg_sigma': 1, - 'drift_sigma': 0.5, - 'overall': 0.1, + nd_orig = {'auto_reg_sigma': 0.6, + 'drift_sigma': 0.4, + 'temporal_noise': 5, 'snr': 30, - 'spatial_sigma': 0.15, - 'system_sigma': 1, + 'max_activity': 1000, + 'spatial_sigma': 4, } # Create the time course for the signal to be generated @@ -297,12 +298,9 @@ def test_calc_noise(): # Calculate the noise nd_calc = sim.calc_noise(noise, mask) - assert abs(nd_calc['overall'] - nd_orig['overall']) < 0.1, 'overall ' \ - 'calculated ' \ - 'incorrectly' - - assert abs(nd_calc['snr'] - nd_orig['snr']) < 10, 'snr calculated ' \ - 'incorrectly' + # How precise are these estimates + precision = abs(nd_calc['temporal_noise'] - nd_orig['temporal_noise']) + assert precision < 1, 'temporal_noise calculated incorrectly' - assert abs(nd_calc['system_sigma'] - nd_orig['system_sigma']) < 1, \ - 'snr calculated incorrectly' + precision = abs(nd_calc['snr'] - nd_orig['snr']) + assert precision < 10, 'snr calculated incorrectly' From 17347640e1fea0beb5d052dd5163b56199cc682d Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 29 Mar 2017 13:15:45 -0400 Subject: [PATCH 2/8] pep8 update --- brainiak/utils/fmrisim.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index df6204556..71c5df42f 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -57,14 +57,12 @@ import logging from itertools import product -import os import nitime.algorithms.autoregressive as ar import math import numpy as np from pkg_resources import resource_stream from scipy import stats import scipy.ndimage as ndimage -import nibabel __all__ = [ "generate_signal", @@ -947,7 +945,7 @@ def calc_noise(volume, # Calculate the fwhm on a subset of volumes - if volume.shape[3]>100: + if volume.shape[3] > 100: # Take only 100 shuffled TRs trs = np.arange(volume.shape[3]) np.random.shuffle(trs) From 9120c62f9468365c86a89416e66149439895c0ff Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 29 Mar 2017 13:30:56 -0400 Subject: [PATCH 3/8] Made _calc_fwhm less complex --- brainiak/utils/fmrisim.py | 183 ++++++++++++++++++-------------------- 1 file changed, 86 insertions(+), 97 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 71c5df42f..5b887183c 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -684,15 +684,15 @@ def _calc_fwhm(volume, voxel_size=[1.0, 1.0, 1.0], ): """ Calculate the FWHM of a volume - Takes in a volume and mask and outputs the FWHM (mm) of each TR of this + Takes in a 3d volume and mask and outputs the FWHM (mm) of this volume for the non-masked voxels Parameters ---------- - volume : multidimensional array - Functional data to have the FWHM measured. Can be 3d or 4d data + volume : 3 dimensional array + Functional data to have the FWHM measured. - mask : 4 dimensional array + mask : 3 dimensional array A mask of the voxels to have the FWHM measured from voxel_size : length 3 list, float @@ -707,97 +707,80 @@ def _calc_fwhm(volume, # What are the dimensions of the volume dimensions = volume.shape - # Identify the number of TRs - if len(dimensions) == 3: - trs = 1 - else: - trs = dimensions[3] - dimensions = dimensions[0:3] - volume_all = volume # Store for later - - # Preset size - fwhm = np.zeros(trs) - # Iterate through the TRs, creating a FWHM for each TR - for tr_counter in list(range(0, trs)): - - # Preset - v_count = 0 - v_sum = 0 - v_sq = 0 - - d_sum = [0.0, 0.0, 0.0] - d_sq = [0.0, 0.0, 0.0] - d_count = [0, 0, 0] - - # If there is more than one TR, just pull out that volume - if trs > 1: - volume = volume_all[:, :, :, tr_counter] - - # Pull out all the voxel coordinates - coordinates = list(product(range(dimensions[0]), - range(dimensions[1]), - range(dimensions[2]))) - - # Find the sum of squared error for the non-masked voxels in the brain - for i in list(range(len(coordinates))): - - # Pull out this coordinate - x, y, z = coordinates[i] - - # Is this within the mask? - if mask[x, y, z, tr_counter] > 0: - - # Find the the volume sum and squared values - v_count += 1 - v_sum += np.abs(volume[x, y, z]) - v_sq += volume[x, y, z] ** 2 - - # Get the volume variance - v_var = (v_sq - ((v_sum ** 2) / v_count)) / (v_count - 1) - - for i in list(range(len(coordinates))): - - # Pull out this coordinate - x, y, z = coordinates[i] - - # Is this within the mask? - if mask[x, y, z, tr_counter] > 0: - # For each xyz dimension calculate the squared - # difference of this voxel and the next - - inrange = (x < dimensions[0] - 1) - inmask = inrange and (mask[x + 1, y, z, tr_counter] > 0) - included = inmask and (~np.isnan(volume[x + 1, y, z])) - if included: - d_sum[0] += volume[x, y, z] - volume[x + 1, y, z] - d_sq[0] += (volume[x, y, z] - volume[x + 1, y, z]) ** 2 - d_count[0] += 1 - - inrange = (y < dimensions[1] - 1) - inmask = inrange and (mask[x, y + 1, z, tr_counter] > 0) - included = inmask and (~np.isnan(volume[x, y + 1, z])) - if included: - d_sum[1] += volume[x, y, z] - volume[x, y + 1, z] - d_sq[1] += (volume[x, y, z] - volume[x, y + 1, z]) ** 2 - d_count[1] += 1 - - inrange = (z < dimensions[2] - 1) - inmask = inrange and (mask[x, y, z + 1, tr_counter] > 0) - included = inmask and (~np.isnan(volume[x, y, z + 1])) - if included: - d_sum[2] += volume[x, y, z] - volume[x, y, z + 1] - d_sq[2] += (volume[x, y, z] - volume[x, y, z + 1]) ** 2 - d_count[2] += 1 - - # Find the variance - d_var = np.divide((d_sq - np.divide(np.power(d_sum, 2), - d_count)), (np.add(d_count, -1))) - - o_var = np.divide(-1, (4 * np.log(1 - (0.5 * d_var / v_var)))) - fwhm3 = np.sqrt(o_var) * 2 * np.sqrt(2 * np.log(2)) - fwhm[tr_counter] = np.prod(np.multiply(fwhm3, voxel_size)) ** (1 / 3) + # Preset + v_count = 0 + v_sum = 0 + v_sq = 0 + + d_sum = [0.0, 0.0, 0.0] + d_sq = [0.0, 0.0, 0.0] + d_count = [0, 0, 0] + + # Pull out all the voxel coordinates + coordinates = list(product(range(dimensions[0]), + range(dimensions[1]), + range(dimensions[2]))) + + # Find the sum of squared error for the non-masked voxels in the brain + for i in list(range(len(coordinates))): + + # Pull out this coordinate + x, y, z = coordinates[i] + + # Is this within the mask? + if mask[x, y, z] > 0: + + # Find the the volume sum and squared values + v_count += 1 + v_sum += np.abs(volume[x, y, z]) + v_sq += volume[x, y, z] ** 2 + + # Get the volume variance + v_var = (v_sq - ((v_sum ** 2) / v_count)) / (v_count - 1) + + for i in list(range(len(coordinates))): + + # Pull out this coordinate + x, y, z = coordinates[i] + + # Is this within the mask? + if mask[x, y, z] > 0: + # For each xyz dimension calculate the squared + # difference of this voxel and the next + + inrange = (x < dimensions[0] - 1) + inmask = inrange and (mask[x + 1, y, z] > 0) + included = inmask and (~np.isnan(volume[x + 1, y, z])) + if included: + d_sum[0] += volume[x, y, z] - volume[x + 1, y, z] + d_sq[0] += (volume[x, y, z] - volume[x + 1, y, z]) ** 2 + d_count[0] += 1 + + inrange = (y < dimensions[1] - 1) + inmask = inrange and (mask[x, y + 1, z] > 0) + included = inmask and (~np.isnan(volume[x, y + 1, z])) + if included: + d_sum[1] += volume[x, y, z] - volume[x, y + 1, z] + d_sq[1] += (volume[x, y, z] - volume[x, y + 1, z]) ** 2 + d_count[1] += 1 + + inrange = (z < dimensions[2] - 1) + inmask = inrange and (mask[x, y, z + 1] > 0) + included = inmask and (~np.isnan(volume[x, y, z + 1])) + if included: + d_sum[2] += volume[x, y, z] - volume[x, y, z + 1] + d_sq[2] += (volume[x, y, z] - volume[x, y, z + 1]) ** 2 + d_count[2] += 1 + + # Find the variance + d_var = np.divide((d_sq - np.divide(np.power(d_sum, 2), + d_count)), (np.add(d_count, -1))) + + o_var = np.divide(-1, (4 * np.log(1 - (0.5 * d_var / v_var)))) + fwhm3 = np.sqrt(o_var) * 2 * np.sqrt(2 * np.log(2)) + fwhm = np.prod(np.multiply(fwhm3, voxel_size)) ** (1 / 3) return fwhm @@ -953,10 +936,16 @@ def calc_noise(volume, else: trs = list(range(0, volume.shape[3])) - noise_dict['fwhm'] = _calc_fwhm(volume[:, :, :, trs], - mask[:, :, :, trs], - noise_dict['voxel_size'], - ).mean() + # Go through the trs and pull out the fwhm + fwhm = [0] * len(trs) + for tr in list(range(0, len(trs))): + fwhm[tr] = _calc_fwhm(volume[:, :, :, trs[tr]], + mask[:, :, :, trs[tr]], + noise_dict['voxel_size'], + ) + + # Keep only the mean + noise_dict['fwhm'] = np.mean(fwhm) # Calculate the autoregressive and drift noise auto_reg_sigma, drift_sigma = _calc_autoregression(volume, mask) From 9d24c09fe69d9763cc7ff85f531757d2faf002be Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Wed, 29 Mar 2017 16:19:33 -0400 Subject: [PATCH 4/8] Update for edits --- brainiak/utils/fmrisim.py | 60 +++++++++++++++++++++++++++---------- tests/utils/test_fmrisim.py | 2 +- 2 files changed, 45 insertions(+), 17 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 5b887183c..62748fa70 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -52,7 +52,7 @@ Authors: Cameron Ellis (Princeton) 2016-2017 - Chris Boldassano (Princeton) 2016-2017 + Chris Baldassano (Princeton) 2016-2017 """ import logging @@ -750,25 +750,25 @@ def _calc_fwhm(volume, # For each xyz dimension calculate the squared # difference of this voxel and the next - inrange = (x < dimensions[0] - 1) - inmask = inrange and (mask[x + 1, y, z] > 0) - included = inmask and (~np.isnan(volume[x + 1, y, z])) + in_range = (x < dimensions[0] - 1) + in_mask = in_range and (mask[x + 1, y, z] > 0) + included = in_mask and (~np.isnan(volume[x + 1, y, z])) if included: d_sum[0] += volume[x, y, z] - volume[x + 1, y, z] d_sq[0] += (volume[x, y, z] - volume[x + 1, y, z]) ** 2 d_count[0] += 1 - inrange = (y < dimensions[1] - 1) - inmask = inrange and (mask[x, y + 1, z] > 0) - included = inmask and (~np.isnan(volume[x, y + 1, z])) + in_range = (y < dimensions[1] - 1) + in_mask = in_range and (mask[x, y + 1, z] > 0) + included = in_mask and (~np.isnan(volume[x, y + 1, z])) if included: d_sum[1] += volume[x, y, z] - volume[x, y + 1, z] d_sq[1] += (volume[x, y, z] - volume[x, y + 1, z]) ** 2 d_count[1] += 1 - inrange = (z < dimensions[2] - 1) - inmask = inrange and (mask[x, y, z + 1] > 0) - included = inmask and (~np.isnan(volume[x, y, z + 1])) + in_range = (z < dimensions[2] - 1) + in_mask = in_range and (mask[x, y, z + 1] > 0) + included = in_mask and (~np.isnan(volume[x, y, z + 1])) if included: d_sum[2] += volume[x, y, z] - volume[x, y, z + 1] d_sq[2] += (volume[x, y, z] - volume[x, y, z + 1]) ** 2 @@ -887,7 +887,7 @@ def _calc_autoregression(volume, def calc_noise(volume, mask=None, - noise_dict={'voxel_size': [1.0, 1.0, 1.0]}, + noise_dict=None, ): """ Calculates the noise properties of the volume supplied. This estimates what noise properties the volume has. For instance it @@ -919,6 +919,12 @@ def calc_noise(volume, if mask is None: mask = np.ones(volume.shape) + # Update noise dict + if noise_dict is None: + noise_dict = {'voxel_size': [1.0, 1.0, 1.0]} + elif 'voxel_size' not in noise_dict: + noise_dict['voxel_size'] = [1.0, 1.0, 1.0] + # Since you are deriving the 'true' values then you want your noise to # be set to that level @@ -1233,13 +1239,35 @@ def _generate_noise_spatial(dimensions, if len(dimensions) == 4: return - # Convert from fwhm to sigma (relationship discovered empirical, only an - # approximation up to sigma=0-5 which corresponds to fwhm=0-8, relies - # on an assumption of brain size). - def logcurve(x, a, b, c): + def logfunc(x, a, b, c): + """Solve for y given x for log function. + + Parameters + ---------- + x : float + x value of log function + + a : float + x shift of function + + b : float + rate of change + + c : float + y shift of function + + Returns + ---------- + + float + y value of log function + """ return (np.log(x + a) / np.log(b)) + c - spatial_sigma = logcurve(fwhm, -0.36778719, 2.10601011, 2.15439247) + # Convert from fwhm to sigma (relationship discovered empirical, only an + # approximation up to sigma = 0 -> 5 which corresponds to fwhm = 0 -> 8, + # relies on an assumption of brain size). + spatial_sigma = logfunc(fwhm, -0.36778719, 2.10601011, 2.15439247) # Set up the input to the fast fourier transform def fftIndgen(n): diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 81a9ffab0..9df305a29 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -276,7 +276,7 @@ def test_calc_noise(): 'temporal_noise': 5, 'snr': 30, 'max_activity': 1000, - 'spatial_sigma': 4, + 'fwhm': 4, } # Create the time course for the signal to be generated From aab55be69116bae58fe0b1b64b86229d005aea07 Mon Sep 17 00:00:00 2001 From: Mingbo Cai Date: Wed, 22 Mar 2017 16:00:59 -0400 Subject: [PATCH 5/8] [reprsimil] Fix BRSA bugs (#197) * a Changes to be committed: new file: brainiak/brsa/test.txt * deleted: brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * deleted: brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * deleted: brsa/test.txt * deleted: brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * a Changes to be committed: new file: brainiak/brsa/test.txt * deleted: brsa/test.txt * deleted: brsa/test.txt * deleted: brsa/test.txt * deleted: brsa/test.txt * Changed all usage of scipy.linalg.block_diag in accordance to its change in last November. Fixed the bug in posterior estimation of betas since the review of the major pull request might take longer. * reduce testing time for brsa * renamed BRSA example * rename brsa example --- brainiak/reprsimil/brsa.py | 26 ++++++------------- ...ample.ipynb => bayesian_rsa_example.ipynb} | 16 +++++------- tests/reprsimil/test_brsa.py | 19 +++++++------- 3 files changed, 24 insertions(+), 37 deletions(-) rename examples/reprsimil/{brsa_representational_similarity_estimate_example.ipynb => bayesian_rsa_example.ipynb} (97%) diff --git a/brainiak/reprsimil/brsa.py b/brainiak/reprsimil/brsa.py index 11eec9fe0..59a466ced 100755 --- a/brainiak/reprsimil/brsa.py +++ b/brainiak/reprsimil/brsa.py @@ -418,10 +418,8 @@ def _prepare_DF(self, n_T, scan_onsets=None): And we denote A = I - rho1*D + rho1**2 * F""" if scan_onsets is None: # assume that all data are acquired within the same scan. - D = np.diag(np.ones(n_T - 1), -1) + np.diag(np.ones(n_T - 1), 1) - F = np.eye(n_T) - F[0, 0] = 0 - F[n_T - 1, n_T - 1] = 0 + D = self._D_gen(n_T) + F = self._F_gen(n_T) n_run = 1 run_TRs = np.array([n_T]) else: @@ -439,12 +437,8 @@ def _prepare_DF(self, n_T, scan_onsets=None): D_ele = map(self._D_gen, run_TRs) F_ele = map(self._F_gen, run_TRs) - D = [] - for d_ele in D_ele: - D = scipy.linalg.block_diag(D, d_ele) - F = [] - for f_ele in F_ele: - F = scipy.linalg.block_diag(F, f_ele) + D = scipy.linalg.block_diag(*D_ele) + F = scipy.linalg.block_diag(*F_ele) # D and F above are templates for constructing # the inverse of temporal covariance matrix of noise return D, F, run_TRs, n_run @@ -480,9 +474,7 @@ def _prepare_data_XYX0(self, X, Y, X0, D, F, run_TRs, no_DC=False): It will only take effect if X0 is not None. """ - X_base = [] - for r_l in run_TRs: - X_base = scipy.linalg.block_diag(X_base, np.ones(r_l)[:, None]) + X_base = scipy.linalg.block_diag(*map(np.ones, run_TRs)).T res = np.linalg.lstsq(X_base, X) if np.any(np.isclose(res[1], 0)): raise ValueError('Your design matrix appears to have ' @@ -860,7 +852,7 @@ def _fit_RSA_UV(self, X, Y, X0, X0TAX0, est_SNR_AR1_UV**2, n_V, n_T, n_run, rank, n_base) est_sigma_AR1_UV = sigma2**0.5 - est_beta_AR1_UV = est_sigma_AR1_UV * est_SNR_AR1_UV**2 \ + est_beta_AR1_UV = est_SNR_AR1_UV**2 \ * np.dot(estU_chlsk_l_AR1_UV, YTAcorrXL_LAMBDA.T) est_beta0_AR1_UV = np.einsum( 'ijk,ki->ji', X0TAX0_i, @@ -1088,8 +1080,7 @@ def _fit_diagV_noGP( = self._calc_LL(rho1, LTXTAcorrXL, LTXTAcorrY, YTAcorrY, X0TAX0, current_SNR2, n_V, n_T, n_run, rank, n_base) - betas = current_sigma2**0.5 * current_SNR2 \ - * np.dot(L, YTAcorrXL_LAMBDA.T) + betas = current_SNR2 * np.dot(L, YTAcorrXL_LAMBDA.T) residuals = Y - np.dot(X, betas) u, s, v = np.linalg.svd(residuals) X0 = u[:, :self.n_nureg] @@ -1213,8 +1204,7 @@ def _fit_diagV_GP( = self._calc_LL(rho1, LTXTAcorrXL, LTXTAcorrY, YTAcorrY, X0TAX0, current_SNR2, n_V, n_T, n_run, rank, n_base) - betas = current_sigma2**0.5 * current_SNR2 \ - * np.dot(L, YTAcorrXL_LAMBDA.T) + betas = current_SNR2 * np.dot(L, YTAcorrXL_LAMBDA.T) residuals = Y - np.dot(X, betas) u, s, v = np.linalg.svd(residuals) X0 = u[:, :self.n_nureg] diff --git a/examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb b/examples/reprsimil/bayesian_rsa_example.ipynb similarity index 97% rename from examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb rename to examples/reprsimil/bayesian_rsa_example.ipynb index f143c45e6..c6c7620ac 100644 --- a/examples/reprsimil/brsa_representational_similarity_estimate_example.ipynb +++ b/examples/reprsimil/bayesian_rsa_example.ipynb @@ -139,7 +139,7 @@ "outputs": [], "source": [ "noise_bot = 0.5\n", - "noise_top = 1.5\n", + "noise_top = 5.0\n", "noise_level = np.random.rand(n_V) * \\\n", " (noise_top - noise_bot) + noise_bot\n", "# The standard deviation of the noise is in the range of [noise_bot, noise_top]\n", @@ -413,7 +413,6 @@ "metadata": {}, "source": [ "# Fit Bayesian RSA to our simulated data\n", - "### The data should be z-scored along time dimension. Otherwise voxels of high variance in noise can dominate the estimation of spatially shared noise component.\n", "\n", "The nuisance regressors in typical fMRI analysis (such as head motion signal) are replaced by principal components estimated from residuals after subtracting task-related response. n_nureg tells the model how many principal components to keep from the residual as nuisance regressors, in order to account for spatial correlation in noise. \n", "If you prefer not using this approach based on principal components of residuals, you can set auto_nuisance=False, and optionally provide your own nuisance regressors as nuisance argument to BRSA.fit()" @@ -427,9 +426,6 @@ }, "outputs": [], "source": [ - "Y_z = scipy.stats.zscore(Y,axis=0)\n", - "Y_std = np.std(Y,axis=0)\n", - "# Z-scoreing the data\n", "\n", "brsa = BRSA(GP_space=True, GP_inten=True,\n", " n_nureg=10)\n", @@ -438,7 +434,7 @@ "# over both space and intensity.\n", "\n", "\n", - "brsa.fit(X=Y_z, design=design.design_task,\n", + "brsa.fit(X=Y, design=design.design_task,\n", " coords=coords_flat, inten=inten, scan_onsets=scan_onsets)\n", "# The data to fit should be given to the argument X.\n", "# Design matrix goes to design. And so on.\n" @@ -500,7 +496,7 @@ "source": [ "regressor = np.insert(design.design_task,\n", " 0, 1, axis=1)\n", - "betas_point = np.linalg.lstsq(regressor, Y_z)[0]\n", + "betas_point = np.linalg.lstsq(regressor, Y)[0]\n", "point_corr = np.corrcoef(betas_point[1:, :])\n", "point_cov = np.cov(betas_point[1:, :])\n", "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n", @@ -634,14 +630,14 @@ }, "outputs": [], "source": [ - "plt.scatter(betas_simulated / Y_std, brsa.beta_)\n", + "plt.scatter(betas_simulated, brsa.beta_)\n", "plt.xlabel('true betas (response amplitudes)')\n", "plt.ylabel('recovered betas by Bayesian RSA')\n", "ax = plt.gca()\n", "plt.show()\n", "\n", "\n", - "plt.scatter(betas_simulated / Y_std, betas_point[1:, :])\n", + "plt.scatter(betas_simulated, betas_point[1:, :])\n", "plt.xlabel('true betas (response amplitudes)')\n", "plt.ylabel('recovered betas by simple regression')\n", "ax = plt.gca()\n", @@ -664,7 +660,7 @@ }, "outputs": [], "source": [ - "u, s, v = np.linalg.svd(Y_z - signal / Y_std)\n", + "u, s, v = np.linalg.svd(Y - signal)\n", "plt.plot(s)\n", "plt.xlabel('principal component')\n", "plt.ylabel('singular value of simulated noise')\n", diff --git a/tests/reprsimil/test_brsa.py b/tests/reprsimil/test_brsa.py index ad25e2d53..1d5a9a46e 100755 --- a/tests/reprsimil/test_brsa.py +++ b/tests/reprsimil/test_brsa.py @@ -42,17 +42,18 @@ def test_fit(): design = utils.ReadDesign(fname=file_path) - # concatenate it by 4 times, mimicking 4 runs of itenditcal timing - design.design_task = np.tile(design.design_task[:,:-1],[4,1]) - design.n_TR = design.n_TR * 4 + # concatenate it by 2 times, mimicking 2 runs of itenditcal timing + n_run = 2 + design.design_task = np.tile(design.design_task[:,:-1],[n_run,1]) + design.n_TR = design.n_TR * n_run # start simulating some data - n_V = 200 + n_V = 50 n_C = np.size(design.design_task,axis=1) n_T = design.n_TR noise_bot = 0.5 - noise_top = 1.5 + noise_top = 5.0 noise_level = np.random.rand(n_V)*(noise_top-noise_bot)+noise_bot # noise level is random. @@ -110,14 +111,14 @@ def test_fit(): signal = np.dot(design.design_task,betas_simulated) # Adding noise to signal as data - Y = signal + noise + Y = signal + noise + inten scan_onsets = np.linspace(0,design.n_TR,num=5) # Test fitting with GP prior. - brsa = BRSA(GP_space=True,GP_inten=True,verbose=False,n_iter = 200,auto_nuisance=False) + brsa = BRSA(GP_space=True,GP_inten=True,verbose=False,n_iter = 5,auto_nuisance=False) # We also test that it can detect baseline regressor included in the design matrix for task conditions wrong_design = np.insert(design.design_task, 0, 1, axis=1) @@ -167,7 +168,7 @@ def test_fit(): 'the BRSA object should not have parameters of GP if GP is not requested.' # GP parameters are not set if not requested assert brsa.beta0_.shape[0] == n_nureg, 'Shape of beta0 incorrect' - p = scipy.stats.pearsonr(brsa.beta0_[0,:],np.mean(noise,axis=0))[1] + p = scipy.stats.pearsonr(brsa.beta0_[0,:],inten)[1] assert p < 0.05, 'recovered beta0 does not correlate with the baseline of voxels.' # Test fitting with GP over just spatial coordinates. @@ -209,7 +210,7 @@ def test_gradient(): design.n_TR = design.n_TR * n_run # start simulating some data - n_V = 200 + n_V = 30 n_C = np.size(design.design_task,axis=1) n_T = design.n_TR From 34a8b022c612dc1f12299bd33ab2e5e6d0dee426 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mihai=20Capot=C4=83?= Date: Thu, 23 Mar 2017 10:30:16 -0700 Subject: [PATCH 6/8] [dev] Fix LDFLAGS for LLVM on MacOS --- .travis.yml | 3 ++- README.rst | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 934a533e2..7cb828e67 100644 --- a/.travis.yml +++ b/.travis.yml @@ -27,7 +27,8 @@ matrix: env: &macos_env - CC=/usr/local/opt/llvm/bin/clang - CXX=/usr/local/opt/llvm/bin/clang++ - - LDFLAGS="-L/usr/local/opt/llvm/lib $LDFLAGS" + - LDFLAGS="-L/usr/local/opt/llvm/lib + -Wl,-rpath,/usr/local/opt/llvm/lib $LDFLAGS" - CPPFLAGS="-I/usr/local/opt/llvm/include $CPPFLAGS" before_install: &macos_before_install - brew install llvm mpich python3 diff --git a/README.rst b/README.rst index 75c511b12..751800f87 100644 --- a/README.rst +++ b/README.rst @@ -77,7 +77,8 @@ shell, e.g., ``bash -l``:: In addition, you also need to specify the directories that the newly installed `clang` will seek for compiling and linking:: - export LDFLAGS="-L/usr/local/opt/llvm/lib $LDFLAGS" + export LDFLAGS="-L/usr/local/opt/llvm/lib "\ + "-Wl,-rpath,/usr/local/opt/llvm/lib $LDFLAGS" export CPPFLAGS="-I/usr/local/opt/llvm/include $CPPFLAGS" Install updated versions of the following Python packages:: From 5c73f97ad3f5a8a642ee3970431e1ba5f836a45d Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Fri, 12 May 2017 17:05:25 -0400 Subject: [PATCH 7/8] _calc_snr turned into _calc_sfnr with the necessary changes associated with it. Improves accuracy of noise modelling --- brainiak/utils/fmrisim.py | 82 ++++++++++++++++++++++++------------- tests/utils/test_fmrisim.py | 8 ++-- 2 files changed, 57 insertions(+), 33 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 634f854f6..0e7b37ef2 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -785,12 +785,13 @@ def _calc_fwhm(volume, return fwhm -def _calc_snr(volume, +def _calc_sfnr(volume, mask, ): - """ Calculate the SNR of a volume + """ Calculate the SFNR of a volume This takes the middle of the volume and averages the signal within the - brain and compares to the std in the voxels outside the brain. + brain and compares to the temporal standard deviation in the voxels + outside the brain. Parameters ---------- @@ -803,9 +804,15 @@ def _calc_snr(volume, Returns ------- - float - The snr of the volume (how much greater in activity the brain is from - the background) + float
 + The sd of the temporal variability of brain voxels.

 + + float
 + The sfnr of the volume (mean brain activity divided by
 temporal + variability in the average non brain voxels)

 + + float
 + What is the max activity measured here, a point of reference for masking """ @@ -821,14 +828,21 @@ def _calc_snr(volume, # then you would have a lower mean signal by default) mean_signal = (slice_volume[slice_mask > 0]).mean() - # What is the noise in the non brain - background_noise = float(slice_volume[slice_mask == 0].std()) + # What are the brain and non-brain voxels + brain_voxels = volume[mask[:, :, :, 0] > 0] + mask_voxels = volume[mask[:, :, :, 0] == 0] + + # What is the noise in the masked voxels + mask_noise = np.std(mask_voxels, 1).mean() - # What is the max activation of this volume - max_activity = volume.max() + # Assume the mask noise is a combination of random + drift. Thus average + # all masked voxels and find the variation over time, this is the + # variance due to drift. + drift_noise = np.mean(mask_voxels, 0).std() - # What is the brain voxel by time activity. - brain_voxels = volume[mask[:, :, :, 0] > 0] + # Subtract the drift variance from total mask noise to find the system + # noise + background_noise = np.sqrt(mask_noise ** 2 - drift_noise ** 2) # Find the noise to brain voxels temporal_noise = np.std(brain_voxels, 1).mean() @@ -836,20 +850,23 @@ def _calc_snr(volume, # Convert temporal noise into percent signal change temporal_noise = temporal_noise / mean_signal * 100 - # Calculate snr - snr = mean_signal / background_noise + # Calculate sfnr + sfnr = mean_signal / background_noise # Convert from memmap - snr = float(snr) + sfnr = float(sfnr) - return temporal_noise, snr, max_activity + # What is the max activation of this volume + max_activity = volume.max() + return temporal_noise, sfnr, max_activity -def _calc_autoregression(volume, + +def _calc_temporal_noise(volume, mask, auto_reg_order=1, ): - """ Calculate the autoregressive sigma of the data. + """ Calculate the mix of autoregressive and drift noise. Parameters ---------- @@ -928,9 +945,9 @@ def calc_noise(volume, # Since you are deriving the 'true' values then you want your noise to # be set to that level - # Calculate the temporal_noise noise and SNR of the volume - noise_dict['temporal_noise'], noise_dict['snr'], noise_dict[ - 'max_activity'] = _calc_snr(volume, mask) + # Calculate the temporal_noise noise and SFNR of the volume + noise_dict['temporal_noise'], noise_dict['sfnr'], noise_dict[ + 'max_activity'] = _calc_sfnr(volume, mask) # Calculate the fwhm on a subset of volumes @@ -954,7 +971,7 @@ def calc_noise(volume, noise_dict['fwhm'] = np.mean(fwhm) # Calculate the autoregressive and drift noise - auto_reg_sigma, drift_sigma = _calc_autoregression(volume, mask) + auto_reg_sigma, drift_sigma = _calc_temporal_noise(volume, mask) # Calibrate for how sigma is originally calculated auto_reg_sigma = auto_reg_sigma / noise_dict['temporal_noise'] @@ -1557,8 +1574,8 @@ def _noise_dict_update(noise_dict): noise_dict['auto_reg_sigma'] = 0.45 if 'physiological_sigma' not in noise_dict: noise_dict['physiological_sigma'] = 0.1 - if 'snr' not in noise_dict: - noise_dict['snr'] = 30 + if 'sfnr' not in noise_dict: + noise_dict['sfnr'] = 30 if 'max_activity' not in noise_dict: noise_dict['max_activity'] = 1000 if 'voxel_size' not in noise_dict: @@ -1645,7 +1662,7 @@ def generate_noise(dimensions, # Create the base (this inverts the process to make the mask) base = mask * noise_dict['max_activity'] - # Set the amount of background based on the snr value + # Set the amount of background based on the SFNR value # What are the midpoints to be extracted mid_x_idx = int(np.ceil(base.shape[0] / 2)) @@ -1658,12 +1675,16 @@ def generate_noise(dimensions, # What is the mean signal of the non masked voxels in this slice? mean_signal = (slice_volume[slice_mask > 0]).mean() + # Set up the machine noise + noise_system = _generate_noise_system(dimensions_tr=dimensions_tr) + # What is the standard deviation of the background activity - system_sigma = mean_signal / noise_dict['snr'] + # (N.B. You need to subtract the mean system noise from this number + # since system_sigma * noise_system.mean() will later be added to the base) + system_sigma = mean_signal / (noise_dict['sfnr'] - noise_system.mean()) - # Set up the machine noise - noise_system = (_generate_noise_system(dimensions_tr=dimensions_tr) * - system_sigma) + # Increase the size of the system noise based on the SFNR + noise_system *= system_sigma # Convert temporal noise (in percent) to real numbers abs_change = noise_dict['temporal_noise'] * mean_signal / 100 @@ -1671,6 +1692,9 @@ def generate_noise(dimensions, # Sum up the noise of the brain noise = base + (noise_temporal * abs_change) + noise_system + # Reject negative values (only happens outside of the brain + noise[noise < 0] = 0 + return noise diff --git a/tests/utils/test_fmrisim.py b/tests/utils/test_fmrisim.py index 9df305a29..a8464c229 100644 --- a/tests/utils/test_fmrisim.py +++ b/tests/utils/test_fmrisim.py @@ -209,7 +209,7 @@ def test_generate_noise(): stimfunction_tr=stimfunction_tr, tr_duration=tr_duration, mask=mask, - noise_dict={'temporal_noise': 0, 'snr': 10000}, + noise_dict={'temporal_noise': 0, 'sfnr': 10000}, ) temporal_noise = np.std(noise[mask[:,:,:,0]>0],1).mean() @@ -274,7 +274,7 @@ def test_calc_noise(): nd_orig = {'auto_reg_sigma': 0.6, 'drift_sigma': 0.4, 'temporal_noise': 5, - 'snr': 30, + 'sfnr': 30, 'max_activity': 1000, 'fwhm': 4, } @@ -302,5 +302,5 @@ def test_calc_noise(): precision = abs(nd_calc['temporal_noise'] - nd_orig['temporal_noise']) assert precision < 1, 'temporal_noise calculated incorrectly' - precision = abs(nd_calc['snr'] - nd_orig['snr']) - assert precision < 10, 'snr calculated incorrectly' + precision = abs(nd_calc['sfnr'] - nd_orig['sfnr']) + assert precision < 5, 'sfnr calculated incorrectly' From 3b282260b2fa236d0ee63b7f75751eeeb519fea5 Mon Sep 17 00:00:00 2001 From: CameronTEllis Date: Fri, 12 May 2017 17:23:02 -0400 Subject: [PATCH 8/8] Pep8 errors corrected --- brainiak/utils/fmrisim.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/brainiak/utils/fmrisim.py b/brainiak/utils/fmrisim.py index 0e7b37ef2..c93cbcee7 100644 --- a/brainiak/utils/fmrisim.py +++ b/brainiak/utils/fmrisim.py @@ -786,8 +786,8 @@ def _calc_fwhm(volume, def _calc_sfnr(volume, - mask, - ): + mask, + ): """ Calculate the SFNR of a volume This takes the middle of the volume and averages the signal within the brain and compares to the temporal standard deviation in the voxels @@ -1684,7 +1684,7 @@ def generate_noise(dimensions, system_sigma = mean_signal / (noise_dict['sfnr'] - noise_system.mean()) # Increase the size of the system noise based on the SFNR - noise_system *= system_sigma + noise_system *= system_sigma # Convert temporal noise (in percent) to real numbers abs_change = noise_dict['temporal_noise'] * mean_signal / 100