Merge remote-tracking branch 'upstream/Tests' into soazig-Tests

code/utils/scripts/convolution_high_res_script.py

@@ -19,15 +19,13 @@
 
 from __future__ import absolute_import, division, print_function
 import sys, os
-#TODO : i'm gonna fix the location when i finish the test for stimuli.py
-sys.path.append(os.path.join(os.path.dirname(__file__), "../"))
+sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
 import numpy as np
 import matplotlib.pyplot as plt
 import nibabel as nib
 from stimuli import *
 from scipy.stats import gamma
 from organize_behavior_data import *
-from load_BOLD import *
 
 # Create the necessary directories if they do not exist
 dirs = ['../../../txt_output', '../../../txt_output/conv_high_res',\

code/utils/scripts/convolution_normal_script.py

@@ -18,15 +18,13 @@
 
 from __future__ import absolute_import, division, print_function
 import sys, os
-#TODO : i'm gonna fix the location when i finish the test for stimuli.py
-sys.path.append(os.path.join(os.path.dirname(__file__), "../"))
+sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
 import numpy as np
 import matplotlib.pyplot as plt
 import nibabel as nib
 from stimuli import *
 from scipy.stats import gamma
 from organize_behavior_data import *
-from load_BOLD import *
 
 # Create the necessary directories if they do not exist
 dirs = ['../../../txt_output', '../../../txt_output/conv_normal',\

code/utils/scripts/log_regression_script.py

@@ -12,7 +12,7 @@
 from __future__ import absolute_import, division, print_function
 import sys, os
 #TODO : later change this
-sys.path.append(os.path.join(os.path.dirname(__file__), "../"))
+sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
 import numpy as np
 import matplotlib.pyplot as plt
 import nibabel as nib
@@ -22,7 +22,7 @@
 
 
 # Create the necessary directories if they do not exist
-dirs = ['../../fig','../../fig/log_reg_behav']
+dirs = ['../../../fig','../../../fig/log_reg_behav']
 for d in dirs:
     if not os.path.exists(d):
             os.makedirs(d)

code/utils/scripts/multi_beta_script.py

@@ -0,0 +1,102 @@
+
+"""
+Purpose:
+-----------------------------------------------------------------------------------
+We generate beta values for each single voxels for each subject and save them to 
+files for multi-comparison test.
+-----------------------------------------------------------------------------------
+
+"""
+
+
+import sys, os
+##sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
+sys.path.append(os.path.join(os.path.dirname('__file__'), "../functions/"))
+import numpy as np
+from glm import *
+#from convolution_normal_script import X_matrix
+#from convolution_high_res_script import X_matrix_high_res
+import nibabel as nib
+import matplotlib.pyplot as plt
+from smoothing import *
+
+# Create the necessary directories if they do not exist
+dirs = ['../../../txt_output/multi_beta']
+for d in dirs:
+    if not os.path.exists(d):
+            os.makedirs(d)
+
+# Locate the different paths
+project_path = '../../../'
+# TODO: change it to relevant path
+conv_path = project_path + 'txt_output/conv_normal/'
+conv_high_res_path = project_path + 'txt_output/conv_high_res/'
+
+# select your own subject
+subject_list = [str(i) for i in range(1,17)]
+
+conv_list = [str(i) for i in range(1,5)]
+
+txt_paths = [('ds005_sub' + s.zfill(3) + '_t1r1' +'_cond'+ c.zfill(3),\
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv001_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv002_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv003_canonical.txt', \
+              conv_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv004_canonical.txt', \
+              '../../../data/ds005_2/sub' + s.zfill(3) + '/model/model001/task001_run002' \
+              #for me (Min) it's ds005_2... the filtered data set
+              + '.feat/filtered_func_data_mni.nii.gz',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv_001_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv_002_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv_003_high_res.txt',\
+              conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r1' +'_conv_004_high_res.txt') \
+                for s in subject_list \
+                for c in conv_list]
+
+for txt_path in txt_paths:
+# get 4_d image data
+    name = txt_path[0]
+
+    img = nib.load(txt_path[5])
+    data = img.get_data()
+
+    p = 5
+    # p is the number of columns in our design matrix
+    # it is the number of convolved column plus 1 (a column of 1's)
+
+    X_matrix1 = np.loadtxt(txt_path[1])
+    X_matrix2 = np.loadtxt(txt_path[2])
+    X_matrix3 = np.loadtxt(txt_path[3])
+    X_matrix4 = np.loadtxt(txt_path[4])
+    X_matrix = np.ones((len(X_matrix1),p))
+    X_matrix[...,1] = X_matrix1
+    X_matrix[...,2] = X_matrix2
+    X_matrix[...,3] = X_matrix3
+    X_matrix[...,4] = X_matrix4
+
+    X_matrix_high_res1 = np.loadtxt(txt_path[6])
+    X_matrix_high_res2 = np.loadtxt(txt_path[7])
+    X_matrix_high_res3 = np.loadtxt(txt_path[8])
+    X_matrix_high_res4 = np.loadtxt(txt_path[9])
+    X_matrix_high_res = np.ones((len(X_matrix1),p))
+    X_matrix_high_res[...,1] = X_matrix_high_res1
+    X_matrix_high_res[...,2] = X_matrix_high_res2
+    X_matrix_high_res[...,3] = X_matrix_high_res3
+    X_matrix_high_res[...,4] = X_matrix_high_res4
+
+    beta_4d = glm_beta(data,X_matrix)
+
+    # smooth the data
+    # use high resolution matrix and re-run the regression
+    data_smooth = smoothing(data,1,range(data.shape[-1]))
+    beta_4d_smooth_high_res = glm_beta(data_smooth,X_matrix_high_res)
+    beta_4d_smooth_high_res_task = beta_4d_smooth_high_res[...,1]
+    beta_4d_smooth_high_res_gain = beta_4d_smooth_high_res[...,2]
+    beta_4d_smooth_high_res_loss = beta_4d_smooth_high_res[...,3]
+    beta_4d_smooth_high_res_dist = beta_4d_smooth_high_res[...,4]
+
+    location_of_txt= dirs[0]
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_task.txt",beta_4d_smooth_high_res_task.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_gain.txt",beta_4d_smooth_high_res_gain.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_loss.txt",beta_4d_smooth_high_res_loss.ravel())
+    np.savetxt(location_of_txt + '/' +name[0:17]+ "_beta_dist.txt",beta_4d_smooth_high_res_dist.ravel())
+

code/utils/scripts/multi_comparison_script.py

@@ -0,0 +1,86 @@
+
+"""
+Purpose:
+-----------------------------------------------------------------------------------
+We seek the activated voxel positionsi through multi-comparison of beta values across
+subjects
+
+Step
+-----------------------------------------------------------------------------------
+1. calculate the mean of each single beta values across subject and plot them
+2. calculate the variance of each single beta values across subject and plot them
+3. calculate the t-stat of each single beta values across subject and plot them
+4. calculate the p-value of each single betav values across subject and plot them
+"""
+
+
+import sys, os
+##sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
+sys.path.append(os.path.join(os.path.dirname('__file__'), "../functions/"))
+import numpy as np
+from glm import *
+import nibabel as nib
+import matplotlib.pyplot as plt
+#from scipy.stats import sem
+from smoothing import *
+from scipy.stats import t as t_dist
+#from visualization import *
+
+
+dirs = ['../../../txt_output/multi_beta']
+
+task = dict()
+gain = dict()
+loss = dict()
+dist = dict()
+
+#load all of them
+for x in range(1,17):
+	task[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_task.txt')
+
+for x in range(1,17):
+	gain[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_gain.txt')
+
+for x in range(1,17):
+	loss[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_loss.txt')
+
+for x in range(1,17):
+	dist[x] = np.loadtxt(dirs[0]+'/ds005_sub'+str(x).zfill(3)+'_t1r1_beta_dist.txt')
+
+#calculate mean and plot (let's try for task)
+task_sum = task[1]
+for x in range(2,17):
+	task_sum +=task[x]
+
+task_mean = task_sum/16
+#task_mean.shape -> (902629,0)
+
+#calculate variance and plot
+stdlst = []
+for x in range(1,17):
+	stdlst.append(task[x])
+stdarray = np.array(stdlst)
+task_std = stdarray.std(axis=0)
+#task_std.shape -> (902629,0)
+task_std_reshape = task_std.reshape(91,109,91)
+
+
+#calculate t-stat and plot
+task_tstat = task_mean/task_std
+#task_tstat -> (902629,0)
+task_tstat_reshape = task_tstat.reshape(91,109,91)
+	#plot task_tstat
+
+
+#calculate p-value and plot
+task_pval = t_dist.cdf(abs(task_tstat), 15)
+task_pval_reshape = task_pval.reshape(91,109,91)
+1-task_pval
+	#plot task
+
+
+
+
+
+
+

code/utils/tests/test_logistic_reg.py

@@ -0,0 +1,83 @@
+"""test_diagnostics.py
+Tests for the functions in logistic_reg.py module 
+
+Run with:
+
+    nosetests test_logistic_reg.py
+"""
+
+
+import os, sys
+sys.path.append(os.path.join(os.path.dirname(__file__), "../functions/"))
+from logistic_reg import *
+from organize_behavior_data import *
+from nose.tools import assert_equal
+from numpy.testing import assert_almost_equal, assert_array_equal
+import pandas as pd
+import statsmodels.api as sm
+import pylab as pl
+import numpy as np
+from scipy import stats
+from statsmodels.formula.api import logit, ols
+
+def test_add_gainlossratio():
+	"""Tests whether the gain/loss ratio is properly added in the data frame
+
+	"""
+
+	#load the subject 2's combined runs on the dataframe (use organize_behav_data.py)
+	run = load_in_dataframe(2)
+	gain = run.ix[:,1]
+	loss = run.ix[:,2]
+	#add the column for gain/loss ratio
+	run['ratio'] = gain/loss
+	run_mat=run.as_matrix()
+	run_ratio=run_mat[:,7]
+	#compare with the output from the actual output
+	test_run = load_in_dataframe(2)
+	test_run_added = add_gainlossratio(test_run).as_matrix()
+	test_run_added_ratio = test_run_added[:,7]
+	assert_array_equal(run_ratio,test_run_added_ratio)
+
+
+def test_organize_columns():
+	"""Tests whether columns in the data frame are organized or not for logistic regression
+
+	"""
+
+	run = load_in_dataframe(2)
+	run_added = add_gainlossratio(run)
+	a = run_added.drop('onset', 1) 
+	# drop PTval column
+	run_organized = a.drop('PTval', 1) 
+	# get the column names
+	cols = run_organized.columns.tolist() 
+	# put respcat column into front
+	cols.insert(0, cols.pop(cols.index('respcat'))) 
+	cols.insert(3, cols.pop(cols.index('ratio')))
+	# reorganize
+	run_organized = run_organized.reindex(columns= cols) 
+	# drop error(rescap=-1) in experiment
+	run_final = run_organized.drop(run_organized[run_organized.respcat == -1].index) 
+	test_run_final = organize_columns(run).as_matrix()
+	assert_array_equal(run_final.as_matrix(),test_run_final)
+
+
+
+def test_log_regression():
+	"""Tests the results of logistic regression.
+	Explore on the beta coefficient
+
+	"""
+	run = load_in_dataframe(2)
+	run_added = add_gainlossratio(run)
+	run_final = organize_columns(run_added)
+
+	#fit the logistic regression line
+	fitted = logit("respcat ~ gain + loss", run_final).fit()
+	#get the parameters
+	fitted_params = fitted.params.as_matrix()
+	test_fitted_params = log_regression(run_final).as_matrix()
+	assert_array_equal(fitted_params,test_fitted_params)
+
+