pretty_plot.py

# simply plotting suite

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import networkx as nx
import mpl_toolkits.axisartist.floating_axes as floating_axes
from matplotlib.projections import PolarAxes
import pandas as pd
import seaborn as sns
sns.set_style('white')
import matplotlib as mpl
import scipy.stats as stats



########################################################################
# Helper functions
########################################################################


def lims(V):
  # get upper and lower bounds on a LIST OF LISTS (or matrix) D=2
  # returns min, max
  mi, ma = np.inf, 0
  for m in V:
    for n in m:
      if n > ma:
        ma = n
      if n < mi:
        mi = n
  return mi, ma


def condition_by_name(labels, arr, arr2=None, arr3=None):
  # sort by common labels in order so same types show up next to one another
  unique_labels = np.unique([i for i in labels])
  order = []
  for i in unique_labels:
    for x in range(len(labels)):
      if labels[x] == i:
        order.append(x)
  new_labels = [labels[j] for j in order]
  new_arr = [arr[j] for j in order]
  if arr2:
    new_arr2 = [arr2[j] for j in order]
  if arr3:
    new_arr3 = [arr3[j] for j in order]
    return new_labels, new_arr, new_arr2, new_arr3
  if arr2:
    return new_labels, new_arr, new_arr2
  return new_labels, new_arr



def group_by_name(labels, arr, metric='mean'):
  # Group data by name for bar/summary plotting, metric=('mean','median')
  llist = list(np.unique(labels))
  vals = {}
  for l in range(len(labels)):
    if labels[l] not in vals.keys():
      vals[labels[l]] = []
    if metric == 'mean':
      vals[labels[l]].append(np.mean(arr[l]))
    elif metric == 'median':
      vals[labels[l]].append(np.median(arr[l]))
  lab_data, arr_data = [k for k in vals.keys()], [np.mean(d) for d in vals.values()]
  return lab_data, arr_data, [np.std(d) for d in vals.values()]



def rm_nan(vec):
  # Remove nans (generated by illegally conditioned angles fed to cos-1)
  return [i for i in vec if str(i) != 'nan']



def dict_summary(obj):
  """
  Print a summary of the object values in a dictionary. Only works for
  dicts of keys=['files','cellTypes','obj1, 'obj2', ...]
  """
  grouplist = list(set(obj['cellTypes']))
  for k in obj.keys():
    if k != 'files' and k != 'cellTypes': # Don't do this for files/types
      print('Property: %s' %k)
      groups = [[] for i in grouplist] # New group for each key
      for samp in range(len(obj[k])):
        groups[grouplist.index(obj['cellTypes'][samp])].append(obj[k][samp])
      # Now print group stats for each group
      ind_vars = []
      for g in range(len(groups)):
        ind_means = [np.mean(i) for i in groups[g]]
        ind_meds = [np.median(i) for i in groups[g]]
        ind_std = [np.std(i) for i in groups[g]]
        for sub in groups[g]:
          ind_vars.append(np.var(sub))
        ind_qs = [[np.percentile(i,25), np.percentile(i,75)]
                      for i in groups[g]]
        ic = np.var(ind_means)**2 / \
             (np.var(ind_means)**2 + sum([np.var(i) for i in groups[g]])**2)
        print('For %i %s objects' %(len(groups[g]), grouplist[g]))
        print('Mean: %.5f +/- %.5f, Median: %.5f' 
              %(np.mean(ind_means), np.mean(ind_std), np.mean(ind_meds)))
        print('Means +/- std            Median   IQR (25th - 75th):')
        for u in range(len(groups[g])):
          print('%.5f +/- %.5f        %.5f (%.5f - %.5f)'
                %(ind_means[u], ind_std[u], ind_meds[u], ind_qs[u][0], ind_qs[u][1]))
        print('Intraclass correlation: %.5f' %ic)
      g_means = []
      for g in groups:
        g_means.append(np.mean([np.mean(i) for i in g]))
      master_ic = np.var(g_means)**2 / \
                  (np.var(g_means)**2 + sum(ind_vars)**2)
      print('Master intraclass correlation: %.5f' %master_ic)
  return
  

#group_means = [np.mean(k) for k in v_means] # group_means are the master means (only 4)
#master_ic = np.var(group_means)**2 / \
#            (np.var(group_means)**2 + sum([np.var(i) for i in v_means])**2)


def pyplot_defaults():
  """ After calling pandas or seaborn, the default plt can be messed up.
  If this causes inconsistencies, can call this to revert back to defaults. """
  mpl.rcParams.update(mpl.rcParamsDefault)
  return



########################################################################
# Basic plots: box bar, scatter (multi-dimensional plots next section)
########################################################################



def pretty_boxplot(V, labelsin, title=None, ticks=None, axes=None):
  """
  V is a matrix of arrays to be plotted; labels must be same length.
  labels = ['PD', 'LG' ...], ticks = ['791_233', ...]
  ticks=False to omit x-ticks altogether; None for legend=ticks
  """
  # mi, ma = lims(V)
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fcolors = ['darkkhaki', 'royalblue', 'forestgreen','lavenderblush']
  # plotting
  fig = plt.figure()
  ax = fig.add_subplot(111)
  if ticks:
    box = ax.boxplot(V, labels=ticks, showmeans=True, notch=True, 
                     patch_artist=True)
    #ax.set_xticks([i-1 for i in range(len(ticks))])
    #ax.set_xticklabels(ticks, rotation=45) # rotate ticks
  elif ticks is None:
    box = ax.boxplot(V, labels=[' ' for i in range(len(labelsin))],
                     showmeans=True, notch=True, patch_artist=True,
                     showfliers=False)
  else:
    box = ax.boxplot(V, labels=labelsin, showmeans=True, notch=True, 
                     patch_artist=True, showfliers=False)
  for patch, color in zip(box['boxes'], [fcolors[i] for i in C]):
    patch.set_facecolor(color)
  # set y axis range
  # ax.set_ylim([mi, ma])
  # legend
  khaki_patch = mpatches.Patch(color='darkkhaki', 
                label=labelsin[C.index(fcolors.index('darkkhaki'))])
  royal_patch = mpatches.Patch(color='royalblue',   
                label=labelsin[C.index(fcolors.index('royalblue'))])
  forest_patch = mpatches.Patch(color='forestgreen', 
                label=labelsin[C.index(fcolors.index('forestgreen'))])
  lavender_patch = mpatches.Patch(color='lavenderblush', 
                label=labelsin[C.index(fcolors.index('lavenderblush'))])
  plt.legend(handles=[khaki_patch, royal_patch, forest_patch, lavender_patch])
  # titles
  if axes:
    ax.set_xlabel(axes[0], fontsize=15)
    ax.set_ylabel(axes[1], fontsize=15)
  ax.set_title(title, fontsize=20)
  
  plt.show()
  return



def mean_scatter(V, labelsin, title=None, ticks=None, axes=None, 
                   showboth=True, showleg='best'):
  """
  One point per cell, plotted by type. showboth=Also show median.
  """
  Vmean = [np.mean(i) for i in V] # Take the mean of each element (okay for lists and scalars)
  Vmed = [np.median(i) for i in V]
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fcolors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  # plotting
  fig = plt.figure()
  ax = fig.add_subplot(111)
  for v in range(len(V)):
    meanline = ax.scatter(L.index(labelsin[v]), Vmean[v],  c=fcolors[C[v]], s=100, 
               marker='o', edgecolor='k', alpha=0.6)#fcolors[C[v]])
    if showboth:
      medline = ax.scatter(L.index(labelsin[v])+.2, Vmed[v], c=fcolors[C[v]], s=100,
                 marker='s', edgecolor='k', alpha=0.6)#fcolors[C[v]])
  # legend
  khaki_patch = mpatches.Patch(color='darkkhaki', 
                label=labelsin[C.index(fcolors.index('darkkhaki'))])
  royal_patch = mpatches.Patch(color='royalblue', 
                label=labelsin[C.index(fcolors.index('royalblue'))])
  patches = [khaki_patch, royal_patch]
  if len(L) > 2:
    forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[C.index(fcolors.index('forestgreen'))])
    patches.append(forest_patch)
  if len(L) > 3:
    tomato_patch = mpatches.Patch(color='tomato', 
                  label=labelsin[C.index(fcolors.index('tomato'))])
    patches.append(tomato_patch)
  if showleg is not None:
    ax.legend(handles=patches, loc=showleg, fontsize=15) # I often mess with this
  if ticks is None:
    ax.tick_params(axis='x',which='both',bottom='off',top='off',
                         labelbottom='off')
  # title
  if title:
    ax.set_title(title, fontsize=20)
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=25)
    ax.set_ylabel(axes[1], fontsize=25)
  ax.tick_params(axis='y', which='major', labelsize=20) # Size of y ticks
  ax.locator_params(nbins=4)  # Number of y-tick bins
  plt.show(); return 
  


