barplotwithsamples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bar plot with confidence interval and samples\n",
    "This notebook plots a box plot or bar plot from per-subject values extracted from REX (or another software) and plots the 90% confidence interval and all per-subject values as scatterplots.\n",
    "\n",
    "If in addition you supply the ROIs nifti images (one per ROI), you will in addition get a nice glass brain image with all ROIs, and the bar plot will also include additional infos (ROI center coordinates + atlas regions names covered by ROI + same color of text as the glass brain).\n",
    "\n",
    "To read boxplots, here is an [excellent tutorial](https://stackoverflow.com/questions/38794406/why-is-matplotlibs-notched-boxplot-folding-back-on-itself).\n",
    "\n",
    "Version 1.9.3\n",
    "\n",
    "TODO:\n",
    "* allow to vary boxplot width according to sample size, larger = more samples: https://blog.minitab.com/blog/statistics-and-quality-data-analysis/how-to-think-outside-the-boxplot\n",
    "* implement alternative visualization matplotlib.pyplot.errorbar: https://stackoverflow.com/questions/31644018/python-boxplot-showing-means-and-confidence-intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "# BEWARE: autoreload works on functions and on general code, but NOT on new class methods:\n",
    "# if you add or change the name of a method, you have to reload the kernel!\n",
    "# also it will fail if you use super() calls in the classes you change\n",
    "# ALSO AUTORELOAD SHOULD BE THE FIRST LINE EVER EXECUTED IN YOUR IPYTHON NOTEBOOK!!!\n",
    "\n",
    "# Profilers:\n",
    "# http://pynash.org/2013/03/06/timing-and-profiling/\n",
    "# http://mortada.net/easily-profile-python-code-in-jupyter.html\n",
    "# use %lprun -m module func(*args, **kwargs)\n",
    "try:\n",
    "    %load_ext line_profiler\n",
    "    %load_ext memory_profiler\n",
    "except ImportError as exc:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure inside IPython Notebook (must be called before any import of matplotlib, direct or indirect!)\n",
    "%matplotlib inline\n",
    "\n",
    "import glob\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as pltcol\n",
    "import numpy as np\n",
    "import nibabel as nib\n",
    "import pandas as pd\n",
    "import textwrap\n",
    "from nilearn import image\n",
    "from nilearn import plotting\n",
    "try:\n",
    "    from adjustText import adjust_text  # you will need this lib if you want to plot subjects ids without overlap: pip install adjustText\n",
    "except ImportError as exc:  # else we will just plot the subjects labels without any adjustment (so the position can overlap with the points and with each others)\n",
    "    adjust_text = None\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARAMETERS - EDIT ME\n",
    "groups = [8, 2, 3, 4]  # set the number of items/values for each group\n",
    "groups_labels = ['CTR', 'LIS', 'UWS', 'E/MCS']  # set names for each group\n",
    "groups_order = range(len(groups))  # in which order we plot each bar. Leave as is to plot bars in the original order, or provide list of group nb to shuffle as you want (this is only a display parameter, it does not change the results)\n",
    "#subjects_labels = range(1, sum(groups)+1)  # set ID for each subject (can be any string, all that matters is that it is the same order as the input values: first value here is the first value in imported rex data file, etc)\n",
    "#subjects_labels = range(1, groups[0]+1)*len(groups)  # use this for within-subject analyses, where you want the labels to be the same range in both bars\n",
    "show_subjects_labels = True  # show label for each subject's point?\n",
    "show_subjects_paired = False  # paired/longitudinal analysis, both groups contain in fact the same subjects but in different conditions, enable this option to restart counting subjects ids to 1 for each condition\n",
    "rex_data_filepath = ['testconjunc.cluster00%i.rex.data.txt' % i for i in range(1,6)]  # single path (as string) or list of paths to the txt files containing the values to plot, one value per line for each subject. The path needs to be absolute or relative to this script's path. Alternative: use glob.glob(r'path/to/files/roi*.txt')[:-1]\n",
    "rex_maps_filepath = ['testconjunc.cluster00%i.rex.roi.img' % i for i in range(1,6)]  # specify here the filenames, in the same order as the input text files (can be a single string or a list of strings), to plot the ROIs on a brain using nilearn\n",
    "atlas_cmap = plt.cm.Set1  # colormap to use for ROIs (only if map provided), can also specify a specific set of colors like this, with one color for each cluster (here 3): pltcol.ListedColormap(['r', 'g', 'b'], name='from_list', N=None)\n",
    "atlas_name = 'aal2'  # name of the atlas to draw the labels, can be either aal2, anatomytoolbox or custom (for custom you need to provide a atlas/custom/atlas.nii file and a atlas/custom/atlas.txt file with the labels and the row number being the value for each ROI in the nifti file)\n",
    "show_colorbar = False  # display the colorbar on the brain image?\n",
    "show_centroid = True  # display a black dot in the spatial center of each cluster? (only if maps are provided)\n",
    "show_atlas_labels = False  # display below the plot the region labels covered by each ROI using the specified atlas\n",
    "dpi_resolution = 300  # dpi resolution when saving on-disk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------\n",
    "## Start of core routines, run all from here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_maps(list_imgs, voxel_threshold=None):\n",
    "    if voxel_threshold is None:\n",
    "        voxel_threshold = 0.0001 # minimum threshold to consider as a voxel and not just background noise (because background voxels can be 0.000001 for example), can be float or str ('1%' to give a percentage). TODO: autodetect minimum value (can be -4, 0.02, etc) as the background and use it as the threshold value.\n",
    "    # Load masks and resample to first\n",
    "    imgs = []\n",
    "    for img in list_imgs:\n",
    "        im = image.load_img(img)\n",
    "        if imgs:\n",
    "            if im.shape != imgs[0].shape:\n",
    "                im = image.resample_to_img(im, imgs[0])\n",
    "        im = image.threshold_img(im, voxel_threshold)\n",
    "        imgs.append(im)\n",
    "    return imgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SANITY CHECKS\n",
    "#if show_subjects_labels:\n",
    "    #if len(subjects_labels) != sum(groups):\n",
    "        #raise(ValueError('subjects_labels does not contain the same number of subjects as groups!'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check filepaths are lists and not just string\n",
    "if isinstance(rex_data_filepath, str):\n",
    "    rex_data_filepath = [rex_data_filepath]\n",
    "if isinstance(rex_maps_filepath, str):\n",
    "    rex_maps_filepath = [rex_maps_filepath]\n",
    "# Loading data from Rex csv\n",
    "dfraw = []\n",
    "nb_rois = len(rex_data_filepath)\n",
    "for fpath in rex_data_filepath:\n",
    "    dfraw.append(pd.read_csv(fpath, index_col=False, header=None, squeeze=True))\n",
    "    if len(dfraw[-1]) != sum(groups):\n",
    "        raise(ValueError('the number of values in the provided txt file is not the same as the supplied groups count (ie, you did not specify the correct number of subjects!): %i vs %i, please check your parameters, or the input file: %s' % (sum(groups), len(dfraw[-1]), fpath)))\n",
    "    print(dfraw[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the values for each group in a separate dataframe\n",
    "df_g = []\n",
    "for i in range(nb_rois):\n",
    "    start = 0\n",
    "    for g in groups:\n",
    "        # Extract current subset of values for the current number of subjects (g is the number of subjects in current group)\n",
    "        df_g.append(dfraw[i][start:start+g])\n",
    "        # Next start after where we left here (to extract next group's values)\n",
    "        start = start+g\n",
    "        if show_subjects_paired:\n",
    "            # Reset indices to restart at 0 the subjects counting because it's the same subjects (longitudinal analysis)\n",
    "            df_g[-1].index = list(range(0,len(df_g[-1])))\n",
    "        #else:\n",
    "            # Continue the numbering (different groups of subjects)\n",
    "\n",
    "df_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "def comp_ci(a, conf_interval=0.9):\n",
    "    '''Calculates the 90% confidence interval from a vector (from the mean).\n",
    "    From the excellent SO answer by Ulrich Stern: https://stackoverflow.com/a/34474255/1121352'''\n",
    "    return scipy.stats.t.interval(conf_interval, len(a)-1, loc=np.mean(a), scale=scipy.stats.sem(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_cluster_center(im, mricron=False):\n",
    "    from nilearn.image.resampling import reorder_img, coord_transform\n",
    "    # Project coordinates to reduced space\n",
    "    if mricron:\n",
    "        im2 = reorder_img(im, resample='continuous')  # DO NOT USE: this will convert to MRIcron coordinates space (ie, [0, 100]) but it will mess things up for nilearn!\n",
    "        # Get indices of nonzero values\n",
    "        matches = im2.get_data().nonzero()\n",
    "    else:\n",
    "        # Get indices of nonzero values\n",
    "        matches = im.get_data().nonzero()\n",
    "    # Compute the euclidian middle of the cluster, from the nonzero values indices (= coordinates)\n",
    "    center = np.mean(matches, axis=1)\n",
    "    # Project center coordinates to brain space (ie, instead of [0, 100] range, it will be [-50, 50] -- I picked these numbers randomly, you see the idea)\n",
    "    if not mricron:\n",
    "        center = coord_transform(center[0], center[1], center[2], im.affine)  # disable this to get MRIcron space\n",
    "    return center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot ROIs on glass brain images!\n",
    "if rex_maps_filepath:\n",
    "    imgs = load_maps(rex_maps_filepath)\n",
    "\n",
    "    #plotting.plot_prob_atlas(imgs, view_type=\"filled_contours\",\n",
    "    #                    title=\"lala\", colorbar=True, cut_coords=(0,0,0), draw_cross=True, cmap=pltcol.ListedColormap(['b', 'g', 'r', 'c', 'm'], name='from_list', N=None))\n",
    "    #plotting.plot_roi(imgs[0])\n",
    "\n",
    "    centers = []\n",
    "    fig = plotting.plot_glass_brain(None, cmap=atlas_cmap, alpha=0.5)  # initialize the glass brain images\n",
    "    for c, im in enumerate(imgs):\n",
    "        # For each ROI\n",
    "        # Get the center (to plot the marker)\n",
    "        center = find_cluster_center(im)\n",
    "        centers.append(center)\n",
    "        print('Center found at: ' + str(center))\n",
    "        # Assign a unique value to this cluster's voxels (to get a different color)\n",
    "        imdata = im.get_data()  # Convert to numpy structure\n",
    "        imdata[imdata != 0] = c+1  # Assign unique value\n",
    "        im2 = nib.Nifti1Image(imdata, affine=im.affine)  # convert back to a nifti file in-memory to supply to nilearn\n",
    "        # Show colorbar?\n",
    "        cbar = False\n",
    "        if show_colorbar and c == (nb_rois-1):\n",
    "            # Can only plot the colorbar at the last iteration, else nilearn will spit an error (cannot use multiple colorbars)\n",
    "            cbar = True\n",
    "        # Plot the clusters\n",
    "        fig.add_overlay(im2, vmin=1, vmax=nb_rois, cmap=atlas_cmap, colorbar=cbar)\n",
    "        # Plot the markers (clusters' centers)\n",
    "        if show_centroid:\n",
    "            fig.add_markers([center], marker_color=['k'], marker_size=20)\n",
    "\n",
    "# Save figure\n",
    "fig.savefig('rois_glass_brain.png', dpi=dpi_resolution)\n",
    "print('Image saved in rois_glass_brain.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting, datasets\n",
    "\n",
    "def get_atlas_label(atlas, region_idx):\n",
    "    \"\"\"\n",
    "    Get atlas label for one specific region index\n",
    "    \"\"\"\n",
    "    return atlas['labels'][atlas['indices'].index(str(region_idx))]\n",
    "\n",
    "def get_atlas_labels(imgs, atlas_choice='aal2', verbose=False):\n",
    "    \"\"\"\n",
    "    Get the list of atlas regions covered by clusters, from a list of nifti maps loaded in-memory via nibabel\n",
    "    atlas_choice is optional, can be 'aal2' or 'SPM12' or 'anatomytoolbox' or 'custom'\n",
    "    For custom, you can provide an atlas/custom/atlas.nii file and an atlas/custom/atlas.txt file containing each label (the index will be the row number for the label)\n",
    "    \"\"\"\n",
    "    voxel_threshold = 0.0001 # minimum threshold to consider as a voxel and not just background noise (because background voxels can be 0.000001 for example), can be float or str ('1%' to give a percentage). TODO: autodetect minimum value (can be -4, 0.02, etc) as the background and use it as the threshold value.\n",
    "    # Atlas\n",
    "    if atlas_choice == 'anatomytoolbox':\n",
    "        atlas_path = 'masks\\AnatomyToolbox_Atlas_Map.nii'  # TODO: build atlas variable with all infos and data (labels, indices, nib niftiimage with affine etc)\n",
    "    elif atlas_choice == 'aal2' or atlas_choice == 'SPM12':\n",
    "        atlas = datasets.fetch_atlas_aal(version='SPM12', data_dir='atlas')\n",
    "    elif atlas_choice == 'custom':\n",
    "        import re\n",
    "        from sklearn.datasets.base import Bunch\n",
    "        atlas_file = r'atlas/custom/atlas.nii'\n",
    "        with open(r'atlas/custom/atlas.txt', 'r') as f:\n",
    "            atlas_labels = f.read()\n",
    "            atlas_labels = re.sub(r'\\(.+\\)', '', atlas_labels)\n",
    "            atlas_labels = map(str.strip, atlas_labels.split('\\n'))\n",
    "        atlas_indices = map(str, range(1, len(atlas_labels)+1))\n",
    "        atlas = Bunch(maps=atlas_file, labels=atlas_labels, indices=atlas_indices, description='Custom atlas')\n",
    "    else:\n",
    "        raise ValueError('Supplied atlas %s is not recognized!' % atlas_choice)\n",
    "    # If the indices were not stored, probably it follows a simple numbering scheme from 1, so create it\n",
    "    if not 'indices' in atlas:\n",
    "        atlas['indices'] = map(str, range(1, len(atlas['labels'])+1))  # expects the indices to be str\n",
    "\n",
    "    # Show some infos about atlas\n",
    "    atlas_im = image.load_img(atlas.maps)\n",
    "    if verbose:\n",
    "        print('Atlas shape: %s' % str(atlas_im.shape))\n",
    "    nb_regions = len(np.unique(atlas_im.get_data()))- 1\n",
    "    if verbose:\n",
    "        print('%i regions in this atlas: %s' % (nb_regions, str(np.unique(atlas_im.get_data())))) # 48 regions because 0 is background\n",
    "        print('%i labels' % len(atlas['labels']))\n",
    "        print('%i indices: %s' % (len(atlas['indices']), atlas['indices']))\n",
    "        print(atlas.keys())\n",
    "\n",
    "    # Resample masks to atlas size\n",
    "    imgs2 = []\n",
    "    for img in imgs:\n",
    "        if img.shape != atlas_im.shape:\n",
    "            img = image.resample_to_img(img, atlas_im)\n",
    "        img = image.threshold_img(img, voxel_threshold)\n",
    "        imgs2.append(img)\n",
    "    imgs = imgs2\n",
    "    del imgs2\n",
    "    imgs[0].shape\n",
    "\n",
    "    # Extract activated atlas brain regions for each mask\n",
    "    maps_regions = []\n",
    "    maps_regions_idxs = []\n",
    "    maps_regions_count = []\n",
    "    for img in imgs:\n",
    "        # Extract only non zeros voxels indices from mask\n",
    "        im_data = img.get_data()\n",
    "        #np.extract(im_data>0, im_data)\n",
    "        vox_thres = np.nonzero(im_data)\n",
    "        # Compare with atlas regions to extract region indices\n",
    "        atlas_data = atlas_im.get_data()\n",
    "        region_indices = set()\n",
    "        region_count = {}\n",
    "        for x in zip(*vox_thres): # walk through all non zero voxels of mask\n",
    "            region_idx = atlas_data[x] # get equivalent voxel from atlas\n",
    "            if region_idx != 0: # if not background\n",
    "                # Append region index into the set (so that they are unique)\n",
    "                region_indices.add(region_idx)\n",
    "                # Increase the count of voxels activated in this region\n",
    "                region_label = get_atlas_label(atlas, region_idx)\n",
    "                if region_label not in region_count:\n",
    "                    region_count[region_label] = 0\n",
    "                region_count[region_label] += 1\n",
    "        if 0 in region_indices:\n",
    "            region_indices.remove(0) # remove background, not part of the atlas labels\n",
    "        if verbose:\n",
    "            print('Atlas indices of brain regions activated in current mask: %s' % str(sorted(region_indices)))\n",
    "        # Extract brain region names from atlas that are present in this mask\n",
    "        matching_idxs = [int(idx) in region_indices for idx in atlas['indices']]\n",
    "        map_brain_regions = filter(None, [label if match else None for label, match in zip(atlas['labels'], matching_idxs)])\n",
    "        maps_regions.append(list(map_brain_regions))\n",
    "        maps_regions_idxs.append(region_indices)\n",
    "        maps_regions_count.append(region_count)\n",
    "\n",
    "    return maps_regions, maps_regions_idxs, maps_regions_count\n",
    "\n",
    "if rex_maps_filepath:\n",
    "    maps_regions, _, _ = get_atlas_labels(imgs, atlas_choice=atlas_name)\n",
    "    print('Found the following regions covered by the ROIs clusters:')\n",
    "    for i in range(nb_rois):\n",
    "        print('ROI %i: %s' % (i, ', '.join(maps_regions[i])))\n",
    "else:\n",
    "    print('No map found! Please check your rex_maps_filepath parameter.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot!\n",
    "\n",
    "# Plotting parameters\n",
    "ylim = None  # limit y axis to these values. Set to None to use default limits automatically detected by matplotlib.\n",
    "width = 0.6  # width of the bars - do not change, it is an internal parameter and does not impact the visualization\n",
    "margins = [0.8, 0.64] # margin we will set and account for to consistently scale the bars whatever the figsize is, in inches\n",
    "figsize = None # [width*len(groups)*nb_rois+sum(margins), 3]  # figure size, in inches, set to None to use default (don't forget to add +width because the bars are in the middle of tickes, thus to have room you must account for +1 ghost bar on the outer edges of the fig)\n",
    "colors = ['b', 'r', 'm', 'y', 'c', 'g', 'k']\n",
    "bars_color_by_roi = True  # group bars colors by ROI, or differenciate by group/bar? NOTE: can only work if a map is provided\n",
    "bars_alpha = 0.5  # bars coloring transparency, allows to more easily see the subjects points on top\n",
    "ylabel = 'Effect sizes'\n",
    "roi_labels_yoffset = -0.17  # how much shift to move the roi labels downward?\n",
    "roi_labels_maxlen = 8  # this is the maximum length (nb of characters) per tick/group\n",
    "roi_labels_fontsize = 11  # roi labels fontsize\n",
    "plot_params = {'legend.fontsize': 15,  # various plot parameters\n",
    "         #'figure.figsize': figsize,\n",
    "         'axes.labelsize': 20,\n",
    "         'axes.titlesize':30,\n",
    "         'xtick.labelsize':12,\n",
    "         'ytick.labelsize':15}\n",
    "subjects_fontsize = 12  # fontsize for the labels per subjects points\n",
    "subjects_labels_alpha = 0.0 #0.5  # how much transparency for the subjects labels\n",
    "confidence_interval = 90 # use a value between 0 and 100 to define the confidence interval to use\n",
    "plotkind = 'violin'  # use a 'box' plot or a 'bar' plot or a 'violin' plot?\n",
    "boxplot_filled = True  # fill boxplot with color or only the edges?\n",
    "boxplot_ci_bootstrap = None # if an integer > 0 (advised: 1000 to 10000), don't use the specified confidence interval, but calculate the 95% confidence interval by bootstrapping. To disable, set to None.\n",
    "boxplot_whiskers = 1.5 # whiskers for the boxplot, see https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.axes.Axes.boxplot.html\n",
    "\n",
    "# Prepare the ROI colorbar range normalization\n",
    "if rex_maps_filepath:\n",
    "    norm = pltcol.Normalize(vmin=1,vmax=nb_rois)  # need to normalize the values we will input to the colormap to be onpoint with nilearn\n",
    "# Plotting each bar\n",
    "plt.rcParams.update(plot_params)  # update the plot parameters right away, else they won't be updated until next cell execution\n",
    "ticks = np.arange(width, width+(width*len(groups)*nb_rois), width)  # do not modify this\n",
    "fig, ax = plt.subplots()\n",
    "if figsize:\n",
    "    fig.set_size_inches(figsize[0], figsize[1], forward=True)\n",
    "texts = []\n",
    "for roi_id in range(nb_rois):\n",
    "    last_j = 0\n",
    "    for i, gi in enumerate(groups_order):\n",
    "        # Get the data for the selected group\n",
    "        dg = df_g[len(groups)*roi_id + gi]\n",
    "        # Draw bars with error bar\n",
    "        bar_color = colors[i] if not bars_color_by_roi or not rex_maps_filepath else atlas_cmap(norm(roi_id+1))\n",
    "        xpos = ticks[len(groups)*roi_id+i]+(float(width)/2)\n",
    "        if plotkind == 'box':\n",
    "            # Boxplot\n",
    "            bp = ax.boxplot(dg, positions=[xpos], patch_artist=boxplot_filled, notch=True, vert=True,\n",
    "                            conf_intervals=[comp_ci(dg, float(confidence_interval)/100)], bootstrap=boxplot_ci_bootstrap,\n",
    "                            whis=boxplot_whiskers,\n",
    "                            widths=width*0.9, showmeans=True, showfliers=True, sym='_',\n",
    "                            flierprops={'markeredgecolor': bar_color, 'markeredgewidth': 2})\n",
    "            ax.axhline(y=0, color='grey', linestyle='dotted')  # add a horizontal line at 0 to separate positive from negative effect size\n",
    "            # Change color of outline and fill\n",
    "            bp['boxes'][0].set(color=bar_color)\n",
    "            if boxplot_filled:\n",
    "                bp['boxes'][0].set_facecolor(bar_color)\n",
    "                bp['boxes'][0].set_alpha(bars_alpha)\n",
    "            # Change color of whiskers and caps\n",
    "            for whisk in bp['whiskers']:\n",
    "                whisk.set(color='black', linewidth=1)\n",
    "            for cap in bp['caps']:\n",
    "                cap.set(color='black', linewidth=1)\n",
    "            # Median bar\n",
    "            bp['medians'][0].set(color='black')\n",
    "            # Mean point\n",
    "            bp['means'][0].set(color='black', markerfacecolor='black', markeredgecolor='black')\n",
    "        elif plotkind == 'bar':\n",
    "            # Barplot\n",
    "            ax.bar(xpos, dg.mean(), width=width, yerr=(dg.mean() - comp_ci(dg, float(confidence_interval)/100)[1]),\n",
    "                   alpha=bars_alpha, color=bar_color,\n",
    "                   error_kw={'ecolor': 'k', 'elinewidth': 1, 'capsize': 15, 'capthick': 1, 'barsabove': False})\n",
    "        elif plotkind == 'violin':\n",
    "            # Violinplot\n",
    "            parts = ax.violinplot(dg, positions=[xpos], vert=True, showmeans=True, showextrema=False, bw_method=0.5)\n",
    "            for pc in parts['bodies']:\n",
    "                pc.set_facecolor(bar_color)\n",
    "                pc.set_edgecolor(bar_color)\n",
    "                pc.set_alpha(bars_alpha)\n",
    "            parts['cmeans'].set(color='black')\n",
    "            #parts['cbars'].set(color='black')\n",
    "            #parts['cmins'].set(color='black')\n",
    "            #parts['cmaxes'].set(color='black')\n",
    "        # Add scatter points for each subject\n",
    "        scatter_x = ticks[len(groups)*roi_id+i]+(float(width)/2)\n",
    "        if plotkind != 'box':\n",
    "            # If barplot, plot all subjects points, else if boxplot then the boxplot will take care of plotting the outliers only\n",
    "            ax.scatter([scatter_x] * len(dg), dg, color=bar_color, marker='_', s=30)  # multiply scatter_x to repeat x position for as many values as we have in dg\n",
    "        # Add label for each subject scatter point\n",
    "        if show_subjects_labels:\n",
    "            for j, y in enumerate(dg):\n",
    "                text = dg.index.values[j] + 1\n",
    "                # Add label for all subjects if barplot, else only for outliers if boxplot\n",
    "                if plotkind != 'box' or y in bp['fliers'][0].get_ydata():\n",
    "                    t = ax.text(scatter_x, y, text, alpha=subjects_labels_alpha, fontsize=subjects_fontsize)\n",
    "                    texts.append(t)\n",
    "            last_j += j+1\n",
    "\n",
    "# Change the ticks to set the groups names (and place the labels nicely)\n",
    "ax.set_xticks([t + float(width)/2 for t in ticks])  # place in the middle of each bar (position tick t + half of bar width)\n",
    "groups_labels_reordered = [groups_labels[x] for x in groups_order]  # reorder the groups labels\n",
    "ax.set_xticklabels(groups_labels_reordered * nb_rois)\n",
    "# Add extra x-axis ROI labels: ROIs centers, colors and names if available (ie, if maps are provided)\n",
    "if rex_maps_filepath:\n",
    "    for ri in range(nb_rois):\n",
    "        tick_roi_offset = ri*len(groups)  # plus offset by the number of past ROIs\n",
    "        middletick = math.floor(len(groups)*0.5)  # middle tick of first ROI = half of nb of groups\n",
    "        tick_even_offset = 0.5 if len(groups) % 2 == 1 else 0  # offset to be in the middle of two ticks if the number of ticks/groups is even\n",
    "        axis_scaling = (1.0/len(ticks))  # scale to range [0.0, 1.0], because we use transform ax.transAxes (so we are in the space of the axis, not ticks)\n",
    "        t = ax.text((tick_roi_offset+middletick+tick_even_offset)*axis_scaling, roi_labels_yoffset, '%i,%i,%i\\n%s' % (round(centers[ri][0]), round(centers[ri][1]), round(centers[ri][2]), '\\n'.join(textwrap.wrap(', '.join(maps_regions[ri]), width=roi_labels_maxlen*len(groups), break_long_words=True)) if show_atlas_labels else ''),\n",
    "                rotation=0, verticalalignment='top', horizontalalignment='center',\n",
    "                transform=ax.transAxes,\n",
    "                color=atlas_cmap(norm(ri+1)),\n",
    "                fontsize=roi_labels_fontsize)\n",
    "        # pltcol.rgb2hex(atlas_cmap(norm(ri+1)))  # to get the hex value of the color\n",
    "# Force draw the plot (with tight layout)\n",
    "plt.tight_layout()  # MUST BE DONE BEFORE subplots_adjust!\n",
    "# Ensure proper and consistent scaling by setting the margins\n",
    "if figsize:\n",
    "    fig.subplots_adjust(left=margins[0]/figsize[0], right=1-margins[1]/figsize[0])  # Courtesy of ImportanceOfBeingErnest: https://stackoverflow.com/a/48532569/1121352\n",
    "if ylim:\n",
    "    ax.set_ylim(ylim)\n",
    "ax.set_xlim([ticks[0], ticks[-1]+width])\n",
    "plt.ylabel(ylabel)\n",
    "\n",
    "# Adjust label for each subject text placement to avoid overlapping (this needs to be done as the last step, after plotting the whole graph and adjusting the limits, else it can do funny things...)\n",
    "if show_subjects_labels and adjust_text is not None:\n",
    "    #texts = subjects_labels[last_j:end]\n",
    "    adjust_text(texts,\n",
    "                text_from_points=True,\n",
    "                text_from_text=True,\n",
    "        only_move={'text':'xy', 'objects':'x'}, force_text=0.05, force_objects=0.7, lim=10, precision=0.1) #, arrowprops=dict(arrowstyle=\"->\", color='r', lw=0.5))\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "# Save the figure\n",
    "barfilename = 'rois_%splot.png' % plotkind\n",
    "fig.savefig(barfilename, dpi=dpi_resolution, bbox_inches='tight')\n",
    "print('Image saved in %s' % barfilename)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------------------------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST: inconsistent scaling, bar width change with number of bars\n",
    "nb_bars_list = [2, 10]\n",
    "\n",
    "for i, nb_bars in enumerate(nb_bars_list):\n",
    "    figsize = [1+nb_bars, 5]\n",
    "    points = [np.random.randn(10) for x in range(nb_bars)]\n",
    "    fig, ax = plt.subplots()\n",
    "    if figsize:\n",
    "        fig.set_size_inches(figsize[0], figsize[1], forward=True)\n",
    "    ticks = np.arange(1, 1+nb_bars, 1)\n",
    "    for b in range(nb_bars):\n",
    "        ax.bar(ticks[b], points[b].mean(), 1)\n",
    "    plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)\n",
    "    fig.savefig('test%i' % i, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST: consistent scaling with same bar size whatever is the number of bars\n",
    "# Courtesy of ImportanceOfBeingErnest: https://stackoverflow.com/a/48532569/1121352\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "nb_bars_list = [2, 10]\n",
    "margleft = 0.8 # inch\n",
    "margright= 0.64 # inch\n",
    "barwidth = 1 # inch\n",
    "\n",
    "\n",
    "for i, nb_bars in enumerate(nb_bars_list):\n",
    "    # Resize proportionally to the number of bars\n",
    "    axwidth = nb_bars*barwidth # inch\n",
    "    figsize = [margleft+axwidth+margright, 5]\n",
    "    # Prepare the ticks\n",
    "    ticks = np.arange(1, 1+nb_bars, 1)\n",
    "    # Generate random points\n",
    "    points = [np.random.randn(10) for x in range(nb_bars)]\n",
    "    # Make the plot\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    fig.subplots_adjust(left=margleft/figsize[0], right=1-margright/figsize[0])\n",
    "    for b in range(nb_bars):\n",
    "        ax.bar(ticks[b], points[b].mean())\n",
    "    ax.set_xlim(ticks[0]-0.5,ticks[-1]+0.5)\n",
    "    fig.savefig('test%i' % i, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test boxplot plotting\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.figure()\n",
    "a = [1, 2, 3, 3, 3, 4, 5, 7]\n",
    "b = plt.boxplot(a, positions=[2], patch_artist=True, showmeans=True, sym='_', flierprops={'markeredgecolor': 'r', 'markeredgewidth': 2, 'markersize': 10})\n",
    "c = b['boxes'][0]\n",
    "c.set_facecolor('red')\n",
    "c.set_alpha(0.5)\n",
    "print(b['means'][0])\n",
    "d = b['means'][0]\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_examples(cms):\n",
    "    \"\"\"\n",
    "    helper function to plot two colormaps\n",
    "    \"\"\"\n",
    "    np.random.seed(19680801)\n",
    "    data = np.random.randn(30, 30)\n",
    "\n",
    "    fig, axs = plt.subplots(1, 2, figsize=(6, 3), constrained_layout=True)\n",
    "    for [ax, cmap] in zip(axs, cms):\n",
    "        psm = ax.pcolormesh(data, cmap=cmap, rasterized=True, vmin=-4, vmax=4)\n",
    "        fig.colorbar(psm, ax=ax)\n",
    "    plt.show()\n",
    "\n",
    "plot_examples([plt.cm.Set1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}