From b0be0c9ef09a639bfb43ac85671512083e55fef9 Mon Sep 17 00:00:00 2001 From: Ian Charest Date: Wed, 20 Jan 2021 10:27:43 +0000 Subject: [PATCH 1/5] initial commit for the handling of columns with all zeros in design matrices. --- glmdenoise/fit_runs.py | 7 +++++-- glmdenoise/utils/make_poly_matrix.py | 25 ++++++++++++++----------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/glmdenoise/fit_runs.py b/glmdenoise/fit_runs.py index 56fda73..b424518 100644 --- a/glmdenoise/fit_runs.py +++ b/glmdenoise/fit_runs.py @@ -6,8 +6,9 @@ def fit_runs(data, design): """Fits a least square of combined runs. - The matrix addition is equivalent to concatenating the list of data and the list of - design and fit it all at once. However, this is more memory efficient. + The matrix addition is equivalent to concatenating the list of data and + the list of design and fit it all at once. However, this is more memory + efficient. Arguments: runs {list} -- List of runs. Each run is an TR x voxel sized array @@ -16,6 +17,8 @@ def fit_runs(data, design): Returns: [array] -- betas from fit + + TODO handle zero regressors. """ X = np.vstack(design) diff --git a/glmdenoise/utils/make_poly_matrix.py b/glmdenoise/utils/make_poly_matrix.py index c61c99c..8b15952 100644 --- a/glmdenoise/utils/make_poly_matrix.py +++ b/glmdenoise/utils/make_poly_matrix.py @@ -1,6 +1,7 @@ from sklearn.preprocessing import normalize import numpy as np + def make_project_matrix(X): """ Calculates a projection matrix @@ -9,10 +10,13 @@ def make_project_matrix(X): Returns: array: Projection matrix size of X.shape[0] x X.shape[0] + + TODO handle zero regressors. """ X = np.mat(X) return np.eye(X.shape[0]) - (X*(np.linalg.inv(X.T*X)*X.T)) + def make_poly_matrix(n, degrees): """Calculates a matrix of polynomials used to regress them out of your data @@ -30,12 +34,11 @@ def make_poly_matrix(n, degrees): for i, d in enumerate(degrees): polyvector = np.mat(time_points**d) - if i > 0: # project out the other polynomials + if i > 0: # project out the other polynomials polyvector = make_project_matrix(polys[:, :i]) * polyvector polys[:, i] = normalize(polyvector.T) - return polys # make_project_matrix(polys) - + return polys # make_project_matrix(polys) def select_noise_regressors(r2_nrs, pcstop=1.05): @@ -57,15 +60,15 @@ def select_noise_regressors(r2_nrs, pcstop=1.05): best = -np.Inf for p in range(1, numpcstotry): - # if better than best so far - if curve[p] > best: + # if better than best so far + if curve[p] > best: - # record this number of PCs as the best - chosen = p - best = curve[p] + # record this number of PCs as the best + chosen = p + best = curve[p] - # if we are within opt.pcstop of the max, then we stop. - if (best * pcstop) >= curve.max(): - break + # if we are within opt.pcstop of the max, then we stop. + if (best * pcstop) >= curve.max(): + break return chosen From 9b7152a77635ee0475fb6d6c9964e51baf84b497 Mon Sep 17 00:00:00 2001 From: Ian Charest Date: Wed, 20 Jan 2021 12:18:08 +0000 Subject: [PATCH 2/5] now handling designs with zero regressors --- example_singletrial.py | 75 ++++++++++++++++++++++++++++ glmdenoise/fit_runs.py | 9 ++-- glmdenoise/utils/make_poly_matrix.py | 6 +-- glmdenoise/utils/optimiseHRF.py | 15 +++--- 4 files changed, 92 insertions(+), 13 deletions(-) create mode 100644 example_singletrial.py diff --git a/example_singletrial.py b/example_singletrial.py new file mode 100644 index 0000000..8e643d5 --- /dev/null +++ b/example_singletrial.py @@ -0,0 +1,75 @@ +import os +import glob +import numpy as np +import nibabel as nib +import pandas as pd +from glmdenoise.utils.make_design_matrix import make_design +from glmdenoise.utils.gethrf import getcanonicalhrf +from glmdenoise.utils.normalisemax import normalisemax + +sub = 2 +ses = 1 +stimdur = 0.5 +TR = 2 + +proj_path = os.path.join( + '/home', + 'adf', + 'charesti', + 'data', + 'arsa-fmri', + 'BIDS') + +data_path = os.path.join( + proj_path, + 'derivatives', + 'fmriprep', + 'sub-{}', + 'ses-{}', + 'func') + +design_path = os.path.join( + proj_path, + 'sub-{}', + 'ses-{}', + 'func') + +runs = glob.glob( + os.path.join(data_path.format(sub, ses), '*preproc*nii.gz')) +runs.sort() +runs = runs[:-1] + +eventfs = glob.glob( + os.path.join(design_path.format(sub, ses), '*events.tsv')) +eventfs.sort() + +data = [] +design = [] + +hrf = normalisemax(getcanonicalhrf(stimdur, TR)) + +for i, (run, eventf) in enumerate(zip(runs, eventfs)): + print(f'run {i}') + y = nib.load(run).get_fdata().astype(np.float32) + dims = y.shape + y = np.moveaxis(y, -1, 0) + y = y.reshape([y.shape[0], -1]) + + n_volumes = y.shape[0] + + # Load onsets and item presented + onsets = pd.read_csv(eventf, sep='\t')["onset"].values + items = pd.read_csv(eventf, sep='\t')["stimnumber"].values + n_events = len(onsets) + + # Create design matrix + events = pd.DataFrame() + events["duration"] = [stimdur] * n_events + events["onset"] = onsets + events["trial_type"] = items + + # pass in the events data frame. the convolving of the HRF now + # happens internally + design.append(make_design(events, TR, n_volumes, hrf).astype(np.float32)) + data.append(y) + diff --git a/glmdenoise/fit_runs.py b/glmdenoise/fit_runs.py index b424518..d2f7c16 100644 --- a/glmdenoise/fit_runs.py +++ b/glmdenoise/fit_runs.py @@ -1,4 +1,5 @@ import numpy as np +from glmdenoise.utils.optimiseHRF import olsmatrix import warnings warnings.simplefilter(action="ignore", category=FutureWarning) @@ -12,17 +13,17 @@ def fit_runs(data, design): Arguments: runs {list} -- List of runs. Each run is an TR x voxel sized array - DM {list} -- List of design matrices. Each design matrix - is an TR x predictor sizec array + design {list} -- List of design matrices. Each design matrix + is an TR x predictor sized array Returns: [array] -- betas from fit - TODO handle zero regressors. """ X = np.vstack(design) - X = np.linalg.inv(X.T @ X) @ X.T + # X = np.linalg.inv(X.T @ X) @ X.T + X = olsmatrix(X) betas = 0 start_col = 0 diff --git a/glmdenoise/utils/make_poly_matrix.py b/glmdenoise/utils/make_poly_matrix.py index 8b15952..e6b7158 100644 --- a/glmdenoise/utils/make_poly_matrix.py +++ b/glmdenoise/utils/make_poly_matrix.py @@ -1,3 +1,4 @@ +from glmdenoise.utils.optimiseHRF import olsmatrix from sklearn.preprocessing import normalize import numpy as np @@ -11,10 +12,9 @@ def make_project_matrix(X): Returns: array: Projection matrix size of X.shape[0] x X.shape[0] - TODO handle zero regressors. """ - X = np.mat(X) - return np.eye(X.shape[0]) - (X*(np.linalg.inv(X.T*X)*X.T)) + + return np.eye(X.shape[0]) - (X @ olsmatrix(X)) def make_poly_matrix(n, degrees): diff --git a/glmdenoise/utils/optimiseHRF.py b/glmdenoise/utils/optimiseHRF.py index 933ecd3..92c96dd 100644 --- a/glmdenoise/utils/optimiseHRF.py +++ b/glmdenoise/utils/optimiseHRF.py @@ -236,8 +236,9 @@ def olsmatrix(X): """OLS regression what we want to do is to perform OLS regression using - and obtain the parameter estimates. this is accomplished - by inv(X'\*X)\*X'\*y = f\*y where y is the data (samples x cases). + and obtain the parameter estimates. this is accomplished + by np.linalg.inv(X.T @ X) @ X.T @ y = f @ y where y is the + data (samples x cases). what this function does is to return which has dimensions parameters x samples. @@ -268,7 +269,8 @@ def olsmatrix(X): # report warning if not np.any(good) == 1: print( - "regressors are all zeros; we will estimate a 0 weight for those regressors." + "regressors are all zeros. \n" + "we will estimate a 0 weight for those regressors." ) f = np.zeros((X.shape[1], X.shape[0])) return f @@ -276,15 +278,16 @@ def olsmatrix(X): # do it if np.any(bad): print( - "One or more regressors are all zeros; we will estimate a 0 weight for those regressors." + "One or more regressors are all zeros. \n" + "we will estimate a 0 weight for those regressors." ) f = np.zeros((X.shape[1], X.shape[0])) - X = np.mat(X) + f[good, :] = np.linalg.inv( X[:, good].T @ X[:, good]) @ X[:, good].T else: - X = np.mat(X) + f = np.linalg.inv(X.T @ X) @ X.T return f From b59df78dfb22c54fb163e4b02ee055e8cc10382b Mon Sep 17 00:00:00 2001 From: Ian Charest Date: Wed, 20 Jan 2021 14:56:36 +0000 Subject: [PATCH 3/5] fixes #47, tests pass --- glmdenoise/pyGlmdenoise.py | 23 ++++++++++++++++------- glmdenoise/utils/optimiseHRF.py | 5 +++-- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/glmdenoise/pyGlmdenoise.py b/glmdenoise/pyGlmdenoise.py index 2f57837..c209da6 100644 --- a/glmdenoise/pyGlmdenoise.py +++ b/glmdenoise/pyGlmdenoise.py @@ -6,7 +6,11 @@ from joblib import Parallel, delayed from sklearn.preprocessing import normalize from glmdenoise.utils.make_design_matrix import make_design -from glmdenoise.utils.optimiseHRF import mtimesStack, olsmatrix, optimiseHRF +from glmdenoise.utils.optimiseHRF import (mtimesStack, + olsmatrix, + optimiseHRF, + calccod, + calccodStack) from glmdenoise.select_noise_regressors import select_noise_regressors from glmdenoise.utils.normalisemax import normalisemax from glmdenoise.utils.gethrf import getcanonicalhrf @@ -15,7 +19,8 @@ from glmdenoise.whiten_data import whiten_data from glmdenoise.fit_runs import fit_runs from glmdenoise.cross_validate import cross_validate -from glmdenoise.utils.make_poly_matrix import make_poly_matrix, make_project_matrix +from glmdenoise.utils.make_poly_matrix import (make_poly_matrix, + make_project_matrix) warnings.simplefilter(action="ignore", category=FutureWarning) @@ -255,11 +260,15 @@ def fit(self, design, data, tr): print('Calculating overall R2 of final fit...') # The below does not currently work, see #47 - # stackdesign = np.vstack(whitened_design) - # modelfits = mtimesStack(olsmatrix(stackdesign), whitened_data) - # self.results['R2s'] = calccodStack(whitened_data, modelfits) - # self.results['R2runs'] = [calccod(whitened_data[c_run], modelfits[c_run], 0) - # for c_run in range(self.n_runs)] + modelfits = [((olsmatrix(x) @ y).T @ x.T).T for x, y in zip( + whitened_design, whitened_data)] + # calculate run-wise fit + + self.results['R2s'] = calccodStack(self.data, modelfits) + self.results['R2runs'] = [calccod( + cdata, + mfit, + 0, 0, 0) for cdata, mfit in zip(self.data, modelfits)] print('Done') print('Calculating SnR...') diff --git a/glmdenoise/utils/optimiseHRF.py b/glmdenoise/utils/optimiseHRF.py index 92c96dd..abb7e0c 100644 --- a/glmdenoise/utils/optimiseHRF.py +++ b/glmdenoise/utils/optimiseHRF.py @@ -94,7 +94,7 @@ def constructStimulusMatrix(v, prenumlag, postnumlag, wantwrap=0): return f -def calccod(x, y, wantgain=0, wantmeansub=1): +def calccod(x, y, dim=None, wantgain=0, wantmeansub=1): """Calculate the coefficient of determination Args: @@ -153,7 +153,8 @@ def calccod(x, y, wantgain=0, wantmeansub=1): """ # input - dim = np.argmax(x.shape) + if dim is None: + dim = np.argmax(x.shape) # handle gain if wantgain: From 4868b0e9a28ce5b9bac76ecdbbd2960e6ee1e733 Mon Sep 17 00:00:00 2001 From: Ian Charest Date: Wed, 20 Jan 2021 15:35:37 +0000 Subject: [PATCH 4/5] fix for plotting in example.py --- example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example.py b/example.py index 8449e8d..fa068b7 100644 --- a/example.py +++ b/example.py @@ -85,4 +85,4 @@ def format_time(time): fit_dur = format_time(time.time()-start) print('Fit took {}!'.format(fit_dur) -gd.plot_figures() +gd.plot_figures(spatialdims=dims[:3]) From 46ae36231ec103d4f1c6f35e98ce8fc72359603d Mon Sep 17 00:00:00 2001 From: Ian Charest Date: Wed, 20 Jan 2021 16:20:36 +0000 Subject: [PATCH 5/5] some bug tweaks, and figure reporting updated --- glmdenoise/pyGlmdenoise.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/glmdenoise/pyGlmdenoise.py b/glmdenoise/pyGlmdenoise.py index c209da6..3e8bc99 100644 --- a/glmdenoise/pyGlmdenoise.py +++ b/glmdenoise/pyGlmdenoise.py @@ -264,11 +264,11 @@ def fit(self, design, data, tr): whitened_design, whitened_data)] # calculate run-wise fit - self.results['R2s'] = calccodStack(self.data, modelfits) + self.results['R2s'] = calccodStack(whitened_data, modelfits) self.results['R2runs'] = [calccod( cdata, mfit, - 0, 0, 0) for cdata, mfit in zip(self.data, modelfits)] + 0, 0, 0) for cdata, mfit in zip(whitened_data, modelfits)] print('Done') print('Calculating SnR...') @@ -405,13 +405,23 @@ def plot_figures(self, report=None, spatialdims=None): dtype='scaled' ) - # report.plot_image(self.results['R2s'], 'FinalModel') - """ TODO - for r in range(n_runs): + print('plotting R2 final') + report.plot_image( + self.full_image( + self.results['R2s']), + 'R2' + ) + + print('plotting R2 run-wise') + for r in range(self.n_runs): report.plot_image( - self.results['R2srun'][r], 'FinalModel_run%02d') - """ + self.full_image( + self.results['R2runs'][r]), + 'R2run_{}'.format(r) + ) + # PC weights + print('plotting PC weights') weights_mats = self.results['PCA_weights'].ravel() weights = np.asarray(np.concatenate(weights_mats)).ravel() thresh = np.percentile(np.abs(weights), 99) @@ -426,6 +436,7 @@ def plot_figures(self, report=None, spatialdims=None): drange=[-thresh, thresh] ) + print('saving html report') # stores html report report.save()