def pretty_scatter(V, labelsin, title=None, ticks=None, axes=None, 
                   showleg='best', moreD=None, bounds=None):
  """
  moreD = another set of type V to be plot on right axes, len(axes)=3
  """
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fcolors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  # plotting
  fig = plt.figure(figsize=(2,4))
  ax = fig.add_subplot(111)
  ax2 = None
  if moreD is not None:
    if len(moreD) == len(V):
      ax2 = ax.twinx()
  for v in range(len(V)):
    ax.scatter(L.index(labelsin[v]), V[v],  c=fcolors[C[v]], s=80, 
               marker='o', edgecolor='k', alpha=0.6)#fcolors[C[v]]) 
    if ax2 is not None:
      ax2.scatter(L.index(labelsin[v])+.1, moreD[v],  c=fcolors[C[v]], s=80, 
                 marker='s', edgecolor='k', alpha=0.6)#fcolors[C[v]])
  # legend
  khaki_patch = mpatches.Patch(color='darkkhaki', 
                label=labelsin[C.index(fcolors.index('darkkhaki'))])
  royal_patch = mpatches.Patch(color='royalblue', 
                label=labelsin[C.index(fcolors.index('royalblue'))])
  patches = [khaki_patch, royal_patch]
  if len(L) > 2:
    forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[C.index(fcolors.index('forestgreen'))])
    patches.append(forest_patch)
  if len(L) > 3:
    tomato_patch = mpatches.Patch(color='tomato', 
                  label=labelsin[C.index(fcolors.index('tomato'))])
    patches.append(tomato_patch)
  if showleg is not None:
    ax.legend(handles=patches, loc=showleg, fontsize=15) # I often mess with this
  if ticks is None:
    ax.tick_params(axis='x',which='both',bottom='off',top='off',
                         labelbottom='off')
  # Get rid of side labels
  for pos in ['top', 'right']:
    ax.spines[pos].set_visible(False)
  # title
  if bounds is not None:
    ax.set_ylim([bounds[0], bounds[1]])
  if title:
    ax.set_title(title, fontsize=20)
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=15)
    ax.set_ylabel(axes[1], fontsize=15)
    if ax2 is not None:
      ax2.set_ylabel(axes[2], fontsize=15) # color='r'
  plt.show(); return 
  


def simple_scatter(x, y, fit=False, title=None, axes=None, showtext=True):
  """
  """
  import scipy.stats as stats
  if type(x[0]) is list: # Multiple plots
    fcolors = ['darkkhaki', 'royalblue', 'forestgreen','deeppink']
    if len(x) > 4:
      print('Can only handle 4 data groups, got %i' %len(x))
  else:
    x, y = [x], [y]
    fcolors = ['steelblue']
  fig = plt.figure()
  ax = fig.add_subplot(111)
  for v in range(len(x)):
    ax.scatter(x[v], y[v], c='ivory', edgecolor=fcolors[v], 
               linewidth=2, s=50, alpha=0.5)
    beta, alpha, r, p, _ = stats.linregress(x[v], y[v])
    ax.plot([min(x[v]), max(x[v])], 
            [min(x[v])*beta+alpha, max(x[v])*beta+alpha], c=fcolors[v],
            linewidth=2.5)
    print('Group %i fit: R=%.5f, P=%.5f (alpha=%.f5, beta=%.5f)'
          %(v, r, p, alpha, beta))
  if title is not None:
    ax.set_title(title, fontsize=20)
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=25)
    ax.set_ylabel(axes[1], fontsize=25)
  #ax.set_xlim([min([min(i) for i in x]), 0.7*max([max(i) for i in x])])
  #ax.set_ylim([min([min(i) for i in y]), max([max(i) for i in y])])
  if showtext is True:
    ax.text(min(x[0]), max(y[0])*.8, r'$r=%.3f$' %(r), 
            fontsize=25) # Can also print P-value as well
  plt.show()
  return
              


# Simple bar plot

def pretty_bar(v, labelsin, stderr=None, ticks=None, title=None, axes=None,
               showleg='best', bounds=None):
  """
  showleg can also be None or 'uppermiddle' etc.
  """
  #mi, ma = lims(v)
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fcolors = ['darkkhaki', 'royalblue', 'forestgreen','lavenderblush']
  # plotting
  fig = plt.figure()
  ax = fig.add_subplot(111)
  if len(np.shape(v))>1:
    V = [np.mean(i) for i in v]
    Vstd = [np.std(i) for i in v]
  else:
    V=v
    Vstd=stderr
  num = len(V)
  # plotting
  if stderr is not None:
    for e in range(num):
      ax.bar(e,V[e], yerr=Vstd[e], color=fcolors[C[e]], ecolor='k', width=0.5)
  else:
    for e in range(num):
      ax.bar(e,V[e], color=fcolors[C[e]], ecolor='k', width=0.5)
  # legend
  khaki_patch = mpatches.Patch(color='darkkhaki', 
                label=labelsin[C.index(fcolors.index('darkkhaki'))])
  royal_patch = mpatches.Patch(color='royalblue', 
                label=labelsin[C.index(fcolors.index('royalblue'))])
  patches = [khaki_patch, royal_patch]
  if len(L) > 2:
    forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[C.index(fcolors.index('forestgreen'))])
    patches.append(forest_patch)
  if len(L) > 3:
    lavender_patch = mpatches.Patch(color='tomato', 
                  label=labelsin[C.index(fcolors.index('tomato'))])
    patches.append(lavender_patch)
  if showleg is not None:
    plt.legend(handles=patches, loc=showleg)
  if stderr is not None:
    ax.set_ylim([0, 1.5*(max(V)+np.mean(Vstd))])
  else:
    ax.set_ylim([0, max(V)+0.5*min(V)])
  ax.set_xlim([-.5, len(V)])
  if bounds is not None and type(bounds) is list:
    ax.set_ylim([bounds[0], bounds[1]])
  if ticks is None:
    ax.tick_params(axis='x',which='both',bottom='off',top='off',
                         labelbottom='off')
  for pos in ['bottom', 'left', 'top', 'right']:
    ax.spines[pos].set_visible(False)
  # title
  if title:
    ax.set_title(title, fontsize=20)
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=15)
    ax.set_ylabel(axes[1], fontsize=15)
  plt.show(); return  



def pretty_stackedbar(V, labelsin, title=None, colors=None,
                      axes=None, showxticks=True, counts=True, ):
  """
  Stacked bar plots. Groups are labeled (x-tick) and binned items are
  colored. Colors can be supplied by user; must be >= len(kkeys) (one
  unique color per unique item, NOT per group).
  _counts_ determines whether percentages have already been computed,
  if False, elements of _V_ should be dicts (keys are groups).
  """
  # Colors - must be different, can be supplied by user
  if colors is None:
    colors = ['darkorchid', 'orchid', 'violet', 'pink', 'deeppink']
  # Most abundant is on bottom, and so on
  if counts:
    print('Computing counts....')
    plotgroups = [{i: float(V[k].count(i))/len(V[k]) for i in list(set(V[k]))} 
                  for k in range(len(V))]
  else:
    print('Using supplied values...')
    plotgroups = V
  # Fill in the dicts
  kkeys = []
  for gr in plotgroups:
    for k in gr.keys():
      if k not in kkeys:
        kkeys.append(k)
  for u in range(len(plotgroups)):
    for k in kkeys:
      if k not in plotgroups[u].keys():
        plotgroups[u][k] = 0. # Just make it a zero
  bottoms = [0. for i in plotgroups]
  # Plotting
  plt.figure()
  pbars = []
  for k in kkeys:
    pbars.append(plt.bar(np.arange(len(V)), 
                         [plotgroups[i][k] for i in range(len(V))],
                         width=.5, color=colors[k],
                         edgecolor='white',
                         bottom=bottoms))
    bottoms = [bottoms[b]+plotgroups[b][k] for b in range(len(bottoms))]
  # Plot aesthetics
  if axes is not None:
    plt.xlabel(axes[0])
    plt.ylabel(axes[1])
  if title is not None:
    plt.title(title)
  if showxticks:
    plt.xticks(np.arange(len(V))+0.25, labelsin)
  plt.xlim([-.5, len(V)])
  plt.legend(pbars, kkeys)
  plt.show()
  return



def plot_cum_dist(V, labelsin, title=None):
  """
  Plot lines showing cumulative distribution, i.e. for Sholl.
  Labels must be len(V) (one label per array).
  """
  if max(V[0]) > 1:
    normed = []
    for i in V:
      M = max(i)
      normed.append([a/M for a in i])
    V = normed
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fcolors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  fig = plt.figure(); ax = fig.add_subplot(111)
  
  for i in range(len(V)):
    ax.plot(V[i], [a/len(V[i]) for a in range(len(V[i]))], color=fcolors[C[i]],
            linewidth=2)
    ax.plot(V[i], [a/len(V[i]) for a in range(len(V[i]))], color=fcolors[C[i]],
            linewidth=4, alpha=0.5)
  
  # legend
  khaki_patch = mpatches.Patch(color='darkkhaki', 
                label=labelsin[C.index(fcolors.index('darkkhaki'))])
  royal_patch = mpatches.Patch(color='royalblue', 
                label=labelsin[C.index(fcolors.index('royalblue'))])
  patches = [khaki_patch, royal_patch]
  if len(L) > 2:
    forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[C.index(fcolors.index('forestgreen'))])
    patches.append(forest_patch)
  if len(L) > 3:
    lavender_patch = mpatches.Patch(color='tomato', 
                  label=labelsin[C.index(fcolors.index('tomato'))])
    patches.append(lavender_patch)
  plt.legend(handles=patches)
  # title
  if title:
    ax.set_title(title, fontsize=20)
  plt.show()



########################################################################
# Pretty dimensional scatters
########################################################################

