docs: def plotSpikeStats (include = ['allCells', 'eachPop'], statDataIn = {}, timeRange = None, graphType='boxplot', stats = ['rate', 'isicv'], bins = 50,
                   popColors = [], histlogy = False, histlogx = False, histmin = 0.0, density = False, includeRate0=False, legendLabels = None, normfit = False,
                    histShading=True, xlim = None, dpi = 100, figSize = (6,8), fontSize=12, saveData = None, saveFig = None, showFig = True, **kwargs): 
      ''' 
      Plot spike histogram
          - 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)
          - graphType ('boxplot', 'histogram'): Type of graph to use (default: 'boxplot')
          - stats (['rate', |'isicv'| 'sync'| 'pairsync']): Measure to plot stats on (default: ['rate', 'isicv'])
          - bins (int or list of edges): Number of bins (if integer) of edges (if list) for histogram (default: 50)
          - popColors (dict): Dictionary with color (value) used for each population (key) (default: None)
          - figSize ((width, height)): Size of figure (default: (10,8))
          - saveData (None|True|'fileName'): File name where to save the final data used to generate the figure;
              if set to True uses filename from simConfig (default: None)
          - saveFig (None|True|'fileName'): File name where to save the figure;
              if set to True uses filename from simConfig (default: None)
          - showFig (True|False): Whether to show the figure or not (default: True)
  
          - Returns figure handle and statData
      '''
  
      from .. import sim
      print('Plotting spike stats...')
  
      # Set plot style
      colors = []
      params = {
          'axes.labelsize': fontSize,
          'font.size': fontSize,
          'legend.fontsize': fontSize,
          'xtick.labelsize': fontSize,
          'ytick.labelsize': fontSize,
          'text.usetex': False,
          }
      plt.rcParams.update(params)
  
      xlabels = {'rate': 'Rate (Hz)', 'isicv': 'Irregularity (ISI CV)', 'sync':  'Synchrony', ' pairsync': 'Pairwise synchrony'}
  
      # 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]
  
      for stat in stats:
          # create fig
          fig,ax1 = plt.subplots(figsize=figSize)
          fontsiz = fontSize
          xlabel = xlabels[stat]
  
          statData = []
          gidsData = []
          ynormsData = []
  
          # Calculate data for each entry in include
          for iplot,subset in enumerate(include):
  
              if stat in statDataIn:
                  statData = statDataIn[stat]['statData']
                  gidsData = statDataIn[stat].get('gidsData', [])
                  ynormsData = statDataIn[stat].get('ynormsData', [])
  
              else:
                  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
                  try:
                      spkts,spkinds = list(zip(*[(spkt, spkind) for spkt, spkind in zip(spkts, spkinds) 
                          if timeRange[0] <= spkt <= timeRange[1]]))
                  except:
                      pass
  
                  # if scatter get gids and ynorm
                  if graphType == 'scatter':    
                      if includeRate0:
                          gids = cellGids
                      else:
                          gids = set(spkinds)
                      ynorms = [sim.net.allCells[int(gid)]['tags']['ynorm'] for gid in gids]
  
                      gidsData.insert(0, gids)
                      ynormsData.insert(0, ynorms)
  
                  # rate stats
                  if stat == 'rate':
                      toRate = 1e3/(timeRange[1]-timeRange[0])
                      if includeRate0:
                          rates = [spkinds.count(gid)*toRate for gid in cellGids] \
                              if len(spkinds)>0 else [0]*len(cellGids) #cellGids] #set(spkinds)] 
                      else:
                          rates = [spkinds.count(gid)*toRate for gid in set(spkinds)] \
                              if len(spkinds)>0 else [0] #cellGids] #set(spkinds)] 
                      statData.append(rates)
  
  
                  # Inter-spike interval (ISI) coefficient of variation (CV) stats
                  elif stat == 'isicv':
                      import numpy as np
                      spkmat = [[spkt for spkind,spkt in zip(spkinds,spkts) if spkind==gid] 
                          for gid in set(spkinds)]
                      isimat = [[t - s for s, t in zip(spks, spks[1:])] for spks in spkmat if len(spks)>10]
                      isicv = [np.std(x) / np.mean(x) if len(x)>0 else 0 for x in isimat] # if len(x)>0] 
                      statData.append(isicv) 
  
  
                  # synchrony
                  elif stat in ['sync', 'pairsync']:
                      try: 
                          import pyspike  
                      except:
                          print("Error: plotSpikeStats() requires the PySpike python package \
                              to calculate synchrony (try: pip install pyspike)")
                          return 0
                      
                      spkmat = [pyspike.SpikeTrain([spkt for spkind,spkt in zip(spkinds,spkts) 
                          if spkind==gid], timeRange) for gid in set(spkinds)]
                      if stat == 'sync':
                          # (SPIKE-Sync measure)' # see http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony
                          syncMat = [pyspike.spike_sync(spkmat)]
                          #graphType = 'bar'
                      elif stat == 'pairsync':
                          # (SPIKE-Sync measure)' # see http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony
                          syncMat = np.mean(pyspike.spike_sync_matrix(spkmat), 0)
                          
  
                      statData.append(syncMat)
  
              colors.insert(0, popColors[subset] if subset in popColors 
                  else colorList[iplot%len(colorList)])  # colors in inverse order
  
          # if 'allCells' included make it black
          if include[0] == 'allCells':
              #if graphType == 'boxplot':
              colors.insert(len(include), (0.5,0.5,0.5))  # 
              del colors[0]
  
          # boxplot
          if graphType == 'boxplot':
              meanpointprops = dict(marker=(5,1,0), markeredgecolor='black', markerfacecolor='white')
              labels = legendLabels if legendLabels else include
              bp=plt.boxplot(statData, labels=labels[::-1], notch=False, sym='k+', meanprops=meanpointprops, 
                          whis=1.5, widths=0.6, vert=False, showmeans=True, patch_artist=True)
              plt.xlabel(xlabel, fontsize=fontsiz)
              plt.ylabel('Population', fontsize=fontsiz) 
  
              icolor=0
              borderColor = 'k'
              for i in range(0, len(bp['boxes'])):
                  icolor = i
                  bp['boxes'][i].set_facecolor(colors[icolor])
                  bp['boxes'][i].set_linewidth(2)
                  # we have two whiskers!
                  bp['whiskers'][i*2].set_color(borderColor)
                  bp['whiskers'][i*2 + 1].set_color(borderColor)
                  bp['whiskers'][i*2].set_linewidth(2)
                  bp['whiskers'][i*2 + 1].set_linewidth(2)
                  bp['medians'][i].set_color(borderColor)
                  bp['medians'][i].set_linewidth(3)
                  # for f in bp['fliers']:
                  #    f.set_color(colors[icolor])
                  #    print f
                  # and 4 caps to remove
                  for c in bp['caps']:
                      c.set_color(borderColor)
                      c.set_linewidth(2)
  
              ax = plt.gca()
              ax.spines['top'].set_visible(False)
              ax.spines['right'].set_visible(False)
              ax.spines['bottom'].set_visible(False)
              ax.get_xaxis().tick_bottom()
              ax.get_yaxis().tick_left()
              ax.tick_params(axis='x', length=0)
              ax.tick_params(axis='y', direction='out')
              ax.grid(axis='x', color="0.9", linestyle='-', linewidth=1)
              ax.set_axisbelow(True)
              if xlim: ax.set_xlim(xlim)
  
          # histogram
          elif graphType == 'histogram':
              import numpy as np
  
              nmax = 0
              pdfmax = 0
              binmax = 0
              for i,data in enumerate(statData):  # fix 
                  if histlogx:
                      histbins = np.logspace(np.log10(histmin), np.log10(max(data)), bins)
                  else:
                      histbins = bins
  
                  if histmin: # min value 
                      data = np.array(data)
                      data = data[data>histmin]
                  
                  # if histlogy:
                  #     data = [np.log10(x) for x in data]
  
                  if density:
                      weights = np.ones_like(data)/float(len(data)) 
                  else: 
                      weights = np.ones_like(data)
  
                  n, binedges,_ = plt.hist(data,  bins=histbins, histtype='step', color=colors[i], linewidth=2, weights=weights)#, normed=1)#, normed=density)# weights=weights)
                  if histShading:
                      plt.hist(data, bins=histbins, alpha=0.05, color=colors[i], linewidth=0, weights=weights) 
                  label = legendLabels[-i-1] if legendLabels else str(include[-i-1])
                  if histShading:
                      plt.hist([-10], bins=histbins, fc=((colors[i][0], colors[i][1], colors[i][2],0.05)), edgecolor=colors[i], linewidth=2, label=label)
                  else:
                      plt.hist([-10], bins=histbins, fc=((1,1,1),0), edgecolor=colors[i], linewidth=2, label=label)
                  nmax = max(nmax, max(n))
                  binmax = max(binmax, binedges[-1])
                  if histlogx: 
                      plt.xscale('log')
  
                  if normfit:
                      def lognorm(meaninput, stdinput, binedges, n, popLabel, color):
                          from scipy import stats 
                          M = float(meaninput) # Geometric mean == median
                          s = float(stdinput) # Geometric standard deviation
                          mu = np.log10(M) # Mean of log(X)
                          sigma = np.log10(s) # Standard deviation of log(X)                        
                          shape = sigma # Scipy's shape parameter
                          scale = np.power(10, mu) # Scipy's scale parameter
                          x = [(binedges[i]+binedges[i+1])/2.0 for i in range(len(binedges)-1)]  #np.linspace(histmin, 30, num=400) # values for x-axis
                          pdf = stats.lognorm.pdf(x, shape, loc=0, scale=scale) # probability distribution
                          R, p = scipy.stats.pearsonr(n, pdf)
                          print('    Pop %s rate: mean=%f, std=%f, lognorm mu=%f, lognorm sigma=%f, R=%.2f (p-value=%.2f)' % (popLabel, M, s, mu, sigma, R, p))
                          plt.semilogx(x, pdf, color=color, ls='dashed')
                          return pdf
  
  
                      fitmean = np.mean(data)
                      fitstd = np.std(data)
                      pdf=lognorm(fitmean, fitstd, binedges, n, label, colors[i])
                      pdfmax = max(pdfmax, max(pdf))
                      nmax = max(nmax, pdfmax)
  
                      # check normality of distribution
                      #W, p = scipy.stats.shapiro(data)
                      #print 'Pop %s rate: mean = %f, std = %f, normality (Shapiro-Wilk test) = %f, p-value = %f' % (include[i], mu, sigma, W, p)
  
  
              plt.xlabel(xlabel, fontsize=fontsiz)
              plt.ylabel('Probability of occurrence' if density else 'Frequency', fontsize=fontsiz)
              xmax = binmax
              plt.xlim(histmin, xmax)
              plt.ylim(0, 1.1*nmax if density else np.ceil(1.1*nmax)) #min(n[n>=0]), max(n[n>=0]))
              plt.legend(fontsize=fontsiz)
  
              # plt.show() # REMOVE!
  
              # if xlim: ax.set_xlim(xlim)
  
              #from IPython import embed; embed()
  
          # scatter
          elif graphType == 'scatter':
              from scipy import stats
              for i,(ynorms,data) in enumerate(zip(ynormsData, statData)):
                  mean, binedges, _ = stats.binned_statistic(ynorms, data, 'mean', bins=bins)
                  median, binedges, _ = stats.binned_statistic(ynorms, data, 'median', bins=bins)
                  #p25 = lambda x: np.percentile(x, 25)
                  #p75 = lambda x: np.percentile(x, 75)
                  
                  std, binedges, _ = stats.binned_statistic(ynorms, data, 'std', bins=bins)
                  #per25, binedges, _ = stats.binned_statistic(ynorms, data, p25, bins=bins)
                  #per75, binedges, _ = stats.binned_statistic(ynorms, data, p75, bins=bins)
                  
                  label = legendLabels[-i-1] if legendLabels else str(include[-i-1])
                  if kwargs.get('differentColor', None):
                      threshold = kwargs['differentColor'][0]
                      newColor = kwargs['differentColor'][1] #
                      plt.scatter(ynorms[:threshold], data[:threshold], color=[6/255.0,8/255.0,(64+30)/255.0], label=label, s=2) #, [0/255.0,215/255.0,255/255.0], [88/255.0,204/255.0,20/255.0]
                      plt.scatter(ynorms[threshold:], data[threshold:], color=newColor, alpha= 0.2, s=2) #[88/255.0,204/255.0,20/255.0]
                  else:
                      plt.scatter(ynorms, data, color=[0/255.0,215/255.0,255/255.0], label=label, s=2) #[88/255.0,204/255.0,20/255.0]
                  binstep = binedges[1]-binedges[0]
                  bincenters = [b+binstep/2 for b in binedges[:-1]] 
                  plt.errorbar(bincenters, mean, yerr=std, color=[6/255.0,70/255.0,163/255.0], fmt = 'o-',capthick=1, capsize=5) #[44/255.0,53/255.0,127/255.0]
                  #plt.errorbar(bincenters, mean, yerr=[mean-per25,per75-mean], fmt='go-',capthick=1, capsize=5)
              ylims=plt.ylim()
              plt.ylim(0,ylims[1])
              plt.xlabel('normalized y location (um)', fontsize=fontsiz)
              #plt.xlabel('avg rate (Hz)', fontsize=fontsiz)
              plt.ylabel(xlabel, fontsize=fontsiz)
              plt.legend(fontsize=fontsiz)
  
              #from IPython import embed; embed()
  
          # elif graphType == 'bar':
          #     print range(1, len(statData)+1), statData
          #     plt.bar(range(1, len(statData)+1), statData, tick_label=include[::-1], 
          #       orientation='horizontal', colors=colors)
  
          try:
              plt.tight_layout()
          except:
              pass
  
          # save figure data
          if saveData:
              figData = {'include': include, 'statData': statData, 'timeRange': timeRange, 'saveData': saveData, 'saveFig': saveFig, 'showFig': showFig}
  
              _saveFigData(figData, saveData, 'spikeStats_'+stat)
  
          # save figure
          if saveFig: 
              if isinstance(saveFig, basestring):
                  filename = saveFig+'_'+'spikeStat_'+graphType+'_'+stat+'.png'
              else:
                  filename = sim.cfg.filename+'_'+'spikeStat_'+graphType+'_'+stat+'.png'
              plt.savefig(filename, dpi=dpi)
  
          # show fig 
          if showFig: _showFigure()
  
      return fig, {'include': include, 'statData': statData, 'gidsData':gidsData, 'ynormsData':ynormsData}