ephystoolkit.EphysToolkit

   1# class imports
   2import re
   3import math
   4import json
   5import warnings
   6import os.path
   7import csv
   8import mat73
   9
  10import numpy as np
  11import scipy.io
  12import pandas as pd
  13import rhd
  14from scipy import stats
  15from scipy import signal
  16
  17from pathexplorer.PathExplorer import path_explorer
  18
  19
  20class ephys_toolkit:
  21
  22    def __init__(self):
  23        
  24        self._SAMPLING_RATE = 20000
  25        self._modpath = os.path.dirname(__file__)
  26
  27    def _bin_events(self, bin_size, events):
  28        self.frames = bin_size ** -1
  29        self.numerator = self._SAMPLING_RATE / self.frames
  30
  31        return events / self.numerator
  32
  33    def _minmax_norm(self, array):
  34        """
  35        Normalize an array with minmax normalization.
  36        """
  37        return (array - array.min()) / (array.max() - array.min())
  38
  39    def _zscore_norm(self, array):
  40        """
  41        Normalize an array with z score normalization.
  42        """
  43        return stats.zscore(array)
  44
  45    def _average_response(self, array, stim_reps):
  46        """
  47        Return the average frequency of a unit's response
  48        across stimulus repeats.
  49        """
  50        return (array / stim_reps) * self.frames
  51    
  52    def apply_temporal_smoothing(self, x, k, t_axis):
  53        """
  54        Apply temporal smoothing to an array, matrix, or tensor 
  55        with a convolution kernel. 
  56        
  57        Args:
  58        - x: Input data.
  59        - k: Convolution kernel.
  60        - t_axis: The dimension in the input data corresponding to time (0 if one dimensional).
  61        
  62        """
  63        return np.apply_along_axis(
  64                        lambda m: np.convolve(
  65                            m, k, 
  66                            mode='same'
  67                        ), axis=t_axis, arr=x
  68                    )
  69    
  70    def static_grating(
  71            self,
  72            windowSizeX_Pixel, windowSizeY_Pixel, # size of the stimulus window in pixels
  73            windowSizeX_Visual, windowSizeY_Visual, # size of the stimulus window in degrees
  74            spatialFrequency, # spatial frequency in cpd
  75            ang, # orientation angle of the grating in degrees
  76            phi, # phase of the grating in degrees
  77            diameter = None # dimater of the circular patch in degrees
  78    ):    
  79        """
  80        Generate a matrix of pixel intensities   
  81        representing a static grating stimulus
  82        with the given parameters.
  83        
  84        Args:
  85        
  86        - windowSizeX_Pixel: Horizontal size of the stimulus screen in pixels.
  87        - windowSizeY_Pixel: Vertical size of the stimulus screen in pixels.
  88        - windowSizeX_Visual: Horizontal size of the stimulus screen in degrees.
  89        - windowSizeY_Visual: Vertical size of the stimulus screen in degrees.
  90        - spatialFrequency: Spatial frequency in cycles per degre (cpd).
  91        - ang: Orientation angle of the grating in degrees.
  92        - phi: Phase of the grating in degrees
  93        - diameter: Dimater of the circular patch in degrees. 
  94        """
  95        
  96        if ((windowSizeX_Pixel <=10)
  97            or (windowSizeX_Pixel%2 != 0)
  98        ):
  99            print('windowSizeX_Pixel must be greater than ten and even')
 100
 101        if ((windowSizeY_Pixel <=10)
 102            or (windowSizeY_Pixel%2 != 0)
 103        ):
 104            print('windowSizeY_Pixel must be greater than ten and even')
 105
 106        if (math.floor(windowSizeX_Pixel/windowSizeY_Pixel*1000) 
 107            != math.floor(windowSizeX_Visual/windowSizeY_Visual*1000)
 108        ):
 109            print('The ratio of X and Y for Pixel and Visual are different!');
 110
 111        if spatialFrequency < 0:
 112            error('spatialFrequency less than zero')
 113
 114        # x, y 
 115        x_center = windowSizeX_Pixel/2
 116        y_center = windowSizeY_Pixel/2
 117
 118        x_range = np.arange(-x_center, x_center,1)
 119        y_range = np.arange(-y_center, y_center,1)
 120        x, y = np.meshgrid(x_range,y_range)
 121
 122        # theta
 123        theta = -np.radians(ang);
 124        xyTheta = y * np.cos(theta) - x * np.sin(theta);
 125        
 126        # phi
 127        phi = np.radians(phi)
 128
 129        # Spatial Frenquency
 130        degreePerPixel = windowSizeX_Visual / windowSizeX_Pixel
 131        sf = spatialFrequency  * degreePerPixel # cycles / pixel
 132        sw = 2 * np.pi * sf # radian / pixel
 133
 134        # contrast
 135        pixelIntensity = (np.sin(sw * xyTheta - phi)) # range: -1 to 1
 136
 137        # round image
 138        if diameter == None:
 139            pass
 140        else:
 141            diameter = diameter / degreePerPixel
 142            r = diameter/2;
 143            c_mask = ( (x**2+y**2) <= r**2 )
 144            pixelIntensity = pixelIntensity*c_mask 
 145
 146        return pixelIntensity
 147
 148    def drifting_grating(
 149            self,
 150            windowSizeX_Pixel, windowSizeY_Pixel, # size of the stimulus window in pixels
 151            windowSizeX_Visual, windowSizeY_Visual, # size of the stimulus window in degrees
 152            tf,  # Temporal frequency - pass as an integer or floating point value
 153            dt,  # Time step value - pass as a floating point value
 154            t,  # Total duration of the stimulus - pass as an int or float of the appropriate time unit
 155            spatialFrequency, # spatial frequency in cpd
 156            ang, # orientation angle of the grating in degrees
 157            diameter = None # dimater of the circular patch in degrees
 158    ):
 159        """
 160        Returns a list of matricies representing
 161        frames of a drifitng grating stimulus.
 162        
 163        Args:
 164        - windowSizeX_Pixel: Horizontal size of the stimulus screen in pixels.
 165        - windowSizeY_Pixel: Vertical size of the stimulus screen in pixels.
 166        - windowSizeX_Visual: Horizontal size of the stimulus screen in degrees.
 167        - windowSizeY_Visual: Vertical size of the stimulus screen in degrees.    
 168        - tf: Temporal frequency - pass as an integer or floating point value.
 169        - dt: Time step value - pass as a floating point value.
 170        - t: Total duration of the stimulus - pass as an int or float of the appropriate time unit.
 171        - spatialFrequency: Spatial frequency in cycles per degre (cpd).
 172        - ang: Orientation angle of the grating in degrees.
 173        - diameter: Dimater of the circular patch in degrees. 
 174        
 175        """
 176
 177        tensor = []
 178        params = []
 179        deg_step = dt * 360 * tf
 180
 181        d = np.arange(dt, t, dt)
 182
 183        phi = 0
 184        for x in d:
 185            m = self.static_grating(
 186                windowSizeX_Pixel, windowSizeY_Pixel, 
 187                windowSizeX_Visual, windowSizeY_Visual, 
 188                spatialFrequency,
 189                ang,
 190                phi,
 191                diameter = diameter
 192            ) 
 193            tensor.append(m)
 194            params.append([sf, tf, ori, phase])
 195            phi += deg_step
 196
 197        return tensor, params
 198
 199    def _discrete_radius(self, dim=tuple, radius=int):
 200        x = np.linspace(-1, 1, dim[1])
 201        y = np.linspace(-1, 1, dim[0])
 202        
 203
 204        m = []
 205        for i0 in range(len(x)):
 206            row = []
 207            for i1 in range(len(y)):
 208                if ((x[i0] ** 2 + y[i1] ** 2)
 209                        < ((radius / 360) ** 2)):
 210                    row.append(1)
 211                else:
 212                    row.append(0)
 213            m.append(row)
 214
 215        return np.array(m).T
 216
 217    def _gaussian_radius(self, dim:tuple, radius:int):
 218        x, y = np.meshgrid(
 219            np.linspace(-1, 1, dim[0]),
 220            np.linspace(-1, 1, dim[1])
 221        )
 222        
 223        
 224        d = np.sqrt(x * x + y * y)
 225        sigma, mu = radius / 360, 0.0
 226        g = np.exp(-((d - mu) ** 2 / (2.0 * sigma ** 2)))
 227
 228        return g
 229
 230    def spike_sorting_metrics(self, file_path):
 231        """
 232        Returns a dataframe with the spike sorting metrics
 233        of a given recording section.
 234        
 235        Args:
 236        
 237        - file_path: Path to the spike sorting metrics file.
 238        """
 239
 240        with open(file_path, 'r') as sorting_file:
 241            sorting_info = json.load(sorting_file)
 242
 243        spike_sorting_data = [
 244            # [cluster['label'] for cluster in sorting_info['clusters']],
 245            [cluster['metrics']['isolation'] for cluster in sorting_info['clusters']],
 246            [cluster['metrics']['noise_overlap'] for cluster in sorting_info['clusters']]
 247        ]
 248
 249        ss_df = pd.DataFrame(np.array(spike_sorting_data).T)
 250        ss_df.columns = [
 251            # 'cluster', 
 252            'isolation', 'noise_overlap'
 253        ]
 254
 255        return ss_df
 256
 257    # functions to make concatenated across trial and non concatenated across trial rasters
 258    def _concatenated_raster(self, stims, spikes, thresh=tuple):
 259        if thresh == tuple:
 260            r = np.array([spikes - st for st in stims])
 261            raster = np.concatenate(r)
 262        else:
 263            r = np.array([spikes - st for st in stims])
 264            ti = np.where(np.logical_and(r <= thresh[1], r >= thresh[0]))
 265            raster = r[ti]
 266        return raster
 267
 268    def _unconcatenated_raster(self, stims, spikes, thresh=tuple):
 269        if thresh == tuple:
 270            rasters = np.array([spikes - st for st in stims])
 271        else:
 272            rasters = []
 273            for i, st in enumerate(stims):  # enumerate to make an initial array then vstack
 274                unthreshed = spikes - st
 275                i = np.where(np.logical_and(unthreshed <= thresh[1], unthreshed >= thresh[0]))
 276                rasters.append(list(unthreshed[i]))
 277        return rasters
 278
 279    def _linear_extrapolation(self, xd, yd, x):
 280        y2, y1 = yd[-1], yd[-2]
 281        x2, x1 = xd[-1], xd[-2]
 282
 283        return y1+(((x-x1)/(x2-x1))*(y2-y1))
 284
 285    def raster(
 286            self,
 287            stims,  # Array of stimulus onset times
 288            spikes,  # Array of spike onset times
 289            thresh=tuple,  # Bounding threshold around the stimulus onset at t = 0 - pass as a tuple
 290            concatenate=True  # Concatenate rasters across trials - pass False to return unconcatenated rasters
 291    ):
 292        """
 293        Returns an array representing a raster of spike times centered 
 294        around the onset times of a given stimulus.
 295        
 296        Args:
 297        
 298        - stims: Array of stimulus onset times.
 299        - spikes: Array of spike onset times.
 300        - thresh: Bounding threshold around the stimulus onset at t = 0 - pass as a tuple.
 301        - concatenate: Concatenate rasters across trials; pass False to return unconcatenated rasters.
 302        
 303        """
 304
 305        if concatenate:
 306            return self._concatenated_raster(stims, spikes, thresh)
 307        else:
 308            return self._unconcatenated_raster(stims, spikes, thresh)
 309        
 310
 311
 312class load_experiment(ephys_toolkit):
 313    """
 314    Create an experiment object for a given recording block.
 315    Takes the spike data file path as the first argument and the
 316    stimulus data file path as the second argument. Initializing
 317    an experiment object generates some important class attributes:
 318    
 319    .stim_data: A pandas dataframe with the stimulus data.
 320    
 321    .spike_data: A dictionary object with the spiking data
 322                 of all the identified clusters.
 323                 
 324    Args:
 325    
 326    - spikefile: Path to the spike data file.
 327    - stimfile: Path to the stimulus data file.
 328    """
 329
 330    def __init__(self, spikefile, stimfile, logfile, **kwargs):
 331        ephys_toolkit.__init__(self)
 332        
 333        self._stim_m73 = None
 334        self._spike_m73 = None
 335        self._spikefile = spikefile # path to the spike file
 336        self._stimfile = stimfile # path to the stim file
 337        self._depth_data = kwargs['kwargs']['depth_data'] 
 338        
 339        try:
 340            self.spikes_mat = scipy.io.loadmat(spikefile)
 341            self.spikes = self.spikes_mat['Data'][0] # raw spikes dictionary
 342        except NotImplementedError:
 343            self._spike_m73 = True
 344            self.spikes_mat = mat73.loadmat(spikefile)
 345            self.spikes = self.spikes_mat['Data']
 346            
 347        try:
 348            self.stims_mat = scipy.io.loadmat(stimfile)
 349            self.stims = self.stims_mat['StimulusData'][0][0] # raw stims dictionary
 350        except NotImplementedError:
 351            self._stim_m73 = True
 352            self.stims_mat = mat73.loadmat(stimfile)
 353            self.stims = self.stims_mat['StimulusData']
 354            
 355        self._logfile = logfile # stimulus log data
 356        self._nodf = False
 357        self._init_stim_data()
 358        self._init_spike_data()
 359        self.set_time_unit()
 360
 361        self._parameters_matched = False
 362        
 363        if self._nodf:
 364            self.stim_conditions = np.unique(
 365                self.stim_data['stim_condition_ids']
 366            )
 367        else:
 368            self.stim_conditions = self.stim_data[
 369                'stim_condition_ids'
 370            ].unique()
 371        
 372        if self._logfile != None:
 373            self._get_conditions_from_log()
 374        else:
 375            r = r'([A-Z]{1,2}_M\d{1,}_Section_\d{1,}_BLK\d{1,})'
 376            experiment_id = re.search(r, spikefile).group(1)
 377            warn_text = f"""
 378            No stimulus log file found for experiment: {experiment_id}.
 379            If this experiment has more that 255 stimulus conditions, the stimulus 
 380            condition IDs will be incorrect in the self.stim_data attribute.
 381            """
 382            warnings.warn(warn_text, stacklevel = 4)
 383
 384    # Generates the .stim_data attribute.
 385    def _init_stim_data(self):
 386        
 387        if self._stim_m73:
 388            stim_data = {
 389                'stim_start_indicies': self.stims['stimulusOnsets'],
 390                'stim_stop_indicies': self.stims['stimulusOffsets'],
 391                'stim_condition_ids': self.stims['conditionNumbers']
 392            }
 393        else:
 394            stims_starts = self.stims[0]
 395            stims_stops = self.stims[1]
 396            stims_condition_ids = self.stims[2]
 397            stim_data = {
 398                'stim_start_indicies': stims_starts[0],
 399                'stim_stop_indicies': stims_stops[0],
 400                'stim_condition_ids': stims_condition_ids[0]
 401            }
 402        
 403        try:
 404            self.stim_data = pd.DataFrame(stim_data)
 405        except ValueError as e:
 406            if str(e) == "All arrays must be of the same length":
 407                warn_text = f"""
 408                Stimulus start time, stop time, and condition 
 409                arrays are unequal in length. This typically
 410                only happens with natural image recordings.
 411                If you include a stimulus log file in the data
 412                folder, ephystoolkit will attempt to fix this.
 413                Otherwise, the .stim_data attribute will return
 414                a dictionary instead of a dataframe.
 415                """
 416                warnings.warn(warn_text, stacklevel = 4)
 417                self.stim_data = stim_data
 418                self._nodf = True
 419
 420    # Generates the .spike_data attribute.
 421    def _init_spike_data(self):
 422        if self._spike_m73:
 423            ci = [int(i) for i in self.spikes['ChannelID']]
 424            si = self.spikes['SpikePoints']
 425            st = self.spikes['SpikeTimes']
 426            self.spike_data = [
 427                {
 428                    'channel_id': ci[i],
 429                    'spike_index': si[i],
 430                    'spike_time': st[i]
 431                }
 432                for i in range(len(st))
 433            ]
 434        else:
 435            
 436            self.spike_data = [
 437                {
 438                    'channel_id': unit[1][0][0],
 439                    'depth': self._depth_data.loc[self._depth_data.channel == unit[1][0][0]]['distance'].values[0],
 440                    'spike_index': unit[2][0],
 441                    'spike_time': unit[3][0]
 442                }
 443                for unit in
 444                self.spikes
 445            ]
 446        """
 447        A dictionary object with the spiking data.
 448        """
 449        
 450    def _get_conditions_from_log(self):
 451        
 452        with open(self._logfile, 'r') as f:
 453            stimlog = f.readlines()
 454
 455        stimlog = stimlog[1:-2]
 456        
 457        # extract conditions for log
 458        r0 = r'[Cc]ondition#: *(\d{1,})'
 459        r1 = r'[Cc]onditions:(\d{1,} +...*)'
 460        condition_ids = []
 461        
 462        if re.search(r0, stimlog[1]): # if the experiment is gratings or checkerboard
 463            for line in stimlog:
 464                c_id = re.search(r0, line).group(1)
 465                condition_ids.append(int(c_id))
 466                
 467            try: 
 468                # replace conditions in experiment with conditions from log
 469                self.stim_data.stim_condition_ids = condition_ids
 470                self.stim_conditions = self.stim_data.stim_condition_ids.unique()
 471
 472            except ValueError:
 473                pass
 474            
 475        elif re.search(r1, stimlog[1]): # if the experiment is natural images
 476            condition_ids = []
 477            for x in stimlog:
 478                cindex = x.index(':')+1 # index where condition numbers begin
 479                numbers = x[cindex:].strip() # remove whitespace at the end
 480                numbers = numbers.replace('  ', ' ') # remove double spaces
 481                numbers = numbers.split(' ') # split into a list of conditon numbers
 482                numbers = [int(x) for x in numbers] # format str into int
 483                condition_ids += numbers #collect condition ids
 484            
 485            # Attempt to fix unequal stimulus arrays/missing values
 486            if self._nodf:
 487                stim_start_indicies = self.stim_data['stim_start_indicies']
 488                stim_stop_indicies = self.stim_data['stim_stop_indicies']
 489                stim_start_times = self.stim_data['stim_start_times']
 490                stim_stop_times = self.stim_data['stim_stop_times']
 491            else:      
 492                stim_start_indicies = self.stim_data['stim_start_indicies'].values
 493                stim_stop_indicies = self.stim_data['stim_stop_indicies'].values
 494                stim_start_times = self.stim_data['stim_start_times'].values
 495                stim_stop_times = self.stim_data['stim_stop_times'].values
 496
 497            stim_data_fixed = {
 498                'stim_condition_ids': condition_ids,
 499                'stim_start_indicies': stim_start_indicies,
 500                'stim_stop_indicies': stim_stop_indicies,
 501                'stim_start_times': stim_start_times,
 502                'stim_stop_times': stim_stop_times
 503            }
 504            keys = list(stim_data_fixed.keys())
 505
 506            # extrapolate the missing values
 507            for key in keys:
 508                xd = np.arange(len(stim_data_fixed[key]))+1
 509                yd = stim_data_fixed[key]
 510                x = np.arange(len(yd), len(condition_ids))+1
 511                y_pred = self._linear_extrapolation(xd, yd, x)
 512                if 'indicies' in key:
 513                    y_pred = np.round(y_pred)
 514                stim_data_fixed[key] = np.concatenate((yd, y_pred))
 515            
 516            self.stim_data = pd.DataFrame(stim_data_fixed)
 517            self.stim_conditions = self.stim_data.stim_condition_ids.unique()
 518            
 519    def set_time_unit(self, bin_size=0.001):
 520        """
 521        Change the time unit of the relative spike times.
 522        Give a bin size relative to 1 second.
 523        IE: If you want a 1 ms bin size, enter 0.001;
 524        if you want a 10 ms bin size, enter 0.01 etc.
 525        
 526        Args:
 527        
 528        - bin_size: Time unit given relative to 1 second. The default unit is 1ms/0.001s.
 529        """
 530        if self._nodf:
 531            self.stim_data['stim_start_times'] = self._bin_events(bin_size,
 532                self.stim_data['stim_start_indicies'])
 533
 534            self.stim_data['stim_stop_times'] = self._bin_events(bin_size,
 535                self.stim_data['stim_stop_indicies'])  
 536        else:
 537            self.stim_data['stim_start_times'] = self._bin_events(bin_size,
 538                self.stim_data['stim_start_indicies'].values)
 539
 540            self.stim_data['stim_stop_times'] = self._bin_events(bin_size,
 541                self.stim_data['stim_stop_indicies'].values)
 542
 543            self.stim_data = self.stim_data[
 544                ['stim_condition_ids',
 545                 'stim_start_indicies',
 546                 'stim_stop_indicies',
 547                 'stim_start_times',
 548                 'stim_stop_times']
 549            ]
 550
 551        for i, unit in enumerate(self.spike_data):
 552            try:
 553                unit.update({
 554                    'rel_spike_time': self._bin_events(bin_size,
 555                                                       unit['spike_index'])
 556                })
 557            except TypeError as e:
 558                if 'NoneType' in str(e):
 559                    print(f"No spikes recorded for unit {i}, skipping...")
 560                else:
 561                    raise TypeError
 562
 563    def condition_times(self, condition):
 564        """
 565        Returns a dictionary object of the start times and
 566        stop times of a particular stimulus condition.
 567        
 568        Args:
 569        
 570        -condition: Condition id for the chosen stimulus condition.
 571        """
 572
 573        condition_starts = self.stim_data.groupby(
 574            self.stim_data.stim_condition_ids
 575        ).get_group(condition)['stim_start_times'].values
 576
 577        condition_stops = self.stim_data.groupby(
 578            self.stim_data.stim_condition_ids
 579        ).get_group(condition)['stim_stop_times'].values
 580
 581        return {
 582            'start': condition_starts,
 583            'stop': condition_stops
 584        }
 585
 586    def match_condition_parameters(self, params_file):
 587        """
 588        Takes a parameters file and alters the .stim_data dataframe
 589        to match stimulus condition parameters to their corresponding 
 590        condition id.
 591        
 592        Args:
 593        
 594        - params_file: Path to the stimulus parameters file.
 595        """
 596        
 597        # if parameters are given in a log file
 598        if os.path.splitext(params_file)[1] == '.log':
 599            df = self.stim_data
 600            with open(params_file, 'r') as f:
 601                lines = f.readlines()[1:-2]
 602            
 603            # read the lines to make a series of columns to insert
 604            insert = []
 605            for line in lines:
 606                items = re.split(r'\t+', line)[2:-1]
 607                dic = {
 608                    item.split(':')[0]:float(item.split(':')[1].strip(' ')) 
 609                    for item in items
 610                }
 611                insert.append(dic)
 612            
 613            # add the parameters to the stim_data dataframe
 614            insert = pd.DataFrame(insert)
 615            newcol = [list(df.columns)[0]] + list(insert.columns) + list(df.columns)[1:]
 616            self.stim_data = df.join(insert)[newcol]
 617            
 618            # make the parameter map
 619            stop_index = list(self.stim_data.columns).index('stim_start_indicies')
 620            df_new = self.stim_data[list(self.stim_data.columns)[:stop_index]]
 621            df_new.columns = ['condition']+list(df_new.columns)[1:]
 622            df_new = df_new.drop_duplicates().sort_values(by = 'condition')
 623            self.parameter_map = df_new.reset_index().drop('index', axis = 1)
 624            
 625            self._parameters_matched = True
 626            
 627        # if parameters are given in a matlab file        
 628        elif os.path.splitext(params_file)[1] == '.mat':
 629            
 630            # load the parameter file and extract the keys
 631            params = scipy.io.loadmat(params_file)
 632            param_keys = list(params.keys())
 633
 634            # regex to get the parameter names + values
 635            name_identifier = r'Var\d{1,}Val'
 636            val_identifier = r'var\d{1,}value'
 637
 638            if 'fnames' in param_keys: # if experiment is natural images
 639
 640                # get the stim file names
 641                stim_file_names = [
 642                    name[0][0] for 
 643                    name in scipy.io.loadmat(params_file)['fnames']
 644                ]
 645
 646                # make a dictionary for the parameter mapping
 647                ni_filemap = {
 648                    'condition': [],
 649                    'size': [],
 650                    'filter': [],
 651                    'file': []
 652                }
 653
 654                # match condition to parameters
 655                freg = r'^[A-Z]{2,3}'
 656                for i in self.stim_conditions:
 657                    if i%2 == 0: # if condition is even
 658
 659                        #append condition & size
 660                        ni_filemap['condition'].append(i)
 661                        ni_filemap['size'].append(60)
 662
 663                        #append filename
 664                        filename = stim_file_names[int(i/2)-1]
 665                        ni_filemap['file'].append(filename) 
 666
 667                        #append filtering condition
 668                        filt = re.search(freg, filename).group(0)
 669                        ni_filemap['filter'].append(filt)
 670
 671                    else: # if condition is odd
 672                        #append condition & size
 673                        ni_filemap['condition'].append(i)
 674                        ni_filemap['size'].append(30)
 675
 676                        #append filename
 677                        filename = stim_file_names[int((i+1)/2)-1]
 678                        ni_filemap['file'].append(filename)
 679
 680                        #append filtering condition
 681                        filt = re.search(freg, filename).group(0)
 682                        ni_filemap['filter'].append(filt)
 683
 684                self.parameter_map = pd.DataFrame(ni_filemap) 
 685                df = self.stim_data
 686
 687                # get the parameters to insert into the original dataframe
 688                insert = []
 689                for cond in df.stim_condition_ids.values:
 690                    arr = self.parameter_map.loc[self.parameter_map.condition == cond].values[0]
 691                    insert.append(arr)
 692                insert = pd.DataFrame(insert, columns=self.parameter_map.columns.values)
 693
 694                # insert the parameter values into the original dataframe
 695                df1 = df.join(insert)
 696
 697                # reset the stim_data attribute
 698                self.stim_data = df1[
 699                    list([df.columns.values[0]])
 700                    + list(self.parameter_map.columns.values[1:])
 701                    + list(
 702                        df.columns.values[1:]
 703                    )
 704                    ]
 705
 706            else: # if experiment is gratings or checkerboard
 707                # use regex to get the param names and values
 708                name_keys = [key for key in param_keys if re.search(name_identifier, key)]
 709                val_keys = [key for key in param_keys if re.search(val_identifier, key)]
 710
 711                # get the condition numbers
 712                condition_range = len(params[val_keys[0]][0])
 713                condition_ids = [i + 1 for i in range(condition_range)]
 714
 715                # map conditions to parameter values
 716                parameter_map = {'condition': condition_ids}
 717                for i in range(len(name_keys)):
 718                    parameter_name = str(params[name_keys[i]][0])
 719                    parameter_values = np.array(params[val_keys[i]][0])
 720
 721                    parameter_map[parameter_name] = parameter_values
 722
 723                    # parameter dataframe + the original stim_data dataframe
 724                self.parameter_map = pd.DataFrame(parameter_map)
 725                df = self.stim_data
 726
 727                # get the parameters to insert into the original dataframe
 728                insert = []
 729                for cond in df.stim_condition_ids.values:
 730                    arr = self.parameter_map.loc[self.parameter_map.condition == cond].values[0]
 731                    insert.append(arr)
 732                insert = pd.DataFrame(insert, columns=self.parameter_map.columns.values)
 733
 734                # insert the parameter values into the original dataframe
 735                df1 = df.join(insert)
 736
 737                # reset the stim_data attribute
 738                self.stim_data = df1[
 739                    list([df.columns.values[0]])
 740                    + list(self.parameter_map.columns.values[1:])
 741                    + list(
 742                        df.columns.values[1:]
 743                    )
 744                    ]
 745
 746            self._parameters_matched = True
 747
 748    def _pr_unitcols(
 749            self,
 750            include_units,
 751            stim_condition=None, 
 752            columns='cluster_id',
 753            thresh=None,
 754            norm=None 
 755    ):
 756
 757        if thresh is None:
 758            return "Please provide a response boundary threshhold as a tuple (ex: (-50,500))"
 759        else:
 760            thresh_min = thresh[0] - 0.5
 761            thresh_max = thresh[1] + 0.5
 762
 763        if stim_condition is 'all':
 764            stim_condition = self.stim_conditions
 765
 766        population = {}
 767        cond_col = []
 768        param_col = []
 769
 770        for condition in stim_condition:
 771            condition_start_times = self.condition_times(condition)['start']
 772
 773            for i, cluster_id in enumerate(include_units):
 774                if type(cluster_id) == np.ndarray:
 775                    unit_spike_times = cluster_id
 776                    cluster_id = f"n{i}"
 777                    
 778                else:
 779                    unit_spike_times = self.spike_data[cluster_id]['rel_spike_time']
 780
 781                # Raster
 782                x = self.raster(
 783                    condition_start_times, 
 784                    unit_spike_times,
 785                    thresh=thresh)
 786
 787                # Raster to histogram 
 788                h, bins = np.histogram(
 789                    x, 
 790                    bins=np.arange(thresh_min,thresh_max,1)
 791                )
 792
 793                # Nonresponsive unit histogram
 794                if h.min() == 0 and h.max() == 0:
 795                    h = np.zeros(len(h))
 796
 797                # Normalize the response
 798                else:
 799                    pass
 800                    norm_method = {
 801                        'minmax': self._minmax_norm(h),
 802                        'zscore': self._zscore_norm(h),
 803                        'average': self._average_response(h, len(condition_start_times))
 804                    }
 805
 806                    if norm is None: continue
 807
 808                    elif (norm is not None
 809                          and norm in list(norm_method.keys())):
 810                        h = norm_method[norm]
 811                    else:
 812                        raise _UnrecognizedNorm(
 813                            f"Normalization method is not recognized,\
 814                            please choose from the following:\
 815                            {list(norm_method.keys())}"
 816                        )
 817
 818                # fill the population dictionary
 819                if cluster_id not in population:
 820                    population[cluster_id] = h
 821                else:
 822                    population[cluster_id] = np.concatenate(
 823                        (population[cluster_id], h)
 824                    )
 825
 826            # append condition ids to the conditions column
 827            cond_col += [int(condition)]*len(h)
 828
 829            # append condition parameters
 830            if self._parameters_matched:
 831                params = list(
 832                    self.parameter_map.loc[self.parameter_map.condition == condition].values[0][1:]
 833                )
 834                param_col += [params] * len(h)
 835
 836        # Rearrange the dataframe
 837        if self._parameters_matched:
 838
 839            # Make the population & stimulus parameters dataframes
 840            population = pd.DataFrame(population)
 841            param_col = pd.DataFrame(
 842                param_col, columns=self.parameter_map.columns.values[1:]
 843            )
 844
 845            # Retrieve their column labels
 846            pcol = list(population.columns.values)
 847            pacol = list(param_col.columns.values)
 848
 849            # Add the condition id column & the stimulus parameters columns
 850            population['stimulus_condition'] = cond_col
 851            population = population.join(param_col)
 852
 853            # Rearrange column order
 854            new_columns = ['stimulus_condition'] + pacol + pcol
 855            population = population[new_columns]
 856
 857        else:
 858            # Make the population dataframe
 859            population = pd.DataFrame(population)
 860
 861            # Retrieve its column labels
 862            pcol = list(population.columns.values)
 863
 864            # Add the condition id column
 865            population['stimulus_condition'] = cond_col
 866
 867            # Rearrange column order
 868            population = population[['stimulus_condition']+pcol]
 869
 870        return population       
 871
 872    def _pr_stimcols(
 873            self,
 874            include_units, 
 875            stim_condition=None, 
 876            columns='cluster_id',
 877            thresh=None, 
 878            norm=None
 879    ):
 880
 881        if thresh is None:
 882            return "Please provide a response boundary threshhold as a tuple (ex: (-50,500))"
 883        else:
 884            thresh_min = thresh[0] - 0.5
 885            thresh_max = thresh[1] + 0.5
 886
 887        if stim_condition is 'all':
 888            stim_condition = np.sort(self.stim_conditions)
 889
 890        population = {}
 891        unit_col = []
 892
 893        for cluster_id in include_units:
 894            unit_spike_times = self.spike_data[cluster_id]['rel_spike_time']
 895
 896            for condition in stim_condition:
 897                condition_start_times = self.condition_times(condition)['start']
 898
 899                # Raster
 900                x = self.raster(
 901                    condition_start_times, 
 902                    unit_spike_times,
 903                    thresh=thresh)
 904
 905                # Raster to histogram 
 906                h, bins = np.histogram(
 907                    x, 
 908                    bins=np.arange(thresh_min,thresh_max,1)
 909                )
 910
 911                # Nonresponsive unit histogram
 912                if h.min() == 0 and h.max() == 0:
 913                    h = np.zeros(len(h))
 914
 915                # Normalize the response
 916                else:
 917                    pass
 918                    norm_method = {
 919                        'minmax': self._minmax_norm(h),
 920                        'zscore': self._zscore_norm(h),
 921                        'average': self._average_response(h, len(condition_start_times))
 922                    }
 923
 924                    if norm is None: continue
 925
 926                    elif (norm is not None
 927                          and norm in list(norm_method.keys())):
 928                        h = norm_method[norm]
 929                    else:
 930                        raise _UnrecognizedNorm(
 931                            f"Normalization method is not recognized,\
 932                            please choose from the following:\
 933                            {list(norm_method.keys())}"
 934                        )
 935
 936                # fill the population dictionary
 937                if condition not in population:
 938                    population[int(condition)] = h
 939                else:
 940                    population[condition] = np.concatenate(
 941                        (population[int(condition)], h)
 942                    )
 943
 944            # append condition ids to the conditions column
 945            unit_col += [cluster_id]*len(h)
 946
 947        # Rearrange the dataframe
 948        population = pd.DataFrame(population)
 949        pcol = list(population.columns.values)
 950        population['cluster_id'] = unit_col
 951        population = population[['cluster_id']+pcol]
 952
 953        return population       
 954
 955    def population_response(
 956            self,
 957            include_units,  # Units to include in the dataframe
 958            stim_condition=None,  # Stimulus condition(s) to include in the dataframe
 959            columns='cluster_id',  # Set column label arrangement
 960            thresh=None,  # Bounding threshold around the stimulus onset at t = 0 - pass as a tuple
 961            norm=None  # Normalization method - choose from: 'minmax', 'zscore', or 'average'
 962    ):
 963        """
 964        Returns a dataframe of the population response PSTH.
 965        By default, each column label represents the
 966        included units. A single column identifies
 967        the stimulus condition at each row. Each row is the
 968        average response at each time step.
 969
 970        Setting columns = "stimulus_condition" will
 971        return a data frame where each column label
 972        represents a stimulus condition. A single
 973        column identifes the included unit at each row.
 974        Each row is the average response at each time step.
 975
 976        Args:
 977
 978        - include_units: Units to include in the dataframe - pass as a 1d array or list like object.
 979        - stim_condition: Stimulus condition(s) to include in the dataframe - 
 980          pass a list of condition ids or 'all' to include all conditions.
 981        - columns: Set column label arrangement - pass either 'cluster_id' or 'stimulus_condition'. 
 982          Default argument is 'stimulus_condition'.
 983        - thresh: Bounding threshold around the stimulus onset at t = 0 - pass as a tuple.
 984        - norm: Normalization method - pass either 'minmax', 'zscore', or 'average'.
 985        """
 986
 987        # Case checks:
 988        if stim_condition is None:
 989            raise _NoStimulusCondition()
 990        elif stim_condition is 'all': pass
 991        elif set(stim_condition).issubset(set(self.stim_conditions)): pass
 992        else:
 993            raise _UnrecognizedStimulusCondition()
 994
 995        col_formats = ['cluster_id', 'stimulus_condition'] 
 996        if columns not in col_formats:
 997            raise _UnrecognizedColumnsInput(list(columns_dic.keys()))
 998
 999        elif columns == 'cluster_id':