#  -- colors are subjects but all axes are dimensions (not groups)
# Data do not need to be conditioned for these

def pretty_3d(v1, labelsin, v2, v3, shadows='z', title=None, axes=None,
              showleg=None, lines=None, figsize=[5,5], ellipses=True):
  """
  A 3-D scatter plot of various features (must have 3). shadows='x', 'y', 'z' 
  (points are projected onto the XY, XZ and YZ planes), len(axes)==3
  lines = ['x', 'y', 'z']
  """
  # Set up colors, check data
  if len(v1) == len(v2) and len(v2) == len(v3) and len(v3) == len(labelsin):
    pass
  else:
    print('Data are of different lengths!'); return
  L = list(np.unique(labelsin)) # Get color data
  C = [L.index(i) for i in labelsin]
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  # 
  fig = plt.figure(figsize=figsize, dpi=150)
  ax = fig.add_subplot(111, projection='3d')
  # For each data point plot it and its projection
  for u in range(len(v1)):
    ax.scatter(v1[u], v2[u], v3[u], color=colors[C[u]], edgecolor='k',
               alpha=0.7, s=30) # Main scatter point
    if shadows is not None:
      if type(shadows) is not list:
        shadows = list(shadows)
      for sha in shadows:
        shpts = {'x': v1[u], 'y': v2[u], 'z': v3[u]}
        shpts[sha] = 0
        ax.scatter(shpts['x'], shpts['y'], shpts['z'], color=colors[C[u]], 
                   edgecolor='none', alpha=0.2, s=20) # XY plane (Z=0)
    if lines is not None:
      if type(lines) is not list:
        lines = list(lines)
      for la in lines:
        lapts = {'x': v1[u], 'y': v2[u], 'z': v3[u]}
        lapts[la] = 0
        ax.plot([v1[u], lapts['x']], [v2[u],lapts['y']], [v3[u],lapts['z']],
                color=colors[C[u]], alpha=0.2, linewidth=1.5) # XY plane (Z=0)
  if ellipses:
    from matplotlib.patches import Ellipse
    for u in range(len(L)):
      dat1 = [v1[j] for j in range(len(C)) if C[j]==u]
      dat2 = [v2[j] for j in range(len(C)) if C[j]==u]
      cov = np.cov(dat1, dat2)
      lambda_, v = np.linalg.eig(cov)
      lambda_ = np.sqrt(lambda_)
      ell = Ellipse(xy=(np.mean(dat1), np.mean(dat2)), # Initial offset
                    width=lambda_[0]*2, height=lambda_[1]*2,
                    angle=np.rad2deg(np.arccos(v[0,0])), alpha=0.1)
      ell.set_facecolor(colors[u])
      ell.set_edgecolor(colors[u])
      ell.set_linewidth(1)
      # ell.set_linestyle('--')
      ax.add_patch(ell)
      pathpatch_2d_to_3d(ell, z=0, normal='z')
      pathpatch_translate(ell, (0, 0, 0)) # Here can add _additional_ offset in x,y,z
      #ax.add_artist(ell)
  # Labels, legends
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=15)
    ax.set_ylabel(axes[1], fontsize=15)
    ax.set_zlabel(axes[2], fontsize=15)
  if showleg is not None:
    plt.legend(colors, L, loc=showleg)
  ax.grid(False) # Remove the grids
  ax.xaxis.pane.set_edgecolor('black')
  ax.yaxis.pane.set_edgecolor('black')
  ax.w_xaxis.set_pane_color((0.,0.,0.,0.1))
  ax.w_yaxis.set_pane_color((0.,0.,0.,0.15))
  ax.w_zaxis.set_pane_color((1.,1.,1.,0.))
  for spax in ['x', 'y', 'z']:
    plt.locator_params(axis=spax, nbins=3)
  ax.set_zlim([0, (int(max(v3) / 10**int(math.log10(max(v3))))+1 ) * 10**int(math.log10(max(v3)))])
  ax.set_zticks([0, (int(max(v3) / 10**int(math.log10(max(v3))))+1 ) * 10**int(math.log10(max(v3)))])
  #ax.set_xlim([0,12000]) # Mess around with stuff as needed
  plt.show()
  return





def pretty_2d(v1, labelsin, v2, title=None, axes=None, showleg=None,
              ellipses=True, figsize=[5,4], alteig=[]):
  """
  A 2-D scatter plot of various features (must have 2). 
  len(axes)==2. ellipses = 1 std dev. If any elipses look back, enter
  that order into alteig (i.e.: alteig=[1,2] if the 2nd and 3rd ells look off).
  """
  # Set up colors, check data
  if len(v1) == len(v2) and len(v2) == len(labelsin):
    pass
  else:
    print('Data are of different lengths!'); return
  L = list(np.unique(labelsin)) # Get color data
  C = [L.index(i) for i in labelsin]
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  print('Order is: darkkhaki, royalblue, forestgreen, tomato') ## Tell user the order
  fig = plt.figure(figsize=figsize, dpi=100)
  ax = fig.add_subplot(111)
  # For each data point plot it and its projection
  for u in range(len(v1)):
    ax.scatter(v1[u], v2[u], color=colors[C[u]], edgecolor='k',
               alpha=1, s=40) # Main scatter point
  # Add std ellipses via eigenvectors etc
  if ellipses: 
    from matplotlib.patches import Ellipse
    for u in range(len(L)):
      dat1 = [v1[j] for j in range(len(C)) if C[j]==u]
      dat2 = [v2[j] for j in range(len(C)) if C[j]==u]
      cov = np.cov([i for i in dat1 if not pd.isnull(i)], 
                   [i for i in dat2 if not pd.isnull(i)])
      if u in alteig:
        vals, vecs = eigsorted(cov)
        theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
        w, h = np.sqrt(vals) * np.array([0.5,2])
        ell = Ellipse(xy=(np.mean(dat1), np.mean(dat2)),
                      width=2, height=h, angle=theta, alpha=0.1)
      else:
        lambda_, v = np.linalg.eig(cov)
        lambda_ = np.sqrt(lambda_)
        ell = Ellipse(xy=(np.mean(dat1), np.mean(dat2)),
                      width=lambda_[0]*2, height=lambda_[1]*2,
                      angle=np.rad2deg(np.arccos(v[0,0])), alpha=0.1)
      ell.set_facecolor(colors[u])
      ell.set_linewidth(1)
      ell.set_edgecolor(colors[u])
      ax.add_artist(ell)
  # Labels, legends
  if axes is not None:
    ax.set_xlabel(axes[0], fontsize=15)
    ax.set_ylabel(axes[1], fontsize=15)
  if showleg is not None:
    plt.legend(colors, L, loc=showleg)
  ax.grid(False) # Remove the grids
  plt.locator_params(axis='x', nbins=4)
  plt.locator_params(axis='y', nbins=4)
  for pos in ['top', 'right']:
    ax.spines[pos].set_visible(False)
  plt.tight_layout() # OR plt.gca().tight_layout() # gca = get current axes
  plt.show()
  return



def hist_2d(x, y, labelsin=None, axes=None, showstats=False, 
            tint=False, forcebins=10):
  """
  Scatter-plot 2-D data with histograms on either axes.
  Colors codes x-y pairs by groups and then can allow the histograms 
  to be tinted (tint=True) to show composition of each group.
  x,y can be pd.Series from DataFrame.
  """
  # Colors etc
  if labelsin is not None:
    uniq_cols = [list(np.random.random(3)) for i in list(set(labelsin))]
    cols = [uniq_cols[list(set(labelsin)).index(i)] for i in labelsin]
    col_dict = {}
    for l in range(len(labelsin)):
      if cols[l] not in col_dict.values():
        col_dict[len(col_dict.keys())] = cols[l]
  else:
    cols = ['blue' for i in range(len(x))]
  ###   Scatter
  plt.subplot2grid( (4,4), (1,0), rowspan=3, colspan=3)
  for p in range(len(x)):
    plt.plot(x[p], y[p], color=cols[p], marker='o', markeredgecolor='none',
             alpha=0.9, ) # markersize=2
  Xrange, Yrange = plt.xlim(), plt.ylim()
  if axes is not None:
    plt.xlabel(axes[0], fontsize=15); plt.ylabel(axes[1], fontsize=15)
  # Stats
  if showstats:
    from scipy.stats import spearmanr as spear
    plt.text(Xrange[1]/2., Yrange[1]-Yrange[1]*.1, 'spearman rank: %.3f' %spear(x,y)[0])
  ###   Horizontal histogram
  plt.subplot2grid( (4,4), (0,0), colspan=3)
  bin_edges = np.linspace(Xrange[0], Xrange[1], forcebins+1)
  bin_colors = [ [] for i in range(forcebins)]
  for d in range(len(x)): # X values only for horizontal
    bin_colors[len([b for b in bin_edges if x[d]>=b])-1].append(cols[d])
  if tint: # Add a tint to the color -- average colors
    bin_col = [[np.mean([k[i] for k in bin_colors[bin_]]) for i in range(3)]
               for bin_ in range(len(bin_colors))]
    for b in range(len(bin_col)): # len(bin_col) == len(bin_edges)-1
      plt.bar((bin_edges[b]+bin_edges[b+1])/2., len(bin_colors[b]),
              color=bin_col[b], edgecolor='white', 
              width=(bin_edges[b+1]-bin_edges[b])*.95)
  else: # No tint, stack the bars (hori)
    for b in range(len(bin_colors)):
      try:
        bar_colors = [col_ for col_ in col_dict.values() if col_ in bin_colors[b]]
      except:
        print(bin_colors[b], col_dict.values())
      bottom = 0
      for color_ in bar_colors:
        height = bin_colors[b].count(color_)
        plt.bar((bin_edges[b]+bin_edges[b+1])/2., height, bottom=bottom,
                color=color_, edgecolor='white', alpha=0.6,
                width=(bin_edges[b+1]-bin_edges[b])*.95)
        bottom += height
  plt.axis('off')
  ###   Vertical histogram
  plt.subplot2grid( (4,4), (1,3), rowspan=3)
  bin_edges = np.linspace(Yrange[0], Yrange[1], forcebins+1)
  bin_colors = [ [] for i in range(forcebins)]
  for d in range(len(y)): # Y values only for vertical
    bin_colors[len([b for b in bin_edges if y[d]>=b])-1].append(cols[d])
  if tint: # Add a tint to the color -- average colors
    bin_col = [[np.mean([k[i] for k in bin_colors[bin_]]) for i in range(3)]
               for bin_ in range(len(bin_colors))]
    for b in range(len(bin_col)): # len(bin_col) == len(bin_edges)-1
      plt.barh((bin_edges[b]+bin_edges[b+1])/2., len(bin_colors[b]),
              color=bin_col[b], edgecolor='white', 
              height=(bin_edges[b+1]-bin_edges[b])*.95)
  else: # No tint, stack the bars (vert)
    for b in range(len(bin_colors)):
      bar_colors = [col_ for col_ in col_dict.values() if col_ in bin_colors[b]]
      bottom = 0
      for color_ in bar_colors:
        height = bin_colors[b].count(color_)
        plt.barh((bin_edges[b]+bin_edges[b+1])/2., height, left=bottom,
                color=color_, edgecolor='white', alpha=0.6,
                height=(bin_edges[b+1]-bin_edges[b])*.95)
        bottom += height
  plt.axis('off')
  plt.show(); return




