improved commenting & code layout

Author: mario

code/python/cni_toolbox/IRM.py

@@ -1,6 +1,6 @@
 '''
 -----------------------------------------------------------------------------
-                                   LICENSE                             
+                                   LICENSE
 
 Copyright 2020 Mario Senden
 
@@ -27,41 +27,41 @@
 class IRM:
     '''
     Input-referred model (IRM) mapping tool.
-    
+
     irm = IRM(parameters) creates an instance of the IRM class.
     parameters is a dictionary with 5 required keys
       - f_sampling: sampling frequency (1/TR)
       - n_samples : number of samples (volumes)
       - n_rows    : number of rows (in-plane resolution)
       - n_cols    : number of columns (in-plance resolution)
       - n_slices  : number of slices
-    
+
     optional inputs are
       - hrf       : either a column vector containing a single hemodynamic
                     response used for every voxel;
                     or a matrix with a unique hemodynamic response along
                     its columns for each voxel.
                     By default the canonical two-gamma hemodynamic response
                     function is generated internally based on the scan parameters.
-    
+
     This class has the following functions
-    
+
       - hrf = IRM.get_hrf()
       - stimulus = IRM.get_stimulus()
       - tc = IRM.get_timecourses()
       - IRM.set_hrf(hrf)
       - IRM.set_stimulus(stimulus)
       - IRM.create_timecourses()
       - results = IRM.mapping(data)
-    
+
 
     Typical workflow:
     1. irm = IRM(params)
     2. irm.set_stimulus()
     3. irm.create_timecourse(FUN,xdata)
     4. results = irm.mapping(data)
     '''
-    
+
     def __init__(self, parameters, hrf = None):
         self.f_sampling = parameters['f_sampling']
         self.p_sampling = 1 / self.f_sampling
@@ -70,7 +70,7 @@ def __init__(self, parameters, hrf = None):
         self.n_cols = parameters['n_cols']
         self.n_slices = parameters['n_slices']
         self.n_total = self.n_rows * self.n_cols * self.n_slices
-        
+
         if hrf != None:
             self.l_hrf = hrf.shape[0]
             if hrf.ndim>2:
@@ -80,23 +80,23 @@ def __init__(self, parameters, hrf = None):
                                                 self.n_total)))),
                                    axis = 0)
             else:
-                self.hrf_fft = fft(np.append(hrf, 
+                self.hrf_fft = fft(np.append(hrf,
                                      np.zeros(self.n_samples)),
                                    axis = 0 )
         else:
             self.l_hrf = int(32 * self.f_sampling)
-            timepoints = np.arange(0, 
+            timepoints = np.arange(0,
                                    self.p_sampling * (self.n_samples +
-                                                      self.l_hrf) - 1, 
+                                                      self.l_hrf) - 1,
                                    self.p_sampling)
             self.hrf_fft = fft(two_gamma(timepoints), axis = 0)
-            
-            
+
+
     def get_hrf(self):
         '''
         Returns
         -------
-        hrf : floating point array
+        hrf: floating point array
             hemodynamic response function(s) used by the class
         '''
         if self.hrf_fft.ndim>1:
@@ -109,34 +109,34 @@ def get_hrf(self):
                     axis = 0)[0:self.l_hrf, :]
         else:
             hrf = ifft(self.hrf_fft, axis = 0)[0:self.l_hrf]
-        
+
         return np.abs(hrf)
-    
+
     def get_stimulus(self):
         '''
         Returns
         -------
         floating point array (1D)
             stimulus used by the class
         '''
-        
+
         return self.stimulus
-    
+
     def get_timecourses(self):
         '''
         Returns
         -------
         floating point array (time-by-gridsize)
             predicted timecourses
         '''
-        
+
         return np.abs(ifft(self.tc_fft, axis = 0)[0:self.n_samples, :])
-    
+
     def set_hrf(self, hrf):
         '''
         Parameters
         ----------
-        hrf : floating point array
+        hrf: floating point array
             hemodynamic response function
         '''
         self.l_hrf = hrf.shape[0]
@@ -147,100 +147,100 @@ def set_hrf(self, hrf):
                                                 self.n_total)))),
                                axis = 0)
         else:
-            self.hrf_fft = fft(np.append(hrf, 
+            self.hrf_fft = fft(np.append(hrf,
                                      np.zeros(self.n_samples)),
                                axis = 0)
-            
+
     def set_stimulus(self, stimulus):
         '''
-    
         Parameters
         ----------
-        stimulus : floating point array (1D)
+        stimulus: floating point array (1D)
             stimulus used by the class
