Merge 1b59ab6 into 296d2a8

code/Makefile

@@ -13,3 +13,6 @@ test:
 coverage:
 	nosetests utils --with-coverage --cover-package=utils
 	nosetests model --with-coverage --cover-package=model
+
+clean:
+	rm -f *.{nii.gz,txt}

code/model/construct_design_matrix.py

@@ -0,0 +1,55 @@
+from __future__ import division, print_function, absolute_import
+import numpy as np
+
+def read_files(filename):
+	"""read from behav.txt and covert data to array format
+
+	Parameter:
+	---------
+	filename: behavdata.txt file
+
+	Return:
+	------
+	val_arr: 
+		values read from behav.txt and sorted as an array
+	"""
+	with open(filename,"r") as infile:
+		lines = infile.readlines()
+	lines = lines[1:]
+	val = [line.split() for line in lines]
+	val_arr = np.array(val, dtype=float)
+	return val_arr
+
+def construct_mat(array, includeED=True):
+	"""Construct the design matrix using the array return from the
+	read_files function.
+
+	Parameter:
+	----------
+	array: 2-D array
+		return by read_files
+	includeED: boolean (True or False)
+		True (default): include Euclidean distance as the 3rd regressor (4th column of design_matrix)
+		False: not include. Only two regressors (gain and loss) are included.
+
+	Return:
+	------
+	design_matrix: 2-D array 
+		The first column is always intercept, and the last column is the response (0 and 1)
+		The column names are (intercept, gain, loss, distance(if includeED=True), response) 
+	"""
+	array = array[array[:,-2]!=-1,:]
+	if includeED:
+		mat = array[:,[1,2,4,5]]
+		diagOfGambleMat = np.array([np.sum(mat[:,2]==1), np.sum(mat[:,2]==3)])
+		distance = np.sqrt(np.sum((mat[:,:2]-diagOfGambleMat)**2, axis=1))
+		design_matrix = np.ones((mat.shape[0],mat.shape[1]+1))
+		design_matrix[:,1:] = mat
+		design_matrix[:,3] = distance
+	else:
+		mat = array[:,[1,2,5]]
+		design_matrix = np.ones((mat.shape[0],mat.shape[1]+1))
+		design_matrix[:,1:] = mat
+	return design_matrix
+
+

code/model/corr_per_voxel.py

@@ -0,0 +1,38 @@
+from __future__ import division, print_function, absolute_import
+import nibabel as nib
+import matplotlib.pyplot as plt
+import numpy as np
+from stimuli import events2neural
+
+def corr(nii_file, cond_file):
+	"""Find the correlations between time course and each voxel
+	
+	Parameters:
+	----------
+	nii_file: bold.nii file 
+	cond_file: condition file
+
+	Return:
+	-------
+	correlations: an array
+	"""
+	img = nib.load(nii_file)
+	n_trs = img.shape[-1]
+	TR = 2 #The TR (time between scans) is 2 seconds
+	time_course = events2neural(cond_file, 2, n_trs)
+	# Call the events2neural function to generate the on-off values for each volume
+	data = img.get_data() 
+	# Using slicing, drop the first 4 volumes, and the first 4 on-off values
+	data = data[..., 4:]
+	time_course = time_course[4:]
+	n_voxels = np.prod(data.shape[:-1]) # Calculate the number of voxels (number of elements in one volume)
+	data_2d = np.reshape(data, (n_voxels, data.shape[-1])) # Reshape 4D array to 2D array n_voxels by n_volumes
+	correlations_1d = np.zeros((n_voxels,)) # Make a 1D array of size (n_voxels,) to hold the correlation values
+	for i in range(n_voxels): # Loop over voxels filling in correlation at this voxel
+		correlations_1d[i] = np.corrcoef(time_course, data_2d[i, :])[0, 1]
+	correlations = np.reshape(correlations_1d, data.shape[:-1]) # Reshape the correlations array back to 3D
+	plt.imshow(correlations[:, :, 14]) # Plot the middle slice of the third axis from the correlations array
+	return correlations #get correlations of two value
+
+
+