########################################################################
# Horizontal plots: histograms, violins and scatters
########################################################################


def hori_histogram(xdata, labelsin, title=None, axes=None, norm=False,
                       showmean=True, switch=False, llog=False, rrange=None,
                       forcebins=100, shade=True, eps=False, xcnt=True):
  """xcnt is the max of the y axis (on bottom)
  """
  if switch:
    for i in range(len(xdata)-1):
      xdata.append(xdata.pop(0))
      labelsin.append(labelsin.pop(0))
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato', 'darkorchid']
  altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink', 'plum']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  #print(L,C)
  fig = plt.figure(dpi=200) # Give it pub-quality DPI
  plots = [fig.add_subplot(1,len(xdata),i+1) for i in range(len(xdata))]
  if norm is True:
    #tdata = np.linspace(0,100,len(xdata[0]))
    X = []
    for x in xdata:
      X.append([i/max(x) for i in x])
    xdata = X
    minm, maxm = 0, 1.
  elif norm is False and rrange is None:
    minm, maxm = np.inf, 0 # condition the data
    for x in xdata:
      if np.mean(x)-np.std(x) < minm:
        minm = np.mean(x)-2*np.std(x)
      if np.mean(x)+np.std(x) > maxm:
        maxm = np.mean(x)+2*np.std(x)
    if minm < 0:
      minm = 0.
  else:
    minm, maxm = rrange[0], rrange[1]
  for p in range(len(xdata)): # now plot
    if type(forcebins) is list:
      b_e = forcebins
    else:
      b_e = np.linspace(minm, maxm, forcebins) # len/100 bins
    hist, _ = np.histogram(xdata[p], bins=b_e)
    plotbins = [(b_e[i]+b_e[i+1])/2. for i in range(len(b_e)-1)]
    # find the appropriate bar width #print(minm, maxm, p);
    hgt = (maxm-minm)/len([i for i in hist if i != 0]) # as high as there are filled hist elements
    hgt = plotbins[2]-plotbins[1]
    q25, q75 = np.percentile(xdata[p], [25, 75])
    barcols = [colors[C[p]] if q75 > plotbins[t] > q25 # In IQR
               else altcolors[C[p]] for t in range(len(plotbins))]
    for b in range(len(plotbins)):
      plots[p].barh(plotbins[b], hist[b]/max(hist), height=hgt, linewidth=0, alpha=0.9,
                    color=barcols[b], edgecolor=barcols[b]) #=colors[C[p]])
    # show the means:
    if showmean:
      plots[p].plot([0,1], [np.mean(xdata[p]), np.mean(xdata[p])], 
                    linewidth=1., c='purple')
      plots[p].plot([0,1], [np.median(xdata[p]), np.median(xdata[p])],
                    '--', linewidth=1., c='purple', )
    if p == 0: #if first plot, show the axes
      #plots[p].tick_params(axis='x',which='both',bottom='off',top='off',
      #                     labelbottom='off')
      if axes:
        plots[p].set_ylabel(axes[1], fontsize=15)
      plots[p].set_ylim([minm, maxm])
      if llog is True:
        plots[p].set_yscale('log'); plots[p].set_ylim([0, maxm]) ## Log scale
      for pos in ['top', 'right']:
        plots[p].spines[pos].set_visible(False)
    else:
      #plots[p].axis('off')
      #plots[p].tick_params(axis='x',which='both',bottom='off',top='off',
      #                     labelbottom='off')
      plots[p].get_yaxis().set_visible(False)
      if llog is True:
        plots[p].set_yscale('log') ## Log scale
      plots[p].set_ylim([minm,maxm])
      for pos in ['top', 'left', 'right']:
        plots[p].spines[pos].set_visible(False)
    plt.locator_params(nbins=4) #################### Set one x-tick
    plots[p].set_xticks([.5])
    plots[p].set_xticklabels(['%i' %int(max(hist))])
  if title:
    plt.suptitle(title, fontsize=20)
  plt.show()
  return


# A sample legend for horizontal histograms

def hori_bars_legend():
  plt.bar(range(100),np.random.normal(100),facecolor='darkgray', edgecolor='darkgray')
  plt.axvspan(25, 75, 0,1, color='gray', alpha=0.3)
  plt.plot([50,50],[0,1], '-', c='k', linewidth=3)
  plt.plot([55,55],[0,1], '--', c='k', linewidth=3)
  plt.tick_params(axis='x', which='both',bottom='off', top='off', labelbottom='off')
  plt.tick_params(axis='y', which='both',left='off', right='off', labelleft='off')
  plt.show()




def violin_box(xdata, labelsin, title=None, axes=None, norm=False,
                showmean=True, switch=False, llog=False, rrange=None,
                forcebins=100, shade=True, eps=False, xcnt=True):
  """xcnt is the max of the y axis (on bottom)
  """
  if switch:
    for i in range(len(xdata)-1):
      xdata.append(xdata.pop(0))
      labelsin.append(labelsin.pop(0))
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato', 'darkorchid']
  altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink', 'plum']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  #print(L,C)
  fig = plt.figure(dpi=200) # Give it pub-quality DPI
  plots = [fig.add_subplot(1,len(xdata),i+1) for i in range(len(xdata))]
  if norm is True:
    #tdata = np.linspace(0,100,len(xdata[0]))
    X = []
    for x in xdata:
      X.append([i/max(x) for i in x])
    xdata = X
    minm, maxm = 0, 1.
  elif norm is False and rrange is None:
    minm, maxm = np.inf, 0 # condition the data
    for x in xdata:
      if np.mean(x)-np.std(x) < minm:
        minm = np.mean(x)-2*np.std(x)
      if np.mean(x)+np.std(x) > maxm:
        maxm = np.mean(x)+2*np.std(x)
    if minm < 0:
      minm = 0.
  else:
    minm, maxm = rrange[0], rrange[1]
  for p in range(len(xdata)): # now plot
    if type(forcebins) is list:
      b_e = forcebins
    else:
      b_e = np.linspace(minm, maxm, forcebins) # len/100 bins
    hist, _ = np.histogram(xdata[p], bins=b_e)
    plotbins = [(b_e[i]+b_e[i+1])/2. for i in range(len(b_e)-1)]
    # find the appropriate bar width #print(minm, maxm, p);
    hgt = (maxm-minm)/len([i for i in hist if i != 0]) # as high as there are filled hist elements
    # hgt = plotbins[2]-plotbins[1]
    q25, q75 = np.percentile(xdata[p], [25, 75])
    barcols = [colors[C[p]] if q75 > plotbins[t] > q25 # In IQR
               else altcolors[C[p]] for t in range(len(plotbins))]
    # Plot the baseline
    for base in [[q25+hgt,q75,colors[C[p]]], [q25,min(plotbins),altcolors[C[p]]], [max(plotbins),q75+hgt,altcolors[C[p]]]]:
      plots[p].plot([0.,0.], [base[0],base[1]], c=base[2], linewidth=1, alpha=0.9) 
    for b in range(len(plotbins)): # Plot the +bar and -bar
      #plots[p].barh(plotbins[b], 0.5+hist[b]/(2*max(hist)), height=hgt, linewidth=0, alpha=0.9,
      #              color=barcols[b], edgecolor=barcols[b]) #=colors[C[p]])
      plots[p].barh(plotbins[b], hist[b]/(max(hist)), height=hgt, linewidth=0, alpha=0.9,
                    color=barcols[b], edgecolor=barcols[b],
                    left=-(hist[b]/(2*max(hist)))) #=colors[C[p]])
    # show the means:
    if showmean:
      plots[p].plot([-.5,.5], [np.mean(xdata[p]), np.mean(xdata[p])], 
                    linewidth=1., c='purple')
      plots[p].plot([-.5,.5], [np.median(xdata[p]), np.median(xdata[p])],
                    '--', linewidth=1., c='purple', )
    if p == 0: #if first plot, show the axes
      if axes:
        plots[p].set_ylabel(axes[1], fontsize=15)
      plots[p].set_ylim([minm, maxm])
      if llog is True:
        plots[p].set_yscale('log'); plots[p].set_ylim([0, maxm]) ## Log scale
      for pos in ['top', 'right']:
        plots[p].spines[pos].set_visible(False)
    else:
      plots[p].get_yaxis().set_visible(False)
      if llog is True:
        plots[p].set_yscale('log') ## Log scale
      plots[p].set_ylim([minm,maxm])
      for pos in ['top', 'left', 'right']:
        plots[p].spines[pos].set_visible(False)
    plt.locator_params(nbins=4) #################### Set one x-tick
    plots[p].set_xticks([.5])
    plots[p].set_xticklabels(['%i' %int(max(hist))])
    # plots[p].set_xlim([0,1.])
  if title:
    plt.suptitle(title, fontsize=20)
  plt.show()
  return






