-
Notifications
You must be signed in to change notification settings - Fork 6
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #160 from berkeley-stat159/glm
Glm
- Loading branch information
Showing
1 changed file
with
88 additions
and
33 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,51 +1,106 @@ | ||
| import sys | ||
| sys.path.append(".././utils") | ||
| import sys, os | ||
| 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 | ||
| #from convolution_normal_script import X_matrix | ||
| #from convolution_high_res_script import X_matrix_high_res | ||
| from load_BOLD import * | ||
| import nibabel as nib | ||
| import matplotlib.pyplot as plt | ||
| from smoothing import * | ||
|
|
||
| location_of_plot = "../../plots/" | ||
| location_of_data = "../../data/ds005/" | ||
| # Create the necessary directories if they do not exist | ||
| dirs = ['../../../txt_output/mrss',\ | ||
| '../../../fig/glm_fitted'] | ||
| for d in dirs: | ||
| if not os.path.exists(d): | ||
| os.makedirs(d) | ||
|
|
||
| # get 4_d image data | ||
| data = load_img(1,1) | ||
| # 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/' | ||
|
|
||
| beta_4d = glm_multi(data,X_matrix) | ||
| MRSS, fitted, residuals = glm_diagnostics(beta_4d, X_matrix, data) | ||
| # select your own subject | ||
| #subject_list = [str(i) for i in range(1,17)] | ||
| subject_list = ['1','5', '11'] | ||
| run_list = [str(i) for i in range(1,4)] | ||
| conv_list = [str(i) for i in range(1,5)] | ||
|
|
||
| # smooth the data and re-run the regression | ||
| data_smooth = smoothing(data,1,range(data.shape[-1])) | ||
| txt_paths = [('ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv'+ c.zfill(3),\ | ||
| conv_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv001_canonical.txt', \ | ||
| conv_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv002_canonical.txt', \ | ||
| conv_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv003_canonical.txt', \ | ||
| conv_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv004_canonical.txt', \ | ||
| '../../data/ds005/sub' + s.zfill(3) + '/BOLD/task001_run' \ | ||
| + r.zfill(3) + '/bold.nii',\ | ||
| conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv_001_high_res.txt',\ | ||
| conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv_002_high_res.txt',\ | ||
| conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv_003_high_res.txt',\ | ||
| conv_high_res_path + 'ds005_sub' + s.zfill(3) + '_t1r' + r +'_conv_004_high_res.txt') \ | ||
| for r in run_list \ | ||
| for s in subject_list \ | ||
| for c in conv_list] | ||
|
|
||
| beta_4d_smooth = glm_multi(data_smooth,X_matrix) | ||
| MRSS_smooth, fitted_smooth, residuals_smooth = glm_diagnostics(beta_4d_smooth, X_matrix, data_smooth) | ||
| 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) | ||
| MRSS, fitted, residuals = glm_mrss(beta_4d, X_matrix, data) | ||
|
|
||
| # use high resolution to create our design matrix | ||
| beta_4d_high_res = glm_multi(data,X_matrix_high_res) | ||
| MRSS_high_res, fitted_high_res, residuals_high_res = glm_diagnostics(beta_4d_high_res, X_matrix_high_res, data) | ||
| # smooth the data and re-run the regression | ||
| data_smooth = smoothing(data,1,range(data.shape[-1])) | ||
| beta_4d_smooth = glm_beta(data_smooth,X_matrix) | ||
| MRSS_smooth, fitted_smooth, residuals_smooth = glm_mrss(beta_4d_smooth, X_matrix, data_smooth) | ||
|
|
||
| # use high resolution to create our design matrix | ||
| beta_4d_high_res = glm_beta(data,X_matrix_high_res) | ||
| MRSS_high_res, fitted_high_res, residuals_high_res = glm_mrss(beta_4d_high_res, X_matrix_high_res, data) | ||
|
|
||
| print ("MRSS of multiple regression: "+str(np.mean(MRSS))) | ||
| print ("MRSS of multiple regression by using high resoultion design matrix: "+str(np.mean(MRSS_high_res))) | ||
| print ("MRSS of multiple regression using the smoothed data: "+str(np.mean(MRSS_smooth))) | ||
|
|
||
| plt.plot(data[4,22,11], label = "actual") | ||
| plt.plot(fitted[4,22,11], label = "fitted") | ||
| plt.plot(fitted_high_res[4,22,11], label = "fitted_high_res") | ||
| plt.plot(data[4,22,11], label = "actual") | ||
| plt.plot(fitted[4,22,11], label = "fitted") | ||
| plt.plot(fitted_high_res[4,22,11], label = "fitted_high_res") | ||
|
|
||
| plt.title("Subject 001, voxel (4,22,11) actual vs fitted") | ||
| plt.legend(loc = "upper left", fontsize = "smaller") | ||
| plt.savefig(location_of_plot + "glm_fitted.png") | ||
| plt.close() | ||
| plt.title(name[0:17]+"voxel (4,22,11) actual vs fitted") | ||
| plt.legend(loc = "upper left", fontsize = "smaller") | ||
| plt.savefig(dirs[1] + '/'+ name[0:17]+ "_glm_fitted.png") | ||
| plt.close() | ||
|
|
||
| location_of_txt="../txt_files/" | ||
| file = open(location_of_txt+'ds005_mrss_result.txt', "w") | ||
| file.write("MRSS of multiple regression for subject 1, run 1 is: "+str(np.mean(MRSS))+"\n") | ||
| file.write("\n") | ||
| file.write("MRSS of multiple regression for subject 1, run 1, using the smoothed data is: "+str(np.mean(MRSS_smooth))+"\n") | ||
| file.close() | ||
| location_of_txt= dirs[0] | ||
| file = open(location_of_txt+ '/' + name[0:17] +'_mrss_result.txt', "w") | ||
| file.write("MRSS of multiple regression for" +name[0:17]+ "is: "+str(np.mean(MRSS))+"\n") | ||
| file.write("\n") | ||
| file.write("MRSS of multiple regression for" +name[0:17]+ "using the smoothed data is: "+str(np.mean(MRSS_smooth))+"\n") | ||
| file.write("\n") | ||
| file.write("MRSS of multiple regression for" +name[0:17]+ "using high_res design matrix is: "+str(np.mean(MRSS_high_res))+"\n") | ||
| file.close() |