center-tbi-statistics-from-full-rest-dump.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Center-TBI statistics (from full dump via XNAT REST API/pyxnat)\n",
    "Compute statistics from the full database dumped in JSON from XNAT REST API (via pyxnat). Please generate beforehand `xnat_data_extrat.json` by executing `xnat_data_extractor.ipynb` (or get the file somewhere else, since it contains private project data it cannot be stored publicly here).\n",
    "\n",
    "Please edit login.cfg with your credentials before executing this script.\n",
    "\n",
    "Before (re-)running this script, please clear output, shutdown and relaunch kernel, close down and reopen your browser, and then (re-)launch all the cells! Else the memory is not correctly freed (this is a bug in ipywidgets or jupyter notebook).\n"
   ]
  },
  {
   "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",
    "\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",
    "    from fdict import fdict\n",
    "except ImportError as exc:\n",
    "    pass\n",
    "\n",
    "import json\n",
    "import pandas as pd\n",
    "\n",
    "# Setup some display options for pandas\n",
    "pd.set_option('max_columns', 400)\n",
    "pd.set_option('expand_frame_repr', False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Edit the filepath to the json file here\n",
    "json_filepath = 'xnat_data_extract.json'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Auxiliary Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def json_load_dict(json_filepath):\n",
    "    with open(json_filepath, 'rb') as f:\n",
    "        #j = f.read()\n",
    "        jdict = json.load(f)\n",
    "        #del j\n",
    "    return jdict\n",
    "\n",
    "def save_df_as_csv(d, output_file, fields_order=None, csv_order_by=None, encoding='ascii', escapechar='\\\\', index=False, verbose=False, **kwargs):\n",
    "    \"\"\"Save a dataframe in a csv\"\"\"\n",
    "    # Define CSV fields order\n",
    "    # If we were provided a fields_order list, we will show them first, else we create an empty fields_order\n",
    "    if fields_order is None:\n",
    "        fields_order = []\n",
    "    # Then automatically add any other field (which order we don't care, they will be appended in alphabetical order)\n",
    "    fields_order_check = set(fields_order)\n",
    "    for missing_field in sorted(d.columns):\n",
    "        if missing_field not in fields_order_check:\n",
    "            fields_order.append(missing_field)\n",
    "    if verbose:\n",
    "        print('CSV fields order: '+str(fields_order))\n",
    "\n",
    "    # Write the csv\n",
    "    # Note that escapechar is particularly important if you have nested fields (lists, json, etc)\n",
    "    if csv_order_by:\n",
    "        d2 = d.sort_values(csv_order_by)\n",
    "    else:\n",
    "        d2 = d\n",
    "    d2.to_csv(output_file, sep=';', index=index, columns=fields_order, encoding=encoding, escapechar=escapechar, **kwargs)\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting pertinent database fields into a Dataframe\n",
    "The original structure is very nested, here we extract only the relevant fields to get a 2D matrix of only the relevant fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the json dict\n",
    "jdict = json_load_dict(json_filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "\n",
    "# Reorder the json dict so that it is organized per scan session instead of per project\n",
    "def ensure_list(obj):\n",
    "    '''Ensure an object is a collection/list, and not a single object (this simplifies looping)'''\n",
    "    if '@ID' in obj:\n",
    "        # If it has an id attribute, this is a single object, we need to transform it into a list\n",
    "        obj = {obj['@ID']: obj}  # we use dicts here instead of lists (for various technical reasons)\n",
    "    return obj\n",
    "\n",
    "def copy_dict_exclude(d, exclude):\n",
    "    '''Copy a dict but exclude the specified keys'''\n",
    "    return {k:v for k,v in d.items() if k not in exclude}\n",
    "\n",
    "exclude_resources = False  # exclude resources because they use up a lot of columns and are pretty much useless for stats (only useful for download)\n",
    "\n",
    "# Initialization, do not touch\n",
    "scans_excludes = ['@xsi:schemaLocation']\n",
    "if exclude_resources:\n",
    "    scans_excludes.append('resources')\n",
    "\n",
    "# Create a dict per scan\n",
    "jdict2 = {'scans': []}\n",
    "for projectname, project in ensure_list(jdict['projects']).items():\n",
    "    # For each project\n",
    "    for subjectname, subject in ensure_list(project['subjects']).items():\n",
    "        # For each subject\n",
    "        for experimentname, experiment in ensure_list(subject['experiments']).items():\n",
    "            # For each experiment\n",
    "            for scanname, scan in ensure_list(experiment['scans']).items():\n",
    "                # For each scan\n",
    "                # Get the scan infos (excluding the columns/fields we do not want)\n",
    "                scan_temp = copy_dict_exclude(scan, scans_excludes)\n",
    "                if 'resources' in scan_temp:\n",
    "                    # Convert to a list so the list of files does not get flattened as columns\n",
    "                    if 'SNAPSHOTS' in scan_temp['resources'] and 'files' in scan_temp['resources']['SNAPSHOTS']:\n",
    "                        scan_temp['resources']['SNAPSHOTS']['files'] = [deepcopy(scan_temp['resources']['SNAPSHOTS']['files'])]\n",
    "                    if 'files' in scan_temp['resources']:\n",
    "                        scan_temp['resources']['files'] = [deepcopy(scan_temp['resources']['files'])]\n",
    "                # Add the scan infos\n",
    "                jdict2['scans'].append(copy_dict_exclude(scan, scans_excludes))  # add this scan record\n",
    "                currecord = jdict2['scans'][-1]  # keep it at hand so we can add infos of above levels\n",
    "                # Add the experiment infos for this scan (these infos are shared across all scans of this experiment, so we duplicate for each scan)\n",
    "                currecord['experiment'] = copy_dict_exclude(experiment, ['scans', 'xnat:scans', '@xsi:schemaLocation'])\n",
    "                # Same for subject infos\n",
    "                currecord['subject'] = copy_dict_exclude(subject, ['experiments', 'xnat:experiments', '@xsi:schemaLocation'])\n",
    "                # Same for project infos\n",
    "                currecord['project'] = copy_dict_exclude(project, ['subjects', 'xnat:subjects', '@xsi:schemaLocation'])\n",
    "                # Compute agregative scores for QA fields\n",
    "                try:\n",
    "                    a = json.loads(currecord['experiment']['xnat:assessors']['xnat:assessor']['ext:results_dict'], strict=False)  # need strict=False because of non-escaped \\n and \\r (line return + carriage return)\n",
    "\n",
    "                    # Protocol check QA: checks if the machine initialized correctly the parameters of a modality before acquisition\n",
    "                    currecord['qa.protocol_check'] = a['protocol_check']\n",
    "                    # clean line returns and spaces from keys of protocol check\n",
    "                    currecord['qa.protocol_check'] = [[x[0].replace('\\n', '').strip(), x[1]] for x in currecord['qa.protocol_check']]\n",
    "                    # extract per modality\n",
    "                    currecord['qa.protocol_check_t1'] = dict(currecord['qa.protocol_check']).get('T1', 'NA')\n",
    "                    currecord['qa.protocol_check_fmri'] = dict(currecord['qa.protocol_check']).get('rs_fMRI', 'NA')\n",
    "                    # compute aggregated score (all modalities passed or not?)\n",
    "                    currecord['qa.protocol_check_global'] = all(1 if v.lower() == 'pass' else 0 for _, v in a['protocol_check'])\n",
    "\n",
    "                    # Head coverage QA: check if head was not cut in the acquisition, for each scan and convert to list of lists for easy conversion to dataframe after (to stop fdict from automatically nesting elements)\n",
    "                    # The possible values are: GOOD, BAD or CHECK, where CHECK means that it should be manually checked (can be good or bad)\n",
    "                    currecord['qa.head_coverage'] = [[k,v] for k,v in a['head_coverage'].items()]\n",
    "                    # extract head coverage for only this scan\n",
    "                    currecord['qa.head_coverage_scan'] = dict(currecord['qa.head_coverage'])[currecord['id']] if currecord['id'] in dict(currecord['qa.head_coverage']) else None\n",
    "                    # compute aggregated score (all scans passed or not?)\n",
    "                    currecord['qa.head_coverage_global'] = all(1 if v.lower() == 'good' else 0 for v in a['head_coverage'].values())\n",
    "                except Exception:\n",
    "                    pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del jdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Postprocessing: we convert nested dicts into json strings\n",
    "#for scanid in xrange(len(jdict2['scans'])):\n",
    "#    for key, value in jdict2['scans'][scanid].items():\n",
    "#        if isinstance(value, dict):\n",
    "#            for subkey, subvalue in jdict2['scans'][scanid][key].items():\n",
    "#                if isinstance(subvalue, dict):\n",
    "#                    jdict2['scans'][scanid][key][subkey] = json.dumps(subvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You should check these files to see if the structure is correct and has all required fields\n",
    "\n",
    "# Save one record for test\n",
    "with open('center-tbi-statistics-from-full-rest-dump_onerecordtest.json', 'w') as f:\n",
    "    json.dump(currecord, f, indent=4, sort_keys=True)\n",
    "# Save all reordered records\n",
    "with open('center-tbi-statistics-from-full-rest-dump_persession.json', 'w') as f:\n",
    "    json.dump(jdict2, f, indent=4, sort_keys=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flatten keys so that pandas can access all nested dicts\n",
    "jdict3 = [fdict.flatkeys(scan, sep='.') for scan in jdict2['scans']]\n",
    "# Save flattened records\n",
    "with open('center-tbi-statistics-from-full-rest-dump_persession_flattened.json', 'w') as f:\n",
    "    json.dump(jdict3, f, indent=4, sort_keys=True)\n",
    "# Free up memory\n",
    "del jdict2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading final table from json file and showing stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    del jdict3\n",
    "except Exception:\n",
    "    pass\n",
    "with open('center-tbi-statistics-from-full-rest-dump_persession_flattened.json', 'r') as f:\n",
    "    jdict3 = json.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to a pandas dataframe!\n",
    "# TODO: if out of memory error, can try to use json-streamer with a mockup-class to return a dict-like generator to pandas: https://github.com/kashifrazzaqui/json-streamer\n",
    "df = pd.io.json.json_normalize(jdict3)\n",
    "# Free up memory\n",
    "del jdict3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Show structure of df\n",
    "print('The dataframe has %s (lines, columns)' % str(df.shape))\n",
    "print('Columns of the dataframe: ')\n",
    "print(df.columns)\n",
    "df  # to pretty print, just put the dataframe name as the last line, without print or anything"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Count of unique values per columns:')\n",
    "for c in df.columns:\n",
    "    try:\n",
    "        print('* %s: %s' % (c, str(len(df[c].unique()))))\n",
    "    except TypeError as exc:\n",
    "        print('* %s: type error (list), cannot compute length' % c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Count per value for all columns:')\n",
    "#with pd.option_context(\"display.max_rows\", -1, \"display.max_columns\", -1):\n",
    "for c in df.columns:\n",
    "    try:\n",
    "        print \"---- %s ---\" % c\n",
    "        print df[c].value_counts()\n",
    "    except TypeError as exc:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = df.apply(lambda df: (df['project.id'], df['qa.protocol_check_global']), axis=1)\n",
    "print('Number of validated scans per center')\n",
    "a.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Only centers with validated scans:')\n",
    "a.apply(lambda x: x if x[1] == True else None).dropna().value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Total number of validated scans:')\n",
    "len(a.apply(lambda x: x if x[1] == True else None).dropna())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Only centers with unvalidated scans:')\n",
    "a.apply(lambda x: x if x[1] != True else None).dropna().value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df[['project.id', 'qa.protocol_check_global']].groupby(['project.id']).value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = df.apply(lambda df: (df['project.id'], df['qa.head_coverage_global']), axis=1)\n",
    "print('Number of validated head coverage scans per center')\n",
    "a.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phantom_validated_centers = df.where(df['experiment.@visit_id'] == 'Phantom')['project.id'].dropna().unique()\n",
    "print('List of Phantom validated centers:')\n",
    "for center in sorted(phantom_validated_centers):\n",
    "    print('* %s' % center)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### rs_fMRI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Scan types:')\n",
    "for t in sorted(df['@type'].unique()):\n",
    "    print('* '+ str(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show count per scan type\n",
    "with pd.option_context('display.max_rows', None):  # show all rows (remove pandas default limit)\n",
    "    print(df['@type'].value_counts().sort_values(ascending=False))\n",
    "# Alternative way (and equivalent):\n",
    "#print('Scan types:')\n",
    "#for k, t in df['@type'].value_counts().sort_values(ascending=False).iteritems():\n",
    "#    print('* '+k+': '+str(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract only resting state data\n",
    "rsfMRI_filter = ['rs fMRI FE EPI', 'rs fMRI FE_EPI', 'rs fMRI FE_EPI SENSE', 'rs-fMRI', 'rs_fMRI', 'rsfMRI', 'rsfMRI 60CM', 'Resting State fMRI', ' rsfMRI', 'RESTING STATE']\n",
    "df_rest = df.loc[df['@type'].isin(rsfMRI_filter), :]\n",
    "#df_rest = df.loc[df['@type'] == 'rs_fMRI', :]  # deprecated, only for rs_fmri, does not include all variants of the name in db\n",
    "df_rest[['@UID', 'project.id', 'qa.protocol_check_global']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save full database as csv file\n",
    "save_df_as_csv(df, 'xnat_data_extract_fulldb.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save rs_fMRI only data to csv file\n",
    "save_df_as_csv(df_rest, 'xnat_data_extract_rs-fMRI-only.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_rest['experiment.@visit_id'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save rs_fMRI only data with only MR 2 weeks or Healthy Volunteer, which are the ones we can use, to csv file\n",
    "df_rest_csg = df_rest[df_rest['experiment.@visit_id'].isin(['MR 2 weeks', 'Healthy Volunteer'])]\n",
    "save_df_as_csv(df_rest_csg, 'xnat_data_extract_rs-fMRI-MR2weeks-only.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')\n",
    "df_rest_csg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate only validated scans csv list\n",
    "df_validated = df.ix[df['qa.protocol_check_global'] == True]\n",
    "df_rest_validated = df_rest.ix[df_rest['qa.protocol_check_fmri'] == 'PASS']\n",
    "# Save them\n",
    "save_df_as_csv(df_validated, 'xnat_data_extract_validatedonly.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')\n",
    "save_df_as_csv(df_rest_validated, 'xnat_data_extract_rs-fMRI-only-validatedonly.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract all subjects with at least one fmri session AND are ICU (intensive care unit)\n",
    "# IMPORTANT: this is the list of scans we should use from now on for WP15.\n",
    "df_rest_icu = df[(df['subject.@group'] == 'ICU') & df['@type'].isin(rsfMRI_filter)]  # TODO: last step is to filter out GOSE 0-3 to get only real severe TBI (our target population for WP15)\n",
    "save_df_as_csv(df_rest_icu, 'xnat_data_extract_rest-icu.csv', csv_order_by=['project.id', '@UID'], encoding='utf-8')\n",
    "print('ICU subjects with fMRI: total subjects count %i, total sessions count %i.' % (len(df_rest_icu.groupby('subject.@ID')), len(df_rest_icu)))\n",
    "df_rest_icu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stats figures\n",
    "Generate interesting stats figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "matplotlib.style.use('ggplot')  # for nicer plots\n",
    "#plt.xkcd()  # uncomment this line for xkcd style plots!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def df_couples(df, cols):\n",
    "    return df.apply(lambda row: tuple(row[col] for col in cols), axis=1).value_counts()\n",
    "\n",
    "def plot_boolean(df, col1, col2, morecols=None, where=None, groupby=None, truename='validated', falsename='not validated', title=None, nostats=False, nograph=False):\n",
    "    '''Plot the col2 according to col1, counting for the number of unique values (true/false or true/nan)'''\n",
    "    # Prepare stuff\n",
    "    if not morecols:\n",
    "        morecols = []\n",
    "    if not groupby:\n",
    "        groupby = []\n",
    "    elif isinstance(groupby, str):\n",
    "        groupby = [groupby]\n",
    "    # Extract both columns\n",
    "    if where is not None:\n",
    "        a = df.loc[where, [col1, col2]+morecols+groupby]\n",
    "    else:\n",
    "        a = df[[col1, col2]+morecols+groupby]\n",
    "    # Group by the first column (we call this column the categories)\n",
    "    groupby.insert(0, col1)\n",
    "    b = a.groupby(groupby)\n",
    "    # Aggregate using the sum (or count of True values) and size (the total number of entries)\n",
    "    c = b.aggregate(['sum', 'size'])\n",
    "    # Replace nan/none values by 0 (when there is absolutely no true value for one category)\n",
    "    c.fillna(0, inplace=True)\n",
    "    # Because of aggregation we get a hierarchical multiindex, we flatten the indices for easier access\n",
    "    c.columns = c.columns.get_level_values(1)\n",
    "    # Compute the number of false values (true values count - total values count)\n",
    "    c['size'] = c['size'] - c['sum']\n",
    "    # Rename the columns\n",
    "    c = c.rename(columns={'sum': truename, 'size': falsename})\n",
    "    # Group by first level if necessary\n",
    "    if groupby is not None:\n",
    "        d = c.groupby(level=0).sum()\n",
    "    else:\n",
    "        d = c\n",
    "    if not nograph:\n",
    "        # Plot as stacked bars\n",
    "        d.plot.bar(stacked=True, title=title)\n",
    "        plt.show()\n",
    "        # Plot only validated\n",
    "        d[truename].plot.bar(stacked=True, title=(title+' (validated only)'))\n",
    "        plt.show()\n",
    "    if not nostats:\n",
    "        # Print the table\n",
    "        print(c)\n",
    "        # Print stats\n",
    "        print(c.describe())\n",
    "        print('\\nTotal sum:\\n'+str(c.sum()))\n",
    "        print('\\nNumber of nonzeros %s entries:\\n%s' % (str(groupby), str((c != 0).sum())) )\n",
    "    # Return c and d in case user wants to manipulate them further (or debug)\n",
    "    return c, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_categorical(df, col1, col2, morecols=None, where=None, groupby=None, truename='validated', falsename='not validated', title=None, sort=None, nostats=False, nograph=False):\n",
    "    '''Plot the col2 according to col1, counting for the number of unique values (categorical values)'''\n",
    "    # Prepare stuff\n",
    "    if not morecols:\n",
    "        morecols = []\n",
    "    if not groupby:\n",
    "        groupby = []\n",
    "    elif isinstance(groupby, str):\n",
    "        groupby = [groupby]\n",
    "    # Extract both columns\n",
    "    if where is not None:\n",
    "        a = df.loc[where, [col1, col2]+morecols+groupby]\n",
    "    else:\n",
    "        a = df[[col1, col2]+morecols+groupby]\n",
    "    # Replace None/NaN values by a 'NA' string so we can also account for none values\n",
    "    a[col2] = a[col2].astype('str')  # alternative would be to convert fillna('NA'), but then it would replace all nan/none, even though they are different values\n",
    "    # Group by all columns we want (including the categories col2) and then unstack to compute the count for each category\n",
    "    b = a.groupby([col1]+groupby+[col2]).size().unstack(fill_value=0)\n",
    "    # Reorder the columns (will make it easier to read)\n",
    "    if sort is not None:\n",
    "        col_order = list(b.columns)\n",
    "        col_order.sort(sort)\n",
    "        b = b[col_order]\n",
    "    # Group by first level if necessary\n",
    "    if groupby:\n",
    "        # If groupby, count the number of non zeros entries\n",
    "        d = (b != 0).groupby(level=0).sum()\n",
    "    else:\n",
    "        # If no groupby, calculate the sum\n",
    "        d = b.groupby(level=0).sum()\n",
    "    if not nograph:\n",
    "        # Plot as stacked bars\n",
    "        ax = d.plot.bar(stacked=True, title=title, figsize=(12,12))\n",
    "        plt.legend(prop={'size':15}, loc='best')\n",
    "        plt.xticks(size=15)\n",
    "        plt.yticks(size=15)\n",
    "        ax.title.set_size(20)\n",
    "        plt.show()\n",
    "        # Plot only validated\n",
    "        ax = d.plot.bar(subplots=True, title=(title+' (per category)'), figsize=(15,20))\n",
    "        plt.tight_layout()  # tighten the space between plots\n",
    "        plt.xticks(size=15)\n",
    "        for axe in ax:  # resize all titles and legends and labels\n",
    "            axe.title.set_size(15)\n",
    "            axe.legend(prop={'size':15}, loc='best')\n",
    "        plt.suptitle((title+' (per category)'), size=20)  # resize main title\n",
    "        plt.show()\n",
    "    if not nostats:\n",
    "        # Print the table\n",
    "        print(b)\n",
    "        # Print stats\n",
    "        print(b.describe())\n",
    "        print('\\nTotal sum:\\n'+str(b.sum()))\n",
    "        print('\\nNumber of nonzeros %s entries:\\n%s' % (str(groupby), str((b != 0).sum())) )\n",
    "    # Return b and d in case user wants to manipulate them further (or debug)\n",
    "    return b, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_values_repartition(df, col, groupby=None, colrename=None, title='Repartition', limit=30):\n",
    "    '''Plot the values repartition of col'''\n",
    "    # Groupby col1 and extract the first value of col2 (we assume it is the same for each value of col1 eg for each subject)\n",
    "    if groupby is not None:\n",
    "        df_repartition = df.ix[:, (groupby, col)].groupby(groupby).aggregate(lambda x: set(x))\n",
    "        df_repartition = df_repartition.applymap(lambda x: list(x)[0] if x else None)\n",
    "    # Extract only the column we are interested in\n",
    "    df_repartition = df.ix[:, col]\n",
    "    # Rename column to be nicer\n",
    "    if colrename:\n",
    "        df_repartition.columns = [colrename]\n",
    "        col = colrename\n",
    "    # Convert to a series, then convert to str to include nan, then count frequency of each value type, then plot\n",
    "    c = df_repartition.astype('str').value_counts()\n",
    "    # Limit the number of values we plot to the top k most frequent ones (if any limit is set)\n",
    "    if limit > 0 and len(c) > limit:\n",
    "        c = c[:limit]\n",
    "        c.plot(kind='bar', title=title+(' (top %i)' % limit))\n",
    "    else:\n",
    "        c.plot(kind='bar', title=title)\n",
    "    plt.show()\n",
    "    # Display gender table\n",
    "    print(df_repartition)\n",
    "    return df_repartition\n",
    "\n",
    "def plot_values_repartition_multi(df, col, groupby=None, title='Repartition', limit=30):\n",
    "    '''Just like plot_values_repartition() but also works when the groupby will generate multiple different values\n",
    "    (so here instead of using a set and getting the first value because we assume that it is the same for all, we actually compute the existence of every value type for this group/subject)'''\n",
    "    if groupby is not None:\n",
    "        # This is the magic: groupby and then do value_counts() on col (the target column), the value_counts() will account for the groupby\n",
    "        df_repartition = df.ix[:, (groupby, col)].groupby(groupby)[col].value_counts()\n",
    "    else:\n",
    "        # If we do not use groupby, then it reduces down to normal plot_values_repartition()\n",
    "        return plot_values_repartition(df, col, title=title, limit=limit)\n",
    "    # Since we want to know the existence per group rather than the total count per group, we have to replace all values to 1\n",
    "    df_repartition[:] = 1\n",
    "    # Then do the sum per value category of col (level 1 because we are in a Series now)\n",
    "    df_repartition = df_repartition.groupby(level=1).sum()\n",
    "    # Sort the values\n",
    "    df_repartition.sort_values(ascending=False, inplace=True)\n",
    "    # And plot!\n",
    "    if limit > 0 and len(df_repartition) > limit:\n",
    "        # If too many values, we reduce\n",
    "        df_repartition[:limit].plot(kind='bar', title=title+(' (top %i)'%limit))\n",
    "    else:\n",
    "        df_repartition.plot(kind='bar', title=title)\n",
    "    # Force drawing the plot\n",
    "    plt.show()\n",
    "    # Print the repartition table\n",
    "    print(df_repartition)\n",
    "    # Return the repartition table for further processing by user\n",
    "    return df_repartition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_evolution_over_time(df, datecol, col, groupby=None, dateformat='%Y-%m-%d', title='Evolution over time'):\n",
    "    '''Plot the evolution of a BOOLEAN col column over time (as defined by the datecol column)'''\n",
    "    if not title:\n",
    "        title = 'Evolution of %s over time (%s)' % (col, datecol)\n",
    "    # Extract only the columns of interest\n",
    "    if groupby is None:\n",
    "        dtime = df.ix[:, (datecol, col)]\n",
    "    else:\n",
    "        dtime = df.ix[:, (datecol, groupby, col)]\n",
    "    # Convert the datecol to date format\n",
    "    dtime[datecol] = pd.to_datetime(dtime[datecol], format=dateformat)\n",
    "    # Convert to boolean and create the \"false\" column so we can plot both true and false\n",
    "    dtime[col] = (dtime[col] == True)\n",
    "    dtime[col+'_false'] = (~dtime[col])\n",
    "    # Set the date as the index of the df\n",
    "    if groupby is None:\n",
    "        dtime2 = dtime.set_index(datecol)\n",
    "    else:\n",
    "        dtime2 = dtime.set_index([datecol, groupby])\n",
    "    # Sort by date\n",
    "    dtime2.sort_index(inplace=True)\n",
    "    # Compute cumulative sum over time and plot\n",
    "    if groupby is None:  # Plot for all without any grouping\n",
    "        # Compute the cumulative sum over time\n",
    "        dtimecum = dtime2.cumsum()\n",
    "        # Plot as an area graph!\n",
    "        dtimecum.plot(kind='area', title=title)\n",
    "        # Plot also the only validated graph\n",
    "        plt.figure()\n",
    "        dtimecum[col].plot(kind='area', title=title+' (only validated)')\n",
    "    else:  # Groupby first and then plot per group the cumulative sum\n",
    "        dtimecum = dtime2.groupby(level=[1]).cumsum()\n",
    "        dtimecum = dtimecum.reset_index()\n",
    "        dtimecum = dtimecum.set_index(datecol)\n",
    "        #a.groupby('project.id').plot(kind='area')\n",
    "        for subtitle, group in dtimecum.groupby(groupby):\n",
    "            group.plot(kind='area', title=subtitle)\n",
    "            group.plot(kind='area', title=subtitle, subplots=True)\n",
    "    # Force draw the plots\n",
    "    plt.show()\n",
    "    # Return variables for user to manipulate\n",
    "    return dtime2, dtimecum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_boolean(df, 'project.id', 'qa.protocol_check_global', title='Protocol check global (all scan types)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_boolean(df, 'project.id', 'qa.head_coverage_global', truename='Good', falsename='Bad', title='Head coverage global (all scan types)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# One center is validated but has bad head coverage, show which one\n",
    "df_bad_head_coverage = df.ix[(df['qa.head_coverage_global'] != True) & (df['qa.protocol_check_global'] == True)]\n",
    "df_good_head_coverage = df.ix[(df['qa.head_coverage_global'] == True) & (df['qa.protocol_check_global'] == True)]\n",
    "print('One example:')\n",
    "print(df_bad_head_coverage.ix[:, 'qa.head_coverage'].iloc[0])\n",
    "print('List of validated centers but which have some scans with bad head coverage:')\n",
    "print(df_bad_head_coverage['project.id'].unique())\n",
    "print(\"Total number of scans that would be validated but have bad head coverage: %i\" % len(df_bad_head_coverage))\n",
    "df_centers_no_good_head_coverage = set(df_bad_head_coverage['project.id']) - set(df_good_head_coverage['project.id'].unique())\n",
    "print('Validated centers with no good head coverage:')\n",
    "print(df_centers_no_good_head_coverage)\n",
    "print('List of validated scans with bad head coverage:')\n",
    "df_bad_head_coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_boolean(df, 'project.id', 'qa.protocol_check_global', where=(df['@type'].isin(rsfMRI_filter)), title='Protocol check global (rs_fMRI)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the rs_fMRI session that is validated but the type is not the standard \"rs_fMRI\"\n",
    "df.ix[(df['qa.protocol_check_global'] == True) & (df['@type'].isin(rsfMRI_filter)) & (df['@type'] != 'rs_fMRI')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Number of subjects with validated rs_fMRI')\n",
    "plot_boolean(df, 'project.id', 'qa.protocol_check_global', where=(df['@type'].isin(rsfMRI_filter)), groupby=['subject.@ID'], title='Protocol check global (rs_fMRI)', nograph=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Number of subjects with validated rs_fMRI')\n",
    "plot_boolean(df, 'subject.@ID', 'qa.protocol_check_global', where=(df['@type'].isin(rsfMRI_filter)), title='Protocol check global (rs_fMRI)', nograph=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Number of subjects with validated rs_fMRI')\n",
    "plot_boolean(df, 'subject.@ID', 'qa.protocol_check_global', where=(df['@type'].isin(rsfMRI_filter)), title='Protocol check global (rs_fMRI)', nograph=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Number of subjects with validated T1')\n",
    "plot_boolean(df, 'project.id', 'qa.protocol_check_global', where=(df['@type'] == 'T1'), title='Protocol check global (T1)')\n",
    "plot_boolean(df, 'project.id', 'qa.protocol_check_global', where=(df['@type'] == 'T1'), groupby=['subject.@ID'], title='Protocol check global (T1)', nograph=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Find healthy controls\n",
    "df_healthy = df_rest[['project.id', 'subject.id', 'experiment.@visit_id']].ix[(df_rest['experiment.@visit_id'] == 'Healthy Volunteer')]\n",
    "df_healthy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_boolean(df_rest, 'project.id', 'qa.protocol_check_global', where=(df_rest['experiment.@visit_id'] == 'Healthy Volunteer'), groupby=['subject.@ID'], title='Protocol check global (controls with rs_fMRI)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot total number of controls per center (validated or not)\n",
    "df_healthy.groupby(['project.id'])['experiment.@visit_id'].count().plot(kind='bar', title='Total per center of healthy volunteers')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Show number of patients (excluding controls) with validated rs_fMRI\n",
    "plot_boolean(df_rest, 'project.id', 'qa.protocol_check_global', where=(df['experiment.@visit_id'] != 'Healthy Volunteer'), title='Protocol check global (rs_fMRI patients only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot number of validated healthy controls\n",
    "plot_boolean(df_rest, 'project.id', 'qa.protocol_check_global', where=(df['experiment.@visit_id'] == 'Healthy Volunteer'), title='Protocol check global (rs_fMRI controls only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check TR parameter and others among rsfMRI scans\n",
    "print(df_rest['xnat:parameters.xnat:tr'].unique())\n",
    "print(df_couples(df_rest, ['project.id', 'xnat:parameters.xnat:tr']))\n",
    "df_rest['xnat:parameters.xnat:voxelRes@xyz'] = df_rest['xnat:parameters.xnat:voxelRes.@x'].astype('str') + 'x' + df_rest['xnat:parameters.xnat:voxelRes.@y'] + 'x' + df_rest['xnat:parameters.xnat:voxelRes.@z']\n",
    "print(df_rest['xnat:parameters.xnat:voxelRes@xyz'].unique())\n",
    "plot_categorical(df_rest, 'project.id', 'xnat:parameters.xnat:tr', title='TR repartition (rs_fMRI only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Extract all scanner models of validated scans\n",
    "df_rest_validated['experiment.xnat:scannermodel'] = df_rest_validated['experiment.xnat:scanner.@manufacturer'].astype('str') + ' ' + df_rest_validated['experiment.xnat:scanner.@model'].astype('str')\n",
    "# fill the nan ones by the other scanner model field\n",
    "idxs = (df_rest_validated['experiment.xnat:scannermodel'] == 'nan nan')\n",
    "df_rest_validated.ix[idxs, 'experiment.xnat:scannermodel'] = df_rest_validated[idxs]['experiment.xnat:scanner']\n",
    "# Show unique scanner models list\n",
    "print(df_rest_validated['experiment.xnat:scannermodel'].unique())\n",
    "#df_couples(df_rest_validated, ['project.id', 'experiment.xnat:scannermodel'])\n",
    "plot_categorical(df_rest_validated, 'project.id', 'experiment.xnat:scannermodel', title='Scanner model per center (rs_fMRI validated only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df, 'project.id', 'subject.@group', title='Subjects groups (all scans)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df, 'project.id', 'subject.@group', groupby=['subject.id'], title='Subjects groups (all scans types)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df_rest, 'project.id', 'subject.@group', groupby=['subject.id'], title='Subjects groups (rs_fMRI only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df_rest_validated, 'project.id', 'subject.@group', groupby=['subject.id'], title='Subjects groups (validated+rsfMRI)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df_rest_icu, 'project.id', 'qa.head_coverage_scan', title='Head coverage (fmri+ICU)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df_rest_icu, 'project.id', 'qa.protocol_check_fmri', title='Protocol check (fmri+ICU)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def custom_sort(a, b):\n",
    "    '''Sort by MR early first and then by the shortest MR timeframe first'''\n",
    "    if a == b:\n",
    "        return 0\n",
    "    elif a.startswith('MR') and b.startswith('MR'):\n",
    "        if a.split()[1] == 'Early':\n",
    "            return -1\n",
    "        elif int(a.split()[1]) > int(b.split()[1]):\n",
    "            return 1\n",
    "        else:\n",
    "            return -1\n",
    "    else:\n",
    "        return 1\n",
    "\n",
    "plot_categorical(df_rest, 'project.id', 'experiment.@visit_id', groupby=['subject.id'], title='Longitudinal scans availability (rs_fMRI)', sort=custom_sort)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_categorical(df, 'project.id', 'experiment.@visit_id', groupby=['subject.id'], title='Longitudinal scans availability', sort=custom_sort)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# nb of validated rs_fmri subjects\n",
    "len(df_rest_validated.groupby('subject.id'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of all available longitudinal categories\n",
    "df_rest['experiment.@visit_id'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test\n",
    "b = df_rest[['subject.id', 'experiment.@visit_id']].groupby(['experiment.@visit_id'])\n",
    "b.aggregate('sum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Build longitudinal table for every subjects (without project id)\n",
    "def generate_longitudinal_table_noproject(df):\n",
    "    b = df_couples(df, ['subject.id', 'experiment.@visit_id'])\n",
    "    df_longitudinal = pd.DataFrame(columns=df['experiment.@visit_id'].unique())\n",
    "    for k, v in b.iteritems():\n",
    "        df_longitudinal.ix[k] = v\n",
    "    df_longitudinal.fillna('', inplace=True)\n",
    "    return df_longitudinal\n",
    "\n",
    "df_rest_icu_longitudinal = generate_longitudinal_table_noproject(df_rest_icu)\n",
    "save_df_as_csv(df_rest_icu_longitudinal, 'xnat_data_extract_rs-fMRI-icu-longitudinal-noproject.csv', index=True, encoding='utf-8')\n",
    "df_rest_icu_longitudinal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Generate longitudinal data table for ICU rs_fMRI sessions (with project id)\n",
    "def generate_longitudinal_table(df):\n",
    "    b = df_couples(df, ['project.id', 'subject.id', 'experiment.@visit_id'])\n",
    "    index = pd.MultiIndex.from_tuples(df_couples(df, ['project.id', 'subject.id']).index)\n",
    "    df_longitudinal = pd.DataFrame(columns=df['experiment.@visit_id'].unique(), index=index)\n",
    "    for k, v in b.iteritems():\n",
    "        df_longitudinal.ix[(k[0], k[1]), k[2]] = v\n",
    "    return df_longitudinal\n",
    "\n",
    "df_rest_icu_longitudinal = generate_longitudinal_table(df_rest_icu)\n",
    "save_df_as_csv(df_rest_icu_longitudinal.fillna(''), 'xnat_data_extract_rs-fMRI-icu-longitudinal.csv', index=True, encoding='utf-8')\n",
    "df_rest_icu_longitudinal.fillna('')  # fill nan with empty string for easier visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate table of only subjects with at least 2 sessions\n",
    "df_rest_longitudinal_morethan1 = df_rest_icu_longitudinal.dropna(axis=0, thresh=2)\n",
    "print((~df_rest_longitudinal_morethan1.isnull()).sum(axis=1))\n",
    "save_df_as_csv(df_rest_longitudinal_morethan1.fillna(''), 'xnat_data_extract_rs-fMRI-icu-longitudinal_morethan1.csv', index=True, encoding='utf-8')\n",
    "df_rest_longitudinal_morethan1.fillna('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Test: show all sessions for one specific subject\n",
    "df_rest.ix[df_rest['subject.id'] == 'CTBI_S01190', ('subject.id', 'experiment.@visit_id', 'experiment.xnat:date', 'experiment.xnat:time', 'project.@ID', 'qa.head_coverage', 'qa.head_coverage_global', 'qa.protocol_check', 'qa.protocol_check_global', 'xnat:quality')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show all available quality categories\n",
    "df['xnat:quality'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show all usable rs_fmri scans for every centers\n",
    "df_couples(df_rest, ['project.id', 'xnat:quality'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot number of usable rs_fMRI subjects per center\n",
    "df['qa.quality'] = (df['xnat:quality'] == 'usable')\n",
    "df_rest['qa.quality'] = (df_rest['xnat:quality'] == 'usable')\n",
    "c, d = plot_boolean(df_rest, 'project.id', 'qa.quality', groupby=['subject.@ID'], title='Quality (rs_fMRI)', truename='usable', falsename='unusable')\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evolution of the number of rsfMRI sessions acquisition over time\n",
    "plot_evolution_over_time(df_rest, 'experiment.xnat:date', 'qa.protocol_check_global', title='Cumulative number of rs_fMRI sessions over time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evolution of the number of acquired rsfMRI sessions over time - per project\n",
    "plot_evolution_over_time(df_rest, 'experiment.xnat:date', 'qa.protocol_check_global', groupby='project.id', title='Cumulative number of rs_fMRI sessions over time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot evolution of number of acquisitions over time - for all types (including CT)\n",
    "plot_evolution_over_time(df, 'experiment.xnat:date', 'qa.protocol_check_global', title='Cumulative number of sessions over time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot evolution of number of acquisitions over time - for all types (including CT) - per center\n",
    "plot_evolution_over_time(df, 'experiment.xnat:date', 'qa.protocol_check_global', groupby='project.id', title='Cumulative number of sessions over time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Most missing fields\n",
    "df_nan = df.isnull().sum() / len(df)\n",
    "for key, val in df_nan.sort_values(ascending=False).iteritems():\n",
    "    print('%s: %.4f' % (key, val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Most complete fields\n",
    "df_nan = df.count() / len(df)\n",
    "for key, val in df_nan.sort_values(ascending=False).iteritems():\n",
    "    print('%s: %.4f' % (key, val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot distribution/histogram of each column\n",
    "# TODO: turn this into a function, can be useful for other projects\n",
    "def ctbi_add_age(df):\n",
    "    df['subject.age_days'] = pd.to_datetime(df['experiment.xnat:date'], format='%Y-%m-%d') - pd.to_datetime(df['subject.xnat:demographics.xnat:yob'], format='%Y')\n",
    "    df['subject.age'] = df['subject.age_days'].apply(lambda d: d.days / 365 if pd.notnull(d) else d)\n",
    "    return df\n",
    "\n",
    "def plot_repartition_all(df, limit=30, maxlimit=100):\n",
    "    '''Plot the repartition of EVERY columns'''\n",
    "    # Temporarily turn back to matplotlib default style, for nicer plots\n",
    "    #with pd.option_context('mpl_style', 'default'):\n",
    "    for col in df.columns:\n",
    "        try:\n",
    "            if df[col].count() > 0:\n",
    "                if col.startswith('resources'):\n",
    "                    # Skip resources columns (useless info anyway, we don't load them up in df anyway, they are too much bloated)\n",
    "                    continue\n",
    "                print(col)\n",
    "                try:\n",
    "                    plt.figure()\n",
    "                    df[col].fillna(0).plot(kind='hist', title=col)\n",
    "                except TypeError:\n",
    "                    # For categorical type values, we cannot use hist (which expects only continuous values), thus fallback to value_counts().barplot()\n",
    "                    if len(df[col].unique()) < maxlimit:\n",
    "                        df[col].fillna(0).value_counts().plot(kind='bar', title=col)\n",
    "                    else:\n",
    "                        # If too many unique values, plot only the top 30 most frequent ones\n",
    "                        df[col].fillna(0).value_counts()[:limit].plot(kind='bar', title=col+(' (top %i)'%limit))\n",
    "        except TypeError:\n",
    "            pass\n",
    "    plt.show()\n",
    "\n",
    "plot_repartition_all(df = ctbi_add_age(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Find buggy age patients (with age > 100)\n",
    "# TODO: enhance by grouping by subject.id (and remove columns that change like @ID and @type, all the others like yob should not change per subject) - find how to aggregate/groupby when the value is the same in pandas?\n",
    "df_buggy_age = df.ix[df['subject.age'] > 100,:].ix[:, ('@ID', '@type', 'project.id', 'experiment.@ID', 'subject.id', 'subject.age', 'experiment.xnat:date', 'subject.xnat:demographics.xnat:yob')]\n",
    "save_df_as_csv(df_buggy_age.fillna(''), 'xnat_data_extract_buggy-age.csv', index=False, encoding='utf-8')\n",
    "print('Number of patients with wrong yob: %i' % len(df_buggy_age.groupby('subject.id')))\n",
    "df_buggy_age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_categorical(df, 'project.id', 'datatype', title='Acquisitions datatypes (all scans)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Display gender repartition per subjects\n",
    "plot_values_repartition(df, 'subject.xnat:demographics.xnat:gender', groupby='subject.id', colrename='gender', title='Gender repartition (per subject)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display gender repartition per subjects (for rs_fMRI only)\n",
    "plot_values_repartition(df_rest, 'subject.xnat:demographics.xnat:gender', groupby='subject.id', colrename='gender', title='Gender repartition (per subject & rs_fMRI only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display gender repartition per subjects (for rs_fMRI and validated only)\n",
    "plot_values_repartition(df_rest_validated, 'subject.xnat:demographics.xnat:gender', groupby='subject.id', colrename='gender', title='Gender repartition (per subject & rs_fMRI+validated only)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_values_repartition(df, '@type', groupby='subject.id', title='Repartition of type (all scans)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the @type repartition grouped by subject (because DTI needs multiple acquisitions, so it is inflating the real number if we just look at the number of scans, we need to look per subject to compare to T1 and rs_fMRI acquisition numbers)\n",
    "plot_values_repartition_multi(df, '@type', 'subject.id', title='Repartition of @type (per subject)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Build demographics summary table (number of scans, number of subjects, validated or not etc)\n",
    "\n",
    "# Scans\n",
    "nb_scans = [\n",
    "    len(df),\n",
    "    len(df.ix[df['@type'] == 'DTI', :]),\n",
    "    len(df.ix[df['@type'] == 'T1', :]),\n",
    "    len(df.ix[df['@type'].isin(rsfMRI_filter), :])\n",
    "    ]\n",
    "nb_scans_validated = [\n",
    "    len(df_validated),\n",
    "    len(df_validated.ix[df['@type'] == 'DTI', :]),\n",
    "    len(df_validated.ix[df['@type'] == 'T1', :]),\n",
    "    len(df_validated.ix[df['@type'].isin(rsfMRI_filter), :])\n",
    "    ]\n",
    "nb_scans_goodhead = [\n",
    "    len(df.ix[df['qa.head_coverage_global'] == True]),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & (df['@type'] == 'DTI'), :]),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & (df['@type'] == 'T1'), :]),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & df['@type'].isin(rsfMRI_filter), :])\n",
    "    ]\n",
    "nb_scans_useable = [\n",
    "    len(df.ix[df['xnat:quality'] == 'usable']),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & (df['@type'] == 'DTI'), :]),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & (df['@type'] == 'T1'), :]),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & df['@type'].isin(rsfMRI_filter), :])\n",
    "    ]\n",
    "\n",
    "# Subjects\n",
    "nb_subjects = [\n",
    "    len(df.groupby('subject.id')),\n",
    "    len(df.ix[df['@type'] == 'DTI', :].groupby('subject.id')),\n",
    "    len(df.ix[df['@type'] == 'T1', :].groupby('subject.id')),\n",
    "    len(df.ix[df['@type'].isin(rsfMRI_filter), :].groupby('subject.id'))\n",
    "    ]\n",
    "nb_subjects_validated = [\n",
    "    len(df_validated.groupby('subject.id')),\n",
    "    len(df_validated.ix[df['@type'] == 'DTI', :].groupby('subject.id')),\n",
    "    len(df_validated.ix[df['@type'] == 'T1', :].groupby('subject.id')),\n",
    "    len(df_validated.ix[df['@type'].isin(rsfMRI_filter), :].groupby('subject.id'))\n",
    "    ]\n",
    "nb_subjects_goodhead = [\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True)].groupby('subject.id')),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & (df['@type'] == 'DTI'), :].groupby('subject.id')),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & (df['@type'] == 'T1'), :].groupby('subject.id')),\n",
    "    len(df.ix[(df['qa.head_coverage_global'] == True) & df['@type'].isin(rsfMRI_filter), :].groupby('subject.id'))\n",
    "    ]\n",
    "nb_subjects_useable = [\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable')].groupby('subject.id')),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & (df['@type'] == 'DTI'), :].groupby('subject.id')),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & (df['@type'] == 'T1'), :].groupby('subject.id')),\n",
    "    len(df.ix[(df['xnat:quality'] == 'usable') & df['@type'].isin(rsfMRI_filter), :].groupby('subject.id'))\n",
    "    ]\n",
    "\n",
    "# Build the table\n",
    "pd.DataFrame(index=['All', 'DTI', 'T1', 'rs_fMRI'], data=\n",
    "        {\n",
    "        'scans': nb_scans,\n",
    "        'scans_validated': nb_scans_validated,\n",
    "        'scans_goodhead': nb_scans_goodhead,\n",
    "        'scans_useable': nb_scans_useable,\n",
    "        'subjects': nb_subjects,\n",
    "        'subjects_validated': nb_subjects_validated,\n",
    "        'subjects_goodhead': nb_subjects_goodhead,\n",
    "        'subjects_useable': nb_subjects_useable,\n",
    "        }\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_categorical(df_rest_csg, 'project.id', 'subject.@group', groupby='subject.id', title='WP15 Task3 Data Availability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------\n",
    "### Under work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------------------------------\n",
    "## Resources\n",
    "Useful resources:\n",
    "* https://stackoverflow.com/questions/34092808/extract-nested-json-embedded-as-string-in-pandas-dataframe#\n",
    "* https://stackoverflow.com/questions/39899005/how-to-flatten-a-pandas-dataframe-with-some-columns-as-json\n",
    "* http://mindtrove.info/flatten-nested-json-with-pandas/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------------------------------------------\n",
    "## Random code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flatten_json(y):\n",
    "    # From https://medium.com/@amirziai/flattening-json-objects-in-python-f5343c794b10\n",
    "    out = {}\n",
    "\n",
    "    def flatten(x, name=''):\n",
    "        if type(x) is dict:\n",
    "            for a in x:\n",
    "                flatten(x[a], name + a + '_')\n",
    "        elif type(x) is list:\n",
    "            i = 0\n",
    "            for a in x:\n",
    "                flatten(a, name + str(i) + '_')\n",
    "                i += 1\n",
    "        else:\n",
    "            out[name[:-1]] = x\n",
    "\n",
    "    flatten(y)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def JSON2DICTParser(data):\n",
    "    # From https://stackoverflow.com/questions/20680272/parsing-a-json-string-which-was-loaded-from-a-csv-using-pandas\n",
    "    if data and re.sub('[\\[\\]{}\\s]+', '', data):\n",
    "        try:\n",
    "            j1 = json.loads(data)\n",
    "        except ValueError:\n",
    "            print(data)\n",
    "            raise\n",
    "        return j1\n",
    "    else:\n",
    "        return float('nan')\n",
    "\n",
    "# Load the csv into a Pandas DataFrame\n",
    "#df = pd.read_csv(json_filepath, converters={'software_metadata':JSON2DICTParser, 'parameters_dict':JSON2DICTParser, 'results_dict':JSON2DICTParser})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract scans non normalized type\n",
    "df_rest_notnormtype = df_rest[df_rest['@type'] != 'rs_fMRI']\n",
    "save_df_as_csv(df_rest_notnormtype.fillna(''), 'xnat_data_extract_rest-notnormtype.csv', index=False, encoding='utf-8')\n",
    "print(len(df_rest_notnormtype))\n",
    "df_rest_notnormtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract only scans with an exact TR (with a floating value)\n",
    "import numpy as np\n",
    "df_rest_floattr = df_rest[df_rest['xnat:parameters.xnat:tr'].astype('float').apply(lambda x: not x.is_integer() and not pd.isnull(x))]\n",
    "print('List of scanners with floating TR values:')\n",
    "print(df_rest_floattr['experiment.xnat:scanner'].value_counts())\n",
    "print('\\nList of unique floating TR values:')\n",
    "print(df_rest_floattr['xnat:parameters.xnat:tr'].unique())\n",
    "print('\\nList of subjects ids with floating TR value:')\n",
    "print(df_rest_floattr['subject.@ID'].unique())\n",
    "print('\\nList of experiments with floating TR values:')\n",
    "print(df_rest_floattr.ix[:, ('experiment.id', 'subject.id', 'xnat:parameters.xnat:tr')])\n",
    "save_df_as_csv(df_rest_floattr.ix[:, ('experiment.id', 'subject.id', 'xnat:parameters.xnat:tr')].fillna(''), 'xnat_data_extract_rest-floattr.csv', index=False, encoding='utf-8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print all scanners unique names with scans counts (to compare with the scanners with floating TR)\n",
    "df_rest['experiment.xnat:scanner'].value_counts()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}