def violin_spline(xdata, labelsin, title=None, axes=None, norm=False,
                showmean=True, stepfilled=True, llog=False, rrange=None,
                forcebins=100, shade=True, eps=False, xcnt=True):
  """xcnt is the max of the y axis (on bottom)
  """
  from scipy.ndimage.filters import gaussian_filter1d as filt
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato', 'darkorchid']
  altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink', 'plum']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  #print(L,C)
  fig = plt.figure(dpi=200) # Give it pub-quality DPI
  plots = [fig.add_subplot(1,len(xdata),i+1) for i in range(len(xdata))]
  if norm is True:
    #tdata = np.linspace(0,100,len(xdata[0]))
    X = []
    for x in xdata:
      X.append([i/max(x) for i in x])
    xdata = X
    minm, maxm = 0, 1.
  elif norm is False and rrange is None:
    minm, maxm = np.inf, 0 # condition the data
    for x in xdata:
      if np.mean(x)-np.std(x) < minm:
        minm = np.mean(x)-2*np.std(x)
      if np.mean(x)+np.std(x) > maxm:
        maxm = np.mean(x)+2*np.std(x)
    if minm < 0:
      minm = 0.
  else:
    minm, maxm = rrange[0], rrange[1]
  for p in range(len(xdata)): # now plot
    if type(forcebins) is list:
      b_e = forcebins
    else:
      b_e = np.linspace(minm, maxm, forcebins) # len/100 bins
    hist, _ = np.histogram(xdata[p], bins=b_e)
    occupied = len([i for i in hist if i != 0])
    if type(forcebins) is int and occupied < forcebins/2.: # Sparse bins
      hist, b_e = np.histogram(xdata[p], bins=int(occupied*2))
    plotbins = [(b_e[i]+b_e[i+1])/2. for i in range(len(b_e)-1)]
    # find the appropriate bar width #print(minm, maxm, p);
    hgt = (maxm-minm)/occupied # as high as there are filled hist elements
    # hgt = plotbins[2]-plotbins[1]
    q25, q75 = np.percentile(xdata[p], [25, 75])
    barcols = [colors[C[p]] if q75 > plotbins[t] > q25 # In IQR
               else altcolors[C[p]] for t in range(len(plotbins))]
    # Plot the baseline
    for base in [[q25,q75,colors[C[p]]], [q25,min(plotbins),altcolors[C[p]]], [max(plotbins),q75+hgt,altcolors[C[p]]]]:
      plots[p].plot([0.,0.], [base[0],base[1]], c=base[2], linewidth=1, alpha=0.9) 
    fit = filt(hist, 1.5)
    q1_inds = [[plotbins[u], fit[u]] for u in range(len(fit)) if plotbins[u] < q25+hgt]
    q4_inds = [[plotbins[u], fit[u]] for u in range(len(fit)) if plotbins[u] > q75-hgt]
    iqr_inds = [[plotbins[u], fit[u]] for u in range(len(fit)) if q25 < plotbins[u] < q75]
    plots[p].fill_betweenx([q[0] for q in q1_inds],[q[1]/max(fit) for q in q1_inds], 
                          [-q[1]/max(fit) for q in q1_inds], color=altcolors[C[p]], alpha=0.9)
    plots[p].fill_betweenx([q[0] for q in q4_inds],[q[1]/max(fit) for q in q4_inds], 
                          [-q[1]/max(fit) for q in q4_inds], color=altcolors[C[p]], alpha=0.9)
    plots[p].fill_betweenx([q[0] for q in iqr_inds],[q[1]/max(fit) for q in iqr_inds], 
                          [-q[1]/max(fit) for q in iqr_inds], color=colors[C[p]], alpha=0.9)
    if stepfilled:
      xsteps = [plotbins[0], plotbins[0], plotbins[0], plotbins[1]]
      ysteps = [0,            hist[0],      hist[1], hist[1]]   # 
      for u in range(1,len(plotbins)-1):
        ysteps.append(hist[u+1]); ysteps.append(hist[u+1]) # ysteps xsteps-1 -> xsteps+1
        xsteps.append(plotbins[u]); xsteps.append(plotbins[u+1]) # now xsteps -> ysteps+1
      xsteps.append(xsteps[-1]); ysteps.append(plotbins[-1]) # Back to 0
      ysteps = np.array(ysteps)/max(ysteps)
      plots[p].plot(ysteps, xsteps, color='white', linewidth=0.5)
      plots[p].plot(-ysteps, xsteps, color='white', linewidth=0.5)
    # show the means:
    if showmean:
      plots[p].plot([-.5,.5], [np.mean(xdata[p]), np.mean(xdata[p])], 
                    linewidth=1., c='purple')
      plots[p].plot([-.5,.5], [np.median(xdata[p]), np.median(xdata[p])],
                    ':',linewidth=2., c='purple', )
    if p == 0: #if first plot, show the axes
      if axes:
        plots[p].set_ylabel(axes[1], fontsize=15)
      plots[p].set_ylim([minm, maxm])
      if llog is True:
        plots[p].set_yscale('log'); plots[p].set_ylim([0, maxm]) ## Log scale
      for pos in ['top', 'right']:
        plots[p].spines[pos].set_visible(False)
    else:
      plots[p].get_yaxis().set_visible(False)
      if llog is True:
        plots[p].set_yscale('log') ## Log scale
      plots[p].set_ylim([minm,maxm])
      for pos in ['top', 'left', 'right']:
        plots[p].spines[pos].set_visible(False)    
    plt.locator_params(nbins=4) 
    plots[p].set_xticks([.0])
    if xcnt: #################### Set one x-tick
      plots[p].set_xticklabels(['%i' %int(max(hist))])
    else:
      plots[p].set_xticklabels([''])
  if title:
    plt.suptitle(title, fontsize=20)
  plt.show()
  return




# Horizontal scatters


