Merge pull request #15 from mingujo/master

basic setups & summary stats

.gitignore

@@ -56,5 +56,9 @@ coverage.xml
 doc/build/
 doc/source/_static/test.png
 
+#data
+data/ds005
+data/ds005_raw_checksums.txt
+
 \#*
 *~

code/utils/correlation.py

@@ -0,0 +1,32 @@
+import numpy as np
+
+
+"""get respcat only"""
+respcat=behav[:,5]
+"""-1 if respcat is 0"""
+np.place(respcat, respcat==0, -1)
+
+"""get respnum only"""
+respnum=behav[:,4]
+
+"""apply the weight"""
+weighted_respcat=np.multiply(respcat, respnum)
+
+"""get loss only"""
+loss=behav[:,2]
+
+"""get gain only"""
+gain=behav[:,1]
+
+"""check the correlation between loss and weighted_respcat"""
+
+print("subject 1")
+print(np.corrcoef(loss, weighted_respcat))
+
+"""check the correlation between gain and weighted_respcat"""
+print(np.corrcoef(gain, weighted_respcat))
+
+"""Subject 1 is more responsive to gain"""
+
+
+

code/utils/diagnostics.py

@@ -0,0 +1,116 @@
+""" Diagnostics.py
+
+A collection of utility functions for diagnostics on FMRI data
+
+See test_* functions in this directory for nose tests
+"""
+
+import numpy as np
+import pdb
+
+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]``.
+    """
+    std_values = [np.std(np.ravel(data[...,i])) for i in range(0, np.shape(data)[3])]
+    return np.array(std_values)
+
+def iqr_outliers(arr_1d, iqr_scale=1):
+    """ 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')
+    sevenfifth = np.percentile(arr_1d,75)
+    twenfifth = np.percentile(arr_1d,25)
+    iqr = sevenfifth-twenfifth
+    low = twenfifth-iqr_scale*iqr
+    hi = iqr_scale*iqr+sevenfifth
+    outliers1 = np.array(np.nonzero(arr_1d > hi))[0]
+    outliers2 = np.array(np.nonzero((arr_1d < low)))[0]
+    outliers = np.hstack((outliers1,outliers2))
+    outliers.sort()
+    return outliers,(low, hi)
+
+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.
+    """
+
+    rms_values=[]
+    for i in range(arr_4d.shape[-1]-1):
+        diff_vol = arr_4d[...,i+1] - arr_4d[...,i]
+        rms_values.append(np.sqrt(np.mean(diff_vol ** 2)))
+    return np.array(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 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]``.
+    """
+    i=0
+    while (i < len(diff_indices)):
+        if i is (len(diff_indices)-1):
+            diff_indices = np.insert(diff_indices,i+1,diff_indices[i]+1)
+            break
+        if diff_indices[i+1] == (diff_indices[i]+1):
+            i=i+1
+            continue
+        diff_indices = np.insert(diff_indices,i+1,diff_indices[i]+1)
+        i=i+2
+    return diff_indices
+
+

code/utils/get_outlier.py

@@ -0,0 +1,25 @@
+"""take out the volumne outliers in an image"""
+
+import nibabel as nib
+import matplotlib.pyplot as plt
+import numpy as np
+import numpy.linalg as npl
+from itertools import product
+import diagnostics
+reload(diagnostics) 
+
+img=nib.load('ds005/sub001/BOLD/task001_run001/bold.nii')
+data=img.get_data()
+
+"""get std"""
+volstd = diagnostics.vol_std(data)
+outliers_index, thres = diagnostics.iqr_outliers(volstd)
+
+
+
+
+
+
+
+
+

code/utils/organize_behavior_data.py

@@ -0,0 +1,19 @@
+import numpy as np
+
+project_location="../../"
+data_location=project_location+"data/ds005/"
+
+"""subject 1"""
+
+behav1=np.loadtxt(data_location+'sub001/behav/task001_run001/behavdata.txt',skiprows=1)
+behav2=np.loadtxt(data_location+'sub001/behav/task001_run001/behavdata.txt',skiprows=1)
+behav3=np.loadtxt(data_location+'sub001/behav/task001_run001/behavdata.txt',skiprows=1)
+
+
+"""concatenate them to be 1"""
+behav_total_run=np.concatenate((behav1,behav2,behav3),axis=0)
+
+"""delete the rows that contain -1 in respcat (these are errors in experiment so we should take them out"""
+behav_total_run=np.delete(behav_total_run, behav_total_run[:,5]==-1,axis=0)
+
+#your behav_total_run is ready to use!