1000            population = self._pr_unitcols(
1001            include_units,
1002            stim_condition=stim_condition,
1003            columns=columns,
1004            thresh=thresh,
1005            norm=norm)
1006
1007        elif columns == 'stimulus_condition':
1008            population = self._pr_stimcols(
1009            include_units,
1010            stim_condition=stim_condition,
1011            columns=columns,
1012            thresh=thresh,
1013            norm=norm)
1014
1015        return population
1016    
1017    def map_rf(self, rfunc, xpix, ypix, xvis, yvis, radius):
1018        """
1019        Map the spatiotemporal receptive field of a unit with 
1020        a response function across a stimulus space. 
1021        
1022        Args:
1023        
1024        - rfunc: The unit response function. Dim0 must be equal to
1025          the total number of stimulus conditions.  
1026        - xpix: The width of the image in pixels.
1027        - ypix: The height of the image in pixels.
1028        - xvis: The azimuth of the image in degrees.
1029        - yvis: The elevation of the image in degrees.
1030        - radius: The radius of the circular frame around the image in degrees.
1031        """
1032        
1033        
1034        self._stim_frame_map(xpix, ypix, xvis, yvis, radius)
1035        fmap = np.array([x for x in self.frame_map.values()])
1036        if fmap.T.shape[-1] == rfunc.shape[0]:
1037            return np.dot(fmap.T, rfunc)
1038        else:
1039            raise _InvalidRfuncShape((fmap.T.shape[-1], rfunc.shape[0]))
1040                    
1041    def _stim_frame_map(self, xpix, ypix, xvis, yvis, radius):
1042        frame_map = {}
1043        for i in range(len(self.parameter_map)):
1044            con = i+1
1045            dx = self.parameter_map.iloc[i]
1046            
1047            frame = self.static_grating(
1048                xpix, ypix,
1049                xvis, yvis,
1050                float(dx['Spatial Freq']), 
1051                float(dx['Orientation']), 
1052                float(dx['Phase'])*360,
1053                diameter = radius
1054                )
1055
1056            frame_map[con] = frame
1057
1058        self.frame_map = frame_map
1059
1060    def _stim_frames(self):
1061        cond = self.stim_data.stim_condition_ids.values
1062        self.stim_frames = np.array([self.frame_map[i] for i in cond])
1063
1064class lfp_tools(ephys_toolkit):
1065
1066    #### initialize class with the necessary paths ###
1067    def __init__(self):
1068        ephys_toolkit.__init__(self)
1069        self._low = 1
1070        self._high = 120
1071        
1072    ### process the lfp data by running the necssary methods ###
1073    def process_lfp(self, intan_file, probe = None):
1074        """
1075        This method runs the lfp processing pipeline and
1076        generates two important attributes:
1077        
1078        - .lfp_heatmaps: contains a dictionary the heatmaps of the 
1079          lfp data for each channel column and each stimulus contrast.
1080          
1081        - .depth_data: contains the depth in micrometeres of each 
1082          channel normalized to the center of layer 4 as 0. 
1083          Negative values are ventral to layer 4 and positive values
1084          are dorsal to layer 4.
1085        
1086        """
1087        self._probe_fol = os.path.join(self._modpath, 'probes')
1088        valid_probes0 = os.listdir(self._probe_fol)
1089        valid_probes = [os.path.splitext(x)[0] for x in valid_probes0]
1090        if probe == None:
1091            inputtext = fr"""
1092            Please enter the probe used in this project. Available probes include: 
1093            {valid_probes}. 
1094            If the probe used in your recording is not included in the list, 
1095            navigate to the folder where the ephystoolkit lirbary is installed. 
1096            
1097            If you are using Anaconda on windows, this folder should be located at:
1098            C:\Users\(username)\conda\envs\(env_name)\lib\site-packages\ephystoolkit
1099            
1100            If you are using Anaconda on Linux, this folder should be located at:
1101            ~/anaconda3/envs/(env_name)/lib/(Python_version)/site-packages/ephystoolkit'
1102            
1103            Once you are in the module folder, navigate the to folder called 'probes'
1104            and copy a csv file mapping your channel ids to their distance from the tip
1105            of the probe. Each line index should correspond to the channel id. Indexes
1106            for this purpose start at 1. For reference, check out any of the included
1107            csv files."""
1108            
1109            probe = input(inputtext)
1110        while probe not in valid_probes:
1111            inputtext = fr"""
1112            {probe} is not included in the list of available probes. The list of available 
1113            probes includes: 
1114            {valid_probes}
1115            If the probe used in your recording is not included in the list, 
1116            navigate to the folder in which the ephys_toolkit module is installed. 
1117
1118            If you are using Anaconda on windows, this folder should be located at:
1119            C:\Users\(username)\conda\envs\(env_name)\lib\site-packages\ephystoolkit
1120
1121            If you are using Anaconda on Linux, this folder should be located at:
1122            ~/anaconda3/envs/(env_name)/lib/(Python_version)/site-packages/ephystoolkit
1123
1124            Once you are in the module folder, navigate the to folder called 'probes'
1125            and copy a csv file mapping your channel ids to their distance from the tip
1126            of the probe. Each line index should correspond to the channel id. Indexes
1127            for this purpose start at 1. For reference, check out any of the included
1128            csv files.
1129            """
1130            probe = input(inputtext)
1131        else:
1132            pass
1133        
1134        self._probe = probe
1135        self.intan_file = intan_file
1136        self.channel_file = os.path.join(self._probe_fol, f"{probe}.csv")
1137        
1138        self._load_lfp_data()
1139        
1140        print("Sorting channel depth...")
1141        self._sort_depth()
1142        
1143        print("Retrieving stimulus onset times...")
1144        self._get_onsets()
1145        
1146        print("Generating lfp heatmap...")
1147        self._get_lfp_heatmap()
1148        
1149        print("Normalizing channel depth to layer 4 depth...")
1150        self._find_l4_distances()
1151        self._normalize_depth()
1152        print("Done!")
1153    
1154    ### bandpass filter for the raw voltage data ###
1155    def _bandpass_filter(self, signal):
1156        
1157        #SAMPLING_RATE = 20000.0
1158        nyqs = 0.5 * self._SAMPLING_RATE
1159        low = self._low/nyqs
1160        high = self._high/nyqs
1161
1162        order = 2
1163
1164        b, a = scipy.signal.butter(order, [low, high], 'bandpass', analog = False)
1165        y = scipy.signal.filtfilt(b, a, signal, axis=0)
1166
1167        return y
1168
1169    ### load the necessary data ###
1170    def _load_lfp_data(self):
1171        
1172        # load the rhd file
1173        self.intan_data = data = rhd.read_data(self.intan_file)
1174        self.amplifier_data = self.intan_data['amplifier_data']
1175        self.board_data = self.intan_data['board_dig_in_data']
1176
1177        # load the channel file
1178        self.channel_data = []
1179        with open(self.channel_file, newline='') as channels:
1180            for row in csv.reader(channels):
1181                self.channel_data.append(np.array(row).astype(float))
1182        self.channel_data = np.array(self.channel_data)
1183
1184        # bandpass filter each channel
1185        print("Bandpass filtering the data...")
1186        self.amplifier_data_filt = np.array([
1187            self._bandpass_filter(channel) for
1188            channel in self.amplifier_data
1189        ]
1190        )
1191        
1192    ### sort amplifier data by channel depth ###
1193    def _sort_depth(self):
1194        ## 0 = deepest
1195        
1196        # adjust single channels offset on x by a negligable amount
1197        if self._probe == '64D':
1198            self.channel_data[np.where(self.channel_data == 16)] = 20
1199            self.channel_data[np.where(self.channel_data == -16)] = -20
1200        if self._probe == '64H':
1201            self.channel_data[np.where(self.channel_data == 20)] = 22.5
1202            self.channel_data[np.where(self.channel_data == -20)] = -22.5
1203            self.channel_data[np.where(self.channel_data == 180.1)] = 177.6
1204            self.channel_data[np.where(self.channel_data == 220.1)] = 222.6
1205        
1206        # index ordered by distance from tip
1207        self.y_index = np.argsort(self.channel_data[:,1])[::-1]
1208        
1209        # list of channel depths ordered by distance from tip
1210        self.c_depth = self.channel_data[self.y_index][:,1]
1211
1212        # unique x coordinates for the channels
1213        xarr = np.unique(self.channel_data[:,0])
1214        
1215        # get channel depths along specific columns
1216        self.channels_bycol = {}
1217        for x_d in xarr: 
1218            col = np.where(self.channel_data[self.y_index][:,0] == x_d)[0]
1219            self.channels_bycol[x_d] = {
1220                'y_index': self.y_index[col],
1221                'c_depth': self.channel_data[self.y_index[col]][:,1]
1222            }
1223        
1224        # get amplifier data along specific columns
1225        self.ampdata_bycol = {}    
1226        for col, subdict in self.channels_bycol.items():
1227            self.ampdata_bycol[col] = self.amplifier_data_filt[
1228                subdict['y_index']
1229            ]
1230
1231    
1232    ### Get the onset times of each stimuls contrast ###
1233    def _get_onsets(self):
1234        
1235        # get indices of when contrast 0 is on
1236        stim0 = self.board_data[1]
1237        stim0_on_index = (np.where(stim0 == True)[0])
1238
1239        # get indices of when contrast 1 is on
1240        stim1 = self.board_data[2]
1241        stim1_on_index = (np.where(stim1 == True)[0])
1242
1243        # get the onset indices of contrast 0
1244        first0_on = stim0_on_index[0]
1245        self.c0_onsets = [first0_on]
1246        for i0, val in enumerate(stim0_on_index[1:]):
1247            i1 = i0+1
1248            if stim0_on_index[i0]+1 != (stim0_on_index[i1]):
1249                self.c0_onsets.append(val)
1250
1251        # get the onset indices of contrast 1
1252        first1_on = stim1_on_index[0]
1253        self.c1_onsets = [first1_on]
1254        for i0, val in enumerate(stim1_on_index[1:]):
1255            i1 = i0+1
1256            if stim1_on_index[i0]+1 != (stim1_on_index[i1]):
1257                self.c1_onsets.append(val)
1258                
1259    ### get the heatmap of lfp data averaged across stimulus repeats ###
1260    def _get_lfp_heatmap(self):
1261
1262        self.lfp_heatmaps = {}
1263        for col, subdict in self.ampdata_bycol.items():
1264            self.lfp_heatmaps[col] = {
1265                'contrast0': np.array([
1266                    subdict[:,i:i+5000] for i in self.c0_onsets
1267                ]).mean(0),
1268                'contrast1': np.array([
1269                    subdict[:,i:i+5000] for i in self.c1_onsets
1270                ]).mean(0),
1271
1272            }
1273        
1274    ### assign a depth to each channel relative to layer 4 ###
1275    def _find_l4_distances(self):
1276        
1277        self.l4_distances = {}
1278        self.l4_distance_list = []
1279        for col, chandata in self.channels_bycol.items():
1280            lfp_heatmap = self.lfp_heatmaps[col]
1281            self.l4_distances[col] = {}
1282            for c, contrast in lfp_heatmap.items():
1283            
1284                # find the timepoint of the deepest sink
1285                pe_delay = 1000
1286                tmin = np.where(contrast[:,pe_delay:] == contrast[:,pe_delay:].min())[1][0]
1287                tmin_curve = contrast[:,tmin+pe_delay]
1288
1289                ## Find the center of layer 4 via cubic spline interpolation
1290
1291                # flip the order because the interpolation function
1292                # strictly requires an x axis with ascending order
1293                x = chandata['c_depth'][::-1]
1294                y = tmin_curve
1295
1296                # initialize the interpolation function
1297                cs = scipy.interpolate.CubicSpline(x, y)
1298
1299                # interpolate the data and flip it back to the correct order
1300                xnew = np.linspace(0, x[-1], 1575)[::-1]
1301                ynew = cs(xnew)[::-1]
1302
1303                # get the layer 4 distance from tip
1304                l4_distance = xnew[np.where(ynew == ynew.min())[0][0]].round(2)
1305                self.l4_distances[col][c] = l4_distance
1306                self.l4_distance_list.append(l4_distance)
1307        self.l4_distance_list = np.array(self.l4_distance_list)
1308    
1309    def _normalize_depth(self):
1310        
1311        # find the average distance across channel
1312        # columns and stimulus contrasts
1313        avg_distance = self.l4_distance_list.mean()
1314
1315        # normalize channel distances by distance from layer 4
1316        # negative distances = below layer 4
1317        # positive distances = above layer 4
1318        l4_normalization = self.c_depth - avg_distance
1319
1320        # compile the depth data into a dataframe
1321        depth_data = {
1322            'channel': self.y_index+1,
1323            'distance': l4_normalization
1324        }
1325        self._depth_data = pd.DataFrame(depth_data)        
1326
1327class load_project(lfp_tools, load_experiment):
1328    """
1329    Initialize the load_project class with a full path to the
1330    directory containing the project files. If the LFP data rhd
1331    file is included in the project directory, this class will 
1332    automatically process the LFP data for the project. 
1333    
1334    The .workbook attribute contains a list of dictionaries
1335    with the following structure:
1336    
1337
1338      
1339          {
1340          
1341              'section_id': int - Identification number of the recording 
1342               section,
1343              
1344              'spike_sorting_metrics': Pandas dataframe - Contains spike
1345               sorting metrics data,
1346              
1347              'lfp_heatmaps': dictionary - Contains the lfp heatmap
1348               at each channel column and stimulus contrast,
1349              
1350              'depth_data': Pandas dataframe - Contains the distance of each
1351               channel of the probe relative to layer 4,
1352              
1353              'blocks': [
1354              
1355                  {
1356                  
1357                  'block_id': int - Identification number of the recording block, 
1358                  
1359                  'experiment', experiment object - Contains the experiment data for
1360                   the given block
1361                  
1362                  },
1363                  
1364                  ]
1365                  
1366          },
1367    
1368    Args:
1369    
1370    - project_path: Path to the directory containing the project files.
1371    """
1372
1373    def __init__(self, project_path, probe = None, use_lfp_file = 0):
1374        lfp_tools.__init__(self)
1375        self._lfp_index = use_lfp_file
1376        self._probe = probe
1377        self._ppath = project_path
1378        self._init_project_workbook()
1379
1380    # generate the workbook of project data
1381    def _init_project_workbook(self):
1382        explorer = path_explorer()
1383
1384        # find and sort spike files
1385        spike_files = explorer.findext(self._ppath, '.mat', r='firings')
1386        spike_files.sort(key=lambda x: int(re.search(r'BLK(\d{1,})', x).group(1)))
1387
1388        # find and sort stim files
1389        stim_files = explorer.findext(self._ppath, '.mat', r='stimulusData')
1390        stim_files.sort(key=lambda x: int(re.search(r'BLK(\d{1,})', x).group(1)))
1391        
1392        # find and sort log files
1393        log_files = explorer.findext(self._ppath, '.log')
1394        log_files.sort(key=lambda x: int(re.search(r'BLK(\d{1,})', x).group(1)))
1395        
1396        # count how many log files are in the data folder
1397        log_diff = len(spike_files)-len(log_files)
1398        for i in range(log_diff):
1399            log_files.append(None)
1400
1401        # zip matching blocks
1402        matched_block_files = zip(spike_files, stim_files, log_files)
1403
1404        # find metrics files
1405        metrics_files = explorer.findext(self._ppath, '.json', r='metrics_isolation')
1406        metrics_files.sort(key=lambda x: int(re.search(r'Section_(\d{1,})', x).group(1)))
1407        
1408        # find checkerboard rhd files
1409        rhd_files = explorer.findext(self._ppath, '.rhd')
1410        rhd_files.sort(key=lambda x: int(re.search(r'Section_(\d{1,})', x).group(1)))
1411        
1412        # group together multiple checkerboard recordings
1413        # from the same section
1414        i = 0
1415        for rhd_file in rhd_files[1:]:
1416            if type(rhd_files[i])!=list:
1417                r = r"[A-Z][A-Z]_M\d{1,}_Section_(\d)"
1418                id_sec_post = re.search(r, rhd_file).group(0)
1419                id_sec_pre = re.search(r, rhd_files[i]).group(0)
1420                if id_sec_post == id_sec_pre:
1421                    rhd_files[i] = [rhd_files[i], rhd_file]
1422                    post_index = rhd_files.index(rhd_file)
1423                    del rhd_files[post_index]
1424
1425            elif (type(rhd_files[i])==list):
1426                id_sec_post = re.search(r, rhd_file).group(0)
1427                id_sec_pre = re.search(r, rhd_files[i][0]).group(0)
1428                if id_sec_post == id_sec_pre:
1429                    rhd_files[i].append(rhd_file)
1430                    post_index = rhd_files.index(rhd_file)
1431                    del rhd_files[post_index]
1432            else:
1433                i+=1
1434
1435        ################################################################################
1436
1437        # compile the workbook
1438        self.workbook = []
1439        
1440        # compile metrics
1441        for metrics_file in metrics_files:
1442            section_parent = int(re.search(r'Section_(\d{1,})', metrics_file).group(1))
1443            df = self.spike_sorting_metrics(metrics_file)
1444            self.workbook.append(
1445                {
1446                    'section_id': section_parent,
1447                    'spike_sorting_metrics': df,
1448                    'lfp_heatmaps': None,
1449                    'depth_data': None,
1450                    'blocks': []
1451                }
1452            )
1453            
1454        # compile rhd
1455        if rhd_files == []:
1456            warn_text = """
1457            No LFP data files found in this project folder. LFP data will be
1458            unavailable for this project. Please include the LFP data rhd file 
1459            in the project folder if you want to access the LFP data.
1460            """
1461            warnings.warn(warn_text, stacklevel = 4)
1462        else:
1463            for rhd_file in rhd_files:
1464                if (type(rhd_file) == list) & (self._lfp_index == 0):
1465                    warn_text = """
1466                    WARNING: More than one LFP data file found in this section.
1467                    Defaulting to the first LFP file found in the section. If 
1468                    You wish to change which file to be used, pass an integer 
1469                    value to the use_lfp_file parameter of the load_project class
1470                    corresponding to the file you wish to use. 
1471                    (ie: load_project(use_lfp_file = 2) lets you use the 2nd LFP 
1472                    data file.).
1473                    """
1474                    warnings.warn(warn_text, stacklevel = 4)
1475                    
1476                    rhd_file = rhd_file[self._lfp_index]
1477                elif (type(rhd_file) == list) & (self._lfp_index != 0):
1478                    use_lfp_file = use_lfp_file-1
1479                    rhd_file = rhd_file[self._lfp_index]
1480                else:
1481                    pass
1482                
1483                self.process_lfp(rhd_file, self._probe)
1484                section_parent = int(re.search(r'Section_(\d{1,})', rhd_file).group(1))
1485                self.workbook[section_parent-1]['lfp_heatmaps'] = self.lfp_heatmaps
1486                self.workbook[section_parent-1]['depth_data'] = self._depth_data
1487
1488        # match experiment (block) objects to section
1489        for matched_files in list(matched_block_files):
1490            # a regex to get the experiment identity
1491            ex_r = r'[A-Z]{2}_M\d+_Section_\d+_BLK\d+' 
1492            experiment_id = re.search(ex_r, matched_files[0]).group(0)
1493            
1494            section_child = int(re.search(r'Section_(\d{1,})', matched_files[0]).group(1))
1495            block = int(re.search(r'BLK(\d{1,})', matched_files[0]).group(1))
1496            
1497            experiment = load_experiment(*matched_files, kwargs = {'depth_data': self._depth_data})
1498            self.workbook[section_child - 1]['blocks'].append({
1499                'block_id': block,
1500                'experiment': experiment
1501            })
1502            print(f"Sucessfully loaded {experiment_id}.")
1503
1504
1505# Class Errors
1506class _UnrecognizedNorm(Exception):
1507    """
1508    Exception raised for unrecognized user input
1509    for array normalization.
1510    """
1511
1512    def __init__(self, message):
1513        self.message = message
1514        super().__init__(self.message)
1515
1516
1517class _NoStimulusCondition(Exception):
1518    """
1519    Exception raised for no user input in the
1520    stim_condition argument in the 
1521    self.population_response method.
1522    """
1523
1524    def __init__(self):
1525        self.message = """
1526        Please pass a stimulus condition or enter 'all' to generate 
1527        a joint condition population response matrix.
1528        """
1529        super().__init__(self.message)
1530
1531
1532class _UnrecognizedStimulusCondition(Exception):
1533    """
1534    Exception raised for invalid user input in the
1535    stim_condition argument in the 
1536    self.population_response method.
1537    """
1538
1539    def __init__(self):
1540        self.message = """
1541        The stimulus condition does not exist within this experiment.
1542        Please enter a valid stimulus condition.
1543        """
1544        super().__init__(self.message)
1545
1546
1547class _UnrecognizedColumnsInput(Exception):
1548    """
1549    Exception raised for invalid user input in the
1550    columns argument in the 
1551    self.population_response method.
1552    """
1553
1554    def __init__(self, arg):
1555        self.message = f"""
1556        Invalid input for 'columns'. Please select one of: {arg}.
1557        """
1558        super().__init__(self.message)
1559        
1560class _InvalidRfuncShape(Exception):
1561    """
1562    Exception raised for invalid shape of the response
1563    function given for receptive field mapping.
1564    """
1565
1566    def __init__(self, arg):
1567        self.message = f"""
1568        Shape mistmatch: Frame map dim - 1 ({arg[0]}) != Response function dim 0 ({arg[1]}).
1569        """
1570        super().__init__(self.message)