def hori_scatter(xdata, labelsin, title=None, axes=None, bounds=False,
                 showmean=True, switch=False, llog=False, counts=False,
                 jittery=False, shade=True, fill=False, bench=False,
                 eps=False, retplot=False, size=10): # retplot returns plt object
  """
  xdata is list of lists (distribution); jittery = y jitter of +/-0.25 units
  shade=IQR shading (bad for int data); bounds=number to normalize, False for 
  no bounds, true for mean +/- std bounds
  """
  if switch:
    for i in range(len(xdata)-1):
      xdata.append(xdata.pop(0))
      labelsin.append(labelsin.pop(0))
  altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink']
  if eps:
    colors = altcolors
  else:
    colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fig = plt.figure(figsize=(6,4), dpi=150)
  plots = [fig.add_subplot(1,len(xdata),i+1) for i in range(len(xdata))]
  if bounds == 'norm':
    X = []
    for x in xdata:
      X.append([i/max(x) for i in x])
    xdata = X
    minm, maxm = 0, 1.
  elif bounds is True:
    minm, maxm = np.inf, 0 # condition the data
    for x in xdata:
      if np.mean(x)-np.std(x) < minm:
        minm = np.mean(x)-np.std(x)
      if np.mean(x)+np.std(x) > maxm:
        maxm = np.mean(x)+np.std(x)
    if minm < 0:
      minm = 0.
  elif type(bounds) is list:
    minm, maxm = bounds[0], bounds[1]
    xdata = [[bounds[0] if i < bounds[0] else i for i in xdat ] for xdat in xdata]
    xdata = [[bounds[1] if i > bounds[1] else i for i in xdat ] for xdat in xdata]
  else:
    minm, maxm = 1, max([max(u) for u in xdata])
  for p in range(len(xdata)): # now plot
    xd = np.random.random(len(xdata[p]))
    if jittery:
      yd = np.random.randn(len(xdata[p]))/12.
    else:
      yd = np.zeros(len(xdata[p])) # Scatter!
    plots[p].scatter(xd, xdata[p]+yd, color=[colors[C[p]] if fill else 'none' for i in [1]][0],
                     edgecolor=colors[C[p]], alpha=0.3, s=size, linewidth=1,
                     ) # linewidth=3) #colors[C[p]]
    if showmean: # Showmeans and medians
      alph = {True: 1.0, False: 0.5}
      alph = alph[shade] # Set an alpha that changes depending on shading
      plots[p].plot([0,1], [np.mean(xdata[p]), np.mean(xdata[p])],
                    linewidth=1., c='black', alpha=1.) #colors[C[p]]
      plots[p].plot([0,1], [np.median(xdata[p]), np.median(xdata[p])], ':',
                    linewidth=1., c='black', alpha=1.)
      q25, q75 = np.percentile(xdata[p], [25, 75])
      plots[p].plot([0.5, 0.5], [q25, q75], color='black',  #colors[C[p]]
                       alpha=0.6, linewidth=1.)
    if bench:
      plots[p].plot([0,1], [bench, bench], c='purple', linewidth=1., alpha=0.3)
    if p == 0: #if first plot, show the axes
      plots[p].tick_params(axis='x',which='both',bottom='off',top='off',
                           labelbottom='off')
      if axes:
        plots[p].set_ylabel(axes[1], fontsize=15)
      plots[p].set_ylim([minm-.5, maxm+.5])
      if llog is True:
        plots[p].set_yscale('log'); plots[p].set_ylim([0, maxm]) ## Log scale
      if counts is True:
        plots[p].set_title('%i' %len(xdata[p]))
      for pos in ['top', 'right']: # Also hide these borders for all plots
        plots[p].spines[pos].set_visible(False)
    else:
      plots[p].tick_params(axis='x',which='both',bottom='off',top='off',
                           labelbottom='off')
      plots[p].get_yaxis().set_visible(False)
      if llog is True:
        plots[p].set_yscale('log') ## Log scale
      plots[p].set_ylim([minm-.5,maxm+.5])
      if counts is True:
        plots[p].set_title('%i' %len(xdata[p]))
      for pos in ['left', 'top', 'right']: # Also hide these borders for all plots
        plots[p].spines[pos].set_visible(False)
  if title:
    plt.suptitle(title, fontsize=20)
  if retplot:
    return fig
  plt.show()
  return




def scatter_distribution(xdata, labelsin, moreD, moreDmean=True, title=None, 
                 axes=None, bounds=False, showmean=True, switch=False, 
                 llog=False, counts=False, jittery=False,
                 fill=False, bench=False, forcebins=100, histbounds=None):
  """
  Plot a sample (see function below) to see what it looks like. Basically
  same as above except plots two distributions -- one as dist, one scatter.
  Good to compare against expected values. moreD is the scatter data.
  purp: refers to xdata (distribution), not the moreD (scatter)
  """
  if switch:
    for i in range(len(xdata)-1):
      xdata.append(xdata.pop(0))
      labelsin.append(labelsin.pop(0))
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fig = plt.figure(figsize=(6,4), dpi=150)
  plots = [fig.add_subplot(1,len(xdata),i+1) for i in range(len(xdata))]
  if bounds is True:
    minm, maxm = np.inf, 0 # condition the data
    for x in range(len(xdata)):
      if min([min(xdata[x]), min(moreD[x])]) < minm:
        minm = min([min(xdata[x]), min(moreD[x])])
      if max([max(xdata[x]), max(moreD[x])]) > maxm:
        maxm = max([max(xdata[x]), max(moreD[x])])
    if minm < 0:
      minm = 0.
  elif type(bounds) is list:
    minm, maxm = bounds[0], bounds[1]
    xdata = [[bounds[0] if i < bounds[0] else i for i in xdat ] for xdat in xdata]
    xdata = [[bounds[1] if i > bounds[1] else i for i in xdat ] for xdat in xdata]
  else:
    minm, maxm = 1, max([max(u) for u in xdata])
  for p in range(len(xdata)): ########## Plotting details ###########
    if type(forcebins) is list:
      b_e = forcebins
    else:
      if histbounds is None:
        b_e = np.linspace(minm, maxm, forcebins) # len/100 bins
      else:
        b_e = np.linspace(histbounds[0], histbounds[1], forcebins) # len/100 bins
    hist, _ = np.histogram(xdata[p], bins=b_e)
    plotbins = [(b_e[i]+b_e[i+1])/2. for i in range(len(b_e)-1)]
    try:   ############################## find appropriate bar width
      hgt = (maxm-minm)/len([i for i in hist if i != 0]) 
    except: # as high as there are filled hist elements
      print(b_e, min(xdata[p]), max(xdata[p]))
    hgt = plotbins[2]-plotbins[1] ######### histogram #########
    t25, t75 = np.percentile(xdata[p], [25, 75])
    barcols = [colors[C[p]] if t75 > plotbins[t] > t25 # In IQR
               else altcolors[C[p]] for t in range(len(plotbins))]
    for b in range(len(plotbins)):
      plots[p].barh(plotbins[b], hist[b]/max(hist), height=hgt, linewidth=0, alpha=0.8,
                    color=barcols[b], edgecolor=barcols[b]) #=colors[C[p]])
    xxd = np.random.random(len(moreD[p]))+0 # Here can change to 1 to shift right
    if jittery: #### MoreD #####
      yyd = np.random.randn(len(moreD[p]))/12.
    else:
      yyd = np.zeros(len(moreD[p])) ####### scatter #########
    moreD[p] = [maxm if y > maxm else y for y in moreD[p]]
    plots[p].scatter(xxd, moreD[p]+yyd, color=['none' if not fill else colors[C[p]] for i in [1]][0],
                     edgecolor=colors[C[p]], alpha=0.8, s=20) # linewidth=3) #colors[C[p]]
    if moreDmean:
      plots[p].plot([0,1+0], [np.mean(moreD[p]), np.mean(moreD[p])],
                    linewidth=1, c='black', alpha=0.9)
      more_iqr = np.percentile(moreD[p], [25, 75])
      plots[p].plot([0.5, 0.5], more_iqr, linewidth=1, color='black', alpha=0.9)
    if showmean:
      plots[p].plot([0,1+0], [np.mean(xdata[p]), np.mean(xdata[p])],
                    linewidth=1, c=altcolors[C[p]], alpha=1)
    if bench:
      plots[p].plot([0,1], [bench, bench], c='purple', linewidth=2, alpha=0.3)
    if p == 0: #if first plot, show the axes
      if axes:
        plots[p].set_ylabel(axes[1], fontsize=10)
      plots[p].set_ylim([minm, maxm+0.1])
      if llog is True:
        plots[p].set_yscale('log'); plots[p].set_ylim([0, maxm]) ## Log scale
      if counts is True:
        plots[p].set_title('%i' %len(xdata[p]))
      for pos in ['top', 'right']: # Also hide these borders for all plots
        plots[p].spines[pos].set_visible(False)
    else:
      plots[p].get_yaxis().set_visible(False)
      if llog is True:
        plots[p].set_yscale('log') ## Log scale
      plots[p].set_ylim([minm,maxm+0.1])
      if counts is True:
        plots[p].set_title('%i' %len(xdata[p]))
      for pos in ['left', 'top', 'right']: # Also hide these borders for all plots
        plots[p].spines[pos].set_visible(False)
    plt.locator_params(nbins=4) #################### Set one x-tick
    plots[p].set_xticks([0.5])
    plots[p].set_xticklabels(['%i' %int(max(hist))])
  if title:
    plt.suptitle(title, fontsize=20)
  plt.show()
  return






########################################################################
# Angles
########################################################################


