docs: def calculateRate (include = ['allCells', 'eachPop'], peakBin = 5, timeRange = None): 
      ''' 
      Calculate avg and peak rate of different subsets of cells for specific time period
          - include (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): List of data series to include. 
              Note: one line per item, not grouped (default: ['allCells', 'eachPop'])
          - timeRange ([start:stop]): Time range of spikes shown; if None shows all (default: None)
          - peakBin (int): Histogram bin size used to calculate peak firing rate; if None, peak rate not calculated (default: 5)
          - Returns list with rates
      '''
  
      from .. import sim
  
      print('Calculating avg and peak firing rates ...')
  
      # Replace 'eachPop' with list of pops
      if 'eachPop' in include: 
          include.remove('eachPop')
          for pop in sim.net.allPops: include.append(pop)
  
      # time range
      if timeRange is None:
          timeRange = [0,sim.cfg.duration]
  
      avg, peak, histData = [], [], []
  
      # Plot separate line for each entry in include
      for iplot,subset in enumerate(include):
          cells, cellGids, netStimLabels = getCellsInclude([subset])
          numNetStims = 0
  
          # Select cells to include
          if len(cellGids) > 0:
              try:
                  spkinds,spkts = list(zip(*[(spkgid,spkt) for spkgid,spkt in zip(sim.allSimData['spkid'],sim.allSimData['spkt']) if spkgid in cellGids]))
              except:
                  spkinds,spkts = [],[]
          else: 
              spkinds,spkts = [],[]
  
          # Add NetStim spikes
          spkts, spkinds = list(spkts), list(spkinds)
          numNetStims = 0
          if 'stims' in sim.allSimData:
              for netStimLabel in netStimLabels:
                  netStimSpks = [spk for cell,stims in sim.allSimData['stims'].items() \
                  for stimLabel,stimSpks in stims.items() for spk in stimSpks if stimLabel == netStimLabel]
                  if len(netStimSpks) > 0:
                      lastInd = max(spkinds) if len(spkinds)>0 else 0
                      spktsNew = netStimSpks 
                      spkindsNew = [lastInd+1+i for i in range(len(netStimSpks))]
                      spkts.extend(spktsNew)
                      spkinds.extend(spkindsNew)
                      numNetStims += 1
  
          if peakBin:
              histo = np.histogram(spkts, bins = np.arange(timeRange[0], timeRange[1], peakBin))
              histoT = histo[1][:-1]+peakBin/2
              histoCount = histo[0] 
  
              histData.append(histoCount)
  
              histoCount = histoCount * float((1000.0 / peakBin)) / float((len(cellGids)+numNetStims)) # convert to firing rate
              peak.append(float(max(histoCount)))
  
          spktsRange = [spkt for spkt in spkts if timeRange[0] <= spkt <= timeRange[1]]
          avg.append(float(len(spktsRange)) / float((len(cellGids)+numNetStims)) / float((timeRange[1]-timeRange[0])) * 1000.0)
  
      return include, avg, peak