code/model/diagnostics.py

@@ -0,0 +1,98 @@
+""" Diagnostics.py
+
+A collection of utility functions for diagnostics on FMRI data
+
+See test_* functions in this directory for nose tests
+"""
+from __future__ import division, print_function, absolute_import
+import numpy as np
+
+def vol_std(data):
+    """ Return standard deviation across voxels for 4D array `data`
+
+    Parameters
+    ----------
+    data : 4D array
+        4D array from FMRI run with last axis indexing volumes.  Call the shape
+        of this array (M, N, P, T) where T is the number of volumes.
+
+    Returns
+    -------
+    std_values : array shape (T,)
+        One dimensonal array where ``std_values[i]`` gives the standard
+        deviation of all voxels contained in ``data[..., i]``.
+    """
+    data2d = np.reshape(data, (np.prod(data.shape[:-1]),data.shape[-1]))
+    return np.std(data2d, axis=0)
+
+def iqr_outliers(arr_1d, iqr_scale=1.5):
+    """ Return indices of outliers identified by interquartile range
+
+    Parameters
+    ----------
+    arr_1d : 1D array
+        One-dimensional numpy array, from which we will identify outlier
+        values.
+    iqr_scale : float, optional
+        Scaling for IQR to set low and high thresholds.  Low threshold is given
+        by 25th centile value minus ``iqr_scale * IQR``, and high threshold id
+        given by 75 centile value plus ``iqr_scale * IQR``.
+
+    Returns
+    -------
+    outlier_indices : array
+        Array containing indices in `arr_1d` that contain outlier values.
+    lo_hi_thresh : tuple
+        Tuple containing 2 values (low threshold, high thresold) as described
+        above.
+    """
+    # Hint : np.lookfor('centile')
+    # Hint : np.lookfor('nonzero')
+    IQR = np.percentile(arr_1d, 75) - np.percentile(arr_1d, 25)
+    low_threshold = np.percentile(arr_1d, 25) - iqr_scale * IQR
+    high_threshold = np.percentile(arr_1d, 75) + iqr_scale * IQR
+    outlier_indices=np.nonzero((arr_1d < low_threshold) + (arr_1d > high_threshold))[0]
+    return (outlier_indices, (low_threshold, high_threshold))
+
+def vol_rms_diff(arr_4d):
+    """ Return root mean square of differences between sequential volumes
+
+    Parameters
+    ----------
+    data : 4D array
+        4D array from FMRI run with last axis indexing volumes.  Call the shape
+        of this array (M, N, P, T) where T is the number of volumes.
+
+    Returns
+    -------
+    rms_values : array shape (T-1,)
+        One dimensonal array where ``rms_values[i]`` gives the square root of
+        the mean (across voxels) of the squared difference between volume i and
+        volume i + 1.
+    """
+    diff = np.subtract(arr_4d[...,1:], arr_4d[...,:-1])
+    diff2d = np.reshape(diff, (np.prod(diff.shape[:-1]),diff.shape[-1]))
+    rms_values = np.sqrt(np.mean(diff2d**2, axis=0))
+    return rms_values
+
+def extend_diff_outliers(diff_indices):
+    """ Extend difference-based outlier indices `diff_indices` by pairing
+
+    Parameters
+    ----------
+    diff_indices : array
+        Array of indices of differences that have been detected as outliers.  A
+        difference index of ``i`` refers to the difference between volume ``i``
+        and volume ``i + 1``.
+
+    Returns
+    -------
+    extended_indices : array
+        Array where each index ``j`` in `diff_indices has been replaced by two
+        indices, ``j`` and ``j+1``, unless ``j+1`` is present in
+        ``diff_indices``.  For example, if the input was ``[3, 7, 8, 12, 20]``,
+        ``[3, 4, 7, 8, 9, 12, 13, 20, 21]``.
+    """
+    return np.unique(np.append(diff_indices, diff_indices+1))
+
+