def circular_hist(angles, labelsin, title=None, same=None, leg=True,
                  ninety=False, save=None, shade=True, eps=True,
                  only=False):
  """
  # IMPORT DEPENDENCIES FROM TOP. Same indicates same group, should be int.
  I.e.: GM same=0, LG same=1, etc.
  ninety=True if only want x-ticks plotted to 90, else 180 is default
  only: to display a single distribution, pass two distributions [fake, real]
  and set only=True.
  """
  if eps:
    altcolors = ['palegoldenrod', 'lightskyblue', 'lightgreen', 'lightpink']
  else:
    altcolors = ['purple' for i in range(4)]
  def to_radians(angs):
    return [i*np.pi/180. for i in angs]
  def r_bin(bins, target): # Return the target bin value, always start from below
    j = [i for i in bins]
    for i in j:
      if i > target:
        return i
  def angulize(angs, nbins=100): # Do everything for the plotting except plot
    if max(angs) > 2*np.pi:
      angs = to_radians(angs)
    rads, thetas_b = np.histogram(angs, bins=nbins)
    width = np.pi/(nbins)
    # Normalize hist height and center the bins
    rads = [i/max(rads) for i in rads]
    thetas = [(thetas_b[i]+thetas_b[i+1])/2. for i in range(len(thetas_b)-1)]
    q25, q75 = np.percentile(angs, [25, 75])
    b25, b75 = r_bin(thetas, q25), r_bin(thetas, q75)
    return angs, rads, thetas, width, b25, b75
  # If it's just one object, plot it simply; else nest the lists
  if type(angles[0]) is not list:
    angles = [angles]
  # Else, create the nested plots
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  ax = plt.subplot(111, polar=True)
  for A in range(len(angles)):
    angs, rads, thetas, width, t25, t75 = angulize(rm_nan(angles[A]))
    if same is None:
      s = A
    else:
      s = same
    barcols = [colors[s] if t75 > thetas[t] > t25 # In IQR
               else altcolors[s] for t in range(len(thetas))]
    if only and A == 0: # Only 
      bar = ax.bar([np.pi], [0.], color=barcols, width=width, bottom=2.+2*A)
    else:
      for t in range(len(thetas)):
        bar = ax.bar(thetas[t], rads[t], color=barcols[t],
                     edgecolor=barcols[t], width=width, bottom=2.+2*A)
    [b.set_facecolor(colors[s]) for b in bar.patches]
    [b.set_edgecolor(colors[s]) for b in bar.patches]
    if only is False or (only is True and A != 0): # Only stuff again
      #if shade: ### Shade
      #  iqr = ax.bar(np.linspace(t25, t75, 100), np.ones(100)*1.5, 
      #             width=np.pi/(200), bottom=2.+2*A)
      #  [i.set_facecolor(colors[s]) for i in iqr.patches]
      #  [i.set_alpha(0.3) for i in iqr.patches]
      #  [i.set_linewidth(0.) for i in iqr.patches]
      #else:
      #  ax.bar(t25, 0.1, width=t75-t25, bottom=2.+2*A + 0., 
      #         facecolor=altcolors[s], edgecolor='white')
      #  #ax.bar(t75, 1, width=np.pi/200, bottom=2.+2*A, facecolor=altcolors[s], edgecolor='none')
      mean = ax.bar(np.mean(angs), 1.5, width=np.pi/400, bottom=2.+2*A)
      med = ax.bar(np.median(angs), 1.25, width=np.pi/400, bottom=2.+2*A)
      k=['k','orange']
      for m in [med.patches[0], mean.patches[0]]:
        m.set_facecolor(k[[med.patches[0], mean.patches[0]].index(m)])
        m.set_linewidth(0.)
  khaki_patch = mpatches.Patch(color='darkkhaki',   # Legend
                label=labelsin[0])
  patches = [khaki_patch]
  if len(angles) > 1:
    royal_patch = mpatches.Patch(color='royalblue', 
                  label=labelsin[1])
    patches.append(royal_patch)
  if len(angles) > 2:
    forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[2])
    patches.append(forest_patch)
  if len(angles) > 3:
    lavender_patch = mpatches.Patch(color='tomato', 
                  label=labelsin[3])
    patches.append(lavender_patch)
  if type(same) is int:
    patches = [mpatches.Patch(color=colors[same], label=labelsin[0])]
  if leg:
    plt.legend(handles=patches, loc=4)
  if title:                                          # Title
    ax.set_title(title, fontsize=40)
  ax.set_yticklabels([])
  ax.set_yticks([])
  if ninety:
    ax.set_xticks([0,np.pi/6, 2*np.pi/6, np.pi/2])
    ax.set_xticklabels([0,30, 60, 90], fontsize=20)
  else:
    ax.set_xticks([0,np.pi/3, 2*np.pi/3, np.pi])
    ax.set_xticklabels([0,60,120,180], fontsize=20)
  if save is not None:
    plt.savefig(save)
  else:
    plt.show()
  return



##########################################################################
# Stats
##########################################################################


def stats_plots(V, labelsin, title=None, plotOne=None, axes=None):
  """
  4 plots of basic statistical properties. IC = intraclass correlation, 
  or the noise sources between the groups. 
  """
  import scipy.stats as stats
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  var = [np.var(i) for i in V]
  skew = [stats.skew(i) for i in V]
  kurt = [stats.kurtosis(i) for i in V]
  uniq = list(set(labelsin))
  v_sort = [[] for u in uniq] # Make a blank list, preparing for IC
  v_means = [[] for u in uniq] # v_means is a list of list of means for each cell of each type
  v_var, v_skew, v_kurt = [[] for u in uniq], [[] for u in uniq], [[] for u in uniq]
  for v in range(len(V)):
    i = uniq.index(labelsin[v])
    v_sort[i].append(V[v])
    v_means[i].append(np.mean(V[v]))
    v_var[i].append(np.var(V[v]))
    v_skew[i].append(stats.skew(V[v]))
    v_kurt[i].append(stats.kurtosis(V[v]))
  # ic = var_between^2 / (var_between^2 + var_within^2)  
  ic = []
  for v in range(len(uniq)):
    I = np.var(v_means[v])**2 / \
        (np.var(v_means[v])**2 + sum([np.var(i) for i in v_sort[v]])**2)
    ic.append([I])
  print(ic)
  group_means = [np.mean(k) for k in v_means] # group_means are the master means (only 4)
  master_ic = np.var(group_means)**2 / \
              (np.var(group_means)**2 + sum([np.var(i) for i in v_means])**2)
  print('Master IC for this set: %.5f' %master_ic)
  ## Plotting stuff
  fig = plt.figure()
  if plotOne is None:
    axs = [fig.add_subplot(221), fig.add_subplot(222), 
           fig.add_subplot(223), fig.add_subplot(224)]
    tits = ['Variance', 'Skew', 'Kurtosis', 'Intraclass correlation']
    plot_vars = [v_var, v_skew, v_kurt, ic]
  else:
    axs = [fig.add_subplot(221)]
    tits = ['']
    plot_vars = {'kurt': [v_kurt], 'var': [v_var], 'skew': [v_skew], 'ic': [ic]}
    plot_vars = plot_vars[plotOne]
  for a in axs: # For each plot
    for u in range(len(uniq)): # For each cell type
      a.scatter(np.ones(len(plot_vars[axs.index(a)][u]))*u, plot_vars[axs.index(a)][u], 
                c=colors[u], s=100, edgecolor='k', alpha=0.6)
      if axs.index(a) == 3:
        a.set_yticks([0,0.12,0.24])
      else:
        a.locator_params(axis='y', nbins=4)
      a.set_xticks([])
      a.set_title(tits[axs.index(a)])
  # Legend and title
  if plotOne is not None and axes is not None:
    if plotOne == 'ic':
      axs[0].set_yticks([0,0.12,0.24])
    axs[0].set_xlabel(axes[0], fontsize=15)
    axs[0].set_ylabel(axes[1], fontsize=15)
  if title is not None:
    plt.suptitle(title, fontsize=20)
  plt.show()
  




def line_scatter(xdata, ydata, labelsin=None, title=None, 
                 lines=True, groups=None, ax_titles=None):
  # a scatter plot with lines through the data
  # groups are a list of listed indices for which ydata belong together
  # i.e.: groups=[[0,1],[2,3]]
  colorlist = ['forestgreen','limegreen','royalblue','lightskyblue',
                'deeppink','orchid']
  markerlist = ['v','o','*','s']
  fig = plt.figure()
  ax = fig.add_subplot(111) # change this to 121 for legend outside fig
  
  if groups is not None:
    L = len(groups)
    print(L)
    cols = []
    for i in groups:
      cols.append(colorlist[i[0]])
      cols.append(colorlist[i[1]])
    marks = [markerlist[i] for i in range(2*L)]
  else:
    L = len(ydata)
    print(L)
    cols = [colorlist[i] for i in range(L)]
    marks = [markerlist[i] for i in range(L)]
    
  for i in range(len(ydata)):
    ax.scatter(xdata, ydata[i], c=cols[i], edgecolor=cols[i],
               marker=marks[i], s=40)
    ax.plot(xdata, ydata[i], c=cols[i], linewidth=3, alpha=0.2)
  
  forest_patch = mpatches.Patch(color='forestgreen', 
                  label=labelsin[1])
  green_patch = mpatches.Patch(color='limegreen',
                  label=labelsin[0])
  patches = [forest_patch, green_patch]
  if L > 2:
    royal_patch = mpatches.Patch(color='royalblue', 
                  label=labelsin[2])
    blue_patch = mpatches.Patch(color='lightskyblue',
                  label=labelsin[3])
    patches.append(royal_patch)
    patches.append(blue_patch)
  if L > 4:
    pink_patch = mpatches.Patch(color='deeppink', 
                  label=labelsin[4])
    orchid_patch = mpatches.Patch(color='orchid',
                  label=labelsin[5])
    patches.append(pink_patch)
    patches.append(orchid_patch) 
  # This can be useful for putting the legend outside the fig  vvv
  #plt.legend(handles=patches,loc='upper left')# bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.)
  if title:
    ax.set_title(title, fontsize=40)
  ax.spines['left'].set_position('zero')
  ax.spines['bottom'].set_position('zero')
  if ax_titles:
    ax.set_xlabel(ax_titles[0], fontsize=25)
    ax.xaxis.set_label_coords(0.5,0)
    ax.set_ylabel(ax_titles[1], fontsize=25)
    ax.yaxis.set_label_coords(0,0.5)
  plt.xticks([-3,3],fontsize=25)
  plt.yticks([-40,20],fontsize=25)
  plt.show()



