DEV_Atlas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=top></a>\n",
    "# DEV: Atlas Construction\n",
    "\n",
    "This notebook served the purpose of model selection and hyperparameter optimization for atlas prediction across channels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "\n",
    "----\n",
    "\n",
    "1. [Preparations](#prep)\n",
    "    1. [Imports & Data Loading](#prep_imports_loading)\n",
    "    2. [Preprocessing](#prep_preprocessing)\n",
    "2. [Evaluation of Regressors](#regressors)\n",
    "    1. [Baseline: Random Assignment (Dummy)](#reg_random)\n",
    "    2. [k Nearest Neighbors](#reg_knn)\n",
    "    3. [RBF MO-SVR](#reg_svr)\n",
    "    4. [MT-Lasso](#reg_lasso)\n",
    "    5. [MT-ElasticNet](#reg_enet)\n",
    "    6. [Random Forest](#reg_forest)\n",
    "    7. [MultiLayer Perceptron](#reg_mlp)\n",
    "    8. [Bagging Ensemble (SVR)](#reg_bag)\n",
    "    9. [Gradient Boosting Regressor (sklearn)](#sk_reg_boost)\n",
    "    10. [Gradient Boosting Regressor (xgboost)](#xg_reg_boost)\n",
    "3. [Comparative Assessment](#comp)\n",
    "    1. [Current Run](#assess_this_run)\n",
    "    2. [All Runs](#assess_all_runs)\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=prep></a>\n",
    "\n",
    "## 1. Preparations\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports & Data Loading <a id=prep_imports_loading></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import modules\n",
    "\n",
    "# External, general\n",
    "from __future__ import division\n",
    "import os, sys, pickle\n",
    "import numpy as np\n",
    "np.random.seed(42)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# External, specific\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from sklearn import preprocessing as sklearn_prep\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn import model_selection, metrics, base, dummy\n",
    "from sklearn import neighbors, svm, multioutput, tree, linear_model, ensemble, neural_network\n",
    "import xgboost\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "# Internal\n",
    "import katachi.utilities.loading as ld\n",
    "import katachi.utilities.atlas_helpers as atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Configuration\n",
    "\n",
    "# Target channel\n",
    "channel_2  = \"NLStdTomato\"\n",
    "#channel_2 =  \"tagRFPtUtrCH\"\n",
    "#channel_2 = \"mKate2GM130\"\n",
    "\n",
    "# Source feature space\n",
    "space_type_source = \"TFOR\"\n",
    "#space_type_source = \"CFOR\"\n",
    "#space_type_source = \"TFOR+CFOR\"\n",
    "\n",
    "# Target feature space\n",
    "space_type_target = \"TFOR\"\n",
    "#space_type_target = \"CFOR\"\n",
    "\n",
    "# Preprocessing steps\n",
    "#shape_type    = \"raw\"\n",
    "#shape_type    = \"z_normed\"\n",
    "shape_type    = \"pca\"\n",
    "#sec_type      = \"raw\"\n",
    "#sec_type      = \"z_normed\"\n",
    "sec_type      = \"pca\"\n",
    "restd_shape   = False\n",
    "restd_sec     = False\n",
    "\n",
    "# Removals and additions\n",
    "shape_num_PCs = 20\n",
    "sec_num_PCs   = 20\n",
    "\n",
    "# Cross-validation\n",
    "num_CVs = 3\n",
    "\n",
    "# Hyperparameter optimization\n",
    "score_weighted = True  # Whether to weigh the multioutput scores by variance when aggregating them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Function to parse relevant IDs from IDR bulk data\n",
    "\n",
    "def parse_from_IDR(dir_path, target):\n",
    "    \n",
    "    # Get all samples\n",
    "    samples = [d for d in os.listdir(dir_path) if len(d)==10\n",
    "               and os.path.isdir(os.path.join(dir_path, d))]\n",
    "    \n",
    "    # Select relevant samples\n",
    "    relevant_samples = []\n",
    "    for d in samples:\n",
    "        \n",
    "        # Get measured files\n",
    "        files = [i for i in os.listdir(os.path.join(dir_path, d))\n",
    "                 if i.startswith(d) and i.endswith('measured.tsv')]\n",
    "            \n",
    "        # Find target\n",
    "        if any([target in f for f in files]):\n",
    "            relevant_samples.append(d)\n",
    "    \n",
    "    return relevant_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load data\n",
    "\n",
    "# Target dir and target IDs\n",
    "dirpath = r'data/experimentA/extracted_measurements/'\n",
    "\n",
    "# Target IDs\n",
    "prim_IDs = parse_from_IDR(dirpath, channel_2)\n",
    "print \"Found %i training IDs!\" % len(prim_IDs)\n",
    "\n",
    "# Prep loader\n",
    "loader = ld.DataLoaderIDR(dirpath, recurse=True, verbose=True)\n",
    "\n",
    "# Load shape space from IDR\n",
    "if space_type_source == \"TFOR\":\n",
    "    shape_fspace, _, fspace_idx = loader.load_dataset(\"shape_TFOR_raw_measured.tsv\", IDs=prim_IDs)\n",
    "elif space_type_source == \"CFOR\":\n",
    "    shape_fspace, _, fspace_idx = loader.load_dataset(\"shape_CFOR_raw_measured.tsv\", IDs=prim_IDs)\n",
    "elif space_type_source == \"TFOR+CFOR\":\n",
    "    shape_fspace_TFOR, _, fspace_idx = loader.load_dataset(\"shape_TFOR_raw_measured.tsv\", IDs=prim_IDs)\n",
    "    shape_fspace_CFOR, _, _ = loader.load_dataset(\"shape_CFOR_raw_measured.tsv\", IDs=prim_IDs)\n",
    "    shape_fspace = np.concatenate([shape_fspace_TFOR, shape_fspace_CFOR], axis=1)\n",
    "else:\n",
    "    raise IOError(\"Invalid space_type_source.\")\n",
    "\n",
    "# Load target channel space from IDR\n",
    "if space_type_target == \"TFOR\":\n",
    "    sec_fspace, _, _ = loader.load_dataset(channel_2+\"_TFOR_raw_measured.tsv\")\n",
    "elif space_type_target == \"CFOR\":\n",
    "    sec_fspace, _, _ = loader.load_dataset(channel_2+\"_CFOR_raw_measured.tsv\")\n",
    "\n",
    "# Report\n",
    "print \"shape_fspace.shape:\", shape_fspace.shape\n",
    "print \"sec_fspace.shape:\", sec_fspace.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing <a id=prep_preprocessing></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Outlier Removal\n",
    "\n",
    "# NOTE: This is a rather simplistic approach but it gets rid of the few really weird cases.\n",
    "\n",
    "# Identify outliers based on extreme absolute values across the standardized sec_fspace\n",
    "outlier_fspace = np.abs(sklearn_prep.StandardScaler().fit_transform(sec_fspace))\n",
    "outliers_sec = outlier_fspace.sum(axis=1) > np.percentile(outlier_fspace.sum(axis=1), 95)\n",
    "\n",
    "# Identify outliers based on extreme absolute values across the standardized shape_fspace\n",
    "outlier_fspace = np.abs(sklearn_prep.StandardScaler().fit_transform(shape_fspace))\n",
    "outliers_shape = outlier_fspace.sum(axis=1) > np.percentile(outlier_fspace.sum(axis=1), 95)\n",
    "\n",
    "# Merge outlier mask\n",
    "outliers = outliers_sec | outliers_shape\n",
    "\n",
    "# Remove outliers from all imported datasets\n",
    "print \"Removing\", np.sum(outliers), \"outliers!\"\n",
    "shape_fspace = shape_fspace[~outliers]\n",
    "sec_fspace   = sec_fspace[~outliers]\n",
    "fspace_idx   = fspace_idx[~outliers]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Standardization & PCA\n",
    "\n",
    "# Shape space\n",
    "shape_fspace_z = sklearn_prep.StandardScaler().fit_transform(shape_fspace)\n",
    "shape_fspace_pca = PCA().fit_transform(shape_fspace_z)\n",
    "if restd_shape:\n",
    "    shape_fspace_pca = sklearn_prep.StandardScaler().fit_transform(shape_fspace_pca)\n",
    "\n",
    "# Secondary channel space\n",
    "sec_fspace_z   = sklearn_prep.StandardScaler().fit_transform(sec_fspace)\n",
    "sec_fspace_pca = PCA().fit_transform(sec_fspace_z)\n",
    "if restd_sec:\n",
    "    sec_fspace_pca = sklearn_prep.StandardScaler().fit_transform(sec_fspace_pca)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Selection & Splitting <a id=prep_ml></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Select features\n",
    "\n",
    "# Crop\n",
    "shape_fspace_pca = shape_fspace_pca[:, :shape_num_PCs]\n",
    "sec_fspace_pca   = sec_fspace_pca[:, :sec_num_PCs]\n",
    "\n",
    "# Decide which shape space to use\n",
    "if shape_type == \"raw\":\n",
    "    X = shape_fspace\n",
    "elif shape_type == \"z_normed\":\n",
    "    X = shape_fspace_z\n",
    "elif shape_type == \"pca\":\n",
    "    X = shape_fspace_pca\n",
    "else:\n",
    "    raise IOError(\"Invalid shape_type. Must be 'z_normed' or 'pca'.\")\n",
    "    \n",
    "# Decide which secondary space to use\n",
    "if sec_type == \"raw\":\n",
    "    y = sec_fspace\n",
    "elif sec_type == \"z_normed\":\n",
    "    y = sec_fspace_z\n",
    "elif sec_type == \"pca\":\n",
    "    y = sec_fspace_pca\n",
    "else:\n",
    "    raise IOError(\"Invalid sec_type. Must be 'z_normed' or 'pca'.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Split data for cross-validation and evaluation\n",
    "\n",
    "# Shuffle split\n",
    "cv_sets = model_selection.ShuffleSplit(n_splits=num_CVs, test_size=0.3, random_state=42)\n",
    "\n",
    "# Prepare generic cross validation scorers\n",
    "scoring = {'eplained_variance'  : metrics.make_scorer(metrics.explained_variance_score),\n",
    "           'mean_squared_error' : metrics.make_scorer(metrics.mean_squared_error),\n",
    "           'r2_score'           : metrics.make_scorer(metrics.r2_score, multioutput='uniform_average'),\n",
    "           'r2_weighted'        : metrics.make_scorer(metrics.r2_score, multioutput='variance_weighted'),\n",
    "           'spearman_r_score'   : metrics.make_scorer(atlas.spearman_r_score),\n",
    "           'spearman_p_bonf'    : metrics.make_scorer(atlas.spearman_p_bonf)}\n",
    "\n",
    "# Whether to use the relevant-only r2 scorer\n",
    "if score_weighted:\n",
    "    hyperopt_scoring = metrics.make_scorer(metrics.r2_score, multioutput='variance_weighted')\n",
    "else:\n",
    "    hyperopt_scoring = metrics.make_scorer(metrics.r2_score, multioutput='uniform_average')\n",
    "    \n",
    "# Prepare a single train-test split for visualization\n",
    "out = model_selection.train_test_split(X, y, test_size=0.3, random_state=42)\n",
    "X_train, X_test, y_train, y_test = out\n",
    "\n",
    "# Prepare dict for result collection\n",
    "result_dict = {}\n",
    "\n",
    "# Report\n",
    "print \"Preparations complete!\"\n",
    "print \"  Final source fspace (full, train, test):\", X.shape, X_train.shape, X_test.shape\n",
    "print \"  Final target fspace (full, train, test):\", y.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=regressors></a>\n",
    "\n",
    "## 2. Evaluation of Regressors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline: Random Assignment (Dummy) <a id=reg_random></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "### Baseline: random assignment\n",
    "print \"BASELINE: RANDOM ASSIGNMENT\"\n",
    "\n",
    "# Create dummy regressor that randomly assigns stuff\n",
    "class RandomRegressor(base.BaseEstimator, base.RegressorMixin):\n",
    "    \"\"\"Dummy classifier that returns random samples of fit target values as predictions.\"\"\"\n",
    "    \n",
    "    def __init__(self, random_state=None):\n",
    "        if random_state:\n",
    "            np.random.seed(random_state)\n",
    "    \n",
    "    def fit(self, X, y):\n",
    "        self.y = y\n",
    "        return self\n",
    "\n",
    "    def predict(self, X, y=None):\n",
    "        choice = np.random.choice(np.arange(self.y.shape[0]), \n",
    "                                  size=X.shape[0])\n",
    "        return self.y[choice,:]\n",
    "\n",
    "# Prepare \"regressor\"\n",
    "baseline_random = RandomRegressor(random_state=42)\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(baseline_random, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"Random\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(baseline_random, \"Random Baseline\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### k Nearest Neighbors <a id=reg_knn></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Simple kNN regression\n",
    "\n",
    "# Prepare regressor\n",
    "knn_reg = neighbors.KNeighborsRegressor(n_neighbors=5, weights='distance', n_jobs=10)\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(knn_reg, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"kNN\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(knn_reg, \"k Nearest Neighbors\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RBF MO-SVR <a id=reg_svr></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Hyperparam screening for SVR\n",
    "\n",
    "# Param grid\n",
    "gd = 1.0 / X_test.shape[1]\n",
    "print 'gd:', gd, '\\n'\n",
    "param_grid = [{'estimator__C': [0.1, 1.0, 10.0, 20.0], \n",
    "               'estimator__epsilon': [0.01, 0.1, 0.5], \n",
    "               'estimator__gamma': [gd, gd*0.1, gd*0.01]}]\n",
    "\n",
    "# Prep regressor\n",
    "svr = svm.SVR(kernel='rbf')\n",
    "multi_svr = multioutput.MultiOutputRegressor(svr)\n",
    "\n",
    "# Run grid search\n",
    "clf = model_selection.GridSearchCV(multi_svr, param_grid, scoring=hyperopt_scoring, n_jobs=12, verbose=2)\n",
    "clf.fit(X, y)\n",
    "\n",
    "# Available outputs\n",
    "print \"\\nOutputs:\"\n",
    "print sorted(clf.cv_results_.keys())\n",
    "\n",
    "# key results\n",
    "print \"\\nResults:\"\n",
    "print clf.best_estimator_\n",
    "print clf.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multioutput (separate!) regression with RBF SVR\n",
    "\n",
    "## Get optimized regressor\n",
    "multi_svr = clf.best_estimator_\n",
    "\n",
    "# Manually prepare regressor\n",
    "#svr       = svm.SVR(kernel='rbf', C=10.0, epsilon=0.1,   # raw-NLS-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pca-NLS-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=10.0, epsilon=0.5,   # pcaW-NLS-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.01,  # raw-NLS-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=10.0, epsilon=0.5,   # pca-NLS-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pcaW-NLS-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=10.0, epsilon=0.1,   # raw-UtrCH-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pca-UtrCH-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.5,    # pcaW-UtrCH-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.1,   # raw-UtrCH-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pca-UtrCH-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.001)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.01,   # pcaW-UtrCH-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.001)\n",
    "#svr       = svm.SVR(kernel='rbf', C=10.0, epsilon=0.1,   # raw-GM130-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.1,   # pca-GM130-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.5,   # pcaW-GM130-TFOR-TFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=20.0, epsilon=0.1,   # raw-GM130-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pca-GM130-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.01)\n",
    "#svr       = svm.SVR(kernel='rbf', C=1.0, epsilon=0.5,    # pcaW-GM130-CFOR-CFOR\n",
    "#                    gamma=1.0 / X_test.shape[1] * 0.1)\n",
    "#multi_svr = multioutput.MultiOutputRegressor(svr, n_jobs=y_train.shape[1])\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_svr, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True, \n",
    "                                        n_jobs=num_CVs)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"SVR\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_svr, \"MO-SVR (RBF)\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot of predictability over target dimensions for SVR for publication\n",
    "\n",
    "# Fit & predict for scoring\n",
    "multi_svr.fit(X_train, y_train)\n",
    "y_train_pred = multi_svr.predict(X_train)\n",
    "y_test_pred = multi_svr.predict(X_test)\n",
    "\n",
    "# Compute score for each dimension\n",
    "score_train = [metrics.r2_score(y_train[:,dim], y_train_pred[:,dim])\n",
    "               for dim in range(y_train.shape[1])]\n",
    "score_test  = [metrics.r2_score(y_test[:,dim], y_test_pred[:,dim])\n",
    "               for dim in range(y_test.shape[1])]\n",
    "\n",
    "# Prep\n",
    "fig,ax = plt.subplots(1, 1, figsize=(12,3))\n",
    "\n",
    "# Plot scores for training set\n",
    "ax.plot(score_train, 'k.-', alpha=0.25, label='train')\n",
    "\n",
    "# Plot scores for test set\n",
    "ax.plot(score_test, 'k.-', label='test')\n",
    "\n",
    "# Cosmetics\n",
    "ax.legend(loc=1)\n",
    "ax.set_xticks(np.arange(0, y_train.shape[1], 1))\n",
    "ax.set_xlim([-0.2, y_train.shape[1]-0.8])\n",
    "ax.set_ylim([-0.1, 1.1])\n",
    "ax.set_xlabel('target dimensions', fontsize=16)\n",
    "ax.set_ylabel('r-squared score', fontsize=16)\n",
    "ax.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "# Done\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MT-Lasso <a id=reg_lasso></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multivariable Linear Regression by Multi-Task Lasso\n",
    "\n",
    "# Prepare regressor\n",
    "#multi_lasso = linear_model.MultiTaskLasso(alpha=1/(10e4*X_train.shape[0]), random_state=42)\n",
    "multi_lasso = linear_model.MultiTaskLassoCV(random_state=42, n_jobs=10)  # CV determines hyperparam alpha\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_lasso, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"Lasso\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_lasso, \"MT-Lasso CV\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MT-ElasticNec <a id=reg_enet></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multivariable Linear Regression by Multi-Task Elastic Net\n",
    "\n",
    "# Prepare regressor\n",
    "#multi_enet = linear_model.MultiTaskElasticNet(alpha=1/10e4, l1_ratio=0.5, random_state=52)\n",
    "multi_enet = linear_model.MultiTaskElasticNetCV(random_state=52, n_jobs=10)  # CV determines hyperparam alpha\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_enet, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"eNet\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_enet, \"MT-ElasticNet CV\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Forest <a id=reg_forest></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Hyperparam screening for Random Forest\n",
    "\n",
    "# Param grid\n",
    "param_grid = [{'bootstrap' : [True, False],\n",
    "               'n_estimators' : [10, 100, 500, 1000, 1500], \n",
    "               'max_depth' : [10, 50, 100, None],\n",
    "               'max_features' : ['auto', 'sqrt']}]\n",
    "\n",
    "# Prep regressor\n",
    "random_forest = ensemble.RandomForestRegressor(random_state=42)\n",
    "\n",
    "# Run grid search\n",
    "clf = model_selection.GridSearchCV(random_forest, param_grid, n_jobs=12, verbose=2, scoring=hyperopt_scoring)\n",
    "clf.fit(X, y)\n",
    "\n",
    "# Available outputs\n",
    "print \"\\nOutputs:\"\n",
    "print sorted(clf.cv_results_.keys())\n",
    "\n",
    "# key results\n",
    "print \"\\nResults:\"\n",
    "print clf.best_estimator_\n",
    "print clf.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multivariable Random Forest Regression\n",
    "\n",
    "# Prepare regressor\n",
    "random_forest = clf.best_estimator_\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # raw-NLS-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1000, max_depth=50,       # pca-NLS-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1000, max_depth=50,       # pcaW-NLS-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=10,       # raw-NLS-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1000, max_depth=50,       # pca-NLS-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1000, max_depth=50,       # pcaW-NLS-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # raw-UtrCH-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # pca-UtrCH-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # pcaW-UtrCH-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # raw-UtrCH-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=10,       # pca-UtrCH-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=10,       # pcaW-UtrCH-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # raw-GM130-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # pca-GM130-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=50,       # pcaW-GM130-TFOR-TFOR\n",
    "#                                               max_features='sqrt', bootstrap=False,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1000, max_depth=10,       # raw-GM130-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=10,       # pca-GM130-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "#random_forest = ensemble.RandomForestRegressor(n_estimators=1500, max_depth=10,       # pcaW-GM130-CFOR-CFOR\n",
    "#                                               max_features='sqrt', bootstrap=True,\n",
    "#                                               random_state=42, n_jobs=10)\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(random_forest, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"Forest\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(random_forest, \"Random Forest\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MultiLayer Perceptron <a id=reg_mlp></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Hyperparam screening for MLP\n",
    "\n",
    "# Param grid\n",
    "param_grid = [{'activation' : ['logistic', 'tanh', 'relu'],\n",
    "               'solver'     : ['lbfgs', 'adam'],\n",
    "               'alpha'      : [0.001, 0.0001, 0.00001]}]\n",
    "# Note: This does not optimize the  hidden_layer_sizes but some preliminary \n",
    "#       testing showed that it's no use; it just increases overfitting\n",
    "\n",
    "# Prep regressor\n",
    "perceptron = neural_network.MLPRegressor(random_state=42, max_iter=1000)\n",
    "\n",
    "# Run grid search\n",
    "clf = model_selection.GridSearchCV(perceptron, param_grid, n_jobs=12, verbose=1, scoring=hyperopt_scoring)\n",
    "clf.fit(X, y)\n",
    "\n",
    "# Available outputs\n",
    "print \"\\nOutputs:\"\n",
    "print sorted(clf.cv_results_.keys())\n",
    "\n",
    "# key results\n",
    "print \"\\nResults:\"\n",
    "print clf.best_estimator_\n",
    "print clf.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multivariable Regression by Multi-Layer-Perceptron\n",
    "\n",
    "# Prepare regressor\n",
    "perceptron = clf.best_estimator_\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-5,    # raw-NLS-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='relu', alpha=1e-4,        # pca-NLS-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='relu', alpha=1e-4,        # pcaW-NLS-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-4,    # raw-NLS-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pca-NLS-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pcaW-NLS-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='tanh', alpha=1e-3,        # raw-UtrCH-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='relu', alpha=1e-5,        # pca-UtrCH-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='relu', alpha=1e-5,        # pcaW-UtrCH-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-5,    # raw-UtrCH-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pca-UtrCH-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pcaW-UtrCH-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='relu', alpha=1e-3,        # raw-GM130-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pca-GM130-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pcaW-GM130-TFOR-TFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-5,    # raw-GM130-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pca-GM130-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "#perceptron = neural_network.MLPRegressor(activation='logistic', alpha=1e-3,    # pcaW-GM130-CFOR-CFOR\n",
    "#                                         solver='adam', random_state=42)\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(perceptron, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"MLP\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(perceptron, \"MultiLayer Perceptron\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bagging Ensemble (SVR) <a id=reg_bag></a>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Ensemble prediction using Bagging (using best option: SVR)\n",
    "\n",
    "# Prepare regressor\n",
    "svr_bag = ensemble.BaggingRegressor(base_estimator=svr, n_estimators=5, \n",
    "                                    max_samples=0.8, max_features=0.8, bootstrap=False,\n",
    "                                    random_state=43, n_jobs=10)\n",
    "multi_svr_bag = multioutput.MultiOutputRegressor(svr_bag)\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_svr_bag, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"baggingSVR\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_svr_bag, \"(MO) Bagging Ensemble of SVR (RBF)\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Boosting Regressor (sklearn) <a id=sk_reg_boost></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Hyperparam screening for GBR\n",
    "\n",
    "# Param grid\n",
    "param_grid = [{'estimator__learning_rate'  : [0.001, 0.01, 0.1], \n",
    "               'estimator__n_estimators'   : [100, 300, 600],\n",
    "               'estimator__subsample'      : [0.1, 0.5, 1.0],\n",
    "               'estimator__max_depth'      : [3, 5, 7]}]\n",
    "\n",
    "# Prep regressor\n",
    "boost       = ensemble.GradientBoostingRegressor(random_state=42)\n",
    "multi_boost = multioutput.MultiOutputRegressor(boost, n_jobs=y_train.shape[1])\n",
    "\n",
    "# Run grid search\n",
    "clf = model_selection.GridSearchCV(multi_boost, param_grid, n_jobs=14, verbose=2, scoring=hyperopt_scoring)\n",
    "clf.fit(X, y)\n",
    "\n",
    "# Available outputs\n",
    "print \"\\nOutputs:\"\n",
    "print sorted(clf.cv_results_.keys())\n",
    "\n",
    "# key results\n",
    "print \"\\nResults:\"\n",
    "print clf.best_estimator_\n",
    "print clf.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multioutput (separate!) regression with Gradient Boosting Regressor\n",
    "\n",
    "# Grab optimized regressor\n",
    "multi_boost = clf.best_estimator_\n",
    "\n",
    "## Manually prepare regressor\n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'huber', learning_rate=0.01,  # raw-NLS-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=5, random_state = 42)\n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-NLS-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=5, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-NLS-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=5, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # raw-NLS-CFOR-CFOR\n",
    "#                                                 n_estimators=300, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-NLS-CFOR-CFOR\n",
    "#                                                 n_estimators=100, subsample=0.5,\n",
    "#                                                 max_depth=5, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-NLS-CFOR-CFOR\n",
    "#                                                 n_estimators=300, subsample=0.5,\n",
    "#                                                 max_depth=5, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # raw-UtrCH-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-UtrCH-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-UtrCH-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # raw-UtrCH-CFOR-CFOR\n",
    "#                                                 n_estimators=300, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-UtrCH-CFOR-CFOR\n",
    "#                                                 n_estimators=100, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-UtrCH-CFOR-CFOR\n",
    "#                                                 n_estimators=300, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # raw-GM130-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-GM130-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-GM130-TFOR-TFOR\n",
    "#                                                 n_estimators=600, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # raw-GM130-CFOR-CFOR\n",
    "#                                                 n_estimators=100, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pca-GM130-CFOR-CFOR\n",
    "#                                                 n_estimators=100, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#boost       = ensemble.GradientBoostingRegressor(loss = 'ls', learning_rate=0.01,     # pcaW-GM130-CFOR-CFOR\n",
    "#                                                 n_estimators=100, subsample=0.5,\n",
    "#                                                 max_depth=3, random_state = 42) \n",
    "#multi_boost = multioutput.MultiOutputRegressor(boost, n_jobs=y_train.shape[1])\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_boost, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True, \n",
    "                                        n_jobs=num_CVs)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"boost\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_boost, \"MO-GBR\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Boosting Regressor (xgboost) <a id=xg_reg_boost></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Multivariate-Multioutput (separate!) regression with XGBoost\n",
    "\n",
    "# Prepare regressor\n",
    "#xgb = xgboost.XGBRegressor(random_state = 42)      # Default hyperparams\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-NLS-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 5)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-NLS-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 5)\n",
    "xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-NLS-TFOR-TFOR\n",
    "                           learning_rate  = 0.01,\n",
    "                           n_estimators   = 600,\n",
    "                           subsample      = 0.5,\n",
    "                           max_depth      = 5)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-NLS-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 300,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-NLS-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 100,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 5)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-NLS-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 300,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 5)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-UtrCH-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-UtrCH-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-UtrCH-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-UtrCH-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 300,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-UtrCH-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 100,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-UtrCH-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 300,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-GM130-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-GM130-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-GM130-TFOR-TFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 600,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # raw-GM130-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 100,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pca-GM130-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 100,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "#xgb = xgboost.XGBRegressor(random_state   = 42,     # pcaW-GM130-CFOR-CFOR\n",
    "#                           learning_rate  = 0.01,\n",
    "#                           n_estimators   = 100,\n",
    "#                           subsample      = 0.5,\n",
    "#                           max_depth      = 3)\n",
    "\n",
    "multi_xgb = multioutput.MultiOutputRegressor(xgb, n_jobs=y_train.shape[1])\n",
    "\n",
    "# Perform cross-validation\n",
    "print \"\\nPerforming cross validation...\"\n",
    "scores = model_selection.cross_validate(multi_xgb, X, y, scoring=scoring,\n",
    "                                        cv=cv_sets, return_train_score=True, \n",
    "                                        n_jobs=num_CVs)\n",
    "atlas.report_cv_scores(scores)\n",
    "result_dict[\"xgb\"] = scores\n",
    "\n",
    "# Plot result of an example fit\n",
    "report = atlas.visualize_regression(multi_boost, \"MO-XGB\",\n",
    "                                    X_train, X_test, y_train, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=comp></a>\n",
    "\n",
    "## 3. Comparative Assessment\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Current Run <a id=assess_this_run></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot test R^2 values\n",
    "\n",
    "# Interactive choice of score\n",
    "from ipywidgets import interact\n",
    "@interact(score=sorted(result_dict.values()[0].keys()))\n",
    "def plot_scores(score=\"test_r2_score\"):\n",
    "    \n",
    "    # Prep\n",
    "    fig   = plt.figure(figsize=(6,2))\n",
    "    width = 0.75\n",
    "\n",
    "    # Get relevant values\n",
    "    score_means = np.array([np.mean(s[score]) for s in result_dict.values()])\n",
    "    score_vals  = np.array([s[score] for s in result_dict.values()])\n",
    "    score_names = np.array(result_dict.keys())\n",
    "\n",
    "    # Sort by mean\n",
    "    score_vals  = score_vals[np.argsort(score_means), :]\n",
    "    score_names = score_names[np.argsort(score_means)]\n",
    "    score_means = np.sort(score_means)\n",
    "\n",
    "    # Plot bars\n",
    "    plt.bar(range(len(score_means)), score_means,\n",
    "            width=width, color='skyblue', edgecolor='', zorder=0)\n",
    "\n",
    "    # Add scatter\n",
    "    if not score==\"score_time\":\n",
    "        for i in range(len(score_means)):  \n",
    "            plt.scatter([i+width/2 for j in score_vals[i,:]], score_vals[i,:], \n",
    "                        color='0.3', alpha=0.5, edgecolor='', s=10)\n",
    "\n",
    "    # Cosmetics\n",
    "    plt.xlim([-width/2, len(score_means)-1+width+width/2])\n",
    "    plt.ylim([0.0, max([1.1, score_vals.flatten().max()+0.1*score_vals.flatten().max()])])\n",
    "    plt.xticks(np.arange(len(score_means))+width/2, score_names, \n",
    "               ha=\"right\", fontsize=8, rotation=45)\n",
    "    plt.ylabel(score)\n",
    "    plt.title(\"Performance Comparison\")\n",
    "\n",
    "    # Show\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Save result\n",
    "\n",
    "# Save\n",
    "with open(\"other/AtlasScores_\" \n",
    "          + shape_type + ('W' if score_weighted else '') + \"_\"\n",
    "          + channel_2 + \"_\" \n",
    "          + space_type_source + '_to_' \n",
    "          + space_type_target \n",
    "          + \".pkl\", \"wb\") as outfile:\n",
    "    pickle.dump(result_dict, outfile, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All Runs <a id=assess_all_runs></a>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Load all existing results\n",
    "\n",
    "# Get existing result files\n",
    "fnames = [fname for fname in os.listdir(\"other/\")\n",
    "          if fname.startswith(\"AtlasScores_\") \n",
    "          and fname.endswith(\".pkl\")]\n",
    "\n",
    "# Load existing results\n",
    "all_results = {}\n",
    "for fname in fnames:\n",
    "    with open(os.path.join(\"other/\", fname), \"rb\") as infile:\n",
    "        all_results[fname[7:-4]] = pickle.load(infile)\n",
    "\n",
    "# Report\n",
    "incl = len(all_results.keys())\n",
    "print \"Loaded\", incl, \"results:\"\n",
    "for key in all_results.keys(): print \" \", key"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "### Plot test R^2 values\n",
    "\n",
    "# Interactive choice of score\n",
    "from ipywidgets import interact\n",
    "@interact(score=sorted(all_results.values()[0].values()[0].keys()))\n",
    "def plot_scores(score=\"test_r2_weighted\"):\n",
    "    \n",
    "    # Prep\n",
    "    fig, axs = plt.subplots(int(np.ceil(len(all_results.keys())/2)), 2,\n",
    "                           figsize=(12,int(2*np.ceil(len(all_results.keys())/2))))\n",
    "    width = 0.75\n",
    "\n",
    "    # For each result...\n",
    "    for ax, key in zip(axs.flatten(), sorted(all_results.keys())):  \n",
    "        result_dict = all_results[key]\n",
    "            \n",
    "        # Get relevant values\n",
    "        try:\n",
    "            score_means = np.array([np.mean(s[score]) for s in result_dict.values()])\n",
    "            score_vals  = np.array([s[score] for s in result_dict.values()])\n",
    "            score_names = np.array(result_dict.keys())\n",
    "        except KeyError:\n",
    "            if 'r2_weighted' in score:\n",
    "                score = score.replace('weighted', 'score')\n",
    "                score_means = np.array([np.mean(s[score]) for s in result_dict.values()])\n",
    "                score_vals  = np.array([s[score] for s in result_dict.values()])\n",
    "                score_names = np.array(result_dict.keys())\n",
    "            else:\n",
    "                raise\n",
    "\n",
    "        # Sort by mean\n",
    "        score_vals  = score_vals[np.argsort(score_means), :]\n",
    "        score_names = score_names[np.argsort(score_means)]\n",
    "        score_means = np.sort(score_means)\n",
    "\n",
    "        # Plot bars\n",
    "        ax.bar(range(len(score_means)), score_means,\n",
    "               width=width, color='skyblue', edgecolor='', zorder=0)\n",
    "\n",
    "        # Add scatter\n",
    "        if not score==\"score_time\":\n",
    "            for i in range(len(score_means)):  \n",
    "                ax.scatter([i+width/2 for j in score_vals[i,:]], score_vals[i,:], \n",
    "                            color='0.3', alpha=0.5, edgecolor='', s=10)\n",
    "\n",
    "        # Cosmetics\n",
    "        ax.set_xlim([-width/2, len(score_means)-1+width+width/2])\n",
    "        ax.set_ylim([0.0, max([1.1, score_vals.flatten().max()+0.1*score_vals.flatten().max()])])\n",
    "        ax.set_xticks(np.arange(len(score_means))+width/2)\n",
    "        ax.set_xticklabels(score_names, ha=\"right\", fontsize=8, rotation=45)\n",
    "        ax.set_ylabel(score)\n",
    "        ax.set_title(key)\n",
    "\n",
    "    # Show & save\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "[Back to Top](#top)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}