-
         '''
         self.stimulus = stimulus
-        
+
     def create_timecourses(self, FUN, xdata):
-        ''' 
+        '''
         creates predicted timecourses based on the stimulus protocol
         and a range of parameters for an input referred model.
 
         Required inputs are
-        - FUN  : function handle 
-                defining the input referred model
+        - FUN: function handle
+            defining the input referred model
         - xdata: dictionary with p elements (p = number of parameters).
-                 Each cell element needs to contain a column vector
-                 of variable length with parameter values to be explored.
-        ''' 
-        
+            Each cell element needs to contain a column vector
+            of variable length with parameter values to be explored.
+        '''
+        print('\ncreating timecourses')
+
         self.xdata = xdata
         self.n_predictors = len(xdata)
+
         n_observations = np.zeros(self.n_predictors)
         for idx, key in enumerate(self.xdata):
             n_observations[idx] = np.size(self.xdata[key])
-            
         self.n_points = np.prod(n_observations).astype(int)
-        
+
         idx_all = np.arange(self.n_points)
         self.idx = np.zeros((self.n_points, self.n_predictors))
-        
         for p in range(self.n_predictors):
             self.idx[:, p] = (idx_all // (np.prod(n_observations[p+1::]))) \
                 % n_observations[p]
-        self.idx = self.idx.astype(int)       
+        self.idx = self.idx.astype(int)
+
         tc = np.zeros((self.n_samples + self.l_hrf,
                        self.n_points))
+
         x = np.zeros(self.n_predictors)
-        print('\ncreating timecourses')
+
         for j in range(self.n_points):
             for p, key in enumerate(self.xdata):
                 x[p] = self.xdata[key][self.idx[j, p]]
             tc[0:self.n_samples, j] = FUN(self.stimulus, x)
             i = int(j / self.n_points * 21)
             sys.stdout.write('\r')
-            sys.stdout.write("[%-20s] %d%%" 
+            sys.stdout.write("[%-20s] %d%%"
                    % ('='*i, 5*i))
 
         self.tc_fft = fft(tc, axis = 0)
-        
+
     def mapping(self, data, threshold = 100, mask = []):
         '''
         identifies the best fitting timecourse for each voxel and
-        returns a dictionary with keyes corresponding to the
+        returns a dictionary with keys corresponding to the
         parameters specified in xdata plus a key 'corr_fit'
         storing the fitness of the solution.
 
         Required inputs are
-        - data : floating point array
+        - data: floating point array
             empirically observed BOLD timecourses
             whose rows correspond to time (volumes).
 
- 
         Optional inputs are
          - threshold: float
              minimum voxel intensity (default = 100.0)
-         - mask: boolean array 
+         - mask: boolean array
              binary mask for selecting voxels (default = None)
         '''
-        data = np.reshape(data.astype(float), 
+        print('\nmapping model parameters')
+
+        data = np.reshape(data.astype(float),
                         (self.n_samples,
                         self.n_total))
-        
+
         mean_signal = np.mean(data, axis = 0)
         data = zscore(data, axis = 0)
-        
+
         if np.size(mask)==0:
                 mask = mean_signal >= threshold
 
         mask = np.reshape(mask,self.n_total)
         voxel_index = np.where(mask)[0]
         n_voxels = voxel_index.size
-        
+
         mag_d = np.sqrt(np.sum(data[:, mask]**2, axis = 0))
-        
-            
+
         results = {'corr_fit': np.zeros(self.n_total)}
         for key in self.xdata:
             results[key] = np.zeros(self.n_total)
@@ -254,60 +254,60 @@ def mapping(self, data, threshold = 100, mask = []):
                                             axis = 1), axis = 0)), axis = 0))
             tc = tc[:, 0:self.n_samples]
             mag_tc = np.sqrt(np.sum(tc**2, axis = 1))
-            print('\nmapping model parameters')
+
             for m in range(n_voxels):
                 v = voxel_index[m]
-                
+
                 CS = np.matmul(tc, data[:, v]) / (mag_tc * mag_d[m])
                 idx_remove = (CS == np.Inf)| (CS == np.NaN);
                 CS[idx_remove] = 0
-                
+
                 results['corr_fit'][v] = np.max(CS)
                 idx_best = np.argmax(CS)
-                
+
                 for pos, key in enumerate(self.xdata):
                     results[key][v] = self.xdata[key][self.idx[idx_best, pos]]
 
                 i = int(m / n_voxels * 21)
                 sys.stdout.write('\r')
-                sys.stdout.write("[%-20s] %d%%" 
+                sys.stdout.write("[%-20s] %d%%"
                     % ('='*i, 5*i))
-        
+
         else:
             for m in range(n_voxels):
                 v = voxel_index[m]
-                    
+
                 tc = np.transpose(
                     zscore(
                         np.abs(
                             ifft(self.tc_fft *
                                  np.expand_dims(self.hrf_fft[:, v],
                                                 axis = 1), axis = 0)), axis = 0))
-                
+
                 tc = tc[:, 0:self.n_samples]
                 mag_tc = np.sqrt(np.sum(tc**2, axis = 1))
-                    
+
                 CS = np.matmul(tc, data[:, v]) / (mag_tc * mag_d[m])
                 idx_remove = (CS == np.Inf)| (CS == np.NaN);
                 CS[idx_remove] = 0
-                
+
                 results['corr_fit'][v] = np.max(CS)
                 idx_best = np.argmax(CS)
-                
+
                 for pos, key in enumerate(self.xdata):
                     results[key][v] = self.xdata[key][self.idx[idx_best, pos]]
 
                 i = int(m / n_voxels * 21)
                 sys.stdout.write('\r')
-                sys.stdout.write("[%-20s] %d%%" 
+                sys.stdout.write("[%-20s] %d%%"
                     % ('='*i, 5*i))
-        
-    
+
+
         for key in results:
             results[key] = np.squeeze(
             np.reshape(results[key],
                        (self.n_rows,
                         self.n_cols,
                         self.n_slices)))
-            
+
         return results