def pretty_distribution(data, benchmark=None, bins=None, bars=False, title=None):
  # plots a pretty distribution of the Data, with or without a benchmark
  from scipy.interpolate import spline
  if bins is None:
    bins = len(data)/20
  hist, bin_edges = np.histogram(data, bins=bins)
  bin_e = [(bin_edges[i]+bin_edges[i+1])/2 for i in range(len(bin_edges)-1)]
  # interpolate for smoothed spline curve
  xnew = np.linspace(bin_e[0], bin_e[-1], 300)
  h_down = [np.mean(hist[i*4:i*4+4]) for i in range(int(len(hist)/4))] # maybe need -1 sometimes
  b_down = [np.mean(bin_e[i*4:i*4+4]) for i in range(int(len(bin_e)/4))]
  power_smooth = spline(b_down, h_down, xnew)
  fig = plt.figure()
  ax = fig.add_subplot(111)
  # draw curved line
  ax.plot(xnew, power_smooth, linewidth=2, c='royalblue')
  # fill under smoothed line or bars
  if not bars:
    ax.fill(xnew, power_smooth, facecolor='royalblue',alpha=0.4)
  else:
    ax.hist(data, bins=bins, edgecolor='lightskyblue', color='royalblue', alpha=0.4)
  if benchmark:
    ax.plot([benchmark, benchmark],[0,max(hist)*.8], color='deeppink', linewidth=3)
  if title:
    ax.set_title(title, fontsize=40)
  plt.show()



def pretty_line(xdata, ydata, labels=None, axlabels=None, title=None):
  if axlabels:
    if len(axlabels) != 2:
      print('Axlabels should have 2 items')
      axlabels=None
  colorlist = ['forestgreen','royalblue','deeppink','darkkhaki']
  fig = plt.figure()
  ax = fig.add_subplot(111)
  # use simple heuristics for determining if there are multiple y's
  if len(ydata) > 1 and len(ydata) <= 4 and len(ydata) != len(xdata):
    for y in range(len(ydata)):
      ax.plot(xdata, ydata[y], linewidth=2, c=colorlist[y])
      ax.plot(xdata, ydata[y], linewidth=3, c=colorlist[y], alpha=0.2)
  # if multiple x's and y's
  elif len(ydata) > 1 and len(ydata) <= 4 and len(ydata) == len(xdata):
    for y in range(len(ydata)):
      ax.plot(xdata[y], ydata[y], linewidth=2, c=colorlist[y])
      ax.plot(xdata[y], ydata[y], linewidth=3, c=colorlist[y], alpha=0.2)
  else:
    ax.plot(xdata, ydata, linewidth=2, c=colorlist[0])
    ax.plot(xdata, ydata, linewidth=3, c=colorlist[0], alpha=0.2)
  # simple legend !
  if labels:
    if len(labels) != len(ydata):
      print('Num labels must equal num ydata!')
    else:
      patches = []
      for y in range(len(labels)):
        patch = mpatches.Patch(color=colorlist[y], label=labels[y])
        patches.append(patch)
      ax.legend(handles=patches)
  if axlabels:
    ax.set_xlabel(axlabels[0], fontsize=15)
    ax.set_ylabel(axlabels[1], fontsize=15)
  if title:
    ax.set_title(title, fontsize=25)
  ax.set_ylim([-0.5,3.5])
  ax.set_xlim([-0.5,5.5])
  plt.show()



def grid_fits(xvals, yvals, labelsin):
  # A grid of plots, order = GM, LG, LP, PD
  colors = ['darkkhaki', 'royalblue', 'forestgreen','tomato']
  L = list(np.unique(labelsin))
  C = [L.index(i) for i in labelsin]
  fig = plt.figure()
  # subs = [0,4,1,2,3,8,12,9,13,5,6,7,14,10,15,11] # row major
  # subs = [0,2,1,4,5,8,10,9,11,3,6,7,14,12,15,13]
  subs = [0,1,4,8,12,2,3,6,7,5,9,13,11,10,15,14]
  for u in range(16):
    # First plot the data
    ax = fig.add_subplot(4,4,subs[u]+1)
    ax.scatter([i+np.random.random(1)*.5 for i in xvals[u]],
                yvals[u], facecolor='none', 
                edgecolor=colors[C[u]], s=20, lw=1, alpha=0.4)
    # Also plot the fit
    b, a, r, p, err = stats.linregress(xvals[u], yvals[u])
    meany = [a+b*y for y in yvals[u]]
    low_bound = np.linspace(min(meany)-np.sqrt(err), max(meany)-np.sqrt(err), 10)
    upp_bound = np.linspace(min(meany)+np.sqrt(err), max(meany)+np.sqrt(err), 10)
    ax.plot(xvals[u], meany, color=colors[C[u]], linewidth=1, alpha=0.8)
    ax.fill_between(np.linspace(0,max(xvals[u]),10), low_bound, upp_bound, facecolor=colors[C[u]],
                    color=colors[C[u]], alpha=0.4)
    ax.set_xlim([-1,max(xvals[u])+1])
    ax.set_ylim([-.1, max(yvals[u])+.1])
    ax.set_xticks([max(xvals[u])])
    ax.set_yticks([1])
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
  plt.tight_layout()
  plt.show()
  return
  
    

########################################################################
# Other stuff
########################################################################


def pretty_dendrogram(nodes):
  import numpy as np
  from scipy.cluster.hierarchy import linkage
  import matplotlib.pyplot as plt
  #
  linkage_matrix = linkage(nodes, "single")
  #
  plt.figure()
  plt.clf()
  #
  show_leaf_counts = True
  ddata = augmented_dendrogram(linkage_matrix,
                 color_threshold=1,
                 p=600,
                 truncate_mode='lastp',
                 show_leaf_counts=show_leaf_counts,
                 )
  plt.title("show_leaf_counts = %s" % show_leaf_counts)

  plt.show()



def pretty_skeleton(geo, color=None):
  """
  Draw a pretty skeleton, ignores z-values.
  """
  # Generate the pairs to be connected
  branchpts = []
  for b in geo.branches: 
    branchpts.append([[n.x, n.y] for n in b.nodes])
  #
  # Plotting
  plt.figure(figsize=(3,6), dpi=175)
  if color is None:
    color='black'
  for b in branchpts: # Plot the background skeleton first
    for s in range(len(b)-1):
      plt.plot([b[s][0], b[s+1][0]],
               [b[s][1], b[s+1][1]], color=color, alpha=0.3)
  plt.show()
  
  
  
########################################################################
# Helper functions
########################################################################


from scipy.cluster.hierarchy import dendrogram
import matplotlib.pyplot as plt
import networkx as nx


def augmented_dendrogram(*args, **kwargs):

    ddata = dendrogram(*args, **kwargs)

    if not kwargs.get('no_plot', False):
        for i, d in zip(ddata['icoord'], ddata['dcoord']):
            x = 0.5 * sum(i[1:3])
            y = d[1]
            plt.plot(x, y, 'ro')
            plt.annotate("%.3g" % y, (x, y), xytext=(0, -8),
                         textcoords='offset points',
                         va='top', ha='center')

    return ddata


## Ellipse helpers! ##

from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import art3d
from mpl_toolkits.mplot3d import proj3d
import numpy as np

def rotation_matrix(v1,v2):
  """
  Calculates the rotation matrix that changes v1 into v2.
  """
  v1/=np.linalg.norm(v1)
  v2/=np.linalg.norm(v2)
  cos_angle=np.dot(v1,v2)
  d=np.cross(v1,v2)
  sin_angle=np.linalg.norm(d)
  if sin_angle == 0:
    M = np.identity(3) if cos_angle>0. else -np.identity(3)
  else:
    d/=sin_angle
    eye = np.eye(3)
    ddt = np.outer(d, d)
    skew = np.array([[    0,  d[2],  -d[1]],
                  [-d[2],     0,  d[0]],
                  [d[1], -d[0],    0]], dtype=np.float64)
    M = ddt + cos_angle * (eye - ddt) + sin_angle * skew
  return M


def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'):
  """
  Transforms a 2D Patch to a 3D patch using the given normal vector.
  The patch is projected into they XY plane, rotated about the origin
  and finally translated by z.
  """
  if type(normal) is str: #Translate strings to normal vectors
    index = "xyz".index(normal)
    normal = np.roll((1,0,0), index)
  path = pathpatch.get_path() #Get the path and the associated transform
  trans = pathpatch.get_patch_transform()
  path = trans.transform_path(path) #Apply the transform
  pathpatch.__class__ = art3d.PathPatch3D #Change the class
  pathpatch._code3d = path.codes #Copy the codes
  pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color    
  verts = path.vertices #Get the vertices in 2D
  M = rotation_matrix(normal,(0, 0, 1)) #Get the rotation matrix
  pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts])


def pathpatch_translate(pathpatch, delta):
  """
  Translates the 3D pathpatch by the amount delta.
  """
  pathpatch._segment3d += delta


def eigsorted(cov):
  """
  Necessary for the 2d ellipses created in pretty_2d scatter plot.
  """
  vals, vecs = np.linalg.eigh(cov)
  order = vals.argsort()[::-1]
  return vals[order], vecs[:,order]



def dist3(pt0, pt1):
  """
  In theory this can handle node objects.
  """
  if type(pt0) is not list:
    try:
      pt0 = [pt0.x, pt0.y, pt0.z] # Try for a node
    except:
      print('what do I do with ' + pt0)
  if type(pt1) is not list:
    try:
      pt1 = [pt1.x, pt1.y, pt1.z] # Try for a node
    except:
      print('what do I do with ' + pt1)
  try:
    return math.sqrt(sum([(pt0[i]-pt1[i])**2 for i in range(3)]))
  except:
    print('dimension mismatch')
    print(pt0, pt1)
  return None






########################################################################

if __name__ == "__main__":
  print('Module is used interactively.')