ols inv when zero regressors

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])

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)
+

glmdenoise/fit_runs.py

@@ -1,25 +1,29 @@
 import numpy as np
+from glmdenoise.utils.optimiseHRF import olsmatrix
 import warnings
 warnings.simplefilter(action="ignore", category=FutureWarning)
 
 
 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
-        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
+
     """
 
     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

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(whitened_data, modelfits)
+        self.results['R2runs'] = [calccod(
+            cdata,
+            mfit,
+            0, 0, 0) for cdata, mfit in zip(whitened_data, modelfits)]
         print('Done')
 
         print('Calculating SnR...')
@@ -396,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)
@@ -417,6 +436,7 @@ def plot_figures(self, report=None, spatialdims=None):
                     drange=[-thresh, thresh]
                 )
 
+        print('saving html report')
         # stores html report
         report.save()
 

glmdenoise/utils/make_poly_matrix.py

@@ -1,6 +1,8 @@
+from glmdenoise.utils.optimiseHRF import olsmatrix
 from sklearn.preprocessing import normalize
 import numpy as np
 
+
 def make_project_matrix(X):
     """ Calculates a projection matrix
 
@@ -9,9 +11,11 @@ def make_project_matrix(X):
 
     Returns:
         array: Projection matrix size of X.shape[0] x X.shape[0]
+
     """
-    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):
     """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

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:
@@ -236,8 +237,9 @@ def olsmatrix(X):
     """OLS regression
 
     what we want to do is to perform OLS regression using <X>
-    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 <f> which has dimensions
     parameters x samples.
@@ -268,23 +270,25 @@ 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
 
     # 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