docs: def hierarchical_timecourse(
          df,
          cluster_cols=True,
          cluster_rows=False,
          n_col_clusters=False,
          n_row_clusters=False,
          fcol=None,
          z_score=0,
          method='ward',
          cmap=cm.PuOr_r,
          return_clusters=False,
          rdistance_fn=distance.pdist,
          cdistance_fn=distance.pdist,
          xlabel='Timepoint',
          ylabel='log$_2$ Fold Change'
      ):
      """
      Hierarchical clustering of samples across timecourse experiment.
  
      Peform a hiearchical clustering on a pandas DataFrame and display the resulting clustering as a
      timecourse density plot.
  
      Samples are z-scored along the 0-axis (y) by default. To override this use the `z_score` param with the axis to `z_score`
      or alternatively, `None`, to turn it off.
  
      If a `n_col_clusters` or `n_row_clusters` is specified, this defines the number of clusters to identify and highlight
      in the resulting heatmap. At *least* this number of clusters will be selected, in some instances there will be more
      if 2 clusters rank equally at the determined cutoff.
  
      If specified `fcol` will be used to colour the axes for matching samples.
  
      :param df: Pandas ``DataFrame`` to cluster
      :param n_col_clusters: ``int`` the ideal number of highlighted clusters in cols
      :param n_row_clusters: ``int`` the ideal number of highlighted clusters in rows
      :param fcol: ``dict`` of label:colors to be applied along the axes
      :param z_score: ``int`` to specify the axis to Z score or `None` to disable
      :param method: ``str`` describing cluster method, default ward
      :param cmap: matplotlib colourmap for heatmap
      :param return_clusters: ``bool`` return clusters in addition to axis
      :return: matplotlib axis, or axis and cluster data
      """
  
      dfc, row_clusters, row_denD, col_clusters, col_denD, edges = _cluster(df,
         cluster_cols=cluster_cols, cluster_rows=cluster_rows, n_col_clusters=n_col_clusters,
         n_row_clusters=n_row_clusters, z_score=z_score, method='ward',
         rdistance_fn=rdistance_fn, cdistance_fn=cdistance_fn
                                                 )
  
      # FIXME: Need to apply a sort function to the DataFrame to order by the clustering
      # so we can slice the edges.
      dfh = dfc.iloc[row_denD['leaves'], col_denD['leaves']]
      dfh = dfh.mean(axis=0, level=[0, 1])
  
      vmax = np.max(dfh.values)
      color = ScalarMappable(norm=Normalize(vmin=0, vmax=vmax), cmap=viridis)
  
      fig = plt.figure(figsize=(12, 6))
  
      edges = [0] + edges + [dfh.shape[1]]
  
      for n in range(len(edges) - 1):
          ax = fig.add_subplot(2, 4, n + 1)
  
          dfhf = dfh.iloc[:, edges[n]:edges[n + 1]]
          xpos = dfhf.index.get_level_values(1)
  
          mv = dfhf.mean(axis=1)
          distances = [distance.euclidean(mv, dfhf.values[:, n]) for n in range(dfhf.shape[1])]
          colors = [color.to_rgba(v) for v in distances]
          order = np.argsort(distances)[::-1]
          for y in order:
              ax.plot(xpos, dfhf.values[:, y], c=colors[y], alpha=0.5, lw=1)  # dfhf.index.get_level_values(1),
              ax.set_xticks(xpos)
  
          if n > 3:
              ax.set_xticklabels(xpos)
              ax.set_xlabel(xlabel)
          else:
              ax.set_xticklabels([])
  
          if n % 4 != 0:
              ax.set_yticklabels([])
          else:
              ax.set_ylabel(ylabel)
  
          ax.set_ylim((-3, +3))
  
      fig.subplots_adjust(hspace=0.15, wspace=0.15)
  
      if return_clusters:
          return fig, dfh, edges
  
      else:
          return fig