Pyrcca_usage_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Usage example of pyrcca\n",
    "\n",
    "In this example, Pyrcca is used to find linear relationships between two artificially constructed interdependent datasets. Both datasets have 1000 samples. One dataset has 4 dimensions, and the other has 5 dimensions. These datasets depend on two latent variables. Each dimension in both datasets is constructed as a sum of an independent component (25%) and one of the two latent variables (75%). Each dataset is split in two halves: a training set and a test set. \n",
    "\n",
    "Pyrcca is used to find a CCA mapping between the training sets. The quality of the mapping is assessed using cross-dataset prediction with test data. Each test set is projected onto the canonical space, and the other test set is predicted using the CCA mapping. The prediction accuracy is quantified by correlating predicted test data with actual test data.\n",
    "\n",
    "The analysis is first run with predefined hyperparameters. Then, grid search with cross-validation is used to find optimal hyperparameters.\n",
    "\n",
    "This example demonstrates that Pyrcca can be used to recover linear relationships for two artificially generated datasets. It is designed to be easily adaptable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pyrcca with predefined hyperparameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import rcca\n",
    "\n",
    "# Initialize number of samples\n",
    "nSamples = 1000\n",
    "\n",
    "# Define two latent variables (number of samples x 1)\n",
    "latvar1 = np.random.randn(nSamples,)\n",
    "latvar2 = np.random.randn(nSamples,)\n",
    "\n",
    "# Define independent components for each dataset (number of observations x dataset dimensions)\n",
    "indep1 = np.random.randn(nSamples, 4)\n",
    "indep2 = np.random.randn(nSamples, 5)\n",
    "\n",
    "# Create two datasets, with each dimension composed as a sum of 75% one of the latent variables and 25% independent component\n",
    "data1 = 0.25*indep1 + 0.75*np.vstack((latvar1, latvar2, latvar1, latvar2)).T\n",
    "data2 = 0.25*indep2 + 0.75*np.vstack((latvar1, latvar2, latvar1, latvar2, latvar1)).T\n",
    "\n",
    "# Split each dataset into two halves: training set and test set\n",
    "train1 = data1[:nSamples/2]\n",
    "train2 = data2[:nSamples/2]\n",
    "test1 = data1[nSamples/2:]\n",
    "test2 = data2[nSamples/2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis\n",
    "There are many more samples than dimensions in the dataset, therefore kernelization and regularization are unnecessary.\n",
    "\n",
    "There are two latent variables, and CCA is expected to find two underlying canonical components that capture the relationships between the datasets. However, it is possible to compute up to four canonical components, so the analysis is first used to find four components. This is done in order to test whether the first two components capture the relationships present in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training CCA, regularization = 0.0000, 4 components\n"
     ]
    }
   ],
   "source": [
    "# Create a cca object as an instantiation of the CCA object class. \n",
    "nComponents = 4\n",
    "cca = rcca.CCA(kernelcca = False, reg = 0., numCC = nComponents)\n",
    "\n",
    "# Use the train() method to find a CCA mapping between the two training sets.\n",
    "cca.train([train1, train2])\n",
    "\n",
    "# Use the validate() method to test how well the CCA mapping generalizes to the test data.\n",
    "# For each dimension in the test data, correlations between predicted and actual data are computed.\n",
    "testcorrs = cca.validate([test1, test2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import matplotlib.pyplot as plt\n",
    "from brewer2mpl import qualitative\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Canonical correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The canonical correlations are:\n",
      "\n",
      "Component 1: 0.95\n",
      "\n",
      "Component 2: 0.95\n",
      "\n",
      "Component 3: 0.10\n",
      "\n",
      "Component 4: 0.00\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEZCAYAAACQK04eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGglJREFUeJzt3XucnHV96PHPkkBWrlGpoBBcDZdq5SJQRIHjoBzDpYeU\neqGIFrAvpedIEns5UqyWrVUr5xwBEw9CLSL2UGNPiRa6lkgxq8hNCGG5JkIUuaMgEm4JCdn+8f0N\nOzuZnX1mZ56ZnXk+79drX/vc5pnvTuD7/Ob7/J7fDyRJkiRJkiRJkiRJkiRJkiSppc4CvtbkOQaA\nzcBWTUfTWiXgwSZe/1Xg060JRVJRfRC4BXgGeAT4HnBYRyNqjQG6P/GfClybWyTqWtPtP2p1lz8D\nzgM+B7wGmAP8X+D4TgbVRWbW2Daj7VFIUkY7Ea3899Y55hDgBuAp4tvAEmDriv2bgdOBn6ZjvlKx\nr48oSdwPPA5cCuyY9g2k1/4R8AvgV8CnKl47CPxjxfrhwPXpPR4ATknbjwNWAU+n7WdXvKb8HhM1\njuYAy4BfAk+kv410/GRxfyTF/cMUy3XAuek8nwW2Af5POuYxojzTn85RYnyL/y+B+4B1wF3A76ft\nbwJeADYR/06/Ttu/Afxtxes/CtwLPAn8K/Dain31/n32TPH/hvj8l9b4jCT1mKOBjdT/1nggkfy3\nAl4P3A0sqti/GbiCSIxziCQ6L+37CJGQBoDtgMuBb6Z9A+m1FwGzgP2A9cA+af/ZjCX+1xNJ8USi\nNf0qYP+0753A76TlfYkkO7/qPWr9fTOAEeBLwCtSDO9oIO5vpNf1E+WYjcDH03v1E9+ivgvMBrZP\nn9EX0jlKjE/87wN2TcsfAJ4Fdknrp7BlqecS4uIC8C4iaR9AXGwWE8m8rNa/z3vSvm8R91JIr30H\nknreycCjDb7mE0QruWwz4xPGt4FPpuVrgD+p2Lc38CKRHAfSa19Xsf8mIvHB+Bb/WUTyzeJ8ouUN\n9RP/24kkWGtflrgHKvafSrTsy/qI5P3Gqvf7WVouUb/Gv4qxUtup1E/8FwNfrNi3XYp1j7Re79/n\nUuLCu1udWDRNWePXVD0J7Ez9/4b2Bv6NuEA8DXweeHXVMY9VLD9PtHAhSg6VCfEBoia+S8W2iV5b\naQ5jSbPa24AVRBL/DVHWqI6vljkpts019mWJuzpxV67/FrAtsJIorzwF/DvxWdfyR0SyLx/7lox/\nQ61YnyP+XSuTefVnvENa/iRxkfoJcCdwWsb31DRg4tdU3QBsAE6oc8xXifLOnsQ9gb8i+39zjzC+\nZbwHUa9+vME4HwDmTrDvn4iSyu5EWeXCjPE9mOKpdSM2S9yjVa+pXH+CqM2/GXhl+pnN2H2CSq8H\n/p4oE70qHXsnkZBrvc9ksW5HXDQenuR1EH/Px4iLxOnABYz/lqJpzMSvqXoa+GuiF898opW6NXAM\ncE46ZnvixuLzwG8D/32Sc/YxlrS+BfwpkZi2J2rcS6ndyq7nn4CjgPcTLe9XM1bj355oJb9I3Iv4\nIJMnS4iy0qNEmWRboi5fLok0G/dm4hmE84nWP0RyfU+NY7dL8T5B/L98GtHiL3ucuKhV3lCv/oxP\nIz6PWSnWG4mLZS19FcvvT+eG+LY0SuP/NuoQE7+acS7RpfPTRLnkAeB/AN9J+/+CSKbriJbpUsYn\n1lot3/K2rxN1+h8RpZrngQV1XjvReR4AjgX+nChjrCJuBpNi/WyK7zNEDbv6PLVsBv4b8U3mAeIb\nQPn+QqNxj9bYdibRU+dG4gJ7NVE2qz7H3cQN5huIksxbgB9XHHcN0dPnMeLfp/r9riH+7suJ1v8b\ngD/MGOvBKb5niN5AC4meTBJfJ1odd9Q5ZjHRC2IEeGs7gpIk5ecIIplPlPiPJZ70hLjRdmM7gpIk\n5WuAiRP/hUT/6rLVjO/9IElqsU7X+HdjfFe2hxi7YSRJykGnEz+M7ykA2XpVSJKmqNYgUe30MPEw\nTNnu1OhDvP/++4+OjIy0LShJ6hEjxJAc43S6xX8F8eQhwKFEf+AtHtAZGRlhdHR02v+cffbZHY+h\nl378PP0sp+tPt3yejD2zMk7eLf5vEQNh7UzU8s9m7GGSi4gePccSfZafw8e+JSl3eSf+kzIcc0bO\nMUiSKnS61NNTSqVSp0PoKX6ereNn2Vrd/nlW96iZrkZTvUqSlFFfXx/UyPO2+CWpYEz8LTA0NMS8\nefMolUrMmzePoaGhTockSRPqdD/+rjc0NMSiRYtYu3bty9vKy8cdd1ynwupqQ0NDLF68mA0bNjBr\n1iwWLlzoZym1kIm/SYsXLx6X9CES/5IlS0xWU+CFVMqfpZ4mbdiwoeb29evXtzmS3lDvQiqpNUz8\nTZo1a1bN7f39/W2OpDd4IZXyZ+Jv0sKFC5k7d/yUrnPnzmXBggUTvEL1eCGV8meNv0nluvOSJUtY\nv349/f39LFiwwHr0FC1cuJC1a9eOK/d4IZVaywe4NO0MDQ15IZVaYKIHuEz8ktSjfHJXkgSY+CWp\ncEz8klQwJn5JKhgTvyQVjIlfkgrGxC9JBWPil6SCMfFLUsGY+CWpYEz8klQwJn5JKhgTvyQVjIlf\nkgrGxC9JBWPil6SCMfFLUsGY+CWpYEz8klQwJn5JKhgTvyQVjIlfkgrGxC9JBWPil6SCyTvxHw2s\nBu4Fzqyxf2fgKuA24E7g1JzjkaTC68vx3DOANcBRwMPAzcBJwD0VxwwCs4CziIvAGmAXYFPVuUZH\nR0dzDFWSek9fXx/UyPN5tvgPAe4D7gc2AkuB+VXHPArsmJZ3BJ5ky6QvSWqhmTmeezfgwYr1h4C3\nVR3zNeAHwCPADsAHcoxHkkS+iT9LbeZTRH2/BMwFrgb2B56pPnBwcPDl5VKpRKlUakGIktQ7hoeH\nGR4envS4PGv8hxI1/KPT+lnAZuCcimO+B3weuC6tX0PcBL6l6lzW+CWpQZ2o8d8C7AUMANsAJwJX\nVB2zmrj5C3FTdx/gZznGJEmFl2epZxNwBrCc6OFzMdGj5/S0/yLgC8AlwAhxEfok8OscY5Kkwsuz\n1NNKlnokqUGdKPVIkqYhE78kFYyJX5IKxsQvSQVj4pekgjHxS1LBmPglqWBM/JJUMCZ+SSoYE78k\nFYyJX5IKJkvify8xZ+46Ypz8Z9KyJKkLZRmkbS3we4yfK7fdHKRNkhrUzCBtj9HZpC9JaqEs4/Hf\nAnwb+C7wYto2CizLKyhJUn6yJP6dgBeA91RtN/FLUhdyIhZJ6lHN1PjnAN8BfpV+Lgd2b2VwkqT2\nyZL4LyEmSX9d+rkybZMkdaEspZ4RYP8M2/JkqUeSGtRMqedJ4MPADOJm8IeAJ1oZnCSpfbK0+AeA\nJcChaf16YAHwQE4x1WKLX5IaNFGL3149ktSjJkr89frxnwmcQ7T2q40CC1sSmSSpreol/rvT75VE\noi/rq1qXJHWReon/yvT7eeCfq/Z9IJ9wJEl5y1LjXwW8NcO2PFnjl6QGTaXGfwxwLLAbsLjixTsA\nG1scnySpTeol/keI+v789Luc+NcBf5pzXJKknGQp9WzD2HDMnWKpR5IaNJVST9kA8AXgzcAr0rZR\n4I0tik2S1EZZB2m7ENgElIBLgctyjEmSlKMspZ5bgQOBO4B9q7a1i6UeSWpQM6We9cQAbfcBZxA3\nfbdrZXCSpPbJ0uI/hJhsfTbwt8COwP8Cbswxrmq2+CWpQZ0apO1o4HziG8M/EGP/VCsB5wFbE8M9\nl2ocY+KXpAZNJfFfWWffKHD8JO85A1gDHAU8DNwMnER8eyibDVwHzAMeAnam9lj/Jn5JatBUavxf\nqrMvSxY+hLgvcH9aX0o8DFaZ+D9IzOH7UFp3ghdJylm9xD9csbwtMen6mgbOvRvwYMX6Q8Dbqo7Z\niyjxrCCGgvgy8I8NvIckqUFZ+vEfTwzKtjytv5WYfH0yWb4VbE10Cz2WKPd8hrgYSJJykqU75yDR\nUl+R1leR7andh4lvCWVzGCvplD1IlHdeSD8/IiZxv3eLIAYHX14ulUqUSqUMIUhScQwPDzM8PDzp\ncVl69dxEJP7KoZhvB/ab5HUzidLQu4m+/z9hy5u7vw18hWjtz0rvdSJjk8CUeXNXkhrUzANcdwEn\np2P3IqZcvD7D6zYRD3wtJ3r4XEwk/dPT/ouA1cBVxIVkM/A1tkz6kqQWytLi3xb4NPCetL6ceJBr\nfV5B1WCLX5IaNNUHuGYCVwNH5hBTI0z8ktSgiRL/ZL16NhElmNk5xCRJ6oAsNf7niJE5r07LEF01\nF+YVlCQpP1kS/+XAMsb65feRrY++JGkaylLj/w9qD5zWTtb4JalBzdT4X8IavyT1DGv8klQwWRL/\nMqzxS1LPyDoRyyxg77S8GtiYTzgTssYvSQ1qZsiGEnAp8Iu0vgdwCvDDFsUmSWqjLC3+W4nB1cpj\n8e9NTKpyYF5B1WCLX5IaNNVePTA2ymbZT8n2TUGSNA1lSeAriYnS/x9x5TgZuCXPoCRJ+clS6ukH\nPg4cltavBS4ANuQVVA2WeiSpQVMdnRNgO2II5pfS+gyil8/zrQouAxO/JDWomRr/D4BXVKxvSwzj\nIEnqQlkS/yzg2Yr1Z4jkL0nqQlkS/3PAQRXrBxMTo0uSulCWXj2fAP4ZeDStv5aYEF2S1IWyDtmw\nDbBPWl4DvJhPOBPy5q4kNaiZXj3TgYlfkhrUTK8eSVIPMfFLUsHUu7l7EPXH3b+1xbFIktqgXo1/\nmPqJ/8jWhlKXNX5JapA3dyWpYJqZiAVgX+BNxIBtZd9sPixJUrtlafEPAu8EfgcYAo4Bfgy8L7+w\ntmCLX5Ia1Ex3zvcBRxFP7p4G7A/MbmVwkqT2yZL4XyCGZN4E7AT8EpiTZ1CSpPxkqfHfDLwS+Box\n89ZzwPV5BiVJyk+jvXreAOwIjOQQSz3W+CWpQc3U+E9grKb/c+AXwO+3LDJJUltlafGPEDd0K90G\nHND6cCZki1+SGtRMi7/WxWFGswFJkjojS+JfCZwLzAX2BM5L2yRJXShL4l8AbAS+DSwF1gMfz3j+\no4HVwL3AmXWO+12iu+gfZDyvJGmK8hyrZwYxW9dRwMNEt9CTgHtqHHc18DxwCXB5jXNZ45ekBk1l\nrJ4vA4uAK2vsGwWOn+Q9DwHuA+5P60uB+WyZ+BcA/0K0+iVJOauX+MuDsH1piufeDXiwYv0h4G01\njpkPvItI/DbrJSln9RJ/+Qbu8BTPnSWJnw/8ZTq2j+4ZJlqSulaWIRsOB84GBiqOHwXeOMnrHmb8\nmD5ziFZ/pYOIEhDAzsTInxuBK6pPNjg4+PJyqVSiVCplCF2SimN4eJjh4eFJj8vSwl4DfIKYavGl\niu1PTPK6mem17wYeAX5C7Zu7ZZcQ9xOW1djnzV1JalAzE7H8Bvj3KbznJuAMYDnRc+diIumfnvZf\nNIVzSpKalKXF/0UicS8DNlRsb+dk67b4JalBzcy5O0ztG7VOti5J05iTrUtSwTQzSNtsxsbnWUn0\n69+plcFJktonS+L/OrAOeD/wAeAZogeOJKkLTXU8/lrb8mSpR5Ia1Eyp5wXgiIr1w4kB1SRJXShL\ni/8AYtyecl3/KeAU2jvvri1+SWpQK3r17Jh+r2tFQA0y8UtSg5p5crcfeC8xVs+MdJJR4LOtC0+S\n1C5ZEv+/EsM2rCRm35IkdbEspZ47gbfkHcgkLPVIUoOa6dVzPbBfqwOSJHVGlhb/PcCewM8ZG6Rt\nlPZeDGzxS1KDmunVMzDB9vunHk7DTPyS1KBmevXcn36/hujhI0nqYllq/McD9xKlnh8SF4KpTMwi\nSZoGsiT+zwFvB34KvIGYSvGmPIOSJOUnS+LfSMyvuxXxANcK4OA8g5Ik5SdLjf8pYAfgWuAy4JfA\ns3kGJUnKT5ZePdsRT+xuBZxMjNlzGfBkjnFVs1ePJDVoKt059wJ2AX5ctf1w4FFgbauCy8DEL0kN\nmsqTu+dTeyTOdWmfJKkL1Uv8uwC319h+O9G7R5LUheol/tl19vkglyR1qXqJ/xbgYzW2f5QYolmS\n1IXq3dzdFfgO8CJjif4gYBZwAnGDt128uStJDZrqIG19wJHEePyjwF3AD1odXAYmfklqUCvm3O0k\nE78kNaiZiVgkST3ExC9JBWPil6SCMfFLUsGY+CWpYEz8klQwJn5JKhgTvyQVTDsS/9HAamLC9jNr\n7D8ZGCFG/bwO2K8NMUlSYeX95O4MYA1wFPAwcDNwEnBPxTFvB+4GniYuEoPAoVXn8cldSWpQp57c\nPQS4D7ifmLR9KTC/6pgbiKQPcBOwe84xSVKh5Z34dwMerFh/KG2byB8D38s1IkkquJk5n7+R+syR\nwEeAw2rtHBwcfHm5VCpRKpWaiUuSes7w8DDDw8OTHpd3jf9QomZ/dFo/C9gMnFN13H7AsnTcfTXO\nY41fkhrUqRr/LcBewACwDXAicEXVMXsQSf9D1E76kqQWyrvUswk4A1hO9PC5mOjRc3rafxHw18Ar\nga+mbRuJm8KSpBw4EYsk9SgnYpEkASZ+SSocE78kFYyJX5IKxsQvSQVj4pekgjHxS1LBmPglqWBM\n/JJUMCZ+SSoYE78kFYyJX5IKxsQvSQVj4pekgjHxS1LBmPilHjY0NMS8efMolUrMmzePoaGhToek\naSDvGbgkdcjQ0BCLFi1i7dq1L28rLx933HGdCkvTgDNwST1q3rx5fP/736+5/aqrrupARGo3Z+CS\nCmbDhg01t69fv77NkWi6MfFLPWrWrFk1t/f397c5Ek03Jn6pRy1cuJC5c+eO2zZ37lwWLFjQoYg0\nXXhzV+pR5Ru4S5YsYf369fT397NgwQJv7Mqbu5LUq7y5K0kCTPySVDgmfkkqGBO/JBWMiV+SCsbE\nL0kFY+KXpIIx8UtSwZj4JalgTPySVDAmfkkqGBO/JBVM3on/aGA1cC9w5gTHLE77R4C35hyPJBVe\nnol/BvAVIvm/GTgJeFPVMccCewJ7AR8DvppjPLkbHh7udAg9xc+zdfwsW6M8ef0BBxzQ1ZPX5zke\n/yHAfcD9aX0pMB+4p+KY44FL0/JNwGxgF+DxHOPKzfDwMKVSqdNh9Aw/z9bxs2xeL01en2eLfzfg\nwYr1h9K2yY7ZPceYJGlKFi9ePC7pQyT+JUuWdCiiqcsz8WedOaV6kgBnXJE07Th5fTaHAldVrJ/F\nljd4LwT+sGJ9NVHqqXYbcUHwxx9//PEn+89ttNlMYC0wAGyTAqh1c/d7aflQ4MZ2BSdJyscxwBri\nJu9Zadvp6afsK2n/CHBgW6OTJEmSpCy+TnRBvaPTgfSAOcAK4C7gTmBhZ8Ppev1EV+nbgLuBv+ts\nOD1hBrAKuLLTgaizjiCeOjbxN29X4IC0vD1RKqy+N6TGbJt+zyTuox3ewVh6wZ8BlwFXdDqQqXKs\nnta4Fniq00H0iMcY64nwLPHA3+s6F05PeD793oZorf66g7F0u92JTin/wJZd0buGiV/T2QDxTeqm\nDsfR7bYiLqaPE2W0uzsbTlc7D/ifwOZOB9IME7+mq+2BfwEWES1/Td1mony2O/BfgFJHo+levwf8\nkqjvd21rX601gDX+VtkaWA58otOB9KDPAH/R6SC61BeIIWZ+DjwKPAd8s6MRqeMGMPG3Qh/xP9N5\nnQ6kR+xMDH4I8ArgR8C7OxdOz3gn9uopvG8BjwAbiBbBaZ0Np6sdTpQmbiO+Uq8ihvbW1OwL3Ep8\nnrcT9Wk17510ca8eSZIkSZIkSZIkSZIkSZIkSYWxK7CUmHjnFmAI2KtN730Q8OUpvnY4vb7XzceR\nUCW1UB9wA/Cxim370R1DBK+gGLPDfQN4b6eDkNQ73gX8cIJ92wH/Aawkni49Pm0fIIZl/ntiUpbl\nxMQiEIOO3UhM2bmMsWEJhoEvEiN6rmHswlJi7BH77YFL0nuNACek7RcAN6f3GqyIb6LE/7vAdcRT\nsTelv6O/4ty3MjYg2qnAd4HvE2O9nEGMmXMrcUF8ZUX85xNPLN+R3gPgVen1I+n4fdP2QWIioBXE\nHNcLKuL7UIprFXAhY4MyPgt8LsV9A/Aa4B3Ak8DP0vFvrPH3SlJDFgLnTrBvBrBDWt4ZuDctDwAb\niW8GAN8GTk7LtxMT4AD8DWNj+6wA/ndaPga4Oi2XGEv851TFUr5olJPvjHSecnKtlfi3IRJtuQS0\nfXrdnxNjtQPsA/wCmEUk/nuJi8POwNOMffs5lxiFtPxeF6XlIxgb/2kJMcAawJFEcoZI/D8mBrZ7\nNfBEiuNNxJACM9JxFwAfTsubgeMqPou/SsuXAH+Aet7MTgegwhits28rYkrAI4ik9DqiFQrROr49\nLa8kLgY7AjsRE+AAXAr8/4rzLUu/b03HV3s3cGLF+m/S7xOBjxL/X7yWSJ4TDby3DzFC48q0Xh46\n+jBgcVpeQyT+vYm/fwUxouNz6T3LF6I7GLu4QYz9RPr7yn/rYYwl5RVEkt8hnXeIuEA+SQwbvGv6\nGw8i7qVADND2WFp+Mb2GFP9/rXhvhxsuABO/2uUu4H0T7DuZaAUfCLxEJPtySWdDxXEvVWyvVJ2s\nyq95iYn/G69+zRuI1vrBRGv8kgneK4uJkmfl37K5Yn0z9f9fLF80JzrvixXLlX/zpcCnahy/sSqO\nyveud4FWj3AiFrXLD4iSx0crtpVv7u5ItFRfIsoYr69znj5gHTHVZbl+/2GiNp7V1cDHK9Znpxie\nS+fehSgT1bOG+FZwcFrfgSirXMtYOWpvYA9gNfVb0n1Vy+VvI4cT3wzWVZ23BPwKeGaC844C1xAX\n2t9K216VYqnnGeJzUI8z8audTgCOIrpz3gl8niiXXEYk0NuJJH5PxWuqW6Dl9VOIWv4IcQH57ATv\nOVpj+XNEPf8O4gZnKZ1nFZGkLyPq5vW8SCToJekcy4kL2wXE/1e3E11XTyFa2KMTxELVvlFgPVGm\nugD447R9kCjdjBATgpxS47WV7gE+TdxMHkm/d53kvZcSwzavxJu7ktQ2Rek6qg6yxS9JkiRJkiRJ\nkiRJkiRJkiRJUvv8J/k6BBvRxpqVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcbe3b75290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot canonical correlations (cca.cancorrs)\n",
    "plt.plot(np.arange(nComponents)+1, cca.cancorrs, 'ko')\n",
    "plt.xlim(0.5, 0.5+nComponents)\n",
    "plt.xticks(np.arange(nComponents)+1)\n",
    "plt.xlabel('Canonical component')\n",
    "plt.ylabel('Canonical correlation')\n",
    "plt.title('Canonical correlations')\n",
    "print '''The canonical correlations are:\\n\n",
    "Component 1: %.02f\\n\n",
    "Component 2: %.02f\\n\n",
    "Component 3: %.02f\\n\n",
    "Component 4: %.02f\\n\n",
    "''' % tuple(cca.cancorrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, two canonical components capture correlational relationships between the datasets, and the first two canonical correlations are both 0.95. The third and the fourth canonical components do not capture meaningful relationships between the datasets. The third and the fourth canonical correlations are very low, 0.10 and 0.00."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cross-dataset predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The prediction accuracy for the first dataset is:\n",
      "\n",
      "Dimension 1: 0.92\n",
      "\n",
      "Dimension 2: 0.91\n",
      "\n",
      "Dimension 3: 0.91\n",
      "\n",
      "Dimension 4: 0.91\n",
      "\n",
      "\n",
      "The prediction accuracy for the second dataset is:\n",
      "\n",
      "Dimension 1: 0.91\n",
      "\n",
      "Dimension 2: 0.93\n",
      "\n",
      "Dimension 3: 0.90\n",
      "\n",
      "Dimension 4: 0.91\n",
      "\n",
      "Dimension 5: 0.90\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEZCAYAAACQK04eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FPX9x/FXEsJlDsKlHIFoQC1ohapIRWXRFqgcWtty\nStXaIv2JYqsWgV9L9NdflZ9KEamIB0pF61VbkVQQwQheQaRAiYgFxXIEPICAilyZ3x/fSbLZ7G5m\nszt7ZN7Px2MfmZmdnfnsKp+Z+cx3vl8QERERERERERERERERERERERERERERidjjwP/Y0xcCHzRw\nO3OB/45FQCIiAtuAr4GDwG7gMeCEGG37MeCOCD9zNbAqRvsXSVnpiQ5AGjULGApkA98BziH42XWT\nBm4/rYGf84KMRAcgyUuJX+JlF7AE6GnPVwL/Bfwb2GwvGwqsA/YBbwJn+n2+N7AWOAA8DTT3e88H\nbPebzwdeAD4FPgfuB04HHgS+i7kC2Wuv+zg1JSOAX9gxfQG8CHTwe68SuA740I5xTpjv2wd4215v\nlx1Dpt/7PYFl9n52A1Ps5RnAVGCL/V3XAJ2AAnv//v9mS4Br7emrMb/ZTPs7TwdOAVbY858BC4Fc\nv88H/k6z7Rj3Amf4rdce+ApoE+b7iogA8DFwiT2dD2wEbrfnK4GlQCugGSax7wHOxZzJ/9T+fCbQ\nFPgEmIRJjD8CjlBT6vFRk/gzgPXAvUALe9vn2+9dRd1Sj3/J6GJMguxl73M28LrfupXAIiDH/j6f\nAoNCfPfvYJJ/OtAVeN+OH8wVUDnwK3s/Wfa6ALcCG4Du9vyZQGuCJ/7XgJ/Z01cDR4Hr7XWaA4WY\n3z8TaGt/lz/a64f7nf4E3OW3n0mYg6CISL22Yc6u99nTczAJBkwS8/mtO5e6NfsPgIvs186A994k\neOL/LiYhB7uavZrwif9Raie8EzAHmC5+MZ/v9/4zwOQg+wnmJszZNcBo4L0Q630ADAuyvID6E/8n\n9cRwOeaqCcL/TucFbGsN8ON6ti0ppKG1VREnLOAyTLkhGP/yTFfMWf4NfssyMaWWNOom/lBJLt9+\nrzLSYO19rfGb/wpTiukE/Mdettvv/a8xZ+vBnIopu5wNtMT8W6vadj7wUYjP5QNbIw3ctj1g/kTg\nPuACzFVGOjUlrnC/UylwCHNA3Y25cljUwJgkCanGL4lk+U3/B/hfIM/vlYU5qy7HJF9/XUNsczvm\nDD3YzU0ryDJ/uzBn1lVOwNS1Aw86TszFlHe6Yerq06j59/YfTP09mO32ZwJ9Zf9t6bfspIB1Ar/f\nH4DjmHp9LjDOL4ZwvxPAAuBK+zPPYa58pJFQ4pdk8TAwAVPrTsMk3SGY5P8WcAy4EXMVcAXmXkAw\nqzEHirswSbI5NeWZPUBnat9kTaOmddBfgGuAszAlqT8A71Bzth8oXKuiLEyZ62vMjeVf+r1XjLm6\nmGTvJ5uaGv8jmJvN3eztfxtT4/8McwAah0nWP8OciYeThTlgHMAcOG/1ey/c7wTmRvAVwFjgz/Xs\nR1KMEr8kSuDZ6XuYFjVzMOWIf2NKP2BuWl6BqWN/AYwA/hpie8cxNfJumIS93V4fYDlQhilffOr3\nOcvv/d/a294FnAyMChOzFWRZlVuAMZik+xCmJVLVugeB79txlmNaCfns92YCzwKvABWYA2JVC6Zf\nYJL350APzH2OcLHcjrnJXAG8ZH8vJ78T9vxaTCnojRDfUSSo+ZizrH+FWWc25h/5ekzLDhFJDo8S\n+UNyIlyISeahEv+lwD/s6fMwl9UikngFmNZYoe6liIRVQOjE/yAw0m/+A0xLBBFJnP/BlKOm1Lei\nSCgFhE78L1H7htKrmOZvIiLikmS4uRvYMqK+JnciIhKFRD/AtRPzIEmVzgRpM33WWWdZ69evj1tQ\nIiKNxHpMFyS1JPqMfxE1Tfb6AvsxrYBqWb9+PZZlxfQ1ffr0mG/TjZfi9FaMqRSnXsn/wjyTUofb\nZ/x/AfpjOojajukxsOrhmXmYFj2XYnoi/Arz8IyIiLjI7cQ/2sE6E12OQURE/CS61JMwPp8v0SE4\nojhjJxVihNSJU1JXqoxgZNn1KhERcSgtLQ2C5PlEt+oRkRTTunVr9u3bl+gwxE9eXh579+6tf0Wb\nzvhFJCJpaWno32NyCfXfJNQZv2dr/CIiXqXELyLiMUr8IiIeo8QvIuIxSvwi0mgUFBTQsmVLcnJy\nyMvLo1+/fsybN8/xzeht27aRnp5OZWWwMehjx8l+Nm7cyKBBg2jXrh3p6bFN1WrOmaReW1pM8YL7\nST/+DZUZzRly1Q0MGDQk0WGJhFT86ivMfe5JjlBJU9L55U/GMuR7A+O6jbS0NBYvXszFF1/MwYMH\nKSkpYdKkSZSWljJ//nzH24lXq6Vw+2natCmjRo3i+uuv5/LLL4/pfj13xl/86isMve4qBl43jqHX\nXUXxq68kOqQ6XltazN/unsT4psv4eYtVjG+6jL/dPYnXlhYnOjSRoIpffYXb5s9mXd8OvN+3E+v6\nduC2+bMj+vcVi234y87OZtiwYTzzzDMsWLCAsrIys5/iYnr37k1ubi5dunTh9ttvr/7MRRddBECr\nVq3Izs6mtLSUrVu3cvHFF9O2bVvatWvHlVdeSUVFRfVnZsyYQefOncnJyeH0009nxYoVgEnqd911\nF926daNt27aMHDmy+vmHYPsJdOqpp3LNNdfQo0ePBn3/cDyV+GP9P5Zbihfcz8TOH9daNrHzxxQ/\nMSdBEYmEN/e5J9n3/TNqLdv3/TN48Pmn4rqNYM4991w6d+7MG2+YMeOzsrJYuHAhFRUVFBcXM3fu\nXF588UUAVq1aBUBFRQUHDx7kvPPOA2DatGmUl5ezadMmtm/fTlFREQCbN2/mT3/6E2vWrOHAgQO8\n8sorFBQUADB79mwWLVrEypUrKS8vJy8vj+uvvz7sfuLFU6WecP9jRXpJ6qb0498EX370UJwjqV8s\nLu/dprKZ+44QvFZ92Doe122E0rFjx+onW/v371+9/Mwzz2TUqFG8/vrrXHbZZUFLL4WFhRQWFgLQ\ntm1bfvWrX3HHHWYM+oyMDA4fPkxZWRlt2rShS5cu1Z+bN28ec+bMoWPHjgBMnz6drl27snDhwoQ/\nAOepxH+EStLLtnLimvdobh3jm7Qm7DnnbA5bzRMdWi2VGcHjqcxsEedIwqu6gvI/mN42fzZA0iT/\n15YWs/CO8Uzuvrt62Yw7/gU8pOQfQ01DFA+apWXEdRuh7Ny5k9atWwNQWlrKbbfdRllZGUeOHOHw\n4cOMGDEi5Gf37NnDpEmTeOONNzh48CCVlZXV2+rWrRuzZs2iqKiIsrIyBg0axMyZM+nQoQPbtm3j\nhz/8Ya0bs02aNGHPnjpDjsSdp0o9R3btpNeq5TxZ8AmPnryTJws+odeq5Rwt35Xo0Gppf1Y/fvNe\n7SR/63vNaX/m+SE+kRhzn3uSio4t6LDgWU5+/Ck6LHiWio4tor40j6XH/zi9VtIHmNx9NwtmFSUm\noDBeW1rMLWMG85uRPm4ZMzil7un88idjyVu2sdayvFc2MuHHY+K6jWDeffdddu7cyQUXXADAmDFj\nuPzyy9mxYwf79+9nwoQJ1a1r7C4Oapk6dSoZGRls3LiRiooKnnjiiVqtcUaPHs2qVav45JNPSEtL\nY/LkyQB06dKFJUuWsG/fvurX119/TYcOHYLuJ548dcbf9egBpn2rotaye75VwR/2HEhQRMGt/GgL\nGwb9gLH+VyaDzqby4638JtHB+dm/azu9Pl3HPX6/6S2r9rO3fZ2R3hLm0z074KS6y/fs3h7/YMKo\nuqHvf29nzt1bAFLiyqTqCu/B55/isHWcZmkZTLj2xoiu/GKxDahpKXPgwAFWrlzJTTfdxLhx4+jZ\nsycAX375JXl5eTRt2pTVq1fz1FNPMWjQIIDqppNbt26le/fu1evn5uaSk5PDzp07ufvuu6v39eGH\nH7Jjxw769etHs2bNaN68efX+J0yYwNSpU1mwYAFdunThs88+4+2332b48OFB9xPMN998w5EjRwA4\nfPgwAM2aNYvo9wjGU4n/xNysoMvb55wQ50jCO0IllT0LKe9ZWGv54XfqDEecUE23b+aeXnUPpBPW\nb05QRHUdsoJf1B5KsovdUDf0H3piTkokfjCJO9oSXyy2MWzYMJo0aUJ6ejo9e/bk5ptvZsKECdXv\nP/DAA9x8881MnDiR/v37M3LkSPbv3w9Ay5YtmTZtGv369ePYsWMsWbKE6dOn89Of/pTc3Fy6d+/O\nlVdeyaxZswCTjKdMmcKmTZvIzMykX79+PPTQQwBMmjQJy7IYOHAgu3bton379owaNYrhw4fX2s/R\no0dZunQpffr0qfU9tm3bximnnAKYK5EWLVpQUFDARx99FNXvAx7rnfOWMYMZ33RZneUPHRvIPQtf\njnr7sTL0uqtY17dDneW9S3fz0oOPxz+gEK4deBaTO26ss3zGrjN49JX1CYioLt+wi8kJuCq5+f1c\nvjypN68tWp7AyGr7zUgfP2+xqs7yRw5dyP89UxL/gMJQ75zJR71zhjHkqhuYs+PkWsvu33EyQ8Yl\n1+iPbtU6Yy2vbd2DE0Be+45xjiS0WyfdxpbsQsZu68q1H3di7LaubM0u5JYbJyc6tFpS5Ya+NA6e\nKvVUXTI/9MQc0o8eojKzBVfcOjHpLqVjVet025CrbmDO3VtqlSju33EyV9yaPAfSwN+yVVoGE348\nRr+leJqnSj0Se68tLabY70A6ZFzyHUhTRar8lir1JJ9ISz1K/CISESX+5KMav4iIhKXELyLiMUr8\nIiIeo8QvIuIxSvwiIh6jxC8ijUZjGnpxwYIFnHPOOeTm5pKfn8/kyZM5fjz6LqrBYw9wiYh7YjHu\nQbTbaExDLx46dIj77ruP8847j08//ZThw4dzzz33VPf+6QWWiCSHYP8eVyxZbN1wSaG1+ar06tcN\nlxRaK5YsdrzdWGyjoKDAWr58ea1lq1evttLT062NGzdalmVZixcvtnr16mXl5ORY+fn5VlFRUfW6\n+fn5VlpampWVlWVlZWVZ77zzjrVlyxZrwIABVps2bay2bdtaY8eOtfbv31/9mbvuusvq1KmTlZ2d\nbZ122mnV+6+srLTuvPNOq7Cw0GrTpo01YsQIa+/evSH3U5+ZM2daw4YNC/peqBwJBD2yqNQjIlGL\nxXChbg052liGXnz99dc544wz6l3PCSV+EYlaLIYLdXPI0cChF6v65vcfehGCl14KCwu55JJLyMzM\nrB56sWp9/6EXjx49SpcuXaq7Up43bx6///3v6dixI5mZmUyfPp3nn3+eysrKiEtJ8+fPZ+3atdxy\nyy0N/g38KfGLSNRi0buomz2UBg69OGDAANq3b0+rVq2YN28eX3zxRcjP7tmzh1GjRtG5c2dyc3MZ\nN25c9fr+Qy+eeOKJjB49mvLycoDqoRfz8vLIy8ujR48eDRp68e9//ztTp07l5Zdfrv4O0VLiF5Go\nxaLLc7e6TU/loReXLFnC+PHjWbx4cfVVSiyoVY+IRC0WXZ7Hqtt0q5EMvbhixQrGjh3Liy++yDnn\nnBPRb9BY1HvHW0TiI5n/PRYUFFgtWrSwsrOzrdzcXOv888+3HnjgAauysrJ6neeff97q2rWrlZ2d\nbQ0dOtS64YYbrHHjxlW//7vf/c5q166dlZeXZ5WWllplZWXW2WefbWVlZVm9e/e27r33Xis/P9+y\nLMvasGGD1adPHys7O9tq3bq1NWzYMKu8vNyyLNOqZ+bMmdZpp51mZWdnW4WFhda0adPq7KdVq1ZW\naWlpne8yYMAAKzMzs7rlT1ZWlnXppZcG/d6h/psQolWP290yDwZmARnAI8CMgPfbAgsxw2E3Ae4B\nHg+yHfs7iEiiqVvm5JNM/fFnAJuB7wE7gXeB0cAmv3WKgGbAFMxBYDNwInAsYFtK/CJJQok/+SRT\nf/x9gC3ANuAo8DRwWcA65UCOPZ0DfEHdpC8iIjHk5s3dTsB2v/kdQOBTCg8DK4BdQDYwwsV4REQE\ndxO/k2vBqcA6wAcUAsuAs4CDgStWPSkH4PP58Pl8MQhRRKTxKCkpoaSkpN713Kzx98XU8Afb81OA\nSmrf4P0H8L/Am/b8cmAysCZgW6rxiyQJ1fiTTzLV+NcA3YECoCkwElgUsM4HmJu/YG7qngZ85GJM\nIiKe52ap5xgwEViKaeHzKKZFz3X2+/OAPwCPAesxB6HfAHtdjElEopSXl+f4yVOJj7y8vIjWT5X/\neir1iIhEKBGlHhERSUJK/CIiHqPELyLiMUr8IiIeo8QvIuIxSvwiIh6jxC8i4jFK/CIiHqPELyLi\nMUr8IiIeo8QvIuIxTjppaw/8AtPLZtX6FvAzl2ISEREXOUn8LwIrMYOkVNrL1GOaiEiKctI75zqg\nl9uB1EO9c4qIRCia3jkXA0NiHZCIiCSGkzP+L4GWwBHgqL3MAnLcCioInfGLiEQo1Bm/BmIREWmk\nQiV+p0MvXgZchDnTfx14KWaRiYhIXDk5478LOBd40l5/FGYg9SkuxhVIZ/wiIhGKptTzL0yrnuP2\nfAampc+ZsQrOASV+EZEIRdOqxwJa+c23Qu34RURSlpMa/53AWqDEnu8P3OZWQCIi4i6nrXo6Yur8\nFrAa2O1aRMGp1CMiEqGG1Pi/BWwCzsYk/Kp1qzLw2hjGVx8lfhGRCDUk8T+M6ZythOA1/QGxCMwh\nJX4RkQhF06qnOfCNg2VuUuIXEYlQNK163nK4TEREUkC4Vj0dMDd1WwLfwRw1qvroael+aCIi4oZw\niX8gcDXQCbjXb/lBYKqLMYmIiIuc1Ph/DDzvdiD1UI1fRCRC0fbOORTogbmpW+WO6MNyTIlfRCRC\n0dzcnQeMAG60NzAC6BrL4EREJH6cdtJ2JrAB+DaQBSwBLnAxrkA64xcRiVA0Z/yH7L9fY270HgNO\nillkIiISV046aVsM5AF3A+/Zyx52LSIREXFVpEMvNrdf+x2uPxiYhenD/xFgRpB1fMAfgUzgc3s+\nkEo9IiIRakirnh9R00dPGnX763mhnn1mAJuB7wE7gXeB0ZiO36q0At4EBgE7gLaY5B9IiV9EJEIN\nGXN3GOEHXKkv8fcBtgDb7PmnMWP3+if+McBfMUkfgid9ERGJoXCJ/+oot90J2O43vwM4L2Cd7pgS\nz2tANnAf8ESU+xURkTCctOo5CXgU04QTzINc1zr4nJPaTCamH6BLMeWe32IOBiIi4hInrXoeBx4D\nptnz/waexRwMwtkJ5PvN51NT0qmyHVPeOWS/VgJn2fuopaioqHra5/Ph8/kchC4i4h0lJSWUlJTU\nu56TVj1rgHOAfwK97WXrgF71fK4J5ubuJcAuzJCNgTd3TwfmYM72mwGlwEjg/YBt6eauiEiEGnJz\nt8qXQBu/+b5AhYPPHQMmAksxLXwexST96+z35wEfYEpIG4BKzPMBgUlfRERiyMkZ/9nA/UBPoAxo\nh+mxc72LcQXSGb+ISIQaesafAVxkv063N7AZOBLj+EREJE6cnPG/C5zrdiD10Bm/iEiEoumPv6o7\nhWeAr6h5indtDOOrjxK/iEiEokn8JQRvkz8gupAiosQvIhKhaGr8i4CZLsQkIiIJUN+Tu8cxbe9F\nRKSRUI1fRKSRUo1fRMRjokn8yUCJX0QkQtGMudsKU+55z37dC+TGMjgREYkfJ4l/PnAA+AkwAjiI\n6a1TRERSkJNSz3pMV8n1LXOTSj0iIhGKptRzCLjQb/4C4OvYhCUiIvHm5Iy/F/Bnaur6+4CrUO+c\nIiJJLRateqoSv5O++GNNiV9EJELRlHruxLTsqbBfecDvYxmciIjEj5PE/wNgv9/8PmCIO+GIiIjb\nnCT+dKC533wLoKk74YiIiNucjLn7JLAc054/DbgGc7NXRERSkNObuz8ALrGnl2EGUI8n3dwVEYmQ\n+uoREfGYaFr1iIhII6LELyLiMUr8IiIe46RVzwXAdKDAb30LOMWlmERExEVObu5uBm7CDLV43G/5\n565EFJxu7oqIRCjUzV0nZ/z7gZdjHZCIiCSGkzP+u4AM4AXgsN9yDbYuIpLENNi6iIjH6AEuERGP\n0WDrIiICaLB1ERHP0WDrIiKNlAZbFxERQIOti4g0WrFo1ZNj/z0Qi4AipMQvIhKhhjy5Ow54AriZ\n2u340+z5mQ72OxiYhXkA7BFgRoj1zgXextw8fsHBdkVEpIHCJf6W9t9sgj/AVZ8MYA7wPWAn8C6w\nCNgUZL0ZwBJS57kCEZGUFS7xz7P/vgq8EfDeBQ623QfYAmyz558GLqNu4r8BeB5z1i8iIi5z0qrn\n/iDLZjv4XCdgu9/8DntZ4DqXAXPteRXyRURcFu6M/7vA+UA74NfUlGGyMeWZ+jhJ4rOA2+x101Cp\nR0TEdeESf1Nqkny23/IDwI8dbHsnkO83n4856/d3NqYEBNAW+AFwFHMvoJaioqLqaZ/Ph8/ncxCC\niIh3lJSUUFJSUu96Ts6wuwKfNCCGJphBXC4BdgGrgdHUrfFXeQx4ieCtetScU0QkQtE8ufsIpqO2\nKq2BpQ4+dwyYaK/7PvAMJulfZ79ERCQBnJzxr8M8vVvfMjfpjF9EJELRnPEfx5R7qhQAlTGJSkRE\n4s7JmLvTgFXASnv+ImC8axGJiIirnDafbAf0xTS7fAf43LWIglOpR0QkQg3ppO1bmJuxZ1PTzh5q\n2udrsHURkSTWkMT/MPALNNi6iEhK0mDrIiIe05BumX9E+G4X1H2yiEgKCpf4h2ESf3tMnz0r7OUD\ngLdQ4hcRSUnhEv/V9t9lQA+g3J7vACxwMSYREXGRkwe48oHdfvN7gC7uhCMiIm5z8gDXq5j+dp7C\n3CQYibkKEBGRFOSkVU8a8EPgQnt+JfA31yIKTq16REQi1JBWPVUszMNaBzFn+i0x/fMfjGF8IiIS\nJ05q/OOB54AH7fnOwN9di0hERFzlJPFfjxlc/YA9/yGmiaeIiKQgJ4n/sP2q0gQNii4ikrKcJP7X\nMV0ztwS+jyn7vORmUCIi4h4nrXrSgZ8DA+35pZjhGON51q9WPSIiEWpoJ21NgI3A6S7EFAklfhGR\nCDV06MVjwGZqD70oIiIpzEk7/tZAGbAa+MpeZgHD3QpKRETc4yTx/7f91/9yQXUXEZEUFS7xtwAm\nAN2ADcB84Gg8ghIREfeEu7n7LHAEWAVcCmwDJsUhpmB0c1dEJEINadXzL+BMe7oJ8C7QO+aROaPE\nLyISoYa06jkWYlpERFJYuDP+48DXfvMtgEP2tAXkuBVUEDrjFxGJUEO6Zc5wLRoREUkYJ331iIhI\nI6LELyLiMUr8IiIeo8QvIuIxSvwiIh6jxC8i4jFK/CIiHqPELyLiMfFI/IOBD4B/A5ODvD8WWI/p\nAfRN4NtxiElExLOcjLkbjQzMCF7fA3ZiOnobDWzyW+e7wPtABeYgUQT0DdiOumwQEYlQQ4dejFYf\nYAumS+ejwNPAZQHrvI1J+gClQGeXYxIR8TS3E38nYLvf/A57WSjXAv9wNSIREY9zMvRiNCKpzwwA\nfgb0C/ZmUVFR9bTP58Pn80UTl4hIo1NSUkJJSUm967ld4++LqdkPtuenAJXAjID1vg28YK+3Jch2\nVOMXEYlQomr8a4DuQAHQFBgJLApYpwsm6V9J8KQvIiIx5Hap5xgwEViKaeHzKKZFz3X2+/OA3wF5\nwFx72VHMTWEREXGB26WeWFGpR0QkQokq9YiISJJR4hcR8RglfhERj1HiFxHxGCV+ERGPUeIXEfEY\nJX4REY9R4hcR8RglfhERj1HiFxHxGCV+ERGPUeIXEfEYJX4REY9R4hcR8RglfhERj1HiFxHxGCV+\nERGPUeIXEfEYJX4REY9R4hcR8RglfhERj1HiFxHxGCV+ERGPUeIXEfEYJX4REY9R4hcR8RglfhER\nj1HiFxHxGCV+ERGPUeIXEfEYJX4REY9R4hcR8RglfhERj1HiFxHxGCV+ERGPcTvxDwY+AP4NTA6x\nzmz7/fVAb5fjERHxPDcTfwYwB5P8ewCjgW8FrHMp0A3oDowH5roYTy0lJSXx2lVUFGfspEKMkDpx\nSupyM/H3AbYA24CjwNPAZQHrDAcW2NOlQCvgRBdjqpYq/7gUZ+ykQoyQOnFK6nIz8XcCtvvN77CX\n1bdOZxdjEhHxPDcTv+VwvbQGfk5ERJJMX2CJ3/wU6t7gfRAY5Tf/AcFLPeswBwS99NJLL72cv9YR\nZ02ArUAB0NQOINjN3X/Y032Bd+IVnIiIuOMHwGbMTd4p9rLr7FeVOfb764HvxDU6ERERERFpfOYD\ne4B/JTqQMPKB14AyYCNwY2LDCak5phnuOuB94M7EhlOvDOCfwEuJDiSMbcAGTJyrExuKSONxIeYJ\n4WRO/CcBvezpLEy5LPD+SLJoaf9tgrlHc0ECY6nPr4EngUWJDiSMj4HWiQ5CGjcv9tWzCtiX6CDq\nsZuau/FfApuAjokLJ6yv7b9NMWfUexMYSzidMY0JHqFuE+Jkk+zxSYrzYuJPNQWYK5TSBMcRSjrm\nILUHU556P7HhhPRH4FagMtGB1MMCXgXWAL9IcCwijUoByV3qqZKFSQCXJzoQB3IxpR5fguMIZijw\nJ3vaR3LX+DvYf9thDqgXJjAWaaR0xp+8MoG/AguBvyc4FicqgGLgnEQHEsT5mH6hPgb+AlwM/Dmh\nEYVWbv/9DPgbps8rEYmBApL7jD8Nk5j+mOhA6tEW07EeQAtgJXBJ4sJxpD/Je8bfEsi2p08A3gQG\nJi4ckcbjL8Au4DCmg7hrEhtOUBdgatHrMM36/onp3jrZnAmsxcS5AVNDT3b9Sd5WPSdjfst1mGa8\nU8KvLiIiIiIiIiIiIiIiIiIiIiIiIiIiIiISleOYZxY2Ytqy/5r6OyzrCox2IZabMA+l1cdHzQNh\nw6g7vGg8dASeS8B+RUSidtBvuh2wDCiq5zM+3HkS92OgjYP13Nq/iIgnHAyYPxn43J4uwHQJ8Z79\n+q69/B1gP+ZKYRLmCiDYeh3s5f/EdNdRNW7AQOAte91nMV0l3Ih5snsDsDxInIMxXWW/B9xHTeK/\nGrjfnn6uQDO0AAAChUlEQVQceAB4GzP2tA9YgOm99DG/bQXbP5gBWYrs5RuA0+zl/al5mnutvX4B\nNV2QNLe3v8F+3+cX2wvAy8CHwIwg30tEJO4CEz+YsRPaYcouzexl3YF37enAvndCrXczMNWeTsf0\nfNoWeJ2aks5k4Lf2dKgBUZoD/wEK7flnqOkC4mpqJ/6n7OnhwAGgJ6Z0tQY4y8H+r7enfwk8bE8v\nouZg1hIzBkIBNYn/ZsyYA2AOFp9gfo+rMQegbHt+G9ApyPcTD2iS6ABEHGoKzMEkzOOYpA517wGE\nWm81ZtjNTExvp+sxZ8M9MGfcVZ99i/BOxyTlrfb8QmB8kPUsag5IGzGD65TZ82WYZJ1fz/5fsP+u\nBa6wp9/EdN73pP3+zoD99gNm29ObMYn/VDue5dQcXN+3Ywj8vHiAEr8ks1MwyfszTNmjHBiHOcv9\nJsRnfhVivVWYvu2HYs7GZ2KuJpYBYyKIyQqYD3fz+Yj9txJTOsJvvgnmu4Xbf9VnjlPzb3UGsBgY\ngjkIDArYdriY/Nc7jvl9xIPUH78kq3bAg9SUTnIwZ80AP6UmaR2kpivjcOt1wRxAHrFfvTH3B/pR\nU7Y5gZorhIP2tgJtxpwpn2LPN7RFkVXP/kMpxFwx/B+mjHVawPurgLH29KmY7/0BwQ8GGuLRo5T4\nJZm0oKY55zJgCXCH/d4DwFWYZp6nYcYiBlOyOW4vnxRmvQH2srXACMxN2c8xte+/2Nt5i5pE+pC9\n/8Cbu99gSjvFmBuve6i5CrCofUUQarpKuP0Hqvr8JEw9fz3miuLlgPcfwPy73gA8jfktjgaJLVRM\nIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIi0lD/D2Y3y9LBzWikAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcc18423c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot correlations between actual test data and predictions \n",
    "# obtained by projecting the other test dataset via the CCA mapping for each dimension.\n",
    "nTicks = max(testcorrs[0].shape[0], testcorrs[1].shape[0])\n",
    "bmap1 = qualitative.Dark2[3]\n",
    "plt.plot(np.arange(testcorrs[0].shape[0])+1, testcorrs[0], 'o', color = bmap1.mpl_colors[0])\n",
    "plt.plot(np.arange(testcorrs[1].shape[0])+1, testcorrs[1], 'o', color = bmap1.mpl_colors[1])\n",
    "plt.xlim(0.5, 0.5 + nTicks + 3)\n",
    "plt.ylim(0.0, 1.0)\n",
    "plt.xticks(np.arange(nTicks)+1)\n",
    "plt.xlabel('Dataset dimension')\n",
    "plt.ylabel('Prediction correlation')\n",
    "plt.title('Prediction accuracy')\n",
    "plt.legend(['Dataset 1', 'Dataset 2'])\n",
    "print '''The prediction accuracy for the first dataset is:\\n\n",
    "Dimension 1: %.02f\\n\n",
    "Dimension 2: %.02f\\n\n",
    "Dimension 3: %.02f\\n\n",
    "Dimension 4: %.02f\\n\n",
    "'''% tuple(testcorrs[0])\n",
    "print '''The prediction accuracy for the second dataset is:\\n\n",
    "Dimension 1: %.02f\\n\n",
    "Dimension 2: %.02f\\n\n",
    "Dimension 3: %.02f\\n\n",
    "Dimension 4: %.02f\\n\n",
    "Dimension 5: %.02f\\n\n",
    "'''% tuple(testcorrs[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both datasets are predicted with close to perfect accuracy, with correlations ranging from 0.90 to 0.93 for each dimension."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pyrcca with hyperparameter grid search using cross-validation\n",
    "\n",
    "Instead of predefining hyperparameters, grid search with cross-validation can be used to find the optimal hyperparameters. Since there are two latent variables, two canonical components are expected to be optimal. Additionally, since the latent variables are not necessarily orthogonal, the optimal regularization coefficient should be greater than zero. The exact value of the optimal regularization coefficient will depend on the instantiation, since the data are randomly generated.\n",
    "All possible numbers of components are tested: 1, 2, 3, 4. Four possible regularization coefficients are tested: 0, 1e2, 1e4, 1e6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training CV CCA, regularization = 0.0000, 1 components, fold #1\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #2\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #3\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #4\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #5\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #6\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #7\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #8\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #9\n",
      "Training CV CCA, regularization = 0.0000, 1 components, fold #10\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #1\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #2\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #3\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #4\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #5\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #6\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #7\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #8\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #9\n",
      "Training CV CCA, regularization = 0.0000, 2 components, fold #10\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #1\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #2\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #3\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #4\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #5\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #6\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #7\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #8\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #9\n",
      "Training CV CCA, regularization = 0.0000, 3 components, fold #10\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #1\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #2\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #3\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #4\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #5\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #6\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #7\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #8\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #9\n",
      "Training CV CCA, regularization = 0.0000, 4 components, fold #10\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #1\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #2\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #3\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #4\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #5\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #6\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #7\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #8\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #9\n",
      "Training CV CCA, regularization = 100.0000, 1 components, fold #10\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #1\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #2\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #3\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #4\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #5\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #6\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #7\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #8\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #9\n",
      "Training CV CCA, regularization = 100.0000, 2 components, fold #10\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #1\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #2\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #3\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #4\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #5\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #6\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #7\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #8\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #9\n",
      "Training CV CCA, regularization = 100.0000, 3 components, fold #10\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #1\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #2\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #3\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #4\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #5\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #6\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #7\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #8\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #9\n",
      "Training CV CCA, regularization = 100.0000, 4 components, fold #10\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #1\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #2\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #3\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #4\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #5\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #6\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #7\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #8\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #9\n",
      "Training CV CCA, regularization = 10000.0000, 1 components, fold #10\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #1\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #2\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #3\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #4\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #5\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #6\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #7\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #8\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #9\n",
      "Training CV CCA, regularization = 10000.0000, 2 components, fold #10\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #1\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #2\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #3\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #4\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #5\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #6\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #7\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #8\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #9\n",
      "Training CV CCA, regularization = 10000.0000, 3 components, fold #10\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #1\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #2\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #3\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #4\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #5\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #6\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #7\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #8\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #9\n",
      "Training CV CCA, regularization = 10000.0000, 4 components, fold #10\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #1\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #2\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #3\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #4\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #5\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #6\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #7\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #8\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #9\n",
      "Training CV CCA, regularization = 1000000.0000, 1 components, fold #10\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #1\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #2\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #3\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #4\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #5\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #6\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #7\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #8\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #9\n",
      "Training CV CCA, regularization = 1000000.0000, 2 components, fold #10\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #1\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #2\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #3\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #4\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #5\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #6\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #7\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #8\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #9\n",
      "Training CV CCA, regularization = 1000000.0000, 3 components, fold #10\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #1\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #2\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #3\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #4\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #5\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #6\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #7\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #8\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #9\n",
      "Training CV CCA, regularization = 1000000.0000, 4 components, fold #10\n"
     ]
    }
   ],
   "source": [
    "# Initialize a cca object as an instantiation of the CCACrossValidate class.\n",
    "ccaCV = rcca.CCACrossValidate(kernelcca=False, numCCs = [1, 2, 3, 4], regs = [0., 1e2, 1e4, 1e6])\n",
    "\n",
    "# Use the train() and validate() methods to run the analysis and perform cross-dataset prediction.\n",
    "ccaCV.train([train1, train2])\n",
    "testcorrsCV = ccaCV.validate([test1, test2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal number of components: 2\n",
      "Optimal regularization coefficient: 10000\n"
     ]
    }
   ],
   "source": [
    "print 'Optimal number of components: %d\\nOptimal regularization coefficient: %d' % (ccaCV.best_numCC, ccaCV.best_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the optimal number of components is two, corresponding to the number of latent variables. The optimal regularization coefficient is greater than zero: 1e4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Canonical correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The canonical correlations are:\n",
      "\n",
      "Component 1: 0.95\n",
      "\n",
      "Component 2: 0.95\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEZCAYAAACQK04eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGH5JREFUeJzt3XmUXGWZgPGnCZDIGhUFhWArmysgIKDgsRBGAswQURQR\nFdAjzIwk6OiIuPa4jcyMiImDoCKiosGRqDBRGcQUsktYmjURosiOgkDYEhLS88f7FX27Ul19q7tu\nVXff53dOnb57vb2996v3fve7IEmSJEmSJEmSJEmSJEmSJEltdSLw7TEeoxdYA6wz5mjaqwLcNYb9\nvwl8uj2hSCqrdwOLgceAe4FfAnt1NaL26GXiJ/6jgEsKi0QT1nj7o9bE8i/A14AvAi8EZgD/DRzc\nzaAmkHUbLJvS8SgkKadNiVb+25tssztwBfAw8WlgHrBeZv0a4FjgD2mbb2TW9RAliTuAB4CzgE3S\nut607/uAPwN/BT6Z2bcP+EFmfm/g8vQedwJHpuUHAdcBj6bln8vsU3uP4RpHM4AFwF+AB9P3Rtp+\npLjfn+K+OMVyGXByOs7ngfWB/0rb3E+UZ6alY1QY2uL/BHA7sBy4GXhrWv4K4ClgNfF7+lta/j3g\nC5n9PwjcBjwE/AJ4UWZds9/Ptin+R4if//wGPyNJk8xMYBXNPzXuQiT/dYCXALcAx2fWrwHOIxLj\nDCKJ7p/WvZ9ISL3AhsC5wPfTut607+nAVGBHYAWwQ1r/OQYT/0uIpHgY0Zp+HrBTWvcm4FVp+jVE\nkp1V9x6Nvr8pQD/wVeA5KYY3tBD399J+04hyzCrgQ+m9phGfon4OTAc2Sj+jL6djVBia+A8FtkjT\n7wQeBzZP80eydqnnTOLkAvBmImnvTJxs5hLJvKbR7+ctad2PiWsppH3fgKRJ7wjgvhb3+TDRSq5Z\nw9CEcQ7w8TR9EfCPmXXbA08TybE37fvizPqriMQHQ1v8JxLJN49TiJY3NE/8ryeSYKN1eeLuzaw/\nimjZ1/QQyftlde/3xzRdoXmN/zoGS21H0TzxnwF8JbNuwxTr1mm+2e/nLOLEu2WTWDROWePXaD0E\nbEbzv6Htgf8lThCPAl8Cnl+3zf2Z6SeJFi5EySGbEO8kauKbZ5YNt2/WDAaTZr09gEVEEn+EKGvU\nx9fIjBTbmgbr8sRdn7iz8y8ANgCuIcorDwO/In7WjbyPSPa1bV+d83toFOsTxO81m8zrf8Ybp+mP\nEyep3wM3AUfnfE+NAyZ+jdYVwErgkCbbfJMo72xLXBP4FPn/5u5laMt4a6Je/UCLcd4JbDPMuh8R\nJZWtiLLKaTnjuyvF0+hCbJ64B+r2yc4/SNTmXwk8N72mM3idIOslwLeIMtHz0rY3EQm50fuMFOuG\nxEnjnhH2g/h+jiFOEscCpzL0U4rGMRO/RutR4LNEL55ZRCt1PeAA4KS0zUbEhcUngZcD/zTCMXsY\nTFo/Bj5CJKaNiBr3fBq3spv5EbAf8A6i5f18Bmv8GxGt5KeJaxHvZuRkCVFWuo8ok2xA1OVrJZGx\nxr2GuAfhFKL1D5Fc39Jg2w1TvA8S/8tHEy3+mgeIk1r2gnr9z/ho4ucxNcV6JXGybKQnM/2OdGyI\nT0sDtP67UZeY+DUWJxNdOj9NlEvuBP4Z+Fla/zEimS4nWqbzGZpYG7V8a8u+S9Tpf0eUap4EZjfZ\nd7jj3AkcCHyUKGNcR1wMJsX6+RTfZ4gadv1xGlkD/APxSeZO4hNA7fpCq3EPNFh2AtFT50riBHsh\nUTarP8YtxAXmK4iSzKuBSzPbXUT09Lmf+P3Uv99FxPd9LtH6fynwrpyx7pbie4zoDTSH6Mkk8V2i\n1XFjk23mEr0g+oHXdiIoSVJx3kgk8+ES/4HEnZ4QF9qu7ERQkqRi9TJ84j+N6F9ds4ShvR8kSW3W\n7Rr/lgztynY3gxeMJEkF6Hbih6E9BSBfrwpJ0ig1GiSqk+4hboap2YoGfYh32mmngf7+/o4FJUmT\nRD8xJMcQ3U785wHHEd389iT6A691g05/fz8DA34QaJe+vj76+vq6HYa0Fv8226unp2enRsuLTvw/\nJgbC2oyo5X+OwZtJTid69BxI9Fl+Am/7lqTCFZ34D8+xzXEFxyBJyhgPF3fVYZVKpdshSA35t9kZ\n9T1qxqsBa/yS1Jqenh5okOdt8UtSyZj4S2ThwoXsv//+VCoV9t9/fxYuXNjtkCR1Qbe7c6pDFi5c\nyPHHH8+yZcueXVabPuigg7oVlgTE3+fcuXNZuXIlU6dOZc6cOf5dFsjEXxJz584dkvQhEv+8efP8\nB1NX2SjpPEs9JbFy5cqGy1esWNHhSKShmjVKVAwTf0lMnTq14fJp06Z1OBJpKBslnWfiL4k5c+aw\nzTZDHz27zTbbMHv27GH2kDrDRknnWeMviVqtdN68eaxYsYJp06Yxe/Zsa6jqujlz5rBs2bIh5R4b\nJcXyBi5JXbdw4UIbJQUY7gYuE78kTVLeuStJAkz8klQ6Jn5JKhkTvySVjIlfkkrGxC9JJWPil6SS\nMfFLUsmY+CWpZEz8klQyJn5JKhkTvySVjIlfkkrGxC9JJWPil6SSMfFLUsmY+CWpZEz8klQyJn5J\nKhkTvySVjIlfkkrGxC9JJWPil6SSMfFLUskUnfhnAkuA24ATGqzfDPg1cD1wE3BUwfFIUun1FHjs\nKcBSYD/gHuBq4HDg1sw2fcBU4ETiJLAU2BxYXXesgYGBgQJDlaTJp6enBxrk+SJb/LsDtwN3AKuA\n+cCsum3uAzZJ05sAD7F20pcktdG6BR57S+CuzPzdwB5123wb+C1wL7Ax8M4C45EkUWziz1Ob+SRR\n368A2wAXAjsBj9Vv2NfX9+x0pVKhUqm0IURJmjyq1SrVanXE7Yqs8e9J1PBnpvkTgTXASZltfgl8\nCbgszV9EXAReXHcsa/yS1KJu1PgXA9sBvcD6wGHAeXXbLCEu/kJc1N0B+GOBMUlS6RVZ6lkNHAdc\nQPTwOYPo0XNsWn868GXgTKCfOAl9HPhbgTFJUukVWeppJ0s9ktSibpR6JEnjkIlfkkrGxC9JJWPi\nl6SSMfFLUsmY+CWpZEz8klQyJn5JKhkTvySVjIlfkkrGxC9JJZMn8b+deGbucmKc/MfStCRpAsoz\nSNsy4O8Z+qzcTnOQNklq0VgGabuf7iZ9SVIb5RmPfzFwDvBz4Om0bABYUFRQkqTi5En8mwJPAW+p\nW27il6QJyAexSNIkNZYa/wzgZ8Bf0+tcYKt2BidJ6pw8if9M4iHpL06v89MySdIElKfU0w/slGNZ\nkSz1SFKLxlLqeQh4LzCFuBj8HuDBdgYnSeqcPC3+XmAesGeavxyYDdxZUEyN2OKXpBYN1+K3V48k\nTVLDJf5m/fhPAE4iWvv1BoA5bYlMktRRzRL/LenrNUSir+mpm5ckTSDNEv/56euTwE/q1r2zmHAk\nSUXLU+O/DnhtjmVFssYvSS0aTY3/AOBAYEtgbmbnjYFVbY5PktQhzRL/vUR9f1b6Wkv8y4GPFByX\nJKkgeUo96zM4HHO3WOqRpBaNptRT0wt8GXgl8Jy0bAB4WZtikyR1UN5B2k4DVgMV4Czg7AJjkiQV\nKE+p51pgF+BG4DV1yzrFUo8ktWgspZ4VxABttwPHERd9N2xncJKkzsnT4t+deNj6dOALwCbAfwBX\nFhhXPVv8ktSibg3SNhM4hfjE8B1i7J96FeBrwHrEcM+VBtuY+CWpRaNJ/Oc3WTcAHDzCe04BlgL7\nAfcAVwOHE58eaqYDlwH7A3cDm9F4rH8TvyS1aDQ1/q82WZcnC+9OXBe4I83PJ24Gyyb+dxPP8L07\nzfuAF0kqWLPEX81Mb0A8dH1pC8feErgrM383sEfdNtsRJZ5FxFAQXwd+0MJ7SJJalKcf/8HEoGwX\npPnXEg9fH0meTwXrEd1CDyTKPZ8hTgaSpILk6c7ZR7TUF6X568h31+49xKeEmhkMlnRq7iLKO0+l\n1++Ih7jftlYQfX3PTlcqFSqVSo4QJKk8qtUq1Wp1xO3y9Oq5ikj82aGYbwB2HGG/dYnS0L5E3//f\ns/bF3ZcD3yBa+1PTex3G4ENgary4K0ktGssNXDcDR6RttyMeuXh5jv1WEzd8XUD08DmDSPrHpvWn\nA0uAXxMnkjXAt1k76UuS2ihPi38D4NPAW9L8BcSNXCuKCqoBW/yS1KLR3sC1LnAhsE8BMbXCxC9J\nLRou8Y/Uq2c1UYKZXkBMkqQuyFPjf4IYmfPCNA3RVXNOUUFJkoqTJ/GfCyxgsF9+D/n66EuSxqE8\nNf7f0HjgtE6yxi9JLRpLjf8ZrPFL0qRhjV+SSiZP4l+ANX5JmjTyPohlKrB9ml4CrComnGFZ45ek\nFo1lyIYKcBbw5zS/NXAkcHGbYpMkdVCeFv+1xOBqtbH4tyceqrJLUUE1YItfklo02l49MDjKZs0f\nyPdJQZI0DuVJ4NcQD0r/IXHmOAJYXGRQkqTi5Cn1TAM+BOyV5i8BTgVWFhVUA5Z6JKlFox2dE2BD\nYgjmZ9L8FKKXz5PtCi4HE78ktWgsNf7fAs/JzG9ADOMgSZqA8iT+qcDjmfnHiOQvSZqA8iT+J4Bd\nM/O7EQ9GlyRNQHl69XwY+AlwX5p/EfFAdEnSBJR3yIb1gR3S9FLg6WLCGZYXdyWpRWPp1TMemPgl\nqUVj6dUjSZpETPySVDLNLu7uSvNx969tcyySpA5oVuOv0jzx79PeUJqyxi9JLfLiriSVzFgexALw\nGuAVxIBtNd8fe1iSpE7L0+LvA94EvApYCBwAXAocWlxYa7HFL0ktGkt3zkOB/Yg7d48GdgKmtzM4\nSVLn5En8TxFDMq8GNgX+AswoMihJUnHy1PivBp4LfJt48tYTwOVFBiVJKk6rvXpeCmwC9BcQSzPW\n+CWpRWOp8R/CYE3/T8Cfgbe2LTJJUkflafH3Exd0s64Hdm5/OMOyxS9JLRpLi7/RyWHKWAOSJHVH\nnsR/DXAysA2wLfC1tEySNAHlSfyzgVXAOcB8YAXwoZzHnwksAW4DTmiy3euI7qJvy3lcSdIoFTlW\nzxTiaV37AfcQ3UIPB25tsN2FwJPAmcC5DY5ljV+SWjSasXq+DhwPnN9g3QBw8AjvuTtwO3BHmp8P\nzGLtxD8b+CnR6pckFaxZ4q8NwvbVUR57S+CuzPzdwB4NtpkFvJlI/DbrJalgzRJ/7QJudZTHzpPE\nTwE+kbbtYeIMEy1JE1aeIRv2Bj4H9Ga2HwBeNsJ+9zB0TJ8ZRKs/a1eiBASwGTHy5yrgvPqD9fX1\nPTtdqVSoVCo5Qpek8qhWq1Sr1RG3y9PCXgp8mHjU4jOZ5Q+OsN+6ad99gXuB39P44m7NmcT1hAUN\n1nlxV5JaNJYHsTwC/GoU77kaOA64gOi5cwaR9I9N608fxTElSWOUp8X/FSJxLwBWZpZ38mHrtvgl\nqUVjeeZulcYXan3YuiSNYz5sXZJKZiyDtE1ncHyea4h+/Zu2MzhJUufkSfzfBZYD7wDeCTxG9MCR\nJE1Aox2Pv9GyIlnqkaQWjaXU8xTwxsz83sSAapKkCShPi39nYtyeWl3/YeBIOvvcXVv8ktSidvTq\n2SR9Xd6OgFpk4pekFo3lzt1pwNuJsXqmpIMMAJ9vX3iSpE7Jk/h/QQzbcA3x9C1J0gSWp9RzE/Dq\nogMZgaUeSWrRWHr1XA7s2O6AJEndkafFfyuwLfAnBgdpG6CzJwNb/JLUorH06ukdZvkdow+nZSZ+\nSWrRWHr13JG+vpDo4SNJmsDy1PgPBm4jSj0XEyeC0TyYRZI0DuRJ/F8EXg/8AXgp8SjFq4oMSpJU\nnDyJfxXxfN11iBu4FgG7FRmUJKk4eWr8DwMbA5cAZwN/AR4vMihJUnHy9OrZkLhjdx3gCGLMnrOB\nhwqMq569eiSpRaPpzrkdsDlwad3yvYH7gGXtCi4HE78ktWg0d+6eQuOROJendZKkCahZ4t8cuKHB\n8huI3j2SpAmoWeKf3mSdN3JJ0gTVLPEvBo5psPyDxBDNkqQJqNnF3S2AnwFPM5jodwWmAocQF3g7\nxYu7ktSi0Q7S1gPsQ4zHPwDcDPy23cHlYOKXpBa145m73WTil6QWjeVBLJKkScTEL0klY+KXpJIx\n8UtSyZj4JalkTPySVDImfkkqGRO/JJVMJxL/TGAJ8cD2ExqsPwLoJ0b9vAzYsQMxSVJpFX3n7hRg\nKbAfcA9wNXA4cGtmm9cDtwCPEieJPmDPuuN4564ktahbd+7uDtwO3EE8tH0+MKtumyuIpA9wFbBV\nwTFJUqkVnfi3BO7KzN+dlg3nA8AvC41Ikkpu3YKP30p9Zh/g/cBejVb29fU9O12pVKhUKmOJS5Im\nnWq1SrVaHXG7omv8exI1+5lp/kRgDXBS3XY7AgvSdrc3OI41fklqUbdq/IuB7YBeYH3gMOC8um22\nJpL+e2ic9CVJbVR0qWc1cBxwAdHD5wyiR8+xaf3pwGeB5wLfTMtWEReFJUkF8EEskjRJ+SAWSRJg\n4pek0jHxS1LJmPglqWRM/JJUMiZ+SSoZE78klYyJX5JKxsQvSSVj4pekkjHxS1LJmPglqWRM/JJU\nMiZ+SSoZE78klYyJX5JKxsQvSSVj4pekkjHxS1LJmPglqWRM/JJUMiZ+SSoZE78klYyJX5JKxsQv\nSSVj4pekkjHxS1LJmPglqWRM/JJUMiZ+SSoZE78klYyJX5JKxsQvSSVj4pekkjHxS1LJFJ34ZwJL\ngNuAE4bZZm5a3w+8tuB4JKn0ikz8U4BvEMn/lcDhwCvqtjkQ2BbYDjgG+GaB8SipVqvdDkFqyL/N\nzigy8e8O3A7cAawC5gOz6rY5GDgrTV8FTAc2LzAm4T+Xxi//NjujyMS/JXBXZv7utGykbbYqMCZJ\nKr0iE/9Azu16RrmfJGmc2RP4dWb+RNa+wHsa8K7M/BIal3quJ04Ivnz58uUr/+t6OmxdYBnQC6yf\nAmh0cfeXaXpP4MpOBSdJKsYBwFLiIu+Jadmx6VXzjbS+H9ilo9FJkiRJ0mTxXeAB4MZuByI1MANY\nBNwM3ATM6W440uTwRuLOaBO/xqMtgJ3T9EZEibj+mqCkUejFxK+J4efAvt0OYrJykDZJ400v8en0\nqi7HMWmZ+CWNJxsBPwWOBx7vciyTlolf0nixHnAu8EOi1KOCmPgljQc9wBnALcApXY5FmjR+DNwL\nrCQGxju6u+FIQ+wNrCHu8L8uvWZ2NSJJkiRJkiRJkiRJkiRJkiRJmvi2AOYTD95ZDCwEtuvQe+8K\nfH2U+1bT/pPdLBwRU1Ib9QBXAMdklu1I3Lgz3i2iHE+H+x7w9m4HIWnyeDNw8TDrNgR+A1wD3AAc\nnJb3ArcC3yIeznEBMC2t25l4RnM/sACYnpZXga8QIzsuZfDEUgHOT9MbAWem9+oHDknLTwWuTu/V\nl4lvuMT/OuAy4m7Tq9L3MS1z7GvT+wIcRYw/83/An4DjgI+lba4AnpuJ/xTiztUb03sAPC/t35+2\nf01a3kc8ZGcR8Yzr2Zn43pPiug44jcEhWh4HvpjivgJ4IfAG4CHgj2n7lzX4fiWpJXOAk4dZNwXY\nOE1vBtyWpnuBVcQnA4BzgCPS9A3Ew2UA/g34WppeBPxnmj4AuDBNVxhM/CfVxVI7adSS75R0nFpy\nbZT41ycSba0EtFHa76PAd9KyHYA/A1OJxH8bcXLYDHiUwU8/JxOjUdbe6/Q0/UYGn58wD/hMmt6H\nSM4Qif9SYoCz5wMPpjheAZyXpiFOau9N02uAgzI/i0+l6TOBt6FJb91uB6DSGGiybh3g34lEtwZ4\nMdEKhWgd35CmryFOBpsAmwKXpOVnAf+TOd6C9PXatH29fYHDMvOPpK+HAR8k/i9eRCTP4R5cswNw\nX4oJBocQ3guYm6aXEol/e+L7XwQ8kV6PMHgiupHBkxvEuEqk76/2ve7FYFJeRCT5jdNxFxInyIeA\nvxDXUvYlTkqL0z7PAe5P00+nfUjx/13mvXuG+X41iZj41Sk3A4cOs+4IohW8C/AMkexrJZ2Vme2e\nySzPqk9WtX2eYfi/8fp9Xkq01ncjWuNnDvNeeQyXPLPfy5rM/Bqa/y/WTprDHffpzHT2ez4L+GSD\n7VfVxZF972YnaE0SDsusTvktUfL4YGZZ7eLuJkRL9RmijPGSJsfpAZYDDzNYv38vURvP60LgQ5n5\n6SmGJ9KxNyfKRM0sJT4V7JbmNybKKpcwWI7aHtgaWELzlnRP3XTt08jexCeD5XXHrQB/BR4b5rgD\nwEXEifYFadnzUizNPEb8HDTJmfjVSYcA+xHdOW8CvkSUS84mEugNRBK/NbNPfQu0Nn8kUcvvJ04g\nnx/mPQcaTH+RqOffSFzgrKTjXEck6bOJunkzTxMJel46xgXEie1U4v/qBqLr6pFEC3tgmFioWzcA\nrCDKVKcCH0jL+4jSTT/w5XTc+n2zbgU+TVxM7k9ftxjhvecD/0qUf7y4K0kdUpauo+oiW/ySJEmS\nJEmSJEmSJEmSJEmSJHXO/wMeBrPPxF7OwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcc00064a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot canonical correlations (cca.cancorrs)\n",
    "plt.plot(np.arange(ccaCV.best_numCC)+1, ccaCV.cancorrs, 'ko')\n",
    "plt.xlim(0.5, 0.5+ccaCV.best_numCC)\n",
    "plt.xticks(np.arange(ccaCV.best_numCC)+1)\n",
    "plt.ylim(0,1)\n",
    "plt.xlabel('Canonical component')\n",
    "plt.ylabel('Canonical correlation')\n",
    "plt.title('Canonical correlations')\n",
    "print '''The canonical correlations are:\\n\n",
    "Component 1: %.02f\\n\n",
    "Component 2: %.02f\\n\n",
    "''' % tuple(ccaCV.cancorrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The canonical correlations for both components are both 0.95."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cross-dataset predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The prediction accuracy for the first dataset is:\n",
      "\n",
      "Dimension 1: 0.94\n",
      "\n",
      "Dimension 2: 0.93\n",
      "\n",
      "Dimension 3: 0.93\n",
      "\n",
      "Dimension 4: 0.93\n",
      "\n",
      "\n",
      "The prediction accuracy for the second dataset is:\n",
      "\n",
      "Dimension 1: 0.93\n",
      "\n",
      "Dimension 2: 0.93\n",
      "\n",
      "Dimension 3: 0.93\n",
      "\n",
      "Dimension 4: 0.93\n",
      "\n",
      "Dimension 5: 0.93\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEZCAYAAACQK04eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//FXEnazEAQVIRAFRFErFBcqqANYoSJutcoi\nX7F+RfpzQVsVxa81ar8VXFARBVxQKu7UKkIFUQxQlyBSgqCg8AULCLiwKsiW+/vjnGQmw8zkTjI3\nM5P7fj4e85i7nLn3MwP5zJ1zzj0HREREREREREREREREREREREREREREROL2HHCvXT4DWFHN40wA\n/icRAYmICKwFdgE7gU3As8AhCTr2s8A9cb5mKLAgQecXSVuZyQ5A6jQHOA/IAX4JnEzkq+t61Tx+\nRjVf5wdZyQ5AUpcSv9SWb4BZwPF2vQz4f8BXwEq77TxgCbAV+AA4MeT1XYDFwA7gZaBRyL4AsC5k\nvQB4HfgW+B54DDgWmAj8CvMLZIst+xzBKiOAq21MPwBvAi1D9pUB1wBf2hjHx3i/pwIf2XLf2Bjq\nh+w/Hphjz7MJuN1uzwJGAavse10EtAIK7flD/2aLgavs8lDMZzbWvue7gKOBuXb9O2AqkBfy+vDP\naZyNcQtwQki5w4CfgENjvF8REQDWAL3tcgGwDLjbrpcBs4GmQENMYt8MnIK5kv8v+/r6QAPga2AE\nJjH+FthLsKonQDDxZwGlwENAY3vs0+2+Kzi4qie0yqgXJkF2tuccB8wLKVsGTAdy7fv5FugT5b3/\nEpP8M4G2wOc2fjC/gDYCN9nzZNuyALcAS4EOdv1EoBmRE//7wO/t8lBgH3CtLdMIaIf5/OsDze17\nediWj/U5PQ6MDjnPCMyXoIhIldZirq632uXxmAQDJokFQspO4OA6+xXAmfaxIWzfB0RO/L/CJORI\nv2aHEjvxP0PlhHcI5gumTUjMp4fsfwUYGeE8kdyIuboGGAh8GqXcCqB/hO2FVJ34v64ihgsxv5og\n9ud0WtixFgGXVHFsSSPVrVsVccMBLsBUN0QSWj3TFnOVf33ItvqYqpYMDk780ZJcgd1XFm+w9lyL\nQtZ/wlTFtAL+Y7dtCtm/C3O1HskxmGqXrkATzN9a+bELgP+L8roCYHW8gVvrwtYPBx4FemB+ZWQS\nrOKK9TmVALsxX6ibML8cplczJklBquOXZHJClv8D/C+QH/LIxlxVb8Qk31BtoxxzHeYKPVLjphNh\nW6hvMFfW5Q7B1GuHf+m4MQFTvdMeU69+B8G/t/9g6t8jWWdfE+4n+9wkZNsRYWXC399fgQOY+vo8\nYEhIDLE+J4ApwOX2Na9hfvlIHaHEL6niKWA4pq47A5N0+2GS/4fAfuAGzK+AizFtAZEsxHxRjMYk\nyUYEq2c2A62p3MiaQbB30EvAlcBJmCqpvwIfE7zaDxerV1E2ppprF6Zh+Q8h+2Zifl2MsOfJIVjH\n/zSmsbm9Pf4vMHX832G+gIZgkvXvMVfisWRjvjB2YL44bwnZF+tzAtMQfDEwGPhbFeeRNKPEL8kS\nfnX6KaZHzXhMdcRXmKofMI2WF2PqsX8ALgX+HuV4BzB15O0xCXudLQ/wHrAcU33xbcjrnJD9d9pj\nfwMcBQyIEbMTYVu5m4FBmKT7JKYnUnnZncCvbZwbMb2EAnbfWOBV4B1gO+YLsbwH09WY5P090AnT\nzhErlrsxjczbgbfs+3LzOWHXF2Oqgv4V5T2KRDQZc5X1WYwy4zB/5KWYnh0ikhqeIf6b5EQ4A5PM\noyX+c4F/2uXTMD+rRST5CjG9saK1pYjEVEj0xD8RuCxkfQWmJ4KIJM+9mOqo26sqKBJNIdET/1tU\nblB6F9P9TUREPJIKjbvhPSOq6nInIiI1kOwbuDZgbiQp15oIfaZPOukkp7S0tNaCEhGpI0oxQ5BU\nkuwr/ukEu+x1A7ZhegFVUlpaiuM4CX3cddddCT+mFw/F6a8Y0ylOPVL/gbkn5SBeX/G/BJyFGSBq\nHWbEwPKbZyZhevScixmJ8CfMzTMiIuIhrxP/QBdlrvM4BhERCZHsqp6kCQQCyQ7BFcWZOOkQI6RP\nnJK+0mUGI8fWV4mIiEsZGRkQIc8nu1ePiKSZZs2asXXr1mSHISHy8/PZsmVL1QUtXfGLSFwyMjLQ\n32NqifZvEu2K37d1/CIifuW7qp6Z777DhNdeYC9lNCCTP/xuMP3OPifZYYmI1BpfJf6Z777DbZPH\nsfXXJ1Rsu23yOAAlfxHxDV9V9Ux47YVKSR9g669PYOK0F5MUkYhI7fNV4t9LGZnLV9Nyyqsc9dyL\ntJzyKpnLV7PHOZDs0A5y/5h76dW5FX27HE6vzq24f8y9yQ4ponSI8/3ZM7l5UF9uvSzAzYP68v7s\nmckOSTxSWFhIkyZNyM3NJT8/n+7duzNp0iTXjdFr164lMzOTsrJIc9AnjpvzLFu2jD59+tCiRQsy\nMxObqn1V1bP3mw10/nYJDx63vWLbzQu2sfPw1Jr46/4x9/Lh1DFM7Lq7YtutU8dwP3DryDuTF1iY\ndIjz/dkzmXrPMEZ22FSxbcw9nwFP0rNPv+QFFsH7s2cyc8pjZB74mbKsRvS74vqUizGWRLSf1fQY\nGRkZzJgxg169erFz506Ki4sZMWIEJSUlTJ482fVxaqvXUqzzNGjQgAEDBnDttddy4YUXJvS8vrri\nb7tvR6WkD/Dgcdtpu29HkiKKbNZLE7k/JJkC3N91N7NenpikiCJLhzife/iuSkkfYGSHTUx5pCg5\nAUXx/uyZ/OOBEQxrMIf/bryAYQ3m8I8HRqTNr5Py9rMl3VryebdWLOnWktsmj2Pmu+/U6jFC5eTk\n0L9/f1555RWmTJnC8uXLzXlmzqRLly7k5eXRpk0b7r777orXnHnmmQA0bdqUnJwcSkpKWL16Nb16\n9aJ58+a0aNGCyy+/nO3bg3lkzJgxtG7dmtzcXI499ljmzp0LmKQ+evRo2rdvT/Pmzbnssssq7n+I\ndJ5wxxxzDFdeeSWdOnWq1vuPxVeJ//C87IjbD8s9pJYjia1Bxv7I24m8PVnSIc5vN6+PuH3zpnW1\nHElsM6c8xnWt11Tadl3rNcx8fnySIopPItrPvGqDO+WUU2jdujX/+peZMz47O5upU6eyfft2Zs6c\nyYQJE3jzzTcBWLBgAQDbt29n586dnHbaaQDccccdbNy4kS+++IJ169ZRVFQEwMqVK3n88cdZtGgR\nO3bs4J133qGwsBCAcePGMX36dObPn8/GjRvJz8/n2muvjXme2uKrxF+W1Sjy9vqNazmS2PY6kWvg\n9qZYzVw6xLnbifxffHeK/dfPPPBz5O37dkfcnmr2ErmuOp72s0QcI5ojjzyy4s7Ws846i+OPPx6A\nE088kQEDBjBv3jwgctVLu3bt6N27N/Xr16d58+bcdNNNFeWzsrLYs2cPy5cvZ9++fbRp04ajjz4a\ngEmTJvGXv/yFI488kvr163PXXXcxbdo0ysrKkn4DXGr97/dYvyuuZ/z6oypte2z9UfQbkloDhPYd\nOJxbP638ZXTLokb0HTA8SRFFlhZxFhzLzV/kVdr0p8/zyCg4LkkBRZYuFyXRNIiSShpmZNXqMaLZ\nsGEDzZo1A6CkpISePXty2GGH0bRpUyZNmsQPP/wQ9bWbN29mwIABtG7dmry8PIYMGVJRvn379jzy\nyCMUFRVx+OGHM3DgQDZu3AiYBtyLLrqI/Px88vPz6dSpE/Xq1WPz5oOmHKl1vkr8Pfv046JbHuXJ\n/efw9O4zeHL/OVx8y6Mp14B268g7Of3ykQwvPYIbSpszvPQIug+5LWUaTMulQ5y3jLiNVTntGLy2\nLVetacXgtW1ZndOOm28YmezQKkmXi5Jo/vC7weTPWVZpW/47yxh+yaBaPUYkn3zyCRs2bKBHjx4A\nDBo0iAsvvJD169ezbds2hg8fXtG7xg5xUMmoUaPIyspi2bJlbN++neeff75Sb5yBAweyYMECvv76\nazIyMhg50vzfatOmDbNmzWLr1q0Vj127dtGyZcuI56lNqfObvJb07NMv5RJ9JLeOvDOlEmg0qR5n\neY+QidNeZI9zgKYZWQy/ZFDK3bBX/n/yyefHk7lvN2X1G3PxLdelxf9VOPhzbpiRxfCrbojrc07E\nMSBYXbNjxw7mz5/PjTfeyJAhQyqqd3788Ufy8/Np0KABCxcu5MUXX6RPnz4AFV0nV69eTYcOHSrK\n5+XlkZuby4YNG3jggQcqzvXll1+yfv16unfvTsOGDWnUqFHF+YcPH86oUaOYMmUKbdq04bvvvuOj\njz7i/PPPj3ieSH7++Wf27t0LwJ49ewBo2LBhXJ9HOnNEJDWk8t9jYWGh07hxYycnJ8fJy8tzTj/9\ndOeJJ55wysrKKspMmzbNadu2rZOTk+Ocd955zvXXX+8MGTKkYv+f//xnp0WLFk5+fr5TUlLiLF++\n3OnatauTnZ3tdOnSxXnooYecgoICx3EcZ+nSpc6pp57q5OTkOM2aNXP69+/vbNy40XEcxykrK3PG\njh3rdOzY0cnJyXHatWvn3HHHHQedp2nTpk5JSclB72XNmjVORkaGk5GR4WRmZjoZGRnOUUcdFfF9\nR/s3ASI2Jmh0ThGJi0bnTD0anVNERGJS4hcR8RklfhERn1HiFxHxGSV+ERGfUeIXEfEZJX4REZ9R\n4hcR8RklfhERn1HiF5E6oy5NvThlyhROPvlk8vLyKCgoYOTIkRw4kJhpYn03SJuIeCMRU0fW9Bh1\naerF3bt38+ijj3Laaafx7bffcv755/Pggw9WjP7pB3EM0yQiXor09zh31gzn+t7tnJVXZFY8ru/d\nzpk7a4br4ybiGIWFhc57771XadvChQudzMxMZ9myZY7jOM6MGTOczp07O7m5uU5BQYFTVFRUUbag\noMDJyMhwsrOznezsbOfjjz92Vq1a5fTs2dM59NBDnebNmzuDBw92tm3bVvGa0aNHO61atXJycnKc\njh07Vpy/rKzMue+++5x27do5hx56qHPppZc6W7ZsiXqeqowdO9bp379/xH3RciRRBmlTVY+I1Fgi\npo70avrJujL14rx58zjhhBOqLOeGEr+I1Fgipo70cvrJdJ96cfLkySxevJibb7652p9BKCV+Eamx\nREwd6eX0k+k89eIbb7zBqFGjePvttyveQ00p8YtIjSVi6kivpp9M56kXZ82axbBhw5gxY0bFr5RE\nUK8eEamxREwdmajpJ506MvXi3LlzGTx4MG+++SYnn3xyXJ9BXVFli7eI1I5U/nusS1Mv9uzZ06lf\nv35Fz5/s7Gzn3HPPjfi+o/2bkKSpF/sCjwBZwNPAmLD9zYGpwBGYXx8PAs9FOI59DyKSbJp6MfXE\nO/Wil4k/C1gJnA1sAD4BBgJfhJQpAhoCt2O+BFYChwP7w46lxC+SIpT4U08qzbl7KrAKWAvsA14G\nLggrsxHItcu5wA8cnPRFRCSBvGzcbQWsC1lfD4TfpfAUMBf4BsgBLvUwHhERwdvE7+a34ChgCRAA\n2gFzgJOAneEFy++UAwgEAgQCgQSEKCJSdxQXF1NcXFxlOS/r+Lth6vD72vXbgTIqN/D+E/hf4AO7\n/h4wElgUdizV8YukCNXxp55UquNfBHQACoEGwGXA9LAyKzCNv2AadTsC/+dhTCIivudlVc9+4Dpg\nNqaHzzOYHj3X2P2TgL8CzwKlmC+hW4EtHsYkIjWUn5/v+s5TqR35+flxlU+Xfz1V9YiIxCkZVT0i\nIpKClPhFRHxGiV9ExGeU+EVEfEaJX0TEZ5T4RUR8RolfRMRnlPhFRHxGiV9ExGeU+EVEfEaJX0TE\nZ9wM0nYYcDVmlM3y8g7we49iEhERD7lJ/G8C8zGTpJTZbRoxTUQkTbkZnXMJ0NnrQKqg0TlFROJU\nk9E5ZwD9Eh2QiIgkh5sr/h+BJsBeYJ/d5gC5XgUVga74RUTiFO2KXxOxiIjUUdESv9upFy8AzsRc\n6c8D3kpYZCIiUqvcXPGPBk4BXrDlB2AmUr/dw7jC6YpfRCRONanq+QzTq+eAXc/C9PQ5MVHBuaDE\nLyISp5r06nGApiHrTVE/fhGRtOWmjv8+YDFQbNfPAm7zKiAREfGW2149R2Lq+R1gIbDJs4giU1WP\niEicqlPHfxzwBdAVk/DLy5Zn4MUJjK8qSvwiInGqTuJ/CjM4WzGR6/R7JiIwl5T4RUTiVJNePY2A\nn11s85ISv4hInGrSq+dDl9tERCQNxOrV0xLTqNsE+CXmW6N8jJ4m3ocmIiJeiJX4zwGGAq2Ah0K2\n7wRGeRiTiIh4yE0d/yXANK8DqYLq+EVE4lTT0TnPAzphGnXL3VPzsFxT4hcRiVNNGncnAZcCN9gD\nXAq0TWRwIiJSe9wO0nYisBT4BZANzAJ6eBhXOF3xi4jEqSZX/Lvt8y5MQ+9+4IiERSYiIrXKzSBt\nM4B84AHgU7vtKc8iEhERT8U79WIj+9jmsnxf4BHMGP5PA2MilAkADwP1ge/tejhV9YiIxKk6vXp+\nS3CMngwOHq/n9SrOmQWsBM4GNgCfAAMxA7+Vawp8APQB1gPNMck/nBK/iEicqjPnbn9iT7hSVeI/\nFVgFrLXrL2Pm7g1N/IOAv2OSPkRO+iIikkCxEv/QGh67FbAuZH09cFpYmQ6YKp73gRzgUeD5Gp5X\nRERicNOr5wjgGUwXTjA3cl3l4nVu6mbqY8YBOhdT3XMn5stAREQ84qZXz3PAs8Addv0r4FXMl0Es\nG4CCkPUCglU65dZhqnd228d84CR7jkqKiooqlgOBAIFAwEXoIiL+UVxcTHFxcZXl3PTqWQScDPwb\n6GK3LQE6V/G6epjG3d7AN5gpG8Mbd48FxmOu9hsCJcBlwOdhx1LjrohInKrTuFvuR+DQkPVuwHYX\nr9sPXAfMxvTweQaT9K+x+ycBKzBVSEuBMsz9AeFJX0REEsjNFX9X4DHgeGA50AIzYmeph3GF0xW/\niEicqnvFnwWcaR/H2gOsBPYmOD4REaklbq74PwFO8TqQKuiKX0QkTjUZj798OIVXgJ8I3sW7OIHx\nVUWJX0QkTjVJ/MVE7pPfs2YhxUWJX0QkTjWp458OjPUgJhERSYKq7tw9gOl7LyIidYTq+EVE6ijV\n8YuI+ExNEn8qUOIXEYlTTebcbYqp7vnUPh4C8hIZnIiI1B43iX8ysAP4HXApsBMzWqeIiKQhN1U9\npZihkqva5iVV9YiIxKkmVT27gTNC1nsAuxITloiI1DY3V/ydgb8RrNffClyBRucUEUlpiejVU574\n3YzFn2hK/CIicapJVc99mJ492+0jH/hLIoMTEZHa4ybx/wbYFrK+FejnTTgiIuI1N4k/E2gUst4Y\naOBNOCIi4jU3c+6+ALyH6c+fAVyJaewVEZE05LZx9zdAb7s8BzOBem1S466ISJw0Vo+IiM/UpFeP\niIjUIUr8IiI+o8QvIuIzbnr19ADuAgpDyjvA0R7FJCIiHnLTuLsSuBEz1eKBkO3fexJRZGrcFRGJ\nU7TGXTdX/NuAtxMdkIiIJIebK/7RQBbwOrAnZLsmWxcRSWGabF1ExGd0A5eIiM9osnUREQE02bqI\niO9osnURkTpKk62LiAigydZFROqsRPTqybXPOxIRUJyU+EVE4lSdO3eHAM8Df6JyP/4Muz7WxXn7\nAo9gbgB7GhgTpdwpwEeYxuPXXRxXRESqKVbib2Kfc4h8A1dVsoDxwNnABuATYDrwRYRyY4BZpM99\nBSIiaStW4p9kn98F/hW2r4eLY58KrALW2vWXgQs4OPFfD0zDXPWLiIjH3PTqeSzCtnEuXtcKWBey\nvt5uCy9zATDBrqsiX0TEY7Gu+H8FnA60AP5IsBomB1M9UxU3SfwR4DZbNgNV9YiIeC5W4m9AMMnn\nhGzfAVzi4tgbgIKQ9QLMVX+orpgqIIDmwG+AfZi2gEqKiooqlgOBAIFAwEUIIiL+UVxcTHFxcZXl\n3FxhtwW+rkYM9TCTuPQGvgEWAgM5uI6/3LPAW0Tu1aPunCIicarJnbtPYwZqK9cMmO3idfuB62zZ\nz4FXMEn/GvsQEZEkcHPFvwRz925V27ykK34RkTjV5Ir/AKa6p1whUJaQqEREpNa5mXP3DmABMN+u\nnwkM8ywiERHxlNvuky2Abphulx8D33sWUWSq6hERiVN1Bmk7DtMY25VgP3sI9s/XZOsiIimsOon/\nKeBqNNm6iEha0mTrIiI+U51hmX9L7GEXNHyyiEgaipX4+2MS/2GYMXvm2u09gQ9R4hcRSUuxEv9Q\n+zwH6ARstOstgSkexiQiIh5ycwNXAbApZH0z0MabcERExGtubuB6FzPezouYRoLLML8CREQkDbnp\n1ZMBXAScYdfnA//wLKLI1KtHRCRO1enVU87B3Ky1E3Ol3wQzPv/OBMYnIiK1xE0d/zDgNWCiXW8N\nvOFZRCIi4ik3if9azOTqO+z6l5guniIikobcJP499lGuHpoUXUQkbblJ/PMwQzM3AX6NqfZ5y8ug\nRETEO2569WQC/w2cY9dnY6ZjrM2rfvXqERGJU3UHaasHLAOO9SCmeCjxi4jEqbpTL+4HVlJ56kUR\nEUljbvrxNwOWAwuBn+w2Bzjfq6BERMQ7bhL//9jn0J8LqncREUlTsRJ/Y2A40B5YCkwG9tVGUCIi\n4p1YjbuvAnuBBcC5wFpgRC3EFIkad0VE4lSdXj2fASfa5XrAJ0CXhEfmjhK/iEicqtOrZ3+UZRER\nSWOxrvgPALtC1hsDu+2yA+R6FVQEuuIXEYlTdYZlzvIsGhERSRo3Y/WIiEgdosQvIuIzSvwiIj6j\nxC8i4jNK/CIiPqPELyLiM0r8IiI+o8QvIuIztZH4+wIrgK+AkRH2DwZKMSOAfgD8ohZiEhHxLTdz\n7tZEFmYGr7OBDZiB3gYCX4SU+RXwObAd8yVRBHQLO46GbBARiVN1p16sqVOBVZghnfcBLwMXhJX5\nCJP0AUqA1h7HJCLia14n/lbAupD19XZbNFcB//Q0IhERn3Mz9WJNxFM/0xP4PdA90s6ioqKK5UAg\nQCAQqElcIiJ1TnFxMcXFxVWW87qOvxumzr6vXb8dKAPGhJX7BfC6LbcqwnFUxy8iEqdk1fEvAjoA\nhUAD4DJgeliZNpikfzmRk76IiCSQ11U9+4HrgNmYHj7PYHr0XGP3TwL+DOQDE+y2fZhGYRER8YDX\nVT2JoqoeEZE4JauqR0REUowSv4iIzyjxi4j4jBK/iIjPKPGLiPiMEr+IiM8o8YuI+IwSv4iIzyjx\ni4j4jBK/iIjPKPGLiPiMEr+IiM8o8YuI+IwSv4iIzyjxi4j4jBK/iIjPKPGLiPiMEr+IiM8o8YuI\n+IwSv4iIzyjxi4j4jBK/iIjPKPGLiPiMEr+IiM8o8YuI+IwSv4iIzyjxi4j4jBK/iIjPKPGLiPiM\nEr+IiM8o8YuI+IwSv4iIzyjxi4j4jBK/iIjPKPGLiPiM14m/L7AC+AoYGaXMOLu/FOjicTwiIr7n\nZeLPAsZjkn8nYCBwXFiZc4H2QAdgGDDBw3gqKS4urq1T1YjiTJx0iBHSJ05JX14m/lOBVcBaYB/w\nMnBBWJnzgSl2uQRoChzuYUwV0uWPS3EmTjrECOkTp6QvLxN/K2BdyPp6u62qMq09jElExPe8TPyO\ny3IZ1XydiIikmG7ArJD12zm4gXciMCBkfQWRq3qWYL4Q9NBDDz30cP9YQi2rB6wGCoEGNoBIjbv/\ntMvdgI9rKzgREfHGb4CVmEbe2+22a+yj3Hi7vxT4Za1GJyIiIiIidc9kYDPwWbIDiaEAeB9YDiwD\nbkhuOFE1wnTDXQJ8DtyX3HCqlAX8G3gr2YHEsBZYiolzYXJDEak7zsDcIZzKif8IoLNdzsZUl4W3\nj6SKJva5HqaNpkcSY6nKH4EXgOnJDiSGNUCzZAchdZsfx+pZAGxNdhBV2ESwNf5H4AvgyOSFE9Mu\n+9wAc0W9JYmxxNIa05ngaQ7uQpxqUj0+SXN+TPzpphDzC6UkyXFEk4n5ktqMqZ76PLnhRPUwcAtQ\nluxAquAA7wKLgKuTHItInVJIalf1lMvGJIALkx2IC3mYqp5AkuOI5DzgcbscILXr+Fva5xaYL9Qz\nkhiL1FG64k9d9YG/A1OBN5IcixvbgZnAyckOJILTMeNCrQFeAnoBf0tqRNFttM/fAf/AjHklIglQ\nSGpf8WdgEtPDyQ6kCs0xA+sBNAbmA72TF44rZ5G6V/xNgBy7fAjwAXBO8sIRqTteAr4B9mAGiLsy\nueFE1ANTF70E063v35jhrVPNicBiTJxLMXXoqe4sUrdXz1GYz3IJphvv7bGLi4iIiIiIiIiIiIiI\niIiIiIiIiIiIiIjUyAHMPQvLMH3Z/0jVA5a1BQZ6EMuNmJvSqhIgeENYfw6eXrQ2HAm8loTziojU\n2M6Q5RbAHKCoitcE8OZO3DXAoS7KeXV+ERFf2Bm2fhTwvV0uxAwJ8al9/Mpu/xjYhvmlMALzCyBS\nuZZ2+78xw3WUzxtwDvChLfsqZqiEGzB3di8F3osQZ1/MUNmfAo8STPxDgcfs8nPAE8BHmLmnA8AU\nzOilz4YcK9L5wUzIUmS3LwU62u1nEbybe7EtX0hwCJJG9vhL7f5ASGyvA28DXwJjIrwvEZFaF574\nwcyd0AJT7dLQbusAfGKXw8feiVbuT8Aou5yJGfm0OTCPYJXOSOBOuxxtQpRGwH+Adnb9FYJDQAyl\ncuJ/0S6fD+wAjsdUXS0CTnJx/mvt8h+Ap+zydIJfZk0wcyAUEkz8f8LMOQDmy+JrzOcxFPMFlGPX\n1wKtIrw/8YF6yQ5AxKUGwHhMwjyASepwcBtAtHILMdNu1seMdlqKuRruhLniLn/th8R2LCYpr7br\nU4FhEco5BL+QlmEm11lu15djknVBFed/3T4vBi62yx9gBu97we7fEHbe7sA4u7wSk/iPsfG8R/DL\n9XMbQ/jrxQeU+CWVHY1J3t9hqj02AkMwV7k/R3nNTVHKLcCMbX8e5mp8LObXxBxgUBwxOWHrsRqf\n99rnMkzVESHr9TDvLdb5y19zgODf6hhgBtAP8yXQJ+zYsWIKLXcA8/mID2k8fklVLYCJBKtOcjFX\nzQD/RTCf8rJgAAAA4UlEQVRp7SQ4lHGscm0wXyBP20cXTPtAd4LVNocQ/IWw0x4r3ErMlfLRdr26\nPYqcKs4fTTvML4b7MdVYHcP2LwAG2+VjMO97BZG/DDTFo08p8UsqaUywO+ccYBZwj933BHAFpptn\nR8xcxGCqbA7Y7SNilOtpty0GLsU0yn6Pqft+yR7nQ4KJ9El7/vDG3Z8xVTszMQ2vmwn+CnCo/Isg\n2nK5WOcPV/76EZj6/FLML4q3w/Y/gfm7Xgq8jPks9kWILVpMIiIiIiIiIiIiIiIiIiIiIiIiIiIi\nIiIiIiIiUl3/H/X9AOPuLrI0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcc1ecae490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot correlations between actual test data and predictions \n",
    "# obtained by projecting the other test dataset via the CCA mapping for each dimension.\n",
    "nTicks = max(testcorrsCV[0].shape[0], testcorrsCV[1].shape[0])\n",
    "bmap1 = qualitative.Dark2[3]\n",
    "plt.plot(np.arange(testcorrsCV[0].shape[0])+1, testcorrsCV[0], 'o', color = bmap1.mpl_colors[0])\n",
    "plt.plot(np.arange(testcorrsCV[1].shape[0])+1, testcorrsCV[1], 'o', color = bmap1.mpl_colors[1])\n",
    "plt.xlim(0.5, 0.5 + nTicks + 3)\n",
    "plt.ylim(0.0, 1.0)\n",
    "plt.xticks(np.arange(nTicks)+1)\n",
    "plt.xlabel('Dataset dimension')\n",
    "plt.ylabel('Prediction correlation')\n",
    "plt.title('Prediction accuracy')\n",
    "plt.legend(['Dataset 1', 'Dataset 2'])\n",
    "print '''The prediction accuracy for the first dataset is:\\n\n",
    "Dimension 1: %.02f\\n\n",
    "Dimension 2: %.02f\\n\n",
    "Dimension 3: %.02f\\n\n",
    "Dimension 4: %.02f\\n\n",
    "'''% tuple(testcorrsCV[0])\n",
    "print '''The prediction accuracy for the second dataset is:\\n\n",
    "Dimension 1: %.02f\\n\n",
    "Dimension 2: %.02f\\n\n",
    "Dimension 3: %.02f\\n\n",
    "Dimension 4: %.02f\\n\n",
    "Dimension 5: %.02f\\n\n",
    "'''% tuple(testcorrsCV[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cross-dataset predictions range between 0.93 and 0.94 for all dimensions for both datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variation of the optimal hyperparameters over many iterations\n",
    "The optimal regularization coefficient will vary with different instantiations of the data. It's possible that the optimal number of components would vary as well. Below, the data are randomly generated in the same manner for 1000 iterations. Then, the analysis is repeated to determine which optimal number of canonical components and optimal regularization coefficient persist over different instantiations of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nIterations = 1000\n",
    "all_numCC = np.zeros((nIterations,))\n",
    "all_reg = np.zeros((nIterations,))\n",
    "\n",
    "for ii in range(nIterations):\n",
    "    # Initialize number of samples\n",
    "    nSamples = 1000\n",
    "    \n",
    "    # Define two latent variables (number of samples x 1)\n",
    "    latvar1 = np.random.randn(nSamples,)\n",
    "    latvar2 = np.random.randn(nSamples,)\n",
    "    \n",
    "    # Define independent components for each dataset (number of samples x dataset dimensions)\n",
    "    indep1 = np.random.randn(nSamples, 4)\n",
    "    indep2 = np.random.randn(nSamples, 5)\n",
    "    \n",
    "    # Create two datasets, with each dimension composed as a sum of 75% one of the latent variables and 25% independent component\n",
    "    data1 = 0.25*indep1 + 0.75*np.vstack((latvar1, latvar2, latvar1, latvar2)).T\n",
    "    data2 = 0.25*indep2 + 0.75*np.vstack((latvar1, latvar2, latvar1, latvar2, latvar1)).T\n",
    "    \n",
    "    # Split each dataset into two halves: training set and test set\n",
    "    train1 = data1[:nSamples/2]\n",
    "    train2 = data2[:nSamples/2]\n",
    "    test1 = data1[nSamples/2:]\n",
    "    test2 = data2[nSamples/2:]\n",
    "\n",
    "    # Analysis\n",
    "    # Initialize a cca object as an instantiation of the CCACrossValidate class.\n",
    "    ccaCV = rcca.CCACrossValidate(kernelcca=False, numCCs = [1, 2, 3, 4], regs = [0., 1e2, 1e4, 1e6], verbose = False)\n",
    "\n",
    "    # Use the train() and validate() methods to run the analysis and perform cross-dataset prediction.\n",
    "    ccaCV.train([train1, train2])\n",
    "    \n",
    "    all_numCC[ii] = ccaCV.best_numCC\n",
    "    all_reg[ii] = ccaCV.best_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of times each number of components was chosen:\n",
      "\n",
      "1 component: 14\n",
      "\n",
      "2 components: 972\n",
      "\n",
      "3 components: 14\n",
      "\n",
      "4 components: 0\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print '''Number of times each number of components was chosen:\\n\n",
    "1 component: %d\\n\n",
    "2 components: %d\\n\n",
    "3 components: %d\\n\n",
    "4 components: %d\\n\n",
    "''' % tuple([(all_numCC == i).sum() for i in (1, 2, 3, 4)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal number of components is two for 97% of the iterations, corresponding to the number of latent variables. Pyrcca successfully captures the relationships in the data prescribed by the latent variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of times each regularization coefficient was chosen:\n",
      "\n",
      "Regularization coefficient = 0: 86\n",
      "\n",
      "Regularization coefficient = 1e2: 295\n",
      "\n",
      "Regularization coefficient = 1e4: 296\n",
      "\n",
      "Regularization coefficient = 1e6: 323\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print '''Number of times each regularization coefficient was chosen:\\n\n",
    "Regularization coefficient = 0: %d\\n\n",
    "Regularization coefficient = 1e2: %d\\n\n",
    "Regularization coefficient = 1e4: %d\\n\n",
    "Regularization coefficient = 1e6: %d\\n\n",
    "''' % tuple([(all_reg == i).sum() for i in (0, 1e2, 1e4, 1e6)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal regularization coefficient is greater than zero for over 90% of the iterations. The other three regularization coefficients are optimal for approximately equal proportions of iterations. This result corresponds to the random variation in the data that are generated during each iteration of the analysis."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}