MMCRF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Structured Output Prediction of Anti-Cancer Drug Activity [2/2]\n",
    "\n",
    "Anas Atmani, Benoît Choffin, Domitille Coulomb, Paul Roujansky"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## B - Method 2 : MMCRF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import time\n",
    "from sklearn.metrics import f1_score\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2305, 59)"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_reg = np.loadtxt('data/data_clean/ncicancer_targets.txt')\n",
    "target_reg.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2305, 59)"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target = np.loadtxt('data/data_clean/ncicancer_bin_targets.txt')\n",
    "target.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2305, 2305)"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Kx = np.loadtxt('data/data_clean/ncicancer_input_kernel.txt')\n",
    "Kx.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Data preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> On the github page of the project, it is written: \"We notice that the activity score is computed according -10*log(GI50). We still need binary value as the activity outcome for classification task. According to NCBI, a molecule is 'active' if the activity score is over(**=) 60 and 'inactive' otherwise.\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGICAYAAAAedKdVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8XFd9///XrNply7a8xkvWY8cQEhLCks2UQJYWUpbQ\nfqFtAm2gNL8fFGihgdASmhbKD2ibltLWlCaBtt+2adkCSVjCEjsbWciek9iJ7STeZFu2dmmW+/vj\n3jsajUbS3Jk7i0bv5+ORh6U7d+6cuRpbn3w+53xOxHEcRERERKQxRes9ABERERGZmYI1ERERkQam\nYE1ERESkgSlYExEREWlgCtZEREREGpiCNZEFxBgT0RgaUzXui+61SHOI13sAIgLGmJ8CF+QdcoAR\nYCfwn8CXrLVjBee3WmtfU+L124DPAfcD/zbHubuAe621v+l97wB/Za39k9LezazXfhvwVuC3ve+3\nAD8BLrHW3l7p9cNijLkQ+DvgBGCntfbUKr7WeuArwB8Au7xju8j7GZRwjU8Dfwa0WWvHil1TROYv\nBWsijeMJ4Pe8r6PAYmAL7i/hNxpjLrbWjnuP/wEQJGuyHvgg8J4Szn0rMBjg2kF8DBjL+/4h4LXA\nU1V6vXL9DZAE3gwcqfJrXeL9ly/oz+CrwO2A//kodk0RmacUrIk0jiFr7b0Fx75vjLkPuAX4Q+Cv\nAKy1T1ZrENbah6t17SKvNQAUvudGsBT4gbX2B/V48aA/A2vti8CLVRqOiNRZRDsYiNTfXGVNY8z9\nwHJr7YZi5xtjzgC+BJyOmxF6BPhLa+2teaVG325r7QZjzI3ABuCXwJXAEHAq8CjTy6B/D7QB7wTS\nwHeAj1lrD3rnfJq8MlzeuHf51/K+Xp83jtd7f04pgxpjVgB/DrwJWAk8DXzRWvv1gut+HTe7+B7c\n4OpB4CPW2vuK3cO8557hXf/V3nu6B/iUtfZeY8wG4PmCp7zHWnvjDNe6EvgAsBlIAM8B/2Ct/bu8\nczqAzwCXA8u8c75orf3XvPvmu8lae2XBfbsd2GCt3Vjw2p8Hft+7Rx/zrtMG/EnhNYG9wB8Dx1lr\nD+Rd43Lgv4DNxf4HwBizBPgyboa3x7s3XwO+YK115np/edeZ8Z57j2/wrv0R4L3AScBnrLWfNcas\nBj4LXAp04X62P5UfSM/2+S98TyLzkRYYiMwPtwPrjTFrCx8wxnQDdwBZ4F3Ar+OW7r5ljNmIW2r0\n5z5dj1ti850DvAJ4B/BRL9NVzAeAdcD/Aa4B3gLcYYwJkp1/K26p92Hc0udDRd7LcuAB4GLgOu+9\nPAzcbIz5eMHp/y/wGuD9wLuB1bjvOTHTAIwxFwD3Ad3e834HaAd+Zow5D9jnje0I8H3v6+/NcK33\n4wYuP8Qtl74Dd37YDV6AjDEmCtwGXAX8rXfeD4GvGWOuwC1fft675NtwA5pCN7mXMmfmvXYE92fx\n39bakYLzi13zq0DMu0/53osbFM6Uqf133Hv8h7g/k+941/69Et9fKfc83/XAPwC/AXzTCxa3487n\n/CPg7cALuBnnS73rz/X5F5n3VAYVmR/2e3+uwv1llW8T0Av8i7X2Nshl4v4MSFprB4wxj3jn7iwo\nscWB91lrn53j9Z8DLrXWpr3rHwT+F7gM+J9S3oC19mFjzBAwlpdRKTztI7iZopdZa6137HZjTCvw\nZ8aYrdZafw7ZqDemCe9aHbiBzdm4v+CL+Svc+3dh3vO+h5u9+4K19tXAvcaYFNBXpCyd7yTgb621\n1/oHjDHbgcO4WcOfAhcB5wHvstb+h3faj40x64BfsdbeZIzxM3kPW2t3FXmdbwLHcAORB71jFwDH\nATcXnmytfbHYNY0xPwGuwM1AYYxZg5u9fP8s7/EC4OvW2v/0vv+pMWYY8LNzs74/3J/HrPccN9vm\n+5a19iv+N8aYzwBrcD8Pz3iHv2eM+bH3Pr7PHJ//Wd6byLyhYE1k/nsc95fnV40xl+Bm4W6z1n6k\nhOcOlxCoAdziB2qeb+OWQy+gxGCtRK8HHswL1Hw342YHX4P7CxrgF/4vf48/Z6uj2IW9YO5s4PP5\nz7PWjhtj/gv4Y2NMp7V2qJSBWmv/2LtuN3AKbvB2lvewHyScj7uy938Lnvv2Ul7DO3fMG99vGGP+\n2Fqbxc2Q7QJ+Xup1gK3AfxhjTrfW/hI3cBvFXW08kx8DV3nB1/eBW621n8l7fNb3V+o9z3tq4Vy9\nNwJPAs8VZHG/Dfytt+q1ks+/yLygMqjI/HCc9+e0SeTW2mHgXOC/cef1fAM4YIz5b2PMsjmuW1Jg\nwmRmz3/NLG4GqafE55dqSeFrefZ5fy7OO1ZY/st6f87079pi3DluM10/gluqK4kx5nhvPtlRvDlY\nuHOqYHKl7jLgWN4q3nLdiJthusAY04Jbcr3ZnzdWom/i/syu8L6/Evgva+1sq07fhdvy5RTcMudO\nY8y9eSXZud5f0Hte+HlchlumTxX897fe42sq/PyLzAsK1kTmhzcCO6y1e4s9aK3dYa29ArccdBZu\neemteKtHQ7Ak/xsvy7EMOOgd8oOGWMHzugjmCG4ZtNBq789DAa+X7yjuOGe6voMbzMzJmzN2K+6C\niXOBDmvtZuDDRV5zkTEmWfD8jUXma83IWns38CzuXK5LcYOgaSXQOa4xjrso453GmLOBk4F/meM5\nA9baa6y1JwAn4rZ/OR53UQLM/f4qvedHcVcLv2qG/x7zxlntz79IXSlYE2lwxpg34/4C+vIMj7/B\nGHPQGHOGtTZrrX3QWnsNbnnIX32ZqXAYv1rQDf9y3MDsR973/sKE3AIIY8xmCoK8EsbxE+BMM30y\n22/j9hAru82Hl4G5H7g8fxGCl6l6J3BfgAxYL+7K2RuttXfnlfgu9f70/229Czd79JaC538et2kt\nlP6zucm7zm8C2621O2c5d6ZrfhU3SPocYK21M83twxiz2BjznDHmQ7gnP+etcv0PYK33eZj1/YVw\nz3+KG1Q+Z619wP8Pt1z+SSBb4udfZF7TnDWRxtFpjPFbd0RxS4wX4K56/AFu+4xi7gcmgH8zxlyH\nO3/nQtzy0VXeOf3en28wxjw1V3uLIjZ71/8X3CDlL3ADNX/Xge8CXwT+0Rjzl7itND7N9Iay/cAr\njTG/wvT5SQB/jRuY/cBra7EXNzD8Tdx2DTOtVi3VNbj38sfGmL/2jn0Et8R4ZakXsdYe9Cbx/773\nZx/uRPuP42aL/Hlz38Nd7LDVGLMSt/nvr+KumvwN7xz/Z/M2Y8z3rbVPz/CyX8dtkfF23NW5syl6\nTWvtE8aYe3CDnY/N8R6PGmMeB64zxqRxgx+DW0b9L2ut4y0UmOv9VXLPvwj8FvATr1XJXtxFER8D\n/tVaO+wtJpjr8y8yrymzJtI4NuPOfboH9xfg/8Xtb3UN8GsFE/xzvDlHbwQscANuG4O3An9grf2q\nd84h77G3ArfN1t5iBn+BG4R8G3eV3U3AZf6cKWvtDtxJ7ytwA7c/9c77RcF1/gZ3YcJtFOmw7/UA\ne533/r8AfAs4A7jCWnt9wDFPY639Ce4qxTTu3KYbgWHgfGvtTwNe7i24q2S/irvI4ldxW1rchjvx\nHmttBvd9fh34BG7pdAvwDmutX0q8A7gTt5fYF2cZ+x7czGOKyTLkTGa75ne9a5RSRv1t3Pv0MdyA\n61O47UreV+r7q+SeW2v347ZPeRx3ntr3cT/Dn8RbxVrK519kvlNTXBGRBcRrL/KStfad9R6LiJRG\nZVARkSZnjGnH3cHg5bjtT86u74hEJAgFayIizW8U+F3cnQM+aK19cI7zRaSBqAwqIiIi0sC0wEBE\nRESkgSlYExEREWlgTT1nra9vcEqNt6ennf7+wh1qpBy6l+HS/QyP7mV4dC/Do3sZnma+l729XZFi\nxxdUZi0eL9wJR8qlexku3c/w6F6GR/cyPLqX4VmI93JBBWsiIiIi842CNREREZEGVrM5a972NjcB\nG3A3Gb4Kd/uRG3G3sXkcuNpamzXGXIW7lUgauN5ae6sxpg13q5LlwCDu9jN9tRq/iIiISD3UMrN2\nKRC31r4OdzPivwC+BFxrrT0PiACXeZsBfxA4B7gI+KwxpgV34+LHvHNvBq6t4dhFRERE6qKWwdoz\nQNwYEwW6cTcSPhP4mff4bcCFuNugbLfWjltrjwE7gNOAc4HbC84VERERaWq1bN0xhFsCfRpYBvwa\ncL611m+vMQgswg3kjuU9r9hx/9isenrap60a6e3tKvsNyFS6l+HS/QyP7mV4dC/Do3sZnoV2L2sZ\nrH0YuMNae40xZi1wJ5DMe7wLOAoMeF/Pdtw/NqvCPiy9vV309Q2WO37Jo3sZLt3P8Ohehkf3Mjy6\nl+Fp5ns5UxBayzJoP5OZsSNAAnjYGLPFO3YJcBdwP3CeMabVGLMI2IS7+GA77ry3/HNFREREmlot\nM2t/DXzNGHMXbkbtE8ADwFZjTBJ4CrjFWpsxxtyAG4xFgU9aa8eMMV8BbjLGbAMmgHfVcOwiIiIi\ndVGzYM1aOwS8s8hDFxQ5dyuwteDYCHB5dUYnIiIi0pjUFFdERESkgSlYEwnRoaOj7N7fnBNfRUSk\nPhSsiYTopjssn/v3h8hmnblPFhERKYGCNZEQHT42xvhEhqHRVL2HIiIiTULBmkiIBkcmABjw/hQR\nEamUgjWRkGSyWUbG0gAMjiizJiIi4VCwJhKS4dE0ub3TlFkTEZGQKFgTCclg3jy1gWEFayIiEg4F\nayIhGcrLpg2oDCoiIiFRsCYSkvx5aiqDiohIWBSsiYREZVAREakGBWsiIckvgw6qz5qIiIREwZpI\nSKaUQZVZExGRkChYEwmJv2tBe0tcCwxERCQ0CtZEQuKXPlcv62B0PE0qna3ziEREpBkoWBMJyeDI\nBC2JGEsXtea+FxERqZSCNZGQDI2m6GxL0NWeALTllIiIhEPBmkgIHMdhcCRFV3uC7vYkoMyaiIiE\nQ8GaSAgmUllS6Syd7Qm6O9xgbUDBmoiIhEDBmkgI/CxaV1uSrja3DDowrDKoiIhUTsGaSAj8laBd\n7Qm6OlQGFRGR8ChYEwmBv5igsy1Bt7fAQGVQEREJg4I1kRAMjXpl0PYEXbkFBiqDiohI5RSsiYRg\nMrOWpDUZIxGPqgwqIiKhULAmEoKhvDlrkUiE7vaEFhiIiEgoFKyJhCC3GtSbr9bZnmRwZALHceo5\nLBERaQIK1kRC4JdB/flq3e1JJtJZxlOZeg5LRESagII1kRAMjaaIRKC9NQ6QtyJUpVAREamMgjWR\nEAyOuPuCRiMRgMlea8NaZCAiIpVRsCYSAn8Td1+32neIiEhI4rV6IWPMlcCV3retwOnAucDfAA7w\nOHC1tTZrjLkKeD+QBq631t5qjGkDvgEsBwaBK6y1fbUav8hMslmH4dEUq5d15I51qTGuiIiEpGaZ\nNWvtjdbaLdbaLcCDwAeBPwWutdaeB0SAy4wxK73HzgEuAj5rjGkBPgA85p17M3BtrcYuMpuhsRQO\n5PYEBfIa4ypYExGRytS8DGqMOQvYbK39Z+BM4GfeQ7cBFwJnA9uttePW2mPADuA03Czc7QXnitTd\n5ErQvDJohzZzFxGRcNSsDJrnE8B13tcRa63fiGoQWAR0A8fyzi923D82q56eduLx2JRjvb1dZQ9c\nptK9dB0YGAdg+bLOyXsSd/9qTWScku+T7md4dC/Do3sZHt3L8Cy0e1nTYM0Ysxgw1tqfeIeyeQ93\nAUeBAe/r2Y77x2bV3z8y5fve3i76+gbLGrtMpXs56YW97v9DxBwnd09Sabe/Wl//cEn3SfczPLqX\n4dG9DI/uZXia+V7OFITWugx6PvDjvO8fNsZs8b6+BLgLuB84zxjTaoxZBGzCXXywHbi04FyRuvO3\nmurMK4Mm4jHaWmJaDSoiIhWrdbBmgOfyvv8ocJ0x5h4gCdxird0P3IAbjN0JfNJaOwZ8BdhsjNkG\nvI/JUqpIXRVuNeXrak9qNaiIiFSspmVQa+3/V/D9M8AFRc7bCmwtODYCXF7VAYqUYdDfxL0tOeV4\nV3uCw/vGyDpOrlmuiIhIUGqKK1KhoSKrQcFtjJvJOoyMpesxLBERaRIK1kQq5GfWOtqml0FBvdZE\nRKQyCtZEKjQ0kiKZiNKSmNomZrLXmoI1EREpn4I1kQoNjk5Mm68G+Zk1rQgVEZHyKVgTqdDQSGpK\n2w5ft8qgIiISAgVrIhUYn8gwkc5OW1wA0J3bzF2ZNRERKZ+CNZEKDI56PdbapgdrfhlUvdZERKQS\nCtZEKjC5iXuROWsdXhlUCwxERKQCCtakLMOjKu1B3lZTRTJrnW1xIqgMKiIilVGwJoE98PRB/s+n\nvs/z+wbqPZS6m2mrKYBYNEpHW0ILDEREpCIK1iSwvmOjOA4cPjZW76HUnb97QWeR1h0A3R1Jte4Q\nEZGKKFiTwLJZB4CM9+dCltsXtEhmDdwVoUOjKTLZbC2HJSIiTUTBmgSWyfjBmgKQwRn2BfV1egsP\nhpRdExGRMilYk8D8jJoftC1kk3PWipdBO1vjAAxpM3cRESmTgjUJLOuoDOobGk0RiUC7F5QVSnr7\nhU6kMrUcloiINBEFaxJYRnPWcgZHUnS2JYhGIkUfV7AmIiKVUrAmgeXmrGUW9py18YkMRwbHWNRR\nvAQK0JJw/4pNpBf2vRIRkfIpWJPAtBrU9YA9yEQqyxkn9854jjJrIiJSKQVrEljGm7OWXuDB2rZH\n9wFwzmmrZjynJResKbMmIiLlUbAmgWW9lh0LuQx6oH8E+8JRNq5bzPLFbTOel4y7f8XG08qsiYhI\neRSsSWCTfdYWbmZt+2NuVu2801bPel6uDDqhYE1ERMqjYE0Cyyzw1h3ZrMP2x/bT1hLjlWbm+WoA\nyYSfWVu4WUgREamMgjUJLLvAm+I+sesI/YPjnL1pRW5O2kxatMBAREQqpGBNApvss7Yws0X+woJz\nZ1lY4EvGtcBAREQqo2BNAlvIc9aGRlM8/Gwfq5d1cMKq7jnPT+b6rCmzJiIi5VGwJoHltptagGXQ\ne5/YTzrjcO7LVxGZYdeCfCqDiohIpRSsSWALuQy67dF9xKIRXvuylSWd768GHVcZVEREyqRgTQJb\nqDsYHBsaZ8/BIU7dsGTWLaby+X3WlFkTEZFyKViTwPxmuAutDHro2BgAq5e1l/ycRDxKBAVrIiJS\nPgVrEthC7bN2eMAN1pZ0t5b8nEgkQjIRU581EREpW7yWL2aMuQZ4C5AE/gH4GXAj4ACPA1dba7PG\nmKuA9wNp4Hpr7a3GmDbgG8ByYBC4wlrbV8vxi8svg6YX2HZTfrC2LECwBu6KUGXWRESkXDXLrBlj\ntgCvA84BLgDWAl8CrrXWngdEgMuMMSuBD3rnXQR81hjTAnwAeMw792bg2lqNXabKLNA5a0eOjQPB\nMmvg9lpTnzURESlXLcugFwGPAd8EvgvcCpyJm10DuA24EDgb2G6tHbfWHgN2AKcB5wK3F5wrdZAL\n1hZoZm3pouCZtXFl1kREpEy1LIMuA9YDvwYcD3wHiFpr/fTMILAI6AaO5T2v2HH/2Kx6etqJx6du\nB9Tb21X+OxCAXH+xSCxat/s5MpbiT//pHt72+pN43RybqYfl6PAErckYG9b2lNRjzdfRluDI4Pic\n90qfzfDoXoZH9zI8upfhWWj3spbB2mHgaWvtBGCNMWO4pVBfF3AUGPC+nu24f2xW/f0jU77v7e2i\nr2+w3PGLx59/NTaertv93PHSMeyefu566EVOXlWbv7QHj4ywpLuVQ4eGAj0vGokwPpHh4MGBGYM8\nfTbDo3sZHt3L8OhehqeZ7+VMQWgty6DbgIuNMRFjzGqgA/ixN5cN4BLgLuB+4DxjTKsxZhGwCXfx\nwXbg0oJzpQ4aYSP3kbEU4G7/VAuj42lGxtMs6W4J/Fx/y6mUVoSKiEgZahasWWtvBR7GDca+C1wN\nfBS4zhhzD+4K0VustfuBG3CDsTuBT1prx4CvAJuNMduA9wHX1WrsMlUj7GAwPJoGYGisNsFauStB\nAVr8zdwVrImISBlq2rrDWvuxIocvKHLeVmBrwbER4PIqDU0CaITVoMNekDZco8za4WPBe6z5cpu5\npzLQlgh1XCIi0vzUFFcCa4Qy6PCYl1mrUbB2pMyVoJC/P6hWhIqISHAK1iSwhiiD5jJraRyn+kHj\nIT9YKyez5pdB1WtNRETKoGBNAvODtLqWQb05a1nHYXQ8XfXXOzLgNsQtK1jzyqDKrImISDkUrElg\njVAGHclbWFCLUujhY2NEIxEWdyUDP7cl4S8wULAmIiLBKViTwBpjgcFkNm1otPqZtcMDY/R0JYlF\ng/+V8eesqQwqIiLlULAmgWQdB3+KWCPMWYPqZ9bSmSxHh8bLWgkKBatBRUREAlKwJoFk87JpjbAa\nFKrfvuPo4DiOU95KUFCfNRERqYyCNQkkv/TpMDV4qxXHcWo6Z+1wBStBIW+BwYQyayIiEpyCNQmk\nMDirRyl0IpUlnXFIxt2Pb+MHa1pgICIi5VOwJoEULipI16EU6s9XW97TBlR/y6lKdi+AydWg41pg\nICIiZVCwJoEUBmv1WBHqz1db0dPufl/1zJrXY63MOWtaYCAiIpVQsCaBTC+D1iFYGy3IrNWsDNpS\n1vOTWmAgIiIVULAmgRTOUctkah+A+Jm1RR1JWhKxqgdrRwbG6GiN05qMl/V8ZdZERKQSCtYkkIbI\nrHlz1DraEnS2xataBnUch8PHxsougYI2chcRkcooWJNAGmHO2oiXWWtvjdPRlqjqDgZDoykm0tmy\nV4JCXp81LTAQEZEyKFiTQKYFa3Upg3qZtdYEnW0JxlMZUlWaD1Zp2w6AhMqgIiJSAQVrEkhjlEHd\nTFpHa5zOtgRQvUUGh4+5K0HLbdsBEI1ESMaj6rMmIiJlUbAmgTRCGdSfo9bRlqDDC9aqNW/Nz6wt\nq2DOGrjz1lQGFRGRcihYk0Cml0HrMWfNL4PG6WytbmbtyEBlDXF9yURUCwxERKQsCtYkkEbYbmpo\nLE0yHiURj9WgDOrNWas0sxaPac6aiIiURcGaBOJn1mLRCADpuqwGTeXKn7lgrUpbTh0aGCMei9LV\nnqjoOslElHE1xRURkTIoWJNA/Mya3zusHmXQ4dE07a1ug9pqz1k7MjDG0u4WopFIRddpSbiZNcep\n/f0SEZH5TcGaBOKXPf3NyWtdBs1mHUbH03S0FmTWqhCsjacyDI6kKp6vBm5w6zj12fheRETmNwVr\nEohfBk0m65NZGxlP4+AuLgDobHP/rEaw5i8uqHS+GkAy7vVaU/sOEREJSMGaBOKXQVu8Rq+1bt0x\nktcQFyYza8NV2MXgub0DAKxa0l7xtfxMpNp3iIhIUArWJJBM4Zy1GpdBh/O2mgJoa4kTjUSqkll7\n0PYBcMYpvRVfy9/MXe07REQkKAVrEkguWIvXpwyav4k7QCQSoaMtHnqwNjKW5vHnD3NcbwcrQ8is\nJXP7gypYExGRYBSsSSC5Mqg/Z63GZVC/3OnPWQO3FBp2sPbIzkOkMw5nmeWhXM+/XyqDiohIUArW\nJJBMbs6aG3yka7yRe+GcNXCzbMNjKbIhtsV44OmDAJy1MZxgzV9gMK4FBiIiEpCCNQlkWp+1GmfW\nhsaKZNZaEzgOjI6Hs8hgdDzNY88dYfWyDlYv6wjlmsmEyqAiIlKe+NynhMcY8xAw4H37PPAXwI2A\nAzwOXG2tzRpjrgLeD6SB6621txpj2oBvAMuBQeAKa70Z4FIzkwsM6rwatG0ys5bfay0/41auR3ce\nJp3JcpapfGGBL6nVoCIiUqaaZdaMMa1AxFq7xfvvPcCXgGuttecBEeAyY8xK4IPAOcBFwGeNMS3A\nB4DHvHNvBq6t1dhlUmEZNFPjMqg/Z629YM4ahNdr7QEbbgkU8sqgyqyJiEhAtcysvQJoN8b8wHvd\nTwBnAj/zHr8NeBOQAbZba8eBcWPMDuA04Fzg83nnfqqGYxdPvcugw0XnrLkf4zC2nBqfyPDYzsOs\nXNLOmpBKoJDfZ03BmoiIBFPLYG0E+ALwVeBk3IArYq31f9sPAouAbuBY3vOKHfePzaqnp5241zLB\n19vbVf47EFq9LJYffLS0Jmp6TycyDpEIrDuuJ7eZ/Krl7utH4vGKx7LtkZeYSGe54JXHsXx5d8Xj\n9fUeGgEg0TLz/dJnMzy6l+HRvQyP7mV4Ftq9rGWw9gywwwvOnjHGHMbNrPm6gKO4c9q65jjuH5tV\nf//IlO97e7vo6xssd/wCDAy6WzD5mbXBwfGa3tNjg2O0JeMcOTyUO+Z4Kyz3HRyseCx33r8HgE1r\nF4X6vsZGxgHoPzpS9Lr6bIZH9zI8upfh0b0MTzPfy5mC0FquBn0v8EUAY8xq3EzZD4wxW7zHLwHu\nAu4HzjPGtBpjFgGbcBcfbAcuLThXaqywz1q6DjsY+GVPX1hz1sZTGR7ZeYjlPW2sXd5Z0bUKaYGB\niIiUq5aZtX8BbjTGbMNd/fle4BCw1RiTBJ4CbrHWZowxN+AGY1Hgk9baMWPMV4CbvOdPAO+q4djF\nM327qdrPWVu1dOpcso7c/qCVBWtPPH+EiVSWV21cTiQSqehahdRnTUREyhUoWPPaZ2SstRPGmE3A\nrwL3Wmu3zfVca+1MAdYFRc7dCmwtODYCXB5kvBK+aRu513C7qVQ6y0QqS2drdTJrO/e6UyJP3bCk\nousUk8usTShYExGRYEougxpjXg/sBc41xqzBLVdeA/zEGPPuKo1PGkw9N3L3e6y1F/RSCytY27Pf\nnQOxfkW4JVCYvF/jaZVBRUQkmCBz1v4S+A/gXuC3gCPAatz+Z38S/tCkEU3rs1bDMmhu94K2qcFa\nPBalJRmrqAzqOA67DwzRu7h1WjAYBr8MqtYdIiISVJBg7XTg81458mLgVq8X2g+Ak6oxOGk80/qs\n1bAMOrkv6PTqfWdrgqGx8oO1/sFxhkZTrF9RneXg6rMmIiLlChKsHQO6jDHdwOtwgzSA44HDYQ9M\nGlM9M2v+7gXFtpTqbEtUVAbd7ZVA11UpWItGI8RjUSZUBhURkYCCLDC4Dfhn3Ia0g8AdxpgLgS8D\n363C2KSY7NkSAAAgAElEQVQB+XPUknXYbmo4N2etSGatLc5EKksqnSFR0Ai5FLsPePPVVlav0WJL\nIqrMmoiIBBYks3Y1cDfuTgSXWWvHgNcC24CPVmFs0oAK+6zVNLM2NnNmrSO3yCBd1rX3HHCb7FYr\nswZugKu9QUVEJKiSM2veXLWPFhz789BHJA0ttxo0XvtgbdY5a3krQnu6WgJfe/eBQRZ3JlnUkaxs\nkLNIxqOMqXWHiIgEFLTP2vm47To2AluA9wA7rbVfD39o0oj84CwR9/us1bAMOlp8NShU1r5jYHiC\n/sFxXnHi0soGOIdkIsbAyERVX0NERJpPkD5rlwC3Ay8AK4EY7k4EXzPGvKc6w5NG45dB47EIsWik\ntmXQ8Zkza5XsYrCnBvPVwF2Uoe2mREQkqCBz1j4N/JG19n1AGsBaex1uafSPwh+aNCI/OItGI8Ri\nEdINtBoUysus+YsLqjlfDSCZiJLJOqRrmI0UEZH5L0iwthk3s1bou8AJ4QxHGt1kZi1KLBqteZ+1\neCxCMjH9Y+sHcGUFa7mdC6ocrMW1mbuIiAQXJFg7RPGg7CzgQDjDkUbnZ9ZiUb8MWrvAY2gsTXtr\nougm65Vk1vYcGKKzLcGS7uALE4Lwg0ytCK2ObNbha997ioftwXoPRUQkVEGCtX8GvuzNXYsAJxpj\n3gv8PXBjFcYmDaiwDFrr1aDF5quB22cNgs9ZGxlLcfDoKOtWdBYNAsOU28w9rWCtGg4dG2XbY/u4\n88EX6j0UEZFQBVkN+llgEfC/QAtwB5ACvgSohccCkc06xKIRIpEI8WikZmVQx3EYHk2zoqe96ON+\nZu3I4Hig6/r91apdAgVoURm0qlLe7hBp7RIhIk0mSJ81B/i4MeYzwCZgAnjGa44rC0Qm6xCNuhmo\nWDRasyzR2ESGrOPMmFlra4mzfmUXT+3u58ldRzh1w5KSrlurlaAAyaQ2c68mfyuvlII1EWkyQVp3\ndBhjvg582Fr7gLX2UcAaY/7VGNNWvSFKI8lks5PBWg3LoJNbTU1fCQoQiUS48uKNRCMRbrztacZL\nbD5bq5WgkJ9ZU7BWDX6QltJqWxFpMkHmrP0N8ErgR3nH3gecDfxVmIOSxpXNOsQifmatdmXQyYa4\nMyeD16/s4qJXr+XQsTG+eddzJV1394EhWpIxlvdU//83/Dlr48r8VIUfpKkMKiLNJkiw9hbgSmvt\nvf4Ba+0dwO8Bl4c9MGlMmaxDLDZZBq1VZm1yq6nimTXfZeccz4qeNn74wAvs3Hts1nPHUxn2HR5m\n/fJOolVeXACTq0GVWauOlMqgItKkggRrrcBokeMDQPVrSNIQsvlz1mK1a93hb+LePsOcNV8yEePK\nSzbiOHDj95+etQHtiweHcJzalEBhss+aWndUh8qgItKsggRrPwf+3BjT4R8wxrQDfwZsC3tg0pgy\n3mpQqG0Z1O+f1llkX9BCZl0PW85Yw0uHhvnePbtnPG93DRcXQH5mTcFENaS8xS4qg4pIswnSuuPD\nwM+Al4wxT3vHDDAIXBT2wKQxFQZrDlOzbdUSJFgDuHzLifzy2T5uu283b3rVWtpapn/UH3/uCAAn\nrO4Ob6CzaFGftaqaLIPq/opIcyk5s2at3QGcCnwc+AVwN/AxYKO19qnqDE8aTTbrEIu5Hxv/z1qU\nQoMGa20tcV5/xhomUlnue2r6Bhv9g+M8svMQG1Z2sWppR5ErhC/XFFeZtaqYLIPWrlGziEgtBMms\nYa09BvxTlcYi80BhZg0gnXFIBPokBTccMFgDOPe01Xxr2/Pc9chetpy+Zspj2x/bh+PA+a9YHeo4\nZ6MFBtU12RRX91dEmkvJv2K9uWofAl4LJHG3nMqx1r4p3KFJI8pkndzKST9Yq8WK0KCZNYCerhZO\nO2Epj+w8zJ4Dg7mFBFnHYduj+0jGo5y9aUVVxltMixYYVJVWg4pIswq6N+g1uFtM7QVeKvhPFoBs\nfmYtVwatQbA2liIWjdCajAV63vmnu5mzux7Zlztm9xzl4NFRXrVx+ZyrS8OkBQbV5a8CVbAmIs0m\nyG+qtwCXW2tvr9ZgpPHl91mL+5m1GrRKGBpN09GWCLzZ+mknLmVRZ5J7ntjP5a8/kWQixl2P7AXg\nvBqWQEEbuVebHwTP1q5FRGQ+CpJZGwd2VGsgMj9M6bPmz1mrQWZteDQVqATqi0WjnPvyVYyMp3nQ\n9jE8luIB28fKJe2cfNyiKox0Zi1aYFBVuR0MMg5ZR4sMRKR5BAnW/g34kDGm+q3epSE5jvtLMLfd\nVKw2mbVs1nGDtTJLluedtgqAnz+yl3ufOEA6k+W8V6wKnKWrVCLu/nXTnLXqyG/ZUYtsr4hIrQT5\n7dcB/BbwVmPMTtxMW44WGDQ/f27aZGatNnPWRsbTOEBHGZk1gOU97Wxa38NTu/vpOzZKLBrhdS9b\nFe4gSxCPRYlFIyqDVkn+XLVUuvorlEVEaiXIP2cx4D+qNRBpfFkvKMvtDZrLrFU3WCtnJWih816x\niqd293NkYJwzT+llUUcyrOEFkkzEGJ9Q1qca8oO1dI22QRMRqYWSgzVr7XsqfTFjzHLgQeCNQBq4\nEXCAx4GrrbVZY8xVwPu9x6+31t5qjGkDvgEsx90x4QprbV+l45Fg/AyaXwaN1yizFkawduYpvXS0\nxhkeS3PeK2qfVfMlE1Fl1qpkSrCmFaEi0kSCzFnDGNNjjPm4MeZfjTHLjTHvMMZsLPG5CdyGuv5m\n8F8CrrXWnofbs+0yY8xK4IPAObhbWH3WGNMCfAB4zDv3ZuDaIOOWcEwrg/qZtSpnMcII1hLxGG/f\nciKvPnUFLzt+aVhDC6wlHlNT3CqZEqxpzpqINJGSgzVjzCnA08B7gXcDncDbgQeMMa8r4RJfAP4R\nt0cbwJm4e40C3AZcCJwNbLfWjnu7JewATgPOBW4vOFdqLFcGLVgNWu0yqL97Qblz1nxbTl/D+9+y\nuer7mM4mmYhqNWiVpPICNG05JSLNJMictb8GbrHWXm2MGfSOvRv4KvA54PyZnmiMuRLos9beYYy5\nxjscsdb6/6IOAouAbuBY3lOLHfePzamnp514fGoT1d7erlKeKkVEk25StL3dne/V3dUKQGdXa3Xv\na8zd23PNyu55//PraE+y/8ho0fcx399bveVX47u723Q/Q6L7GB7dy/AstHsZJFh7DfDR/APeHLPP\n4c5Dm817AccYcyFwOm4pc3ne413AUWDA+3q24/6xOfX3j0z5vre3i76+wRnOlrkcPjYGQGrCLeON\njbkZryP9w1W9r/v7hgDITKTn/c8vilui23/gWG41LeizGYbR8XTu64OHBulKBprlIUXocxke3cvw\nNPO9nCkIDfKvmQO0FTm+nII2HoWstedbay+w1m4Bfgn8DnCbMWaLd8olwF3A/cB5xphWY8wiYBPu\n4oPtwKUF50qN+XPT6lUGrWTOWqNIxrXlVLXkb+CuBQYi0kyCBGvfAa43xnR63zvGmBOAvwG+V8Zr\nfxS4zhhzD+7G8LdYa/cDN+AGY3cCn7TWjgFfATYbY7YB7wOuK+P1pELT+6zVZiP3MBYYNIrcllNa\nZBC6qQsMNGdNRJpHkDLoR3An9x/xnnc/sAS4D/ijUi/iZdd8FxR5fCuwteDYCHB5gLFKFUzvs+bG\n+tXuaTWUW2Aw/7uc+pu5jyvzE7oJrQYVkSYV5LdfEngd8AbceWcTwBPW2h9XY2DSeAr7rNWqDDo0\nmqKtJT5ljtd8pcxa9ah1h4g0qyDB2gPA26y1PwJ+VKXxSAOrZxm0swmyaqDN3Kslm3WmfA5TCtZE\npIkESVVEmGMhgTS3GbebqmKw5jgOQ6PpppivBvkLDJRZC1OqoKycTmvOmog0jyDpiq8Btxtj/hV4\nnsmdCACw1v57mAOTxpMpaIqb226qilmMiVSWdCZbcUPcRuFn1sYUrIWqMJOmvUFFpJkECdY+5f35\niSKPOYCCtSbnZ9aihXPWqphZa6aVoAA93S0AHOwfneNMCcLPrMVjUdKZrFp3iEhTCbKR+/yf3S0V\nKcys1aIMmgvWWpsjWFu/wm14uHt/czZ0rJeU12OtrSXG4EhWrTtEpKkoAJOS5YI1r2VHrAZl0KGx\n5sqsrVjSTksyxp6DCtbC5LftaGtx//9Tq0FFpJmUnFkzxqRwy51FWWuToYxIGlY9yqC53QvamyNY\ni0YirF3eyc6XjjGeyuTmsEll/DJou4I1EWlCQeasXcXUYC0OnAJcQYCmuDJ/1aMMOjjSXJk1cEuh\nO148xot9Q5y4elG9h9MUUgWZNbXuEJFmEmTO2o3FjhtjHgJ+F/hGSGOSBuXvDTrZZ80vg1Y/s9Ys\nq0EB1q1wd2zbc0DBWlj84CxXBlXrDhFpImHMWbsXODeE60iDm9ZnzQvaqtkmodkWGIAWGVRDKuUH\na25ZWa07RKSZVBSsGWNagT8A9oczHGlk07abqsVq0CZbYACwelkH8ViEPQcUrIXFz6y1t7ifE7Xu\nEJFmUukCg5h37PfDHJQ0puxM201VsQzabH3WwO0Ftqa3kxf7hkhnssRjWpRdqfzWHaAFBiLSXCpZ\nYADuZu73WmufD29I0qimLzDw5qxVseQ0PJoiHouSTDRXQLN+RRe79w+y7/AIa5d31ns48970BQaa\nsyYizaPk34DeAoNbgbuttTdZa28CMsBQlcYmDaawz1q8RjsYdLbFiXil12ax3ltkoHlr4Sjss1bN\n3n8iIrVWcrBmjHk18Cxuhs3358DjxpjTwx6YNJ4Z+6xVtQzaPJu451u30ltkoHlroUgX9FlT6w4R\naSZBaktfxN3/85q8YxuB/wG+FOagpDHNXAatTrCWyWYZHW/OYO243k4iEbTIICSFZVAtMBCRZhIk\nWDsd+JK1NuMfsNY6uIHaq8IemDSe6X3W/MxadX4xDo+mgebqseZrScRYvbSDPQeHyDqaX1Upvwza\nkogRjUa0N6iINJUgwdoR4NQix08GlB5YAGbqs1atzNqgtxK0qwmDNYB1K7oYn8hwsH+03kOZ9/zM\nWiIeJRGPajWoiDSVIKtBbwb+yRjzJ8AvvGNnAX8B/FvYA5PGU+s+a824e0G+9Ss6uecJd5HBy82K\neg9nXvPnqCXiURIxBWsi0lyCBGufBpYCW4EEEAHSwJeBT4Y+Mmk4frnOL4P6Cw2qVQZtxh5r+dZ7\niww0b61yfp+1RDxKPB5V6w4RaSpB9gZNAx8wxvwxYIAUsMNaO1KtwUlj8Vd9+uXPSCRCLBqpWmat\n2YO1tcu1IjQs08qgWmAgIk0kyA4GMdzs2l5r7Ve8Y78wxtwKfMZbbCBNLFMwZ83/Oq0yaFnaW+Ms\nX9zG7v2DOFpkUJEpwVosynAqVecRiYiEJ8gCg78EfhfYnXdsK/A+4M/CHJQ0psI+awCxaLRqfdaa\nPbMGbr+14bE0fUe1yKASfrCWVGZNRJpQkGDt3cC7rLXf9w9Ya/8ZuBJ4T8jjkgaUcaaWQf2vq7Xd\n1EII1vydDHa+eKzOI5nf/GAtHnPnrKl1h4g0kyDB2mJgf5Hje4DecIYjjczPoEWjU8ugyqyVz19k\n8Mye/jqPZH6bSGeJx6JEIhGtBhWRphMkWLsf+ENjTOEmjf8P8FB4Q5JGlS3YGxTc/UGrlVkbHk0R\nYXILoWZ00ppFxKIRHttxqN5DmddS6SyJuPu5TMRjZLKOmg2LSNMI8lvwT4A7gTcYYx70jp0BrAIu\nDntg0ngK+6yBO2dtPJ2Z6SkVGRxN0d4an5LJazatyTjHr+rm2RePMjqezm2XJMGkMlmSuWDN2wYt\nkyUaj9VzWCIioSg5s2atvR94GfDfQDuQBG4BNlpr767O8KSRFPZZg+qWQYdHU01dAvVtXL+YbNbh\nmReO1nso81Y6nckFaXEv85tKK7MmIs2hpGDNGHOuMeb/Aj8HPoy7T2gKuMNau7eK45MG4je/ndK6\no0p91hzHYXisOTdxL7RpXQ8AT+3WvLVyTUwpg7p/at6aiDSLOWsuxphrgeuAZ4DvAYdxFxucD9xp\njPm0tfb6Eq4Tw231YQAH+H1gDLjR+/5x4GprbdYYcxXwftwdEq631t5qjGkDvgEsx92L9AprbV+w\ntyuVyJVBowWtO6owZ210PEMm6yyIYO3ENYtIxKM8rWCtbKl0lkRMwZqINKdZgzVjzMW4PdSustZ+\nrcjj7wb+xRhzn7X2h3O81psBrLXnGGO24O4pGgGutdb+1Bjzj8Blxph7gA/i7jvaCmwzxvwQ+ADw\nmLX208aY3wSuBT4U4L1KhYr2WatSGXRorPlXgvqSiRgb1y/h8Z2HGFogpd+wpdJZEgkFayLSnOYq\ng34Y+EKxQA3AWvtvwOe882Zlrf0WbgNdgPXAUeBM4GfesduAC4Gzge3W2nFr7TFgB3AacC5we8G5\nUkMz91lzQu/A3+y7FxQ67eRlOIBVC4/AslmHTNbJZdZyc9bUa01EmsRcZdBXAh+d45z/BK4u5cWs\ntWljzE3AW4F3AG/M26ZqEFgEdAP5HUKLHfePzaqnp514wWqw3t6uUoYqRfgtO1au6Abce9nW6gZT\nS5d2TmnpUak9h90tZ1cs61wQP7OXn7gMgF0Hh7n43OZ/v2EaG08D0NGepLe3K5dZ6+pqXRCfnWrT\nPQyP7mV4Ftq9nCtYaweG5jhnFHdlaEmstVcYYz4O3Ae05T3UhZttG/C+nu24f2xW/f1T95jv7e2i\nr0+bZpdrbMz9pXj48BDLl3fT1zeYW3Sw/8AAyUR4bRJe2j8AQCSbXRA/s1PW9ZBMRHnYHlwQ7zdM\nfvNkJ+vQ1zeYC9b6Dg/R3aLWHZXQv5nh0b0MTzPfy5mC0LlSIc8Ar5/jnC3A03MNwBjz28aYa7xv\nR4As8IA3fw3gEuAu3Oa75xljWo0xi4BNuIsPtgOXFpwrNZRxHKKRCJHI1DIoEPqK0IWwe0G+RDzK\nKcctZu+hYY4Njdd7OPNK/r6gAHF/zpr2BxWRJjFXsPYN4HpjzPpiDxpjTsRdKPDVEl7rf4EzjDE/\nB+4A/hC3fHqdt6ggCdxird0P3IAbjN0JfNJaOwZ8BdhsjNmGO/ftuhJeU0KUzTrTGtRWK1g7OugG\nLF3tCyNYA9i43mvhoXlrgUx4TZn9IM2fu6b9QUWkWcxVBr0B+DXgl8aYrwH34Lbu6Mad8H8V8FNr\n7da5XshaOwy8s8hDFxQ5dytum4/8YyPA5XO9jlRPJuNM6bEGk/PYMiGvvNtzwE1xr13eGep1G9km\nL1h7evdRXnPqyjqPZv7wM2uFfdZSWg0qIk1i1mDNWpsyxlyE277j95m66vMA8FfA56s3PGkkmawz\nZaspcPcG9R8Li+M47No/yPLFbbS3LpzM2roVnbS1xNRvLaBcsBabWgYN+38gRETqZc6muNbaCeCT\nXnNcA/TgZtd2WGv1r+ECknWKlEG9TFs6xGDt0LExhsfSnLphSWjXnA9i0ShmbQ+/3HGIw8fGWLqo\nNfeY4zgcHZpg/+Fh9h8ZYXFXC6eftGzK/MGFKjdnLTG1DKrMmog0i5J3jfZabMy5kECaVyaTndJj\nDdwAw38sLLv3uyXQDSsX1tJscOet/XLHIb579y662hMc7B/l4NFRDhwZYWwiM+XcM05exhWXbKS7\nveTF2E2pMLOWa4qrvUFFpEmUHKyJZLLF5qyFXwbd5QVr6xdgsHaqN2/t549Mbrkbj0VZ0dPGqqXt\nrFzazoqedrY/to+Hnz3Ezpfu48pLN3H6ScsYGUvx6M7DPPTsIfr6R7n6rS9j2eK2mV6qaUzOWXPb\ndPi9FbWDgYg0CwVrUrKs17ojX241aIgr73Z7PdYWYrB23PJO3v+WzUykMizvaaN3cRuLu1qm3ffX\nvmwlP/zFC/zPz3Zywy2Psn5FFy/2DU0Jmr951/Nc9eZTa/0Was4vd2ojdxFpVgrWpGSZrJPbf9EX\n98ugIWXW8hcXdCygxQX5Xn3qijnPiUYiXHT2OjYfv4St332S3QcG2bCyizNO6eWMk5bxz999knuf\n2M+lr13PmmUdNRh1/aS81h0Jte4QkSalYE1Klsk4RFtmKoOGk8U4vEAXF5TruN5OPv2eVzE2kaGt\nZfKv81vPP56/+5/H+Pa25/mDX39ZHUdYfRMFrTviat0hIk0mvM0cpells05uQYEv7DLorgW8uKBc\nkUhkSqAGcPpJyzh+VRcPPH0wt2CjWc28wEDBmog0BwVrUrKM4xRZDRruAoPdBxbu4oIwRSIR3nr+\nCQB8667n6jya6koXtu7wg7WQsr0iIvWmYE1KVnS7KX8Hg5B+MS7klaBh27xhCacct4hHdh5m50vH\n6j2cqpkobIobU+sOEWkuCtakZJnMLJm1EMqgjuOwa98AvYtbF+zigjDlZ9e+2cTZtcLWHVoNKiLN\nRsGalMRxHLJVLoP6iwvWr+yu+FriMut62Hz8Ep7c1Y9t0g3itTeoiDQ7BWtSkqzjBmMzlUHDmB+k\nxQXV8evnHg/ArffsrvNIqmNan7WYm2HT3qAi0iwUrElJsl7mrJplUH9xgYK1cJ24ZhEb1y3mieeP\nNOXK0FRqap+1eNz9TKbUZ01EmoSCNSmJ32B0WmYtxDKoFhdUzyWvWQ/Abfc1X3Zt+g4G3nZTat0h\nIk1CwZqUxC+DTsushbQ3qOM47N4/qMUFVfKy45ewdnknv3j6IAePjtZ7OKGaac6aFhiISLNQsCYl\nycxQBs1tN1XhL8bDA2MMjaa0uKBKIpEIl7x6HY4Dd9y3p97DCZUfrCXjBa07FKyJSJNQsCYl8ees\nVasMuluLC6ruVZuWs2xRK9se28fA8ES9hxMav8+aH6TFohGikYj2BhWRpqFgTUriLyCYtt1USGVQ\nzVervlg0ykVnryOVzvKjB1+o93BCk0pniceiRCKT/yMRj0fUukNEmoaCNSlJZqY5ayGVQXftGwBg\n/QoFa9V07mmr6GxLcOeDLzE6nq73cEKRSmdz89R8iVhUZVARaRoK1qQk1SyDZh2H5/YNsGJJO51t\nWlxQTS2JGBeedRwj42m2Pbqv3sMJRSqTzc1X88VjUZVBRaRpKFiTksy0wCCMMui+wyOMjmc4cbUW\nF9TCltPXAPDLHYfqPJJwpNKZaZm1eCyi1h0i0jQUrElJ/DLnzGXQ8oO15/a6m4wrWKuN7o4kx/V2\nsOOlY7mVlPNZsTJoXGVQEWkiCtakJDNuN+V9X8kvxuf2uvPVTli9qOxrSDAb1/WQSmdzgfJ8VjRY\niytYE5HmoWBNSlLNMujOlwZIxqOs6e0of4ASyMb1PQA8vedonUdSuZkza5qzJiLNQcGalGTuBQbl\nZTHGJtK8dGiIDSu7cn2ypPrMusVEgKd399d7KBXJZh0yWYdErMicNWXWRKRJ6LejlGSyz1phZq2y\nOWu79g3iOCqB1lpHa4K1KzrZufcYE95G6PPR5FZTsSnHE7EomayTK9+LiMxnCtakJDP1WYtX2Lpj\npzdn6gQtLqi5jet6SGccdnpzBucjv/FtsdYdoM3cRaQ5KFiTksxVBk2XGaz5iwtOXKPMWq1tXOfN\nW5vHpVA/K1hszhqgeWsi0hQUrElJJhcYFG43Vf4OBo7j8NzeAXq6Wujpaql8kBLIKWsXE4nA03vm\nb7DmZ9biRfqsgTZzF5HmEK/FixhjEsDXgA1AC3A98CRwI+AAjwNXW2uzxpirgPcDaeB6a+2txpg2\n4BvAcmAQuMJa21eLsYtrxjlrFZRBDw+McWx4gjNNb+UDlMDaW+OsX9HFc3sHGE9laEnE5n5Sg5mc\nsza9dQcoWBOR5lCrzNpvAYettecBFwN/D3wJuNY7FgEuM8asBD4InANcBHzWGNMCfAB4zDv3ZuDa\nGo1bPHP1WSsnWMuVQLW4oG42ru8hk3XY8eL87LfmB2szzllTsCYiTaBWwdp/A5/yvo7gZs3OBH7m\nHbsNuBA4G9hurR231h4DdgCnAecCtxecKzXkt+aYsc9aGb8UJ5vhanFBveTmrc3TUuiMmTUvWEtp\nzpqINIGalEGttUMAxpgu4BbczNgXrLX+v6SDwCKgG8j/X/xix/1jc+rpaSdesKS/t7ervDexwHV0\nuM1TFy1qy93D3t4uHD/jFosGvrd7Dg4RjUY482WraE3W5KPY0Orx2XxtVys3/M+j7Nw7MC//brx4\nZBSAxd1tU8bf3enOgezqap2X76uR6P6FR/cyPAvtXtbsN6QxZi3wTeAfrLX/boz5fN7DXcBRYMD7\nerbj/rE59fePTPm+t7eLvr7Bssa/0B096t7LkeFx+voGp9zLWDTC2Hg60L1NZ7LsePEYx/V2MHhs\nlIX+U6nnZ3PDyi6efeEoL7zUP++C5kOHhwGYyPv89fZ2kZpIA9B3aIjulvk3F69R6N/M8OhehqeZ\n7+VMQWhNyqDGmBXAD4CPW2u/5h1+2Bizxfv6EuAu4H7gPGNMqzFmEbAJd/HBduDSgnOlhmbqswZu\nKTRoU9wXDg6RzmQ1X60BbFznzlt7dh7OW/NXg87cukNz1kRk/qvVnLVPAD3Ap4wxPzXG/BS3FHqd\nMeYeIAncYq3dD9yAG4zdCXzSWjsGfAXYbIzZBrwPuK5G4xbPTH3WwG3nEXS7qZ0vqRluo9i4fjEw\nP/utzdxnzW/doTlrIjL/1WrO2oeADxV56IIi524FthYcGwEur87opBQzbeTuHwu6GvS5fVpc0ChO\nXrOYeCzKg8/08fYtJxKNTP8ZN6oZM2vx6JTHRUTmMzXFlZJM9lmb/pEJWgY9dGyUh589RHd7ghVL\n2kMbo5SnJRnj1acu52D/KI8/d6Tewwkktxp02kbu2m5KRJqHgjUpyUx91sDdH7TUMqjjONx029OM\nT2S4/PUnzassTjO78My1APzowRfqPJJg/GAsmZj6T5kfvKUDludFRBqRgjUpyexl0GjJe4Pe9eg+\nntjVz8tPWMrrXrYy1DFK+dav7OLk4xbx+HNH2OetsJwPJmbIrPn9/9JpzVkTkflPwZqUZNYFBiWW\nQSy/YHIAACAASURBVI8MjPGfdz5LW0uMKy42RJRVaygXnuVm1+588KU6j6R0k01xp7bnSGg1qIg0\nEQVrUpKZdjDwj821wMBxHG6+wzI6nuE3fuVklnS3VmWcUr4zTl5GT1cL2x7fx+h4ut7DKcncOxgo\nWBOR+U/BmpRkrjLoXHPW7n58P4/uPMypG3o477RVVRmjVCYei/Irr1zD+ESGbY/tq/dwSpJKz9C6\nw/s+aP8/EZFGpGBNSlJJGTSTzfKfd+6gJRHjyos3qvzZwM5/xWrisSg/fvDF3KKSRjZzZs39jCmz\nJiLNQMGalKSUPmvODL/cn983yNBoitduXsGyxW1VHadUpqs9yWs2r/DaeByu93DmNFOftYRad4hI\nE1GwJiWZK1gDZszEPLnL7d116oYlVRqdhOnCM48D4I77X5gxAG8UfmYtqe2mRKSJKViTksxeBp19\nftCTzx8hEoFNG3qqN0AJzboVXWw+fglP7e5v+LlrE3MsMNB2UyLSDBSsSUlmy6zFvWPFVoSOjqfZ\nuXeADSu76WhNVHeQEporL95IW0ucf//RsxzsH6n3cGbkZ9bi03Yw8PcGVWZNROY/BWtSkmx2tu2m\nvMxakWDtmReOksk6nKqs2ryydFErv/WmUxifyLD11idL3qGi1lLpLPFYdNqiFe0NKiLNRMGalCQz\nWxnUz6wV+cX4hDdfbbPmq807rzl1BWdvWs7Olwb4/j27a/rajuNwz+P72frdJ9n+2Mx931Lp7LT5\naqCmuCLSXOL1HoDMD7MuMIjNXAZ9clc/yUSUE9csqu4AJXSRSITfvsjw7IvH+Pa2XWw+fiknrO6u\n+uuOjqf5xg8s9zxxAIB7ntjPzXdYTjtxKee8fBWnn7Qsd24qnZk2Xw0ms72asyYizUCZNSnJrAsM\nvGOF+4P2D46z99Awp6xdXPQXqjS+jtYEv/erm8g6DltvfTI3R6xa9hwY5DM3/oJ7njjACau7+fi7\nzuDXzz2eZYtaedD2ccMtj/LLZw/lzk9lskU/W4nc3qDKrInI/KfMmpQkO8cOBjC9DPqkSqBNYdOG\nJbzhlcfx44de5KcPv8QbX7U21Otnsll2vjTAL3cc4kcPvEA643Dxq9fxtvNPIB6LYtb18OZzNvDE\n80f40n89wt1P7Of0k93sWiqdpbNt+sIVte4QkWaiYE1Kkp5tb9AZyqDqr9Y83nLuBu5+Yh/fvXsX\n57x8Fe2tlf3TMZ7K8JDt4+Fn+3hiV39uTlpnW4Lf+7VNnHbisinnRyIRNh+/hBVL2nl0xyHGJzK0\nJGOk0sUzawrWRKSZKFiTkpRSBs3vs+Y4Dk/u6qe7I8lxvR21GaRUTVd7kktevZ7//flz3H7/Ht52\n/gllXWfPgUF+9she7n3iQC5AW7aoldecuoKXn7CUTet7aEnGij43Eonwqo3LufXuXTyy8xBnb1ox\nY7AWjUaIRiKasyYiTUHBmpRktjJoPDq9dcdLh4Y5NjzBazav0F6gTeKNZ63lxw+9yA9+sYdfeeUa\nFne2FD0vncmy99AwLxwc4vDAGEcGxukfHKfv6Cj7j7g92xZ3JnnDmet57eaVrFzSXvJn5GwvWPvF\n0wc50/SSyTq5lZ+F4vGIWneISFNQsCYlyWQdopFI0V+qk2XQyV+MTz7vlUDXqwTaLFqSMS4793hu\nvt3yne27+J2LTO6xF/uG+NEDL7Jr/wB7Dw0XzWi1tcQ4/aRlnP+K1bz8xCVFe/bNZU1vB6uWtvPo\nzsMMjbqZuUS8eCYuEYuqDCoiTUHBmpQkk3WKlkCheBn0yd39AGw+XsFaMznvtFXccf8L/PyXe3nT\nq9bS1Z7gW3c9z08eeoms45CIR1m7vJN1K7pYt6KL5Yvb6OlqoaerhbaWyv+58Uuh39m+iweePghM\n3xfUF49FtRpURJqCgjUpSSbrFC2BQn7rDvcXYzqTxe45yqql7fR0FS+VyfwUi0Z5xwUn8OVvPs4/\nffsJDg+MMTSaYkVPG7/xhpN5+QnlZcyC8IO1ux939y2dqS1MPKY5ayLSHBSsSUmys2XWCjZyf+Dp\ng4ynMpx24tKajU9q55Wn9HLC6m6e2ztASzLG5a8/kTeetXba/pzVsqa3k9XLOnh+3yAwubVUoXgs\nythEpiZjEhGpJgVrUpJsCZm1TNbBcRzuuP8FIhF4/SuPq+UQpUYikQhXvflU7nviAOefvnrGhQbV\n9KqNy/n2tueBWTJr8Sj/f3t3HidXVed9/HNr631JJ91ZICEbHCEEEgj76jzDovICRRlBo4KOiuPI\nuMzjiuPo6DPq4zgzDOM8vmAcxMGRTUZkQJSdhC1AgkmAE0ISspCkO+n0Xl3bvc8f51ankvSWrauq\n832/Xv1K1V2qzj2prv7d39myycxYFktE5LDQtPIyKtlRBmtrNnXw1vZuTjmumZbGqrEsooyhyROq\nufzcWUUJ1AAWvaNl4PFwfdY0GlRExgMFazIqvu+PohnU5+EXNgFwyekzxqxscuQ5alINR4Xz9w3X\nZy2nPmsiMg4oWJNRGU0z6Ns7e3ll7Q7mTKtnrhZul8PstDC7NtQ8a/FohJwf4AcK2ESkvClYk1EZ\nzWjQx1/eQgBcrKyajIHzT57GibOa9lmaKm9gySlN3yEiZU4DDGRUhp1nLZwUt7c/y8T6Sk45bvA/\nniKHUmNtBV/84IIh9xeuD5qIDz5xrohIOVBmTUZluGbQWMG8WhedNv2wz7MlMhqx8CZCc62JSLkb\n08yaMeYM4AfW2guNMXOB24AAWAV81lrrG2M+CXwayALftdY+YIypAv4TaAG6gY9Za9vGsuxHutGs\nYFBVEeW8k6aOZbFEhpSff01LTolIuRuzFIgx5svArUBluOnHwI3W2vMAD7jCGDMFuAE4B7gE+Htj\nTAXwGWBleOztwI1jVW5xXJ+1wT8utdVxAC5YcNQhWVJI5FDIN4Nq+g4RKXdj2V71JnBlwfNTgSfD\nxw8BfwqcDiy11qastZ3AWuAk4Fzgd3sdK2NouGbQuUc18L+vXsCV588e41KJDC2uAQYiMk6MWRrE\nWnuvMWZmwSbPWpvvTNINNAD1QGfBMYNtz28b0YQJ1cRie3Ysbm6u2++yH+mCICDnB1RUxPaov8LH\nLS31xSjauKLP5qHT3FxHXThhb119ler2IKjuDh3V5aFzpNVlMdusCm9364AOoCt8PNz2/LYR7drV\nt8fz5uY62tq6D7C4Ry7fdzG1n/MH6k91eWipPg+dfF1m0lkA2nb0UF+h0aAHQp/LQ0d1eeiM57oc\nKggt5rC95caYC8PH7wKeBl4AzjPGVBpjGoDjcYMPlgLv3utYGSM538XVQzWDipSiwqk7RETKWTGD\ntS8B3zbGPAskgHustduAm3DB2GPAN6y1/cC/AfOMMUuATwHfLlKZj0i5MLM21GhQkVKUHw2qAQYi\nUu7GtBnUWrsBODN8vAa4YJBjbgFu2WtbH3DVGBRRBpFvBlVmTcqJ5lkTkfFCs5fKiHIK1qQMabkp\nERkvFKzJiNQMKuUorj5rIjJOKFiTEakZVMpRfs1a9VkTkXKnYE1GpMyalKN8Zi2nPmsiUuYUrMmI\ndmfW9HGR8qHlpkRkvNBfXxlRVs2gUoa0kLuIjBcK1mREvppBpQwNTN2h0aAiUuYUrMmINMBAytHu\n0aDqsyYi5U3BmoxI86xJOdJyUyIyXihYkxHl1wZVM6iUEw0wEJHxQsGajEjNoFKORrPcVBAEpDK5\nsSqSiMgBUbAmI9I8a1KORjMa9PaHLV+8eQnrt3aNVbFERPabgjUZkTJrUo5GWm7qrW3dPLnibZKp\nHDf/eiWdPamxLJ6IyKgpWJMRZTUprpShkRZyv/uJtQAsMs3s6k5x830ryWiaDxEpQfrrKyPSPGtS\njobrs7Zq/U5e3bCLebOa+Mx7T+TMEybz5pYufvF7SxBoqg8RKS0K1mREagaVcjTU1B1+EHDP42/i\nAVddOAfP87j2Xe/gmCl1LPnjVh55aXMRSisiMjQFazIizbMm5WioAQbPv7qdja09nDlvMjMm1wGQ\niEf53JXzqa9JcOeja3nhte2DvubG7d187xcvDrl/vPCDgG/ftox/+tXLxS6KiKBgTUZB86xJOYp4\nHtGIx9s7+3jqlbfpSWbIZH3ue2odsajH+86bvcfxTfWVfO7986lIRPjpb1bzxPIte+x/bUM737/j\nZd7c0sV//n4Nff2Z/S5TEATcv2Q937z1edq7+g/q+g6n9Vu7eGtbN08t30JaU5uIFJ2CNRmRMmtS\nrhYeO4mu3jS3PfQ6X/iXJfzdz5exo7OfPznlaCY1Vu1z/JxpDXz5mlOorY5z+8OW/3l2A+CycT++\n6xWyOZ8FcyfRk8zw4HMb96ssQRBwzxNv8t9L1rNlRy9P/3HrIbjCw2PFGzsAyGR91mzuKHJpRETB\nmoxIAwykXP3F++bzg+vP4gMXzuHo5lo2t/VSUxnjsrNnDnnOMVPq+NriU2mqr+DeJ9fx4ztX8NP7\nV5OIR/jCny3g+ivmMaGugj+8uGnU2bEgCLjr8bU89PxGJjdVk4hHWLpyK36JDmbIB2sAq9e3F7Ek\nIgIK1mQUNMBAyllzYxXvPvMYvnXdaXz/02fyt9edTm1VfNhzpjRV8/XFpzKlqZpV69tpqE3w1Q+f\nyvHHTCARj3Ll+bNdk+rT60Z8/yAI+K9H3+DhFzYxdWI1X/nQQhaZFnZ09vPGptLLWrXu6mPLjl5O\nmDmBRCyiYE2kBChYkxFlFazJONEyoZqJDZWjOrapvpKvLj6FK8+fzY0fWcT0ltqBfWfNm8LRzbU8\ns3Ibm1p7hnyNvv4stz7wGo+8uJlpk2r48odOobG2gnPmTwVg6cptB3dBh0E+q3b68ZM5cc4kNrf1\n0qEJg0WKSsGajEjNoHKkqq9OcNnZM/cJ8CIRjz975xwCdk+uu7fVG9r5m589z7Ort3HMlDq+fM1C\nGmoSAJgZjUysr2TZ6630p7OH+zL2y4q1O/CAk+dMZKFpBtQUKlJsCtZkRGoGFdnXvFlNnDBzAqvW\ntfPE8i2se7uLrTt72dGZ5Be/t/zDr1bQ0Z3m8nNm8o2PnEp9GKiBG6l6zvwppDI5XrJtRbyKPfUk\nM6zZ1MnsafU01Faw8LgWwAWeIlI8sWIXQEpfTstNiezD8zyuunAu375tGbc/bPfZP21SDZ94z/HM\nmlo/6PlnnziF+5duYOnKrQPNosW28s2d+EHAgmMnATBjSh0NtQleXd+OHwREPN2wiRSDgjUZUU7N\noCKDOmZKHTd84CQ2be8mmcqRTGdJprIcNamGS8+YQTwWHfLclgnVHDe9kdc3drCjIznoVCJjbfla\n119twVwXrHmex4kzm1i6ahubW3sGJhEWkbGlYE1GpHnWRIa2YO6kgeBmf50zfwprNnXwzKptXH7u\nrENcsv2TyfqsXLeTlsYqpk2qGdg+b5YL1lavb1ewJlIkateSEanPmsjhsci0kIhHWLLXnGs536e1\nI8nq9e08vnwLdz++lt8uXc8Lr23nrW3dh2VQgt24i1Q6x4JjJ+EVNHeeMLMJgFWHYZDB9l199B7A\nShAiRxpl1mREGg0qcnhUVcRYZFp4ZtU2/uqfnybrB2Sz/kA2ezjNjZWY6RMwMxox0xsPuhl1eThl\nx8Jj98wS1tckmDG5ljc2d5DK5KiID920O1qpTI77nlrHH5Ztor42wQ3vP2nIvn0iomCt5CVTWTwP\nKhPF+6/KhmuDKrMmcuhdfNp0Nm7vxg8gFvWIRyPEYxEm1FXQ3FhFy4QqmhurSKZybG/vY/uuPra1\n97FhazdLVm5lyUq3bFVVRYyKeIRELEoi7l4jHo0Qi0WIha9ZlYhRWREd+BfczVgQwEtr2qipjDH3\n6IZ9ynjirIls3N7Dmk0dzJ898aCu127cxX88+DqtHUkm1FXQ0Z3iB3e8zJ9fdgKL3tFyUK8tMl6V\nTbBmjIkAPwFOBlLAn1trB5/gaBzYurOXPyzbxNJV2/A8eOfCo7j09Bk01FbscVxff4at7X1UxKNU\nV8SoTMSoTEQHzYKl0jl2dvWzs6ufbTvdLOVv7+xl645e/ACa6iqYEP5MrK+keUIVLY1V9PW7Jhdl\n1kQOvRmT6/jOJ84Y3cFzdgdKvh+wqbUHu6kDu3EXrR1JMhmfdDZHTzJDOuuTzfn7VZbzTpo66Kjv\nebOaePC5t1i9vp15s5roSWbo6k3T0ZOivSvFjs5+2rv66U1maGqopKXRfXc01VeSTGXp6kvT2ZNm\n4/buge+0S0+fwRXnzeK1Dbv46W9X85P/XsX7zp/NZWcds0czrJSPIAjI+QGx6PjpYZVMZVn2eivb\nd/Xx3nNnE48V59rKJlgD3gtUWmvPMsacCfwDcEWRy3TQgiCgP52jO5mhuy9NR3eKpSu3sSIclTWp\noZKcH/DwC5t47OUtXLBgGifNnsiazR28umEX67d2MdjygolYhIpElIp4lHgsQldvmt7+ffu5RDyP\n5glVxKMe7V0ptuzoHbKsyqyJlI5IxOOYKXUcM6WOi0+bPugx+T+emaxPJufTn87Rn3IjVvvTOcCN\n+Ix47vd79rR9s2oAc49qIBGP8OhLm3nkxc0Htabp1InVfPw9xzMnfK8Fx07i64tP5aZ7XuG+p9ax\n7LXtVFbEiITlisei1NfEqa9J0FCdoLoyTjbnk87kBgLSiOcRiXhEox5Rz9udTQwzi8lUls7eNF29\naTp703hAXXWcuuoEddVxqsL38zz3nRiNetRUxqmtilNTFaemMkYk4uGF9QXgBwG5nE8mGwwExZGI\new0Pj2hkd5kinkcQBGSyOVIZn0zWJ+f7JGLu+7ki7m6w3XX5ZLLu2jLh9WVyPtmsTzYX7H7su+tO\nxKLE4xESsQi+H9CXytLXn6UvlSWdyYXXFf4fRyPUVsUHrr2mMkYqk6M/nSMZntfakWTrzl62tfex\nvT2J5+HqqSpOXU2CRCxCEAT4we6/X129abr6XP1mcwE1lTEa6yporElQV5MgCNxAtVzOJwgYuOaK\nRJTKRJTG2gpawuTAxIbKUQd7QeDqPpsLBsrkBwGBv7t8fuAyx4WfkVjEI5316e3PkuzP0pfKEARQ\nUxmnqjJGdUWM1l19LPnjVpbZVtIZn0Q8wsWLpu+TMBkr5RSsnQv8DsBa+5wxZlGRy8N9T63jyRVb\nBj4UQQD7foXt3h6ED9xjt8/3GfSLb860ei45fQanHNdMzg9YsnIrDz67gUdedF+W4D58c6Y1MHNq\nHblcQDKdpT/lfulSmdzAL2F/X4b6mgSzptbTVF/JxPoKJjdVM21SDZMnVO9xp9CfzrKrO8XOzn5a\nO5K07krS1pEkEY/SMqH4UwuIyOh5nkcs6hGLRqgC6qsP7HXisQgXLZrOstdaqa9JDPw01iTcd0qD\n+16pqYqzs7Ofto4krR1J2jtTVFXGaKhJuJ9a9z209x/j6S213PjRRfz0/tW8sblz4A9sKcoHY/sT\nsOYThcOd4nnD7x9rEc+jubESPI/OnhRvD3MjH4tGaKhJML2ljop4hM7etLv5bxv6nOHeNxGPEOQD\nr3ywFdkdTHsepDOj69t5sJobKzl3/lTOmT+1aIEagBeU0qdjGMaYW4F7rbUPhc83ArOttUMOi8pm\nc0FsmHmODtbdj67h0WWbiETcXYsXfoj25uHhRcAD8PJ3Z+HxQE1VnIbaCnfnWFvBibMn8o5wBFah\nTNbniZc2sbm1h3lzJnLi7IlUVw6/ILWISLnKZ246ulN0dKfY1d1PTzJDIh6lIh6hIh4jHouEWa6A\nrO+Ty/lkswGZXM5lFLM+VRUu0zOhrpLGOvcHt7MnRVePa8pNprJ7ZIvSGZ+eZJruvgzdvWl6kxn3\nHv7u4CE60L/QZcfw8v3/dt+EZ3M+uVxAzvfxPI+KuOtPmIhHiUY80hmfVCZHOuPK6vbvPiYeC/se\nxgoeRyMDj30/IBW+RiqdHcic1YQ/+cEgvu8yTJms77JgPSk6elL0JDNUJmJUh9mk6qo4k5uqOaq5\nlikTa/a4kc/mfLp702Sy7loikXxgFaW6MjZo03V/OktXT3p3RisaIeJ5pLM5+lM5+tMum9fWkWRb\n2CVn285e+tM5Il4+QPMGPgt+mDED9qiPWCwykDkbyJAWPPY8F1z7A5+RgHgsQm2YXaytdquL9CYz\n9CYzYb1Eeeep05k3e+JYdwEa9M3KKVj7MfCctfau8Plma+3Rw53T1ta9x8U1N9fR1tZ9GEt55FBd\nHlqqz0NHdXnoqC4PHdXloTOe67K5uW7QYK2cegEuBd4NEPZZW1nc4oiIiIgcfuXUZ+0+4CJjzDO4\nNOF1RS6PiIiIyGFXNsGatdYHri92OURERETGUjk1g4qIiIgccRSsiYiIiJQwBWsiIiIiJUzBmoiI\niEgJU7AmIiIiUsIUrImIiIiUMAVrIiIiIiVMwZqIiIhICVOwJiIiIlLCFKyJiIiIlDAvCIJil0FE\nREREhqDMmoiIiEgJU7AmIiIiUsIUrImIiIiUMAVrIiIiIiVMwZqIiIhICVOwJiIiIlLCYsUuwOFi\njIkDPwNmAhXAd4FXgduAAFgFfNZa6xepiGXHGNMCvARcBGRRXR4QY8zXgMuBBPAT4ElUlwck/D3/\nOe73PAd8En0295sx5gzgB9baC40xcxmk/owxnwQ+javf71prHyhagUvYXnW5APgX3GczBXzUWrtd\ndTk6hXVZsO1DwOestWeFz4+IuhzPmbXFwE5r7XnApcDNwI+BG8NtHnBFEctXVsI/ij8FkuEm1eUB\nMMZcCJwNnANcAExHdXkw3g3ErLVnA98Bvofqc78YY74M3ApUhpv2qT9jzBTgBtzn9hLg740xFcUo\nbykbpC7/GRdYXAj8GviK6nJ0BqlLjDELgU/gPpccSXU5noO1u4Fvho89XNR9Ki6LAfAQ8KdFKFe5\n+hHw/4C3w+eqywNzCbASuA/4LfAAqsuDsQaIGWMiQD2QQfW5v94Erix4Plj9nQ4stdamrLWdwFrg\npDEtZXnYuy6vttauCB/HgH5Ul6O1R10aYyYC/wf4fMExR0xdjttgzVrbY63tNsbUAfcANwKetTa/\nZEM30FC0ApYRY8y1QJu19uGCzarLAzMJWARcBVwP3AFEVJcHrAfXBPo6cAtwE/ps7hdr7b24IDdv\nsPqrBzoLjlG9DmLvurTWbgUwxpwN/CXwj6guR6WwLo0xUeDfgS/i6ivviKnLcRusARhjpgOPA7+w\n1v4SKOy3Ugd0FKVg5efjwEXGmCeABcDtQEvBftXl6O0EHrbWpq21FnenXfjlorrcP1/A1edxwMm4\n/muJgv2qz/032PdkV/h47+0yAmPMB3GtEu+x1rahujwQpwLHAv8G/Ao4wRjzTxxBdTlugzVjzGTg\n98BXrLU/CzcvD/sMAbwLeLoYZSs31trzrbUXhP0uVgAfBR5SXR6QJcClxhjPGDMNqAEeVV0esF3s\nvrNuB+Lo9/xgDVZ/LwDnGWMqjTENwPG4wQcyDGPMYlxG7UJr7bpws+pyP1lrX7DWzgv/Bl0NvGqt\n/TxHUF2O29GgwNeBCcA3jTH5vmt/BdxkjEkAr+GaR+XAfAm4RXW5f6y1Dxhjzsd9yUSAzwLrUV0e\nqH8EfmaMeRqXUfs68CKqz4Oxz++2tTZnjLkJF7hFgG9Ya/uLWchSFzbd3QRsBH5tjAF40lr7LdXl\noWGt3Xak1KUXBMHIR4mIiIhIUYzbZlARERGR8UDBmoiIiEgJU7AmIiIiUsIUrImIiIiUMAVrIiIi\nIiVMwZqIlBxjzBRjTMYY8+p+nhczxnyh4Pm1xpjsKM/9W2PM2oLnZxljztmf9xcRORwUrIlIKVqM\nm3/ueGPMeftx3gdxC5Hn3QkcNcpzfwScWfD8Kdys6SIiRTWeJ8UVkfL1MdyyMpcBn2L0qxB4hU+s\ntUkgOZoTrbU9uLVGB30tEZFi0aS4IlJSjDGLgGXABcB5wI3ANGvtrnB/HfB94ANANfAMcAMwFbcW\ncN514b+3WmtjxpjbgJnhkjX59zoNt5rEccCHgcXW2rnGmA3AMeFhT+KWWbvIWjuv4Nw5wFpgobV2\nxV7XUAPcDLwbt/brCuDr1trHwv1x4Fu4oLQp3P8la+1z4f5zge8BpwB9uAzhV621fcaYmbis4zeA\nz+PWmz0ZmIxb0eFiXID6OPBFa+3bw9e4iJQ6NYOKSKm5FtiOW0f1LqAStx5t3l3A/wKuARbhsmEP\nA8/i1mEEF7jdudfr/hy3jmBhs+iHgWettW/sdexpQA4XDF0ZnnuCMWZhwTGLgVf2DtRC3wFOAC4J\n/10O3BcGceCWIfpEWN6TccHa74wxzcaYM4DHcAHraWF9XDHI9VwDnB9eQxx4AheknR2+bwJ4LFw2\nSkTKmII1ESkZYWBxDXCvtdYPg6iXcU2hGLfA4qXA9dbax6y1Ntx3L24t4E5wawaGTaCFngA24fq1\n5dduvBoXiO3BWtsWPuy01rZba5cDf8QFaHmLgduGuJS5QDewPlzA+6+B9wO5MDP4ceBr1trfWGvX\n4tYtvhWXZfsS8KK19q+tta9bax8CrgcuM8bMK3iPm8P9L4d1VgNca61dFQaQ1+D6671/iDKKSJlQ\nsCYipeRyXMByd8G2O3FZrXOB+eG2F/I7rbU7rbVfstZuG+6FrbUB8AtcEAMuO9fIvhmrodwGXGOM\niRhjzgJmAr8c4tj/i2vCbDPGPAF8Dng9XGTa4LJehdeQDYMzC5yIa9otlO+zd2LBtnUFjxcCzUCn\nMabHGNODax6tAY4f5fWJSInSAAMRKSXXhv8+4pJowO6O/vkM2sG4HbjRGHMs8CHgfmttxyjPvQP4\nIfBOXNPog9ba1sEOtNYuMcYcjWuOvBj4bPi+ZwCZEd5nsAER+RvrzBDHpYHVYbn2NtrrE5ESpcya\niJQEY8wUXHDzE2BBwc/JuD5pHwDy2bNFBefVG2Naw8zbsCOmwmbVZ3DNn+/FBW9D2eO1wsDsTrPq\neAAAAexJREFUQVxAdAWDNJ8WlOlvgHOstfdZaz+DmwIkDbwHNyghu9c1RIwxa4wxVwOv4vqdFTo3\n/Pe1Id5yNTAL2GmtXRs2rbbipjGZP8Q5IlImlFkTkVKxGHcD+UNr7VuFO4wxP8QFcmcAvwF+Yoz5\nDNCGGzXZieuQ3xQevwh4fYj3+TluTrUk8LthytONa35tKcig3YbLsPUBDwxz7kzgI8aYT+JGbl6E\nGxX6vLW21xjzr8D3jDE7gDeAL+D63D0OrAKWG2N+BNwSvta/4jJ5r4WjQfd2B2506F3GmK8B/bgR\ns6fjAjkRKWPKrIlIqfgY8MDegRpAOOXFK7im0Gtx/b1+AzyP6/91qbU2hQt2nsRlzz41xPvchRs9\n+Utr7XCrG3wf+AtcVi/vAVyQ91/W2vQw596AG9H5S2ANLhi7zlr7ZLj/K2E5foYbCXoCcIm1dru1\ndhVufrkLcIMa/gO4D7hqqDcLB1NchAsiHwOW4m7G/2SoploRKR+aZ01EZJSMMROBt4GzwlGYIiKH\nnYI1EZERhEHahbj53iZZa7VmqIiMGfVZExEZWRz4d1xWbbARlyIih40yayIiIiIlTAMMREREREqY\ngjURERGREqZgTURERKSEKVgTERERKWEK1kRERERKmII1ERERkRL2/wF0UJOhEvxfTgAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113dd3550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_= np.empty(1)\n",
    "for i in np.unique(target_reg):\n",
    "    y_ = np.append(y_, len(target[target_reg == float(i)]))\n",
    "y_= y_[1:]\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.plot(np.unique(target_reg), y_)\n",
    "plt.xlabel('Activity score', fontsize=15)\n",
    "plt.ylabel('Occurences', fontsize=15)\n",
    "plt.title('Distribution of activity scores', fontsize=17)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe 3 peaks of molecular activity: 40, 60, 80. We make the following assumptions (according to the article):\n",
    "\n",
    "> the peak at 40 can be seen as a measure representing the \"inactivity\";\n",
    "\n",
    "> the peak at 80 can be seen as a measure representing the \"activity\";\n",
    "\n",
    "> there are variations around these values, with a limit at 60. This limit is highly questionable. However, this is the threshold used in the article, but we could imagine testing other limits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Justification of the molecular activity limit at 60: maximum limit such that every molecule is active at least on one cancer cell line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All molecule active at least on one cancer cell line with molecular activity limit at 60: True.\n"
     ]
    }
   ],
   "source": [
    "print('All molecule active at least on one cancer cell line with molecular activity limit at 60: %s.' %(np.sum(target_reg>=60, axis=1)>=1).all())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Markov Network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to construct from the target data (cell-lines bioactivity) a correlation matrix between the pairs of cell lines. Then, we can extract the Markov network from the matrix by favoring high-valued pairs (two methods will be explained afterwards)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 - Correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cor = np.corrcoef(np.transpose(target_reg)) # covariance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of edges: 191 out of 1740.\n"
     ]
    }
   ],
   "source": [
    "' we make use of a threshold of 0.9 to build the Markov network (first method) '\n",
    "thres = .9\n",
    "affinity_matrix = cor.copy()\n",
    "affinity_matrix[np.abs(cor) < thres] = 0\n",
    "print('Number of edges: %d out of %d.' %(np.sum(affinity_matrix!=0)/2, cor.size/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Network representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/networkx/drawing/nx_pylab.py:126: MatplotlibDeprecationWarning: pyplot.hold is deprecated.\n",
      "    Future behavior will be consistent with the long-time default:\n",
      "    plot commands add elements without first clearing the\n",
      "    Axes and/or Figure.\n",
      "  b = plt.ishold()\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/networkx/drawing/nx_pylab.py:138: MatplotlibDeprecationWarning: pyplot.hold is deprecated.\n",
      "    Future behavior will be consistent with the long-time default:\n",
      "    plot commands add elements without first clearing the\n",
      "    Axes and/or Figure.\n",
      "  plt.hold(b)\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/__init__.py:917: UserWarning: axes.hold is deprecated. Please remove it from your matplotlibrc and/or style files.\n",
      "  warnings.warn(self.msg_depr_set % key)\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/rcsetup.py:152: UserWarning: axes.hold is deprecated, will be removed in 3.0\n",
      "  warnings.warn(\"axes.hold is deprecated, will be removed in 3.0\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAE+CAYAAADWJPrxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAABSRJREFUeJzt2bFtw0AQRcGloepUgGKVplgNsLxzoMBQZvkFNIEZgNgL\nf/YCbmutNQDAn30dPQAAzk5MASASUwCIxBQAIjEFgEhMASASUwCIxBQAIjEFgEhMASASUwCIxBQA\nIjEFgEhMASASUwCIxBQAIjEFgEhMASASUwCIxBQAIjEFgEhMASASUwCIxBQAIjEFgEhMASASUwCI\nxBQAIjEFgEhMASAS05mZ53Pmen1dAPjQ5egB/8LjMbPvr/ftdugUAM5HTGdm7vf3CwAf2NZa6+gR\nAHBm/pkCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJ\nKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSm\nABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgC\nQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoA\nkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBE\nYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJ\nKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSm\nABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgC\nQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoA\nkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBE\nYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJ\nKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSm\nABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgC\nQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoA\nkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJKQBE\nYgoAkZgCQCSmABCJKQBEYgoAkZgCQCSmABCJ6czMtv18ABzr+Zy5Xl/3JC5HDwCAN4/HzL6/3rfb\noVN+S0wB+F/u9/d7Attaax09AgDOzD9TAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJTAIjE\nFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJT\nAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwBIBJTAIjEFAAiMQWASEwB\nIBJTAIjEFAAiMQWASEwBIBJTAIi+AaGyHqg1LtoNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11434ed30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "\n",
    "G_correl_matrix = nx.from_numpy_matrix(affinity_matrix) # total graph\n",
    "nx.draw_spectral(G_correl_matrix, node_size=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We remark here that it seems that there are only three cancers in the dataset; on the contrary, the thresholding-based graph-building method created here three isolated subgraphs, with no interaction between one cancer from a subgraph and from another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "221"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(G_correl_matrix.edges())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Maximum weight spanning tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the paper, the maximum spanning tree is computed from additional external data. We used data from Juho Rousu's github page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "' Import the data '\n",
    "cancerCL_corr = np.loadtxt('data/data_clean/ncicancer_cancerCL_corr.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/networkx/drawing/nx_pylab.py:126: MatplotlibDeprecationWarning: pyplot.hold is deprecated.\n",
      "    Future behavior will be consistent with the long-time default:\n",
      "    plot commands add elements without first clearing the\n",
      "    Axes and/or Figure.\n",
      "  b = plt.ishold()\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/networkx/drawing/nx_pylab.py:138: MatplotlibDeprecationWarning: pyplot.hold is deprecated.\n",
      "    Future behavior will be consistent with the long-time default:\n",
      "    plot commands add elements without first clearing the\n",
      "    Axes and/or Figure.\n",
      "  plt.hold(b)\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/__init__.py:917: UserWarning: axes.hold is deprecated. Please remove it from your matplotlibrc and/or style files.\n",
      "  warnings.warn(self.msg_depr_set % key)\n",
      "/Users/benoitchoffin/anaconda/envs/py35/lib/python3.5/site-packages/matplotlib/rcsetup.py:152: UserWarning: axes.hold is deprecated, will be removed in 3.0\n",
      "  warnings.warn(\"axes.hold is deprecated, will be removed in 3.0\")\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAE+CAYAAADWJPrxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVfXe/vEbUEER0hxyyqk8mqA4klkqRxBwDk0stTTj\nNBwqy+E4lI3WgzPqkZxKM83STHNAUMEcyNCcQu1gmKaZDRY4D8Bevz/y53PqaVC38F177/fruvrH\nEm67Lrz357O+ay0vy7IsAQCA6+ZtOgAAAK6OMgUAwEmUKQAATqJMAQBwEmUKAICTKFMAAJxEmQIA\n4CTKFAAAJ1GmAAA4iTIFAMBJlCkAAE6iTAEAcBJlCgCAkyhTAACcRJkCAOAkyhQAACdRpgAAOIky\nBQDASZQpAABOokwBAHASZQoAgJMoUwAAnESZAgDgJMoUAAAnUaYAADiJMgUAwEmUKQAATqJMAQBw\nEmUKAICTKFMAAJxEmUr6bto0nW3XTlq82HQUAIALKmE6gB1cmjlT/vv2qcDPTyViY03HAQC4GCZT\nSbeOHq3d1appxvnzsizLdBwAgIuhTCV59e6tegcOaMbPP2v27Nmm4wAAXAxlepm/v7+WLl2qbUOH\n6uQ993D9FABw1bws9pq/crx5c1XduVOXwsNVav1603EAAC6AyfQ3qg4frv21a2viTz/J4XCYjgMA\ncAGU6W/FxqregQNKLltWr732muk0AAAXwJr3D3z77bdq0aKF5s6dq6ioKNNxAAA2Rpn+iY0bN6p3\n797KzMxUrVq1TMcBANgUa94/0a5dOw0dOlS9evXSxYsXTccBANgUk+lfsCxL9913nypXrqw33njD\ndBwAgA0xmf4FLy8vzZ07V6WWL9exJk24/xQA8H8wmV6lU23bKnDzZuWXK6eSM2dKPMMXAHAZk+lV\nCnzySeWWKKGSeXm6NGuW6TgAABuhTK9WbKx83nhDqT4+GvXllyosLDSdCABgE5TpNQiMi1Olbds0\n9fhxPfroo6bjAABsgjK9Rs2aNdPMmTN1/u23dbRxYw4kAQA4gHS99tasqeCjR3Xq7rsVuGWL6TgA\nAIOYTK9T0Pjx2nLTTXpqxw6dPn3adBwAgEFMpk44f/68atasqVtuuUVZWVny8vIyHQkAYACTqRNK\nly6tbdu26cCBA+rfv7/pOAAAQyhTJ9WpU0fLli3TggULNGPGDNNxAAAGUKY3QOfOnTV69GjFx8dr\n+/btpuMAAIoZ10xvoKioKNXYulWz7rpLPo88wiMHAcBDUKY3kMPh0J4aNdT0+HFZ0dHyWrPGdCQA\nQDFgzXsDeXt7q+HYsdoSEKBVFSuajgMAKCZMpkXgyJEjuvPOOzVv3jxFRUWZjgMAKGKUaRHZvHmz\n7rvvPm3ZskX16tUzHQcAUIRY8xaRNm3a6NVXX1W3bt108uRJ03EAAEWIybSIxcfH6+uvv9ZHH30k\nHx8f03EAAEWAybSIJSYm6syZMxo9erTpKACAIkKZFrGSJUtqyZIlWrRokd577z3TcQAARYA1bzHZ\ns2ePIiIilJqaqmbNmpmOAwC4gZhMi0lISIhmzJihmJgYff/996bjAABuICbTYvbCCy8oPT1d6enp\nKlWqlOk4AIAbgDItZg6HQz179lSlSpU0c+ZM3oEKAG6ANW8x8/b21vz58/XJJ58oKSnJdBwAwA3A\nZGrIwYMH1bp1a73//vsKCwszHQcA4AQmU0Nuu+02vfvuu7r//vt16NAh03EAAE6gTA0KDw/XqFGj\n1L17d505c8Z0HADAdWLNa5hlWYqLi1NeXp6WLFkib28+3wCAq+FvbsO8vLyUlJSk48ePa8yYMabj\nAACuQwnTASD5+vrqww8/VGhoqBo1aqSYmBjTkQAA14DJ1CaqVKmitJAQ3d2jh47ff7/pOACAa8A1\nUxsprFhRPj/9pB8k5X/zjapXr246EgDgKjCZ2ohPv34qrFBB7/n46F+1a6swOlpavNh0LADAX2Ay\ntaHc3Fz9p25d3ZWXJys6Wl5r1piOBAD4E0ymNlS+fHk1SUzUloAALS9XznQcAMBfoExtqnT//qr7\nn//omU8+0aJFi0zHAQD8CW6NsbFq1app9erVCg8PV/Xq1dW2bVvTkQAAv4PJ1OaCg4O1cOFC9erV\nS9nZ2abjAAB+B2XqAiIiIpSQkKBOnTrphx9+MB0HAPAbnOZ1IS+88ILWrl2r9PR0lSlTxnQcAMBl\nlKkLsSxL/fv315kzZ7RkyRL5+PiYjgQAEGtel+Ll5aU5c+YoNzdXw4YNMx0HAHAZZepiSpUqpQ8/\n/FApKSmaNm2a6TgAcEMlJydrWe3aUqVK0jPPmI5z1VjzuqjDhw/r7rvvVlJSkrp37246DgA4xeFw\naODAgZo/f75yS5bUTZcuSRUrSj/+aDraVeE+UxdVu3ZtffTRR+rUqZOqVaumli1bmo4EANclLy9P\nszt00P07dyp+6NBfinThQqlvX9PRrhqTqYtbsWKFHn/8cWVkZKhOnTqm4wDANfnkk0/Ut29frXQ4\nFHzkiNSxo5ScbDrWNWMydXHdunXTkSNH1LlzZ2VkZKh8+fKmIwHAXyosLNRrr72mpKQkzZ49W8Hn\nz0vz5kkDBpiOdl2YTN3E4MGDtWvXLqWkpMjX19d0HAD4Q0eOHFG/fv1UqlQpzZ8/X9WqVTMdyWmc\n5nUTEyZM0M0336y4uDjx+QiAXX3wwQdq0aKFOnfurLVr17pFkUpMpm7l3Llzat++vSIjI/XKK6+Y\njgMAV5w9e1bPPvus0tLStGjRIoWGhpqOdEMxmbqRMmXKaMWKFVq4cKHeeust03EAQJJ0aOxYfVa5\nshp98YV27drldkUqcQDJ7VSuXFnJyclq27atbr31VnXo0MF0JAAeyuFwaNKkSWoyerQi8vPVLiBA\nCgw0HatIMJm6ofr16+uDDz5Q3759lZWVZToOAA907NgxRUZGavny5Wo0ceIvt7y46Endq8E1Uze2\naNEiDR8+XJ9++qnbXOQHYH9Lly7VP//5Tz355JMaOXKkSpRw/yWo+/8JPdgDDzygw4cPq3Pnztq0\naZMCAgJMRwLgxs6cOaNBgwZp48aNWrFihe68807TkYoNa143N2LECLVs2VK9e/dWQUGB6TgA3FRm\nZqaaNGkiSdq1a5dHFanEmtcj5Ofnq2vXrqpTp46SkpLk5eVlOhIAN1GwaJEOvfiixn73nTrOnaue\nPXuajmQEZeohTp8+rTZt2qhPnz7617/+ZToOADeQkpKi+h07qo6k/Bo1VPLoUdORjOGaqYcICAjQ\nqlWr1Lp1a9WuXVuxsbGmIwFwUZZl6dVGjdR83z5VvvxrJU+eNJrJNK6ZepAaNWpo1apVevLJJ5WR\nkWE6DgAXdPz4cd17770KO3xYnSX516r1y3tHBw40Hc0oytTDNG7cWO+884569uypL7/80nQcAC7C\nsizNnz9fISEhCgkJ0V0zZ/5y7+i4cb+8wDsx0XREo7hm6qE+j4hQtQ0bVCYuTmVmzjQdB4CNHTt2\nTI8++qiOHTumuXPnqmnTpqYj2Q6TqYcK2r1bFR0OnZs9W+fPnzcdB4ANWZZ1pTxDQ0O1bds2ivQP\nUKYeyqdfP53z99dWy9LuatVUsGiR6UgAbOTo0aPq2LGjpk2bpnXr1unFF19UqVKlTMeyLcrUUyUm\nqsyZM4oMD9ddeXnaN3SoHA6H6VQADLMsS7Nnz1azZs10zz33KDMzUyEhIaZj2R63xng430cfVYGP\njz48ckQzn3xS06dP56EOgIf6+uuv9Y9//EO5ubnasGGDgoODTUdyGUymni42ViVSUzUkM1Pbt2/X\n8OHDxZk0wLM4HA698cYbatGihdq3b6+tW7dSpNeIyRSSpMDAQKWkpCgsLEwBAQEaPXq06UgAisGh\nQ4f0yCOP6Ny5c9q4caMaNmxoOpJLYjLFFRUqVNC6dev0zjvvaPLkyabjAChCDodD//73vxUaGqpO\nnTopIyODInUCkyl+pUqVKlq/fr3atm0rf39/Pfroo6YjAbjBDh48qIcfflh9P/tMh7285P/NN5KP\nj+lYLo2HNuB35eTkKCwsTGPHjlXfvn1NxwFwAzgcDk2cOFGlhg/XA5Ylf0n+0i+PA/zxR8PpXBuT\nKX7X7bffrtTUVEVERMjf31/33nuv6UgAnLB582adbtdOT1+en3wlWX5+UtmyEh+YnUaZ4g8FBQVp\n9erVio6OVunSpRUVFWU6EoBrlJ+frycqVFDM6dOKlFRSkry9pZtvllffvh7/TN0bhQNI+FPNmjXT\nsmXL1K9fP23atMl0HADXYO3atQoJCVGf/Hx1llSyYkWpZEkpKoqH099gXDPFVUlLS9MDDzyg1atX\nq2XLlqbjAPgTOTk5GjJkiPbt26dJkyap6/nz8nr7bWnAAIl3GRcJyhRXbeXKlYqLi9O6devUuHFj\n03EA/Mbp06c1ZswYvfnmmxo2bJieeeYZ+fr6mo7lEVjz4qp17dpVU6dOVXR0tA4cOGA6DoDLHA6H\n5s2bp/r16+uHH35QVlaWhg8fTpEWIw4g4Zr07t1bZ8+eVYcOHbRx40bVrl3bdCTAo3366ad6+umn\n5e3treXLlys0NNR0JI9EmeKaDRw4UGfOnFFERIQ2b96sqlWrmo4EeJxvv/1Ww4cPV3p6uhISEtS3\nb195e7NsNIX/87guTz/9tB555BFFREToxIkTpuMAHuPChQt6/fXX1bhxY9WsWVPZ2dl68MEHKVLD\nmExx3UaOHKnTp08rMjJS6enpKleunOlIgNuyLEvLly/XkCFD1KRJE23btk1169Y1HQuXcZoXTrEs\nS4MGDdKOHTu0du1a+fv7m44EuJ29e/dq0KBB+v777zVlyhSFh4ebjoTfYC8Ap3h5eSkxMVENGjRQ\n9+7ddeHCBdORALfx888/68knn1T79u3Vo0cP7d69myK1KcoUTvP29tasWbNUsWJF9erVS/n5+aYj\nAS6toKBA06dPV4MGDSRJX3zxheLj41WiBFfm7Io1L26Y/Px89ezZU2XKlNHChQvlwyudgGuWnp6u\nQYMGqXLlykpMTFSjRo1MR8JVoExxQ124cEFdunRRzZo1NWfOHE4YAlfpq6++0tChQ7V7925NnDhR\n9957r7y8vEzHwlXibzrcUH5+flq+fLmys7P1zDPPiM9qwJ87c+aMnnvuOYWGhqpFixbav3+/YmJi\nKFIXQ5nihitbtqxWr16tjIwMPffcc6bjALbkcDi0YMECNWjQQEeOHNGePXs0atQo+fn5mY6G68DV\nbBSJcuXKKTU1Ve3atVPZsmU1atQo05EA29i2bZsGDRqkwsJCLVmyRHfddZfpSHASZYoiU7FiRa1b\nt05t27ZV2bJl9fTTT5uOBBh1/PhxjRo1SqmpqXr99df10EMPca7ATVCmKFLVqlVTWlralUIdOHCg\n6UhAsbt48aKmTJmicePG6ZFHHtF//vMfBQYGXtPXuHTpkrKzs/X555+r9MqVujUtTd8GBqp7/fq8\np9QGKFMUuVq1amndunUKCwtTmTJldP/995uOBBQLy7K0cuVKDR48WEFBQdq6davq1av3p7/H4XDo\nyJEjysrKuvJPy0OH1HD7dqVWqaJv77lHr+/Zo9tOnFBBYaG0Zs0vv5EyNYoyRbH429/+ptTUVHXo\n0EH+/v7q2rWr6UhAkdq/f7+effZZHT16VNOnT1dUVNSv/r1lWfruu++0d+9e7du3T/v27VPQ/v26\nIzNTywIDdTg0VI0aNVLHjh3VbdYslXM4FB0SIr3/vrR4sTRvnkr87W/SgQO/TKYwijJFsWnUqJFW\nrlypzp07691331VERITpSMANl5ubq5dfflkLFy7U888/ryeeeEK5ubnasGGDLr7zjm5NS9PiMmU0\n7fvv5e3treDgYAUHB6t58+aKzcrSzYWFimrVSkpO/t8v6ucnBQb+b2nGxjKJ2gwPbUCxy3rhBX33\nP/+j6s8/r4Yvvmg6DnBDFBQUaMKECRo7dqyCgoJ022236auvvtL+/ftlWZaCgoKUdPiwGn3zjU6E\nhqpwxQpVrlz51/eTXp44uQbqeihTFLuzYWHy37hRP3l768cXX1SDF14wHQm4bgUFBXr//fe194UX\n1PbQIW2sXVunoqPVsGFDBQUFqWHDhv9bmpSl26JMUfwWL5bi46UTJ7S+ZEmVWLtWYWFhplMB1yQ/\nP18LFizQ66+/rqpVq2rphQuqtH271LHjr1e08AhcM0Xx+/+fyOfNU5UWLfT3Xr20YMGC/3NAA7Cj\nixcvat68eUpISFDdunU1Z84ctWvX7tdTJzwOkymMy8jIUExMjObMmaNu3bqZjgP8rgsXLmjOnDka\nN26cgoKCNHr0aLVu3dp0LNgEkymMu/vuu7V69Wp16dJFFy9eVK9evUxHAq44e/asZs6cqQkTJqhF\nixZaunSpWrZsaToWbIYyhS20bNlSa9euVXR0tC5evKh+/fqZjgQPd/r0aSUlJWny5Mm65557lJyc\nrCZNmpiOBZuiTGEbISEhSktLU4cOHXThwgXFxcWZjgQPlJeXp2nTpmnatGmKiIhQWlqagoKCTMeC\nzVGmsJWGDRvq448/Vnh4uC5evKj4+HjTkeAhfvrpJyUmJuqNN95Q586dtXnzZtWvX990LLgIyhS2\nU69ePW3cuFHh4eE6f/68hg4dajoS3NgPP/ygSZMmafbs2erRo4cyMzN12223mY4FF8O7f2BLderU\n0aZNmzRr1iyNGTPGdBy4oePHj2vw4MFq0KCBTp06pZ07d2r27NkUKa4LZQrbqlGjhjZu3KhFixbp\nueeeE3dx4UY4evSonnrqKQUFBcnhcCgrK0tJSUmqVauW6WhwYZQpbK1q1ar6+OOPlZycrCFDhlCo\nuG6HDx/WY489ppCQEPn5+Wn//v1KTExU9erVTUeDG6BMYXuVKlVSenq6tmzZovj4eDkcDtOR4EJy\ncnI0cOBANW/eXBUrVtSBAwc0fvx4ValSxXQ0uBHKFC6hfPnyWr9+vT7//HPFxcWpsLDQdCTY3Bdf\nfKF+/fqpVatWqlmzpnJycvTaa6+pYsWKpqPBDVGmcBmBgYFKSUnR4cOH9dBDD6mgoMB0JNjQ559/\nrtjYWLVr104NGzbUwYMH9dJLL6l8+fKmo8GNUaZwKWXLltXq1av1888/q3fv3rp06ZLpSLCJnTt3\nKiYmRlFRUQoNDdVXX32lUaNG6aabbjIdDR6AMoXLKV26tJYvX67CwkL16NFDFy5cMB0JBmVmZqpL\nly7q2rWrwsLCdPDgQQ0dOlRly5Y1HQ0ehDKFS/L19dWSJUvk7++vbt266dy5c6YjoZht3rxZkZGR\nio2NVefOnXXw4EENGjRIZcqUMR0NHohXsMGlFRQUaODAgfr666+1atUqBQQEmI6EImRZljZs2KBX\nXnlFR48e1ciRI/XQQw+pVKlSpqPBw1GmcHkOh0OPP/64srKytGbNGpUrV850JNxglmUpNTVVr776\nqk6cOKHnnntOffr0UYkSPBEV9kCZwi1YlqVBgwYpIyNDa9euVYUKFUxHwg1gWZZWrlypMWPG6OzZ\ns3r++ecVGxsrHx8f09GAX6FM4TYsy9KIESO0Zs0arVu3TrfccovpSLhODodDy5Yt05gxY2RZlkaP\nHq2YmBh5e3PMA/bEjgRuw8vLSwkJCSpdurTCwsKUlpamatWqmY6Fa1BYWKjFixfrtddeU+nSpfXy\nyy+ra9eu8vLyMh0N+FOUKdyKl5eXXnrpJfn5+aldu3ZKS0tTzZo1TcfCXyh49119P26cxn3/vT6r\nW1cTJkxQVFQUJQqXwZoXbisxMVFTpkzR+vXrea2WTZ07d05vvfWWGg4bpvYXLuhEaKgqfPopJQqX\nw2QKt/XMM8/Iz89PYWFhWr9+verXr286Ei7Lzc3V9OnTNW3aNLVu3VqRzz8vZWSo4oABEkUKF0SZ\nwq09/vjj8vPzU/v27ZWamqrg4GDTkTzasWPHNHnyZL311lvq3r27Pv74Y91xxx2mYwFO42gc3N6A\nAQM0ceJERUREaOfOnabjeKTs7GzFxcWpUaNGKiws1J49ezR37lyKFG6DyRQe4f7775evr686duyo\nFStW6M477zQdySN89tlnSkhI0KZNmxQfH68vv/ySe4DhlihTeIyYmBj5+vqqa9euWrp0qdq0aWM6\nkluyLEtpaWlKSEjQgQMHNGTIEL399tvy9/c3HQ0oMpzmhcdZv369+vTpo0WLFik8PNx0HLdRWFio\nZcuWKSEhQefOndPw4cP1wAMP8NxceATKFB5p8+bN6tmzp+bNm6dOnTqZjuPSLl68qAULFmjcuHEq\nV66cRo4cqW7duvG0IngUyhQe69NPP1X37t01Y8YMxcTEmI7jck6fPq1Zs2Zp8uTJCg4O1ogRI9Su\nXTvuEYVH4qMjPFarVq20Zs0aJQ8YoCPBwSpctMh0JJfw448/avTo0apbt662b9+ulStXKiUlRWFh\nYRQpPBaTKTxbp04qXLNGPpLWlSiht+67T9HR0YqKilKVKlVMp7OVw4cPa+LEiVq4cKFiY2M1dOhQ\n3X777aZjAbbAZArPtn69fCQVSmo2daoiIiK0atUq3XHHHWratKlGjhypjRs36tKlS6aTGrN37149\n+OCDat68ufz9/bVv3z7NmDGDIgX+C5MpPNpPrVqpbGamfKKjVWLNmiu/XlBQoMzMTKWmpiolJUXZ\n2dn6+9//fmVqrVOnjsHUxeOTTz5RQkKCtm/frkGDBunxxx/nxevAH6BM4dEiIiLUp08fDRw48E//\nux9//FHr1q1TSkqKUlNTVb58eUVFRSk6Olrt2rVTmTJliilx0bIsS2vWrFFCQoKOHTumYcOGqX//\n/ipdurTpaICtUabwWBkZGXrwwQeVnZ2tkiVLXvXvczgc2rNnj1JSUpSSkqKdO3eqdevWio6OVnR0\ntBo0aOByB3EKCgq0ePFiJSQkyNvbWyNGjNB9992nEiV4rgtwNShTeKzIyEjFxsYqLi7Oqa9z8uRJ\npaenXylXSVem1vDwcN100003Im6ROH/+vObOnavx48erZs2aGjFihKKjo13uwwBgGmUKj7R161Y9\n8MADOnDgwA19Qo9lWcrOzr5SrBkZGWratOmVa61Nmza1xcMM8vLylJSUpKlTp+rOO+/U8OHD1bp1\na9OxAJdFmcIjRUdHKyYmRo899liRfp9z585p06ZNV661/vzzz4qMjFR0dLQiIyNVqVKlIv3+v3X8\n+HElJiZqzpw56tKli/71r38pKCioWDMA7ogyhcfJzMxUr169lJOTU+zPjT18+PCVE8IbNmzQ7bff\nfuVaa6tWrYrsGmVOTo7Gjx+vJUuW6MEHH9TgwYNVq1atIvlegCeiTOFxOnXqpK5du+qJJ54wmuPS\npUvaunXrlan10KFDCg8Pv7ISvvXWW53+Hjt37tTYsWOVnp6uf/7zn3rqqadUsWLFG5AewH+jTOFR\ntm/frh49eignJ0e+vr6m4/zKd999p7Vr1yolJUVr167VLbfccmVqbdOmjfz8/K7q61iWpY8//lgJ\nCQnav3+/Bg8erH/84x8qW7ZsEf8JAM9FmcKjdO3aVdHR0YqPjzcd5U8VFhZqx44dV6bWrKwstWnT\n5srUWq9evf9z4tbhcOijjz5SQkKCTp48qeHDh6tv3768Ag0oBpQpPMaOHTvUvXt35eTkXPWUZxe5\nublav379lVPCvr6+V6bWu+++WytWrNC4ceMUEBCgkSNHqnv37rY4NQx4CsoUHqN79+6KiIjQU089\nZTqKUyzL0t69e69MrdW2bNFjfn4KfPppBb/8MveIAgZQpvAIu3btUpcuXXTw4EGXm0r/ymovL3WW\npI4dpeRk03EAj8QeCB7hlVde0bBhw9yuSH/66SfNk3SxfXtpwADDaQDPRZnC7R0eP17xycl6wg1v\nCTn1yCOaLsknKEiKjTUdB/BYrHnh9nZWrapm333nlmvQs/7+8j93TqpYUfrxR9NxAI/FZAq3duDA\nAU0/e1b5kZFuuQbNLVtWhZJUo4bpKIBHYzKFW4uLi1ONGjX00ksvmY5SJPLLl1fJvDwmU8AwXlYI\nt/XNN9/oww8/1Jdffmk6SpH5LjxcZT/6SOX79jUdBfBolCnc1sSJE/Xwww+rQoUKpqMUmR39+mlu\nfr4+Skw0HQXwaJQp3NKJEyf09ttvKysry3SUIpWXl6dy5cqZjgF4PA4gwS1NnTpVvXr1UvXq1U1H\nKVKUKWAPTKZwO6dOnVJSUpIyMzNNRylylClgD0ymcDszZ85UZGSkbrvtNtNRilxubi5lCtgAkync\nyoULFzRp0iSlpqaajlIs8vLyVL58edMxAI/HZAq3MnfuXLVo0UKNGzc2HaVYsOYF7IHJFG6joKBA\n48aN07vvvms6SrGhTAF7YDKF23jvvfdUq1Yt3XXXXaajFBvKFLAHJlO4BYfDoYSEBE2aNMl0lGLF\nASTAHphM4RZWrlwpX19fdejQwXSUYsUBJMAeKFO4PMuydHTiRK2yLHktWWI6TrEpKCjQ2bNnFRAQ\nYDoK4PEoU7i8zZs3K3jHDlXdtUuaN890nGJz6tQpBQYGytubH2PANH4K4fLGjRun/D59fnn5txu+\ns/SPcPgIsA8OIMGl7d27V5999pnaHD4s+fmZjlOsOHwE2AeTKVzahAkT9NRTT8nPw4pU4vARYCdM\npnBZ33zzjVasWKGcnBzTUYxgzQvYB5MpXNaUKVPUv39/3XzzzaajGEGZAvbBZAqXdPLkSb355pva\ntWuX6SjGUKaAfTCZwiXNnDlTnTp1Uq1atUxHMYYDSIB9MJnC5Vy8eFFTpkxRcnKy6ShG1dm+XeFH\nj0oNGkixsabjAB6NyRQuZ+HChQoODlZISIjpKEaF7t+v2l984VEPqgDsiskULsXhcGj8+PH697//\nbTqKcasrVVLlypV1iwc9qAKwKyZTuJTVq1erTJkyat++vekoxq3291f25MmseAEboEzhUo5OnKhV\nDodHPdD+j3CaF7APyhQuY9u2bWqQmamqu3dznVC/lOlNN91kOgYAcc0ULmTSpEnq16uXdOKERz3Q\n/o+cPHngGpbXAAAEfUlEQVSSyRSwCS/LsizTIYC/8vXXX6tZs2Y6dOiQAgMDTccxrqCgQL6+vsrP\nz+cVbIAN8FMIlzBlyhQNHDiQIr3s1KlTCggIoEgBm2DNC9s7efKk3n77be3evdt0FNtgxQvYCx9r\nYXtz5sxRdHS0br31VtNRbIOTvIC9MJnC1vLz8zVlyhQtW7bMdBRb4SQvYC9MprC1Dz74QHXr1lXz\n5s1NR7EV1ryAvVCmsC3LsjRx4kQNGTLEdBTbYc0L2AtrXthW9quvakpOju46e9Z0FNthzQvYC5Mp\nbOvM9Om6++RJec+fbzqK7bDmBeyFMoUt5eTk6I3z51UQGcnTjn4Ha17AXihT2NLUqVNVOT5eJVJT\neSvK72DNC9gL10xhO3l5eVqwYIGysrJMR7Et1ryAvTCZwnbefPNNdezYUdWrVzcdxbZY8wL2wmQK\nWykoKNDUqVO1dOlS01FsjTUvYC9MprCVnSNG6P3Tp9Xiq69MR7E11ryAvVCmsJXCN99Uq9xcXv79\nF1jzAvZCmcI2tm/frndKlJAjOprbYf6EZVk6efIka17ARihT2MbWZ5/V0AoV5P3ww9wO8yfOnDkj\nX19flSxZ0nQUAJdRprCNmFOnVDc7mxXvX+B6KWA/lCls49bnn5c6dmTF+xe4XgrYD2UK+4iNlZKT\nWfH+hZLLl+vdvDxp8WLTUQBcRpkCLqb88uUK+fZb1uGAjVCmgIs50Lq1dlWtyjocsBHKFHAxnzdo\noJndurEOB2yEMgVcDKd5AfuhTAEXw3N5AfuhTAEXw9OPAPuhTAEXw5oXsB/KFHAxrHkB+6FMARfD\nZArYD2UKuBiumQL2Q5kCLoY1L2A/lCngYljzAvbjZVmWZToEgKuTn58vPz8/FRQUyMvLy3QcAJcx\nmQIu5NSpUwoMDKRIAZuhTAEXwooXsCfKFHAhHD4C7IkyBVwIkylgT5Qp4EKYTAF7okwBF8IDGwB7\nokwBF8KaF7AnyhRwIax5AXsqYToAgKtXe9s2tTl4UGrUSIqNNR0HwGU8AQlwISk+Pop2OKSOHaXk\nZNNxAFzGmhdwIXMdDm0JDJQGDDAdBcB/YTIFXIi3t7ciIyOVkpJiOgqA/8JkCrgQy7JUvXp10zEA\n/AaTKeBCdnh5qakk7yZNpF27TMcBcBllCrgSb2/JsiQvL8nhMJ0GwGWseQFXEhLyS5GGhJhOAuC/\nMJkCAOAkJlMAAJxEmQIA4CTKFAAAJ1GmAAA4iTIFAMBJlCkAAE6iTAEAcBJlCgCAkyhTAACcRJkC\nAOAkyhQAACdRpgAAOIkyBQDASZQpAABOokwBAHASZQoAgJMoUwAAnESZAgDgJMoUAAAnUaYAADiJ\nMgUAwEmUKQAATqJMAQBwEmUKAICTKFMAAJxEmQIA4CTKFAAAJ1GmAAA4iTIFAMBJlCkAAE6iTAEA\ncBJlCgCAkyhTAACcRJkCAOAkyhQAACdRpgAAOIkyBQDASZQpAABOokwBAHASZQoAgJMoUwAAnESZ\nAgDgJMoUAAAnUaYAADiJMgUAwEn/DxCHAqlL0Mf5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113c088d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "from spanning_tree import maximum_spanning_tree # script spanning_tree.py\n",
    "\n",
    "G_max_spanning_tree = nx.from_numpy_matrix(cancerCL_corr) # total graph\n",
    "G_max_spanning_tree = maximum_spanning_tree(G_max_spanning_tree) # graph mst (minimum spanning tree)\n",
    "nx.draw_spectral(G_max_spanning_tree, node_size=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "58"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(G_max_spanning_tree.edges())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Analysis of the graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def analyze_graph(G):\n",
    "\n",
    "    d = list(nx.connected_component_subgraphs(G))\n",
    "    n_nodes_in_subgraphs = [len(g.nodes()) for g in d]\n",
    "    n_disconnected_nodes = np.sum(np.array(n_nodes_in_subgraphs)==1)\n",
    "    \n",
    "    print('Results:\\n-------------------------------------------------------------------')\n",
    "    if n_disconnected_nodes==0:\n",
    "        print('All %d nodes are connected.' %(len(G.nodes())))\n",
    "    else:\n",
    "        print('%d disconnected subgraphs' %len(d))\n",
    "        print('Minimum number of nodes in a subgraph: %d' %np.min(n_nodes_in_subgraphs))\n",
    "        print('Maximum number of nodes in a subgraph: %d' %np.max(n_nodes_in_subgraphs))\n",
    "        print('%d nodes (ouf of %d, i.e. %.2f%%) are singletons (i.e. disconnected from the graph)' %(n_disconnected_nodes, len(G.nodes()), n_disconnected_nodes/len(G.nodes())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results:\n",
      "-------------------------------------------------------------------\n",
      "25 disconnected subgraphs\n",
      "Minimum number of nodes in a subgraph: 1\n",
      "Maximum number of nodes in a subgraph: 35\n",
      "24 nodes (ouf of 59, i.e. 0.41%) are singletons (i.e. disconnected from the graph)\n"
     ]
    }
   ],
   "source": [
    "# Threshold\n",
    "analyze_graph(G_correl_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results:\n",
      "-------------------------------------------------------------------\n",
      "All 59 nodes are connected.\n"
     ]
    }
   ],
   "source": [
    "# Maximum weight spanning tree\n",
    "analyze_graph(G_max_spanning_tree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Comment:\n",
    "\n",
    "The maximum spanning tree method seems more adequate as all nodes of the *cancer cell lines* graph are connected, which is not the case when thresholding correlation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. MMCRF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Theoretical problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See paper **Efficient Algorithms for Max-Margin Structured Classification, Juho Rousu et al.**\n",
    "\n",
    "Explains very well the optimization problem to be resolved in section 1.2.6 **Marginal dual problem**:\n",
    "\n",
    "$$ \\max_{\\mu}\\mu^{T} l_{H} − \\frac{1}{2}\\mu^{T}K_{H}\\mu $$\n",
    "\n",
    "Then the gradient becomes: \n",
    "\n",
    "$$ g = l_{H} − K_{H}\\mu $$\n",
    "\n",
    "Furthermore, the objective can be written as :  \n",
    "$$obj = \\mu^{T} l_{H} − \\frac{1}{2}\\mu^{T}K_{H}\\mu$$\n",
    "$$obj = \\sum_{i=1}^{m}\\sum_{e\\in E}\\sum_{u\\in Y_{E}} \\mu(e,u) l_{H}(i,u) − \\frac{1}{2}\\sum_{e\\in E}\\sum_{i,i^{'}}\\sum_{u,u^{'}} \\mu_{e}(i,u)K_{e}(i,u;i^{'},u^{'})\\mu_{e}(i^{'},u^{'})$$\n",
    "\n",
    "That's why we need to construct the full loss vector (for each of the 4 possible labels), and $K_{mu_X}$ and $K_mu$ for the gradient descent.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def invert(X):\n",
    "    return [not x for x in X]\n",
    "\n",
    "class model_MMCRF:\n",
    "    \n",
    "    \"\"\"\n",
    "    Implementation of the MMCRF model, which makes use of the structure of the output and of the input to obtain\n",
    "    better performances on structured learning tasks.\n",
    "    \n",
    "    Python adaptation of the Matlab implementation of Rousu et al. in their article (cited above).\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self,params):\n",
    "        \"\"\"\n",
    "        Initialization of hyperparameters (there is a class for them)\n",
    "        \"\"\"\n",
    "        self.params = params\n",
    "        \n",
    "\n",
    "    def profile_init(self):\n",
    "        \"\"\"\n",
    "        Initialization of the profile --> useful to monitor the learning phase\n",
    "        \"\"\"\n",
    "        self.start_time = time.time()\n",
    "        self.next_profile_tm = self.start_time\n",
    "        self.n_err = 0\n",
    "        self.p_err = 0\n",
    "        self.n_err_microlbl = 0\n",
    "        self.p_err_microlbl = 0 \n",
    "        self.n_err_microlbl_prev = 0\n",
    "        self.microlabel_errors = []\n",
    "        self.it = 0\n",
    "        self.err_ts = 0\n",
    "        self.metrics_tr = []\n",
    "        print('alg: M3LBP\\ttm:\\titer:\\tobj:\\tmu:max\\tmu:min\\tdgap:\\ttrain\\tacc:\\tstd:\\tf1:')\n",
    "\n",
    "    def profile_update(self):\n",
    "        \"\"\"\n",
    "        Update of the profile\n",
    "        \"\"\"\n",
    "        tm = time.time()\n",
    "        ' We compute the current training loss '    \n",
    "        m_tr, s_tr, f1_tr = self.assess_performance(self.Kx_tr, self.Y_tr)\n",
    "        print('\\t\\t%d\\t%d\\t%.1e\\t%.1e\\t%.1e\\t%.1e\\t\\t%.2f%%\\t%.2f%%\\t%.3f'\n",
    "            %(np.round(tm-self.start_time), self.it, self.obj, np.max(self.mu), np.min(self.mu), self.primal_ub-self.obj, m_tr*100, s_tr*100, f1_tr))\n",
    "        self.metrics_tr.append([self.it, m_tr, s_tr, f1_tr])\n",
    "        \n",
    "        \n",
    "    def var_init(self):\n",
    "        self.MBProp = np.zeros((2*self.e, 2*self.e)) # for edge to edge propagation\n",
    "        self.MBPropEdgeNode = np.zeros((2*self.e, self.l)) # for edge to node propagation\n",
    "        self.Rmu = []\n",
    "        self.Smu = []\n",
    "        self.term12 = np.zeros((1, self.e))\n",
    "        self.term34 = np.zeros((4, self.e))\n",
    "        \n",
    "    \n",
    "    def learn(self, Kx_tr, Y_tr, G):\n",
    "        \"\"\"\n",
    "        Input data assumed by the algorithm\n",
    "\n",
    "        Kx_tr        X-kernel, assume to be positive semidefinite and normalized (Kx_tr(i,i) = 1)\n",
    "        Y_tr         Y-data: assumed to be class labels encoded {-1,+1}\n",
    "        G            graph of the Markov network e_i = [E(i,1),E(i,2)]\n",
    "\n",
    "        params       parameters used by the learning algorithm\n",
    "\n",
    "        loss         losses associated with different edge labelings\n",
    "        mu           marginal dual variables: these are the parameters to be learned\n",
    "\n",
    "        m            number of training instances\n",
    "        l            number of labels\n",
    "        Ye           Denotes the edge-labelings 1 <-- [-1,-1], 2 <-- [-1,+1], 3 <-- [+1,-1], 4 <-- [+1,+1]\n",
    "        IndEdgeVal   IndEdgeVal{u} = [Ye == u] \n",
    "        Kmu          Kx_tr*mu\n",
    "\n",
    "        primal_ub\n",
    "        profile\n",
    "        obj\n",
    "        params\n",
    "        opt_round\n",
    "        \"\"\"\n",
    "\n",
    "        \"\"\"\n",
    "        E            edges of the Markov network e_i = [E(i,1),E(i,2)]\n",
    "        \"\"\"\n",
    "        \n",
    "        self.Kx_tr = Kx_tr\n",
    "        self.Y_tr = Y_tr\n",
    "        self.G = G\n",
    "        self.E = self.G.edges() \n",
    "        self.e = len(self.E)\n",
    "\n",
    "        self.l = self.Y_tr.shape[1]\n",
    "        self.m = self.Kx_tr.shape[0]\n",
    "        self.mu = np.zeros((4*self.e,self.m))\n",
    "        \n",
    "        self.var_init()\n",
    "        self.profile_init()\n",
    "        \n",
    "        self.compute_loss_vector(self.Y_tr) # automatic scaling\n",
    "\n",
    "        ' Matrices for speeding up gradient computations '\n",
    "        self.Ye = np.resize(self.loss == 0, (4,self.e*self.m))\n",
    "\n",
    "        self.IndEdgeVal = []\n",
    "        for u in range(4):\n",
    "            self.IndEdgeVal.append(np.resize(self.Ye[u,:]!=0, (self.e,self.m)))\n",
    "\n",
    "        ' Initialization of variables '\n",
    "\n",
    "        self.Ye = np.resize(self.Ye, (4*self.e,self.m))\n",
    "        self.Kxx_mu_x = np.zeros((4*self.e,self.m))\n",
    "        self.Kmu = np.zeros((4*self.e*self.m,1))\n",
    "        \n",
    "        ' Starting descent... '\n",
    "        self.obj = 0; self.primal_ub = float('inf'); self.it = 0; self.opt_round = 1;\n",
    "        self.profile_update()\n",
    "        self.prev_obj = 0\n",
    "\n",
    "        ' Repeat until working set converged and close to optima '\n",
    "        while (self.primal_ub - self.obj > self.params.epsilon*self.obj):\n",
    "\n",
    "            self.progress_made = 0\n",
    "            ' Conditional gradient optimization... '\n",
    "            for x in range(self.m):\n",
    "                ' obtain initial gradient for x : Kmu_x '\n",
    "                self.Kmu_x = self.compute_Kmu_x(x, self.Kx_tr[:,x])\n",
    "\n",
    "                ' conditional gradient optimization on x '\n",
    "                x_iter = self.optimize_x(x)\n",
    "\n",
    "                self.it += x_iter\n",
    "\n",
    "                if self.params.verbosity==2:\n",
    "                    self.profile_update()\n",
    "\n",
    "\n",
    "            ' Current full gradient... '\n",
    "            self.Kmu = self.compute_Kmu(self.Kx_tr, self.mu)\n",
    "\n",
    "            self.progress_made = (self.obj > self.prev_obj)\n",
    "            if self.progress_made:\n",
    "                self.prev_obj = self.obj.copy()\n",
    "                self.prev_mu = self.mu.copy()\n",
    "            else: # restore previous solution and finish\n",
    "                self.mu = self.prev_mu.copy()\n",
    "                self.Kmu = self.compute_Kmu(self.Kx_tr, self.Smu, self.Rmu, self.mu)\n",
    "\n",
    "                self.obj = self.mu.flatten().T.dot(self.loss.flatten()) - self.mu.flatten().T.dot(self.Kmu.flatten())/2 \n",
    "\n",
    "            ' Duality gap and primal upper bound '\n",
    "            self.compute_duality_gap()\n",
    "\n",
    "            self.profile_update()\n",
    "\n",
    "            self.opt_round = self.opt_round + 1;\n",
    "            if np.logical_or(self.opt_round > self.params.max_iter, self.progress_made == 0):\n",
    "                self.next_profile_tm = 0\n",
    "                self.profile_update()\n",
    "                break\n",
    "                \n",
    "        print('Descent done.')\n",
    "        \n",
    "    \n",
    "    def compute_loss_vector(self, Y):\n",
    "        \"\"\"\n",
    "        LOSS VECTOR\n",
    "        Compute loss vector to be used in the conditional gradient optimization\n",
    "        \"\"\"\n",
    "\n",
    "        labels = [-1,1]\n",
    "\n",
    "        \"\"\"Initialize loss of size (4, M*E):\n",
    "            M: nb of melecules\n",
    "            E: nb of edges to be predicted for each molecule\n",
    "            4: nb of losses to be computed (for each of the four possible labeling we compute a different loss)\"\"\"\n",
    "        self.loss = np.ones((4,self.m*self.e))\n",
    "\n",
    "        self.E_tail = [elt[0] for elt in self.E]\n",
    "        self.E_head = [elt[1] for elt in self.E]\n",
    "        self.Te1 = Y[:,self.E_tail] # the label of edge tail\n",
    "        self.Te2 = Y[:,self.E_head] # the label of edge head\n",
    "\n",
    "        \"\"\"Scaling the loss : rescale to microlabels by dividing node loss among the adjacent edges\"\"\"\n",
    "        NodeDegree = np.ones((Y.shape[1],1))\n",
    "        for v in range(self.Y_tr.shape[1]):\n",
    "            NodeDegree[v] = np.sum(np.array(self.E).flatten()==v)\n",
    "\n",
    "        NodeDegree = np.hstack(np.repeat(np.expand_dims(NodeDegree, axis = 0), self.m, axis=0)).T\n",
    "\n",
    "        u = 0\n",
    "        for u_1 in labels:\n",
    "            for u_2 in labels:\n",
    "                self.loss[u,:] = np.resize(np.divide((self.Te1 != u_1), NodeDegree[:, self.E_tail]) + np.divide((self.Te2 != u_2), NodeDegree[:, self.E_head]), (self.e*self.m,1))[:,0]\n",
    "                u += 1\n",
    "        self.loss = np.resize(self.loss, (4*self.e,self.m))\n",
    "\n",
    "    \n",
    "    def compute_duality_gap(self):\n",
    "        \"\"\"\n",
    "        DUALITY GAP\n",
    "        \"\"\"\n",
    "\n",
    "        self.loss = np.resize(self.loss, (4,self.e*self.m))\n",
    "        self.Kmu = np.resize(self.Kmu, (4,self.e*self.m))\n",
    "        self.mu = np.resize(self.mu, (4,self.e*self.m))\n",
    "\n",
    "        gradient = self.loss - self.Kmu\n",
    "\n",
    "        self.dgap = float('inf'); LBP_iter = 1; Gmax = -float('inf')\n",
    "\n",
    "        while LBP_iter <= self.e:\n",
    "\n",
    "            LBP_iter *= 2 # no of iterations = diameter of the graph\n",
    "            \n",
    "            Ymax, max_gradient = self.max_gradient_labeling(gradient, self.params) # \"max_gradient\" not to confuse with G\n",
    "            ' \"max_gradient\" of size (4*e,m) '\n",
    "            Gmax = np.maximum(Gmax, max_gradient)\n",
    "            ' C: SVM hyperparameter '\n",
    "            duality_gap = self.params.C*np.maximum(Gmax,0) - np.sum(np.resize(np.sum(gradient*self.mu.reshape(gradient.shape), 0), (self.e,self.m)), 0)\n",
    "            self.dgap = np.sum(duality_gap)\n",
    "\n",
    "            if self.obj+self.dgap < self.primal_ub+1e-6:\n",
    "                break\n",
    "\n",
    "        if self.primal_ub == float('inf'): # initial value\n",
    "             self.primal_ub = self.obj+self.dgap;\n",
    "        else:\n",
    "             self.primal_ub = (self.obj+self.dgap)/min(self.opt_round,10) + self.primal_ub*(1-1/min(self.opt_round,10)); # averaging over a few last rounds\n",
    "\n",
    "        self.loss = np.resize(self.loss, (4*self.e,self.m))\n",
    "        self.Kmu = np.resize(self.Kmu, (4*self.e,self.m))\n",
    "        self.mu = np.resize(self.mu, (4*self.e,self.m))\n",
    "\n",
    "    \n",
    "    def Y2U(self, E, Y_tr):\n",
    "        \"\"\"\n",
    "        CONVERT NODE LABELING TO EDGE LABELING\n",
    "        \"\"\"\n",
    "        \n",
    "        m = Y_tr.shape[0]\n",
    "        e = len(E)\n",
    "\n",
    "        U = np.zeros((4,e*m))\n",
    "        \n",
    "        E_tail = [elt[0] for elt in E]\n",
    "        E_head = [elt[1] for elt in E]\n",
    "        \n",
    "        Te1 = Y_tr[:,E_tail] # the label of edge tail\n",
    "        Te2 = Y_tr[:,E_head] # the label of edge head\n",
    "\n",
    "        U[0,:] = np.resize(np.logical_and(Te1 == -1, Te2 == -1), (1,e*m))\n",
    "        U[1,:] = np.resize(np.logical_and(Te1 == -1, Te2 == -1), (1,e*m))\n",
    "        U[2,:] = np.resize(np.logical_and(Te1 == 1, Te2 == 1), (1,e*m))\n",
    "        U[3,:] = np.resize(np.logical_and(Te1 == 1, Te2 == 1), (1,e*m))\n",
    "\n",
    "        return U\n",
    "    \n",
    "    \n",
    "    def max_gradient_labeling_brute_force(self, Y_tr, gradient):\n",
    "        \"\"\"\n",
    "        Brute force inference of max. gradient labeling\n",
    "        \"\"\"\n",
    "        \n",
    "        m = int(gradient.size/(4*self.e))\n",
    "        \n",
    "        Yall = np.zeros((2.**l,self.l))\n",
    "\n",
    "        for j in range(l):\n",
    "            tmp = np.concatenate((-np.ones((2**j,1)),np.ones((2**j,1))))\n",
    "            Yall[:,j] = np.concatenate(np.repeat(np.expand_dims(tmp, 0), 2**(l-(j+1)), 0))[:,0]\n",
    "\n",
    "        Yall = np.fliplr(Yall)\n",
    "\n",
    "        Uall = np.resize(self.Y2U(self.E, Yall), (4*self.e,2**l))\n",
    "        gradient = np.resize(gradient, (4*self.e,m))\n",
    "\n",
    "        tmp = gradient.T.dot(Uall)\n",
    "\n",
    "        Gmax = np.amax(tmp, 0)\n",
    "        Imax = np.argmax(tmp, 0)\n",
    "        Ymax = Yall[Imax,:]\n",
    "\n",
    "        return Ymax, Gmax\n",
    "    \n",
    "    \n",
    "    def max_gradient_labeling(self, gradient, params):\n",
    "        \"\"\"\n",
    "        MAX GRAD. LABELING VIA LOOPY BELIEF PROPAGATION\n",
    "        Find maximum gradient labeling for a single example using Loopy Belief Propagation\n",
    "\n",
    "        gradient is of size (4*e,m) and contains the gradient (m-vector) for each edge-labeling (4*e in total)\n",
    "        MBProp (2*e,2*e) direction-specific adjacency matrix\n",
    "        \"\"\"\n",
    "        \n",
    "        m = int(gradient.size/(4*self.e))\n",
    "\n",
    "        if params.debugging==1:\n",
    "            Ymax, Gmax = self.max_gradient_labeling_brute_force(gradient)\n",
    "            \n",
    "        else:\n",
    "            \n",
    "            ineg = 1; ipos = 2;\n",
    "\n",
    "            if np.all(self.MBProp==np.zeros((2*self.e, 2*self.e))): # MBProp is set to its initial values\n",
    "                self.MBProp, self.MBPropEdgeNode = self.buildBeliefPropagationMatrix(self.l)\n",
    "\n",
    "            gradient = np.resize(gradient, (4,self.e*m))\n",
    "\n",
    "            ' Edge-labeling specific gradient matrices (m,e) '\n",
    "            Gnn = np.resize(gradient[0,:], (self.e,m)).T # edge-gradients for labeling [-1,-1]\n",
    "            Gnp = np.resize(gradient[1,:], (self.e,m)).T # edge-gradients for labeling [-1,+1]\n",
    "            Gpn = np.resize(gradient[2,:], (self.e,m)).T # edge-gradients for labeling [+1,-1]\n",
    "            Gpp = np.resize(gradient[3,:], (self.e,m)).T # edge-gradients for labeling [+1,+1]\n",
    "\n",
    "            ' SumMsg_*_*: mx|E| matrices storing the sums of neighboring messages from '\n",
    "            ' the head and tail of the edge, respectively, on the condition that '\n",
    "            ' the head (resp. tail) is labeled with -1 --> neg or +1 --> pos. '\n",
    "            SumMsg_head_neg = np.zeros((m,self.e))\n",
    "            SumMsg_head_pos = np.zeros((m,self.e))\n",
    "            SumMsg_tail_neg = np.zeros((m,self.e))\n",
    "            SumMsg_tail_pos = np.zeros((m,self.e))\n",
    "\n",
    "            iTail = np.arange(self.e)\n",
    "            iHead = np.arange(self.e,2*self.e)\n",
    "\n",
    "            ' Iterate until messages have had time to go accros the whole graph: at '\n",
    "            ' most this takes O(e) iterations (i.e. when the graph is a chain) '\n",
    "            for it in range(params.max_iter):\n",
    "                ' find max-gradient configuration and propage gradient value over the edge '\n",
    "                Msg_head_neg = np.maximum(SumMsg_tail_pos+Gpn, SumMsg_tail_neg+Gnn)\n",
    "                Msg_head_pos = np.maximum(SumMsg_tail_pos+Gpp, SumMsg_tail_neg+Gnp)\n",
    "                Msg_tail_neg = np.maximum(SumMsg_head_pos+Gnp, SumMsg_head_neg+Gnn)\n",
    "                Msg_tail_pos = np.maximum(SumMsg_head_pos+Gpp, SumMsg_head_neg+Gpn)\n",
    "\n",
    "                ' Sum up gradients of consistent configurations and propage to neighboring edges '\n",
    "                SumMsg_tail_neg = np.hstack((Msg_tail_neg, Msg_head_neg)).dot(self.MBProp[:,iTail])\n",
    "                SumMsg_tail_pos = np.hstack((Msg_tail_pos, Msg_head_pos)).dot(self.MBProp[:,iTail])\n",
    "                SumMsg_head_neg = np.hstack((Msg_tail_neg, Msg_head_neg)).dot(self.MBProp[:,iHead])\n",
    "                SumMsg_head_pos = np.hstack((Msg_tail_pos, Msg_head_pos)).dot(self.MBProp[:,iHead])\n",
    "\n",
    "\n",
    "            ' find out the labeling: sum up the edge messages coming towards each node '\n",
    "            M_max1 = np.hstack((Msg_tail_neg, Msg_head_neg)).dot(self.MBPropEdgeNode)\n",
    "            M_max2 = np.hstack((Msg_tail_pos, Msg_head_pos)).dot(self.MBPropEdgeNode)\n",
    "            ' pick the label of maximum message value '\n",
    "            Ymax = (M_max1 <= M_max2)*2 - 1\n",
    "\n",
    "            ' find out the max gradient for each example: pick out the edge labelings consistent with Ymax '        \n",
    "            E_tail = [elt[0] for elt in self.E]\n",
    "            E_head = [elt[1] for elt in self.E]\n",
    "            \n",
    "            Te1 = Ymax[:,E_tail] # the label of edge tail\n",
    "            Te2 = Ymax[:,E_head]\n",
    "\n",
    "            Umax = np.zeros((4, self.e*m))\n",
    "\n",
    "            Umax[0,:] = np.resize(np.logical_and(Te1 == -1, Te2 == -1).T, (1,self.e*m))\n",
    "            Umax[1,:] = np.resize(np.logical_and(Te1 == -1, Te2 == 1).T, (1,self.e*m))\n",
    "            Umax[2,:] = np.resize(np.logical_and(Te1 == 1, Te2 == -1).T, (1,self.e*m))\n",
    "            Umax[3,:] = np.resize(np.logical_and(Te1 == 1, Te2 == 1).T, (1,self.e*m))\n",
    "            ' sum up the corresponding edge-gradients '\n",
    "            Gmax = np.resize(np.sum(gradient*Umax, 0), (self.e,m))\n",
    "            Gmax = np.resize(np.sum(Gmax, 0), (m,1))\n",
    "\n",
    "            gradient = np.resize(gradient, (4*self.e,m))\n",
    "\n",
    "        return Ymax, Gmax\n",
    "    \n",
    "    \n",
    "    def buildBeliefPropagationMatrix(self, n_nodes):\n",
    "        \"\"\"\n",
    "        Construct a matrix containing the neighborhood information of the edges.\n",
    "        The matrix consists of four blocks, corresponding to the edges that merge \n",
    "        (e(2) = e'(2)), branch (e(1) = e'(1)), form a chain forward (e(2) = e'(1)) or backward (e(1) = e'(2))\n",
    "        \"\"\"\n",
    "    \n",
    "        E_tail = [elt[0] for elt in self.E]\n",
    "        E_head = [elt[1] for elt in self.E]\n",
    "\n",
    "        MBProp = np.zeros((2*self.e, 2*self.e)) # for edge to edge propagation\n",
    "        MBPropEdgeNode = np.zeros((2*self.e, n_nodes)) # for edge to node propagation\n",
    "\n",
    "        iTail = np.arange(self.e)\n",
    "        iHead = np.arange(self.e,2*self.e)\n",
    "\n",
    "        for node in range(n_nodes):\n",
    "            eTail = np.where(np.asarray(E_tail) == node)[0]\n",
    "            eHead = np.where(np.asarray(E_head) == node)[0]\n",
    "\n",
    "            ' Edges that meet node with their tail '\n",
    "            MBPropEdgeNode[iTail[eTail], node] = 1;\n",
    "            ' Edges that meet node with the head '\n",
    "            MBPropEdgeNode[iHead[eHead], node] = 1;\n",
    "\n",
    "            ' Matrix block for progating messages from edges that meet with their tails at node (eTail) '\n",
    "            Link = MBProp[np.ix_(iTail,iTail)]\n",
    "            Link[np.ix_(eTail,eTail)] = 1\n",
    "            ' remove diagonal; we do not propage messages back to self '\n",
    "            MBProp[np.ix_(iTail, iTail)] = Link - np.diag(np.diag(Link)) # check\n",
    "\n",
    "            ' Matrix block for progating messages via a backward chain (eTail meeting eHead) at node '\n",
    "            ' messages will go from iTail to iTail (excluding self loops) '\n",
    "            Link = MBProp[np.ix_(iTail, iHead)]\n",
    "            Link[np.ix_(eTail, eHead)] = 1\n",
    "            ' remove diagonal; in case there are self loops e = (v,v) in the graph '\n",
    "            MBProp[np.ix_(iTail, iHead)] = Link - np.diag(np.diag(Link))\n",
    "\n",
    "            ' Matrix block for progating messages from edges that meet with their heads at node (eHead) '\n",
    "            Link = MBProp[np.ix_(iHead, iHead)]\n",
    "            Link[np.ix_(eHead, eHead)] = 1 \n",
    "            ' remove diagonal; we do not propage messages back to self '\n",
    "            MBProp[np.ix_(iHead, iHead)] = Link - np.diag(np.diag(Link))\n",
    "\n",
    "            ' Matrix block for progating messages  via a forward chain (eHead meeting eTail) at node; '\n",
    "            Link = MBProp[np.ix_(iHead, iTail)]\n",
    "            Link[np.ix_(eHead, eTail)] = 1\n",
    "            ' remove diagonal; in case there are self loops e = (v,v) in the graph '\n",
    "            MBProp[np.ix_(iHead, iTail)] = Link - np.diag(np.diag(Link))\n",
    "\n",
    "        return MBProp, MBPropEdgeNode\n",
    "    \n",
    "    \n",
    "    def optimize_x(self, x):\n",
    "        \"\"\"\n",
    "        Conditional gradient optimizer for a single example\n",
    "        \"\"\"\n",
    "        \n",
    "        it = 0\n",
    "        while it < self.params.max_iter:\n",
    "            gradient = self.loss[:,x] - self.Kmu_x[:,0]\n",
    "            if np.linalg.norm(gradient) < self.params.tolerance:\n",
    "                break\n",
    "\n",
    "            ' find maximum gradient labeling (worst margin violator) '\n",
    "            Ymax, Gmax = self.max_gradient_labeling(gradient, self.params)\n",
    "\n",
    "            ' gradient towards zero '\n",
    "            G0 = -self.mu[:,x].T.dot(gradient)\n",
    "\n",
    "            ' convert to update direction '\n",
    "            Umax_e = 2*(Ymax[:,self.E_tail]>0) + (Ymax[:,self.E_head]>0)\n",
    "\n",
    "            mu_1 = np.zeros(self.mu[:,x].shape)\n",
    "            if np.all(Gmax > np.maximum(self.params.tolerance, G0)):\n",
    "                for u in range(4):\n",
    "                    mu_1[4*(np.arange(1,self.e+1))-4 + u] = self.params.C*(Umax_e == u)\n",
    "\n",
    "                if np.sum(mu_1) > 0:\n",
    "                    smu_1_te = np.sum(np.resize(mu_1*self.Ye[:,x], (4,self.e)), 0)\n",
    "                    tmp = np.repeat(np.expand_dims(smu_1_te, 0), 4, 0)\n",
    "                    smu_1_te = np.resize(tmp, (self.mu[:,x].shape[0],1))\n",
    "                    kxx_mu_1 = np.asarray(invert(self.Ye[:,x]))*self.params.C + mu_1 - smu_1_te[:,0]\n",
    "                    \n",
    "                else:\n",
    "                    kxx_mu_1 = np.zeros(self.mu[:,x].shape)\n",
    "\n",
    "                Kmu_1 = self.Kmu_x[:,0] + kxx_mu_1 - self.Kxx_mu_x[:,x];\n",
    "\n",
    "            else:\n",
    "                if G0 < self.params.tolerance:\n",
    "                    break\n",
    "                else:\n",
    "                    kxx_mu_1 = np.zeros(self.mu[:,x].shape)\n",
    "                    mu_1 = np.zeros(self.mu[:,x].shape)\n",
    "                    Kmu_1 = self.Kmu_x[:,0] + kxx_mu_1 - self.Kxx_mu_x[:,x]\n",
    "\n",
    "            d_x = mu_1 - self.mu[:,x]\n",
    "\n",
    "            Kd_x = Kmu_1 - self.Kmu_x[:,0]\n",
    "            l = gradient.T.dot(d_x)\n",
    "            q = d_x.T.dot(Kd_x)\n",
    "            alpha = min(l/q,1)\n",
    "\n",
    "            delta_obj = gradient.T.dot(d_x)*alpha - alpha**2/2*d_x.T.dot(Kd_x)\n",
    "\n",
    "            if np.logical_or(delta_obj <= 0, alpha <= 0):\n",
    "                break\n",
    "\n",
    "            self.mu[:,x] = self.mu[:,x] + d_x*alpha\n",
    "            self.Kmu_x[:,0] = self.Kmu_x[:,0] + Kd_x*alpha\n",
    "            self.obj = self.obj + delta_obj\n",
    "            self.Kxx_mu_x[:,x] = (1-alpha)*self.Kxx_mu_x[:,x] + alpha*kxx_mu_1\n",
    "\n",
    "            it += 1\n",
    "\n",
    "        \"\"\"\n",
    "        For speeding up gradient computations: \n",
    "        store sums of marginal dual variables, distributed by the true edge values into Smu\n",
    "        store marginal dual variables, distributed by the pseudo edge values into Rmu\n",
    "        \"\"\"\n",
    "        \n",
    "        temp = self.mu[:,x].reshape((4,self.e))\n",
    "\n",
    "        for u in range(4):\n",
    "            self.Smu[u][:,x] = np.sum(temp).T*self.IndEdgeVal[u][:,x]\n",
    "            self.Rmu[u][:,x] = temp[u,:].T\n",
    "\n",
    "        self.mu[:,x] = np.resize(temp, (4*self.e,1))[:,0]\n",
    "        \n",
    "        return it\n",
    "    \n",
    "\n",
    "    def compute_Kmu_x(self, x, Kx):\n",
    "        \"\"\"\n",
    "        GRADIENT COMPUTATIONS\n",
    "        Gradient for x\n",
    "        \"\"\"\n",
    "\n",
    "        \"\"\"\n",
    "        For speeding up gradient computations:\n",
    "        Smu : store sums of marginal dual variables, distributed by the true edge values\n",
    "        Rmu : store marginal dual variables, distributed by the pseudo edge values\n",
    "        \"\"\"\n",
    "        \n",
    "        m = Kx.shape[0]\n",
    "\n",
    "        if self.Rmu == []:\n",
    "            for u in range(4):\n",
    "                self.Smu.append(np.zeros((self.e,m)))\n",
    "                self.Rmu.append(np.zeros((self.e,m)))\n",
    "\n",
    "        for u in range(4):\n",
    "            Ind_te_u = self.IndEdgeVal[u][:,x] # loss of edges for this label, associated to the moledule x (boolean)\n",
    "            H_u = self.Smu[u].dot(Kx)- self.Rmu[u].dot(Kx)\n",
    "            \n",
    "            self.term12[0,Ind_te_u] = H_u[Ind_te_u].T\n",
    "            self.term34[u,:] = - H_u.T\n",
    "\n",
    "        tmp = np.vstack(np.repeat(np.expand_dims(self.term12,axis=0), 4, axis=0))\n",
    "\n",
    "        Kmu_x = np.resize(tmp + self.term34, (4*self.e,1))\n",
    "        \n",
    "        return Kmu_x\n",
    "    \n",
    "    \n",
    "    def compute_Kmu(self, Kx, mu0):\n",
    "        \"\"\"\n",
    "        Complete gradient\n",
    "        \"\"\"\n",
    "\n",
    "        m, m_oup = Kx.shape\n",
    "\n",
    "        if 0:\n",
    "            for x in range(m):\n",
    "                self.Kmu[:,x] = self.compute_Kmu_x(x, Kx[:,x], self.E, IndEdgeVal)\n",
    "\n",
    "            Kmu = np.resize(Kmu, (4,self.e*m))\n",
    "            \n",
    "        else:\n",
    "\n",
    "            mu0 = np.resize(mu0, (4,self.e*m))\n",
    "            Smu = np.resize(np.sum(mu0, 0),(self.e,m))\n",
    "            term12 = np.zeros((1,self.e*m_oup))\n",
    "            Kmu = np.zeros((4,self.e*m_oup))\n",
    "\n",
    "            for u in range(4):\n",
    "                IndEVu = self.IndEdgeVal[u]\n",
    "                Rmu_u = np.resize(mu0[u,:], (self.e,m))\n",
    "                H_u = Smu*IndEVu\n",
    "                H_u = H_u - Rmu_u\n",
    "                Q_u = H_u.dot(Kx)\n",
    "\n",
    "                term12 = term12 + np.resize(Q_u*IndEVu, (1,self.e*m_oup))\n",
    "                Kmu[u,:] = np.resize(-Q_u, (1,self.e*m_oup))\n",
    "\n",
    "            for u in range(4):\n",
    "                Kmu[u,:] = Kmu[u,:] + term12\n",
    "\n",
    "        return Kmu\n",
    "    \n",
    "    \n",
    "    def compute_w_phi_e(self, Kx):\n",
    "\n",
    "        m = Kx.shape[0]\n",
    "        self.Ye = np.resize(self.Ye, (4,self.e*m))\n",
    "        self.mu = np.resize(self.mu, (4,self.e*m))\n",
    "        m_oup = Kx.shape[1]\n",
    "\n",
    "        if not np.any(self.mu>0): # no match\n",
    "            self.w_phi_e = np.zeros((4, self.e*m_oup))\n",
    "            \n",
    "        else:\n",
    "            self.w_phi_e = np.sum(self.mu, 0)\n",
    "            self.w_phi_e = np.repeat(np.expand_dims(self.w_phi_e, 0), 4, 0)\n",
    "            self.w_phi_e = self.Ye*self.w_phi_e\n",
    "            self.w_phi_e = self.w_phi_e - self.mu\n",
    "            self.w_phi_e = np.resize(self.w_phi_e, (4*self.e,m))\n",
    "            self.w_phi_e = self.w_phi_e.dot(Kx)\n",
    "            self.w_phi_e = np.resize(self.w_phi_e, (4,self.e*m_oup))\n",
    "            \n",
    "        self.Ye = np.resize(self.Ye, (4*self.e,m))\n",
    "        self.mu = np.resize(self.mu, (4*self.e,m))\n",
    "\n",
    "    \n",
    "    def predict(self, Kx):\n",
    "        \"\"\"\n",
    "        MODEL ERROR COMPUTATION\n",
    "        \"\"\"\n",
    "        \n",
    "        self.compute_w_phi_e(Kx)\n",
    "        Ypred, _ = self.max_gradient_labeling(self.w_phi_e, self.params)\n",
    "        \n",
    "        return Ypred\n",
    "    \n",
    "    \n",
    "    def assess_performance(self, Kx, y_true, verbose=False):\n",
    "        \"\"\"\n",
    "        Compute the accuracy score of a kernel sample Kx\n",
    "        y_true are the true labels\n",
    "        \"\"\"\n",
    "        \n",
    "        self.Ypred = self.predict(Kx)\n",
    "        accuracy_per_mol = np.mean(self.Ypred == y_true, 1)\n",
    "        f1 = f1_score(y_true.flatten(),self.Ypred.flatten())\n",
    "        if verbose:\n",
    "            print('Average accuracy: %.2f%%' %(np.mean(accuracy_per_mol)*100))\n",
    "            print('Standard deviation: %.2f%%' %(np.std(accuracy_per_mol)*100))\n",
    "            print('F1 score: %.2f%%' %(f1*100))\n",
    "        \n",
    "        return np.mean(accuracy_per_mol), np.std(accuracy_per_mol), f1\n",
    "    \n",
    "    def plot_results(self):\n",
    "        \"\"\"\n",
    "        Plot the evolution of the metrics with the number of iterations\n",
    "        \"\"\"\n",
    "        \n",
    "        tmp = np.array(self.metrics_tr)\n",
    "        \n",
    "        list_iter = tmp[:,0]\n",
    "        list_mean_acc = tmp[:,1]\n",
    "        list_mean_std = tmp[:,2]\n",
    "        list_f1_score = tmp[:,3]\n",
    "\n",
    "        f, ax = plt.subplots(figsize=(12,5), dpi=60)\n",
    "\n",
    "        ax.plot(list_iter, list_mean_acc, label='mean accuracy')\n",
    "        ax.plot(list_iter, list_mean_std, label='accuracy std')\n",
    "        ax.plot(list_iter, list_f1_score, label='f1 score')\n",
    "        ax.set_ylim([0,1])\n",
    "        ax.set_xlim([min(list_iter),max(list_iter)])\n",
    "        box = ax.get_position()\n",
    "        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])\n",
    "        ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "        ax.set_xlabel('iterations', fontsize=15)\n",
    "        ax.set_title('Evolution of metrics with number of iterations.', fontsize=17)\n",
    "        plt.show()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class params():\n",
    "    \n",
    "    \"\"\"\n",
    "    Parameters for the MMCRF model.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, C, epsilon, tolerance, max_iter, max_CGD_iter, max_LBP_iter, profiling, profile_tm_interval, verbosity, debugging):\n",
    "        self.C = C\n",
    "        self.epsilon = epsilon\n",
    "        self.tolerance = tolerance\n",
    "        self.max_iter = max_iter\n",
    "        self.max_CGD_iter = max_CGD_iter\n",
    "        self.max_LBP_iter = max_LBP_iter\n",
    "        self.profiling = profiling\n",
    "        self.profile_tm_interval = profile_tm_interval\n",
    "        self.verbosity = verbosity\n",
    "        self.debugging = debugging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "class profile():\n",
    "    \"\"\"\n",
    "    EXECUTION PROFILING\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.start_time = time.time()\n",
    "        self.next_profile_tm = self.start_time\n",
    "        self.n_err = 0\n",
    "        self.p_err = 0\n",
    "        self.n_err_microlbl = 0\n",
    "        self.p_err_microlbl = 0 \n",
    "        self.n_err_microlbl_prev = 0\n",
    "        self.microlabel_errors = []\n",
    "        self.it = 0\n",
    "        self.err_ts = 0\n",
    "        print('alg: M3LBP\\ttm:\\titer:\\tobj:\\tmu:max\\tmu:min\\tdgap:')\n",
    "        \n",
    "    def profile_update(self, Ye, obj, mu, primal_ub, params):\n",
    "        \n",
    "        m = Ye.shape[1]\n",
    "        tm = time.time()\n",
    "\n",
    "        print('\\t\\t%d\\t%d\\t%.2f\\t%.2f\\t%.2f\\t%.2f'\n",
    "        %(np.round(tm-self.start_time), self.it, obj, np.max(mu), np.min(mu), primal_ub-obj))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parameters = params(C = 100,\n",
    "                    epsilon = 1,\n",
    "                    tolerance = 1e-8,\n",
    "                    max_iter = 10,\n",
    "                    max_CGD_iter = 100,\n",
    "                    max_LBP_iter = 10,\n",
    "                    profiling = True,\n",
    "                    profile_tm_interval = True,\n",
    "                    verbosity = 1,\n",
    "                    debugging = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test\n",
    "\n",
    "We train the algorithm on a sample of 10 molecules against all 59 cancer cell lines and test it on 90 molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_label = target.copy()\n",
    "Kx_tr = Kx[2000:2010,2000:2010]; y_label_tr = y_label[2000:2010,:]\n",
    "Kx_ts = Kx[2000:2010,2010:2100]; y_label_ts = y_label[2010:2100,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All good!\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "We check that no cancer has zero molecules with effects or the contrary\n",
    "\"\"\"\n",
    "if np.any(np.all(y_label_tr==1, 0)):\n",
    "    print(np.where(np.all(y_label_tr==1, 1)))\n",
    "    if np.any(np.all(y_label_tr==-1, 0)):\n",
    "        print(np.where(np.all(y_label_tr==-1, 1)))\n",
    "else:\n",
    "    print('All good!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We test the model for:\n",
    "\n",
    "- The output graph obtained by applying a maximum spanning tree algorithm \n",
    "\n",
    "- The output graph obtained by correlation thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Maximum spanning tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t8.47%\t0.00%\t0.156\n",
      "\t\t0\t2\t4.8e+00\t4.0e-01\t0.0e+00\t5.1e+05\t\t66.44%\t3.20%\t0.100\n",
      "\t\t0\t5\t4.1e+02\t4.2e-01\t0.0e+00\t5.1e+05\t\t64.07%\t3.54%\t0.078\n",
      "\t\t0\t8\t8.6e+02\t3.8e-01\t0.0e+00\t5.1e+05\t\t64.41%\t3.71%\t0.087\n",
      "\t\t0\t10\t1.2e+03\t3.7e-01\t0.0e+00\t5.1e+05\t\t67.46%\t3.46%\t0.150\n",
      "\t\t0\t13\t1.6e+03\t4.3e-01\t0.0e+00\t5.1e+05\t\t69.83%\t3.20%\t0.110\n",
      "\t\t0\t24\t2.1e+03\t4.4e-01\t0.0e+00\t5.1e+05\t\t68.47%\t3.41%\t0.155\n",
      "\t\t0\t35\t2.6e+03\t4.2e-01\t0.0e+00\t5.1e+05\t\t66.10%\t3.39%\t0.091\n",
      "\t\t0\t37\t3.0e+03\t3.8e-01\t0.0e+00\t5.1e+05\t\t60.00%\t3.95%\t0.092\n",
      "\t\t0\t40\t3.3e+03\t4.1e-01\t0.0e+00\t5.1e+05\t\t66.10%\t3.03%\t0.091\n",
      "\t\t0\t44\t3.8e+03\t5.3e-01\t0.0e+00\t5.1e+05\t\t66.27%\t3.43%\t0.091\n",
      "\t\t0\t44\t3.8e+03\t5.3e-01\t0.0e+00\t5.1e+05\t\t66.27%\t3.43%\t0.091\n",
      "Descent done.\n"
     ]
    }
   ],
   "source": [
    "model = model_MMCRF(parameters)\n",
    "model.learn(Kx_tr, y_label_tr, G_max_spanning_tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average accuracy: 66.27%\n",
      "Standard deviation: 3.43%\n",
      "F1 score: 9.13%\n"
     ]
    }
   ],
   "source": [
    "# Results on training set\n",
    "m_tr, s_tr, f1_tr = model.assess_performance(Kx_tr, y_label_tr, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAEZCAYAAACEpG9LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAJOgAACToB8GSSSgAAIABJREFUeJzt3Xd85FW9//HXd0p6spvNZnezjaUedulFOgpiQVDRiwiC\nVxErdq969ad4xauXa0MvlmsBFaxgR0FRLFylr8Iu/VB32ZLspm2yyaTNzPf3x/nOZGYyySabTKbw\nfj4eeczMt54z35l8P3Oq5/s+IiIiIpUgVOwEiIiIiMwVBTYiIiJSMRTYiIiISMVQYCMiIiIVQ4GN\niIiIVIxIsRMge2aM2QTsM8nqN1lrr53l8W8D7rbWfnQa2zYAr7XWfjcjbZ+11n5zNmnYG8aY1wBf\nBRqBk6y1DxTwXFn5zrN+DfAMsNZa+1ih0jGZzOtQStdoT4wx1wI11toLipiG9wCfAsLAGmttb8a6\ny4EzrbUnBK/PxX1XthUoLfsCh1prf1vsz5RIuVKJTfn4CNCW5++GeU7HB4G3Zbx+HnDdPKch5dPA\nb4FDgUcKfK7cfOfagrseTxQ4HZPJvA57SqsEjDEe8DngKuCIzKAm8EXgrGDbfYCf4wLpQvkucHLw\nvNifKZGypBKb8tFvre0odiIAL/OFtbazWAkBFgJ3WWs3zcO5vKlWWmsTQNGuT851mDKtkqUKqAX+\nlu9zZK0dyHg5H+9r+hzF/kyJlCtPA/SVvj1VJRhjzgBuBlqttbuDZYcADwCrrLXbjTEXAf8P2B94\nEvi4tfY3wba3EVRF5asaMMZ0AKlqqu+llltrvZwqEA94H/BuYAWwEfigtfaOjHx8CTgPV8KwCXiP\ntfbWSfLVhvvF/FJcNcGvgfdba/uMMZkf3P+z1p6Ws+/FwFuAGzPS/hngn8A3gdXAH4CLrLXDwT6X\nAB/D/Up+MJX24Fj58v0z4EJgEHgF8BhBtYExpgX4CnA2MAb8FPiAtXbUGPMvQVr2w/0qv8Jamz5+\nRh42AD+01n4xeH09cJy1dr/g9QuC92Qx8BTwWWB4krT+EjgWOA54Gnintfa2POdcg6v+OA+4AlcF\nehdwsbV2kzHmNOCvQG3G+/ZZ4ARr7Wkzfd+Dz1sLMBS8h88CH7HW/jo4thcc513AAuAe4L3W2keC\n9bnX4TBr7UhOng7Cfe5OAUaBH+C+C8uDvKZcZ629OGffywmqonI+c2+y1l5rjHlF8D4dADwO/Ie1\n9sZg32txn9t1wft4JrAN+DLwYlzJz5PAR621vwm2f2Nw/P8DLiajKsoYswD4b+DVwb5/xH1/tu3p\numXk5a3B+30/8G/W2rsQqTCqiqoMtwG9uBtDynm4X6GpoOYa3D/UI4CfAL80xhw5w/PcAFwJ3Ie7\n+ef6GPAJXLXZUbib0C3GmOUZ21yOu+EfirsRfMcYM+FzaIyJAn8GlgAvAl4OHMP4TbsN6ATeD/zL\nJOl9HnAkcAIu718A/gd34z0HeAnBjcQYcxbweeDfce/Rz4FbjTH7TZHvN+ACl/OBrJsp8CtcEPni\nIO0vBT5ljFkCXI+r+jC4m9B3jDEmT/r/AJyW8foFwJqM9/MlwB+DX/Ypk6X1Lbj2SIfiArAfBkHD\nZC4D3gQcD6wC/muKbXNN+30PvBzYHOzzA+BnwfsOcCnuZvwmXGB2P3BbcJNPSV+HPEHNIuDvQB+u\niudfcd+Nz+KCyjXBpufigvKpHBc8ngrcYIw5DHctv4R7X68CfmyMOTljn4twwciZwIYgfw3A84HD\ngPXAd40xVcH57wK+Rv7P9M+BE4O0ngo0A7/O+f7kvW7GmFfhvisXAWtx7+PP9vAZEClLCmzKx/8Y\nYwZy/yBdZP1T4DUZ25+HC2AA/g34hrX2O9bax621V+BKeP59Jgmw1g4BA8BYbrVY8A/y/cCnrLW/\nCBo7vh9XkvDujE1/ZK39mbX2Sdwv+VXAsjynOxNXonGhtXZDUOrzRuDVxpiDgvMngT5rbc8kSY4C\n77LWPgF8Hffr+WvW2juttX8G/ob7NQ2uVOBz1tpfWmufDEpJ/oQr2Zgs3z8J0nZ/zntxCO7Gc7G1\ndr219h5cm5etuJKsKLDVWrs5KKl5EbAjT/pvAU4xxoSNMQcDPu6GdEqw/iXA7zN3mCKt12S8758P\n0tE6yfsG8Glr7R3W2g3At3DBynTN5H0H2Git/bB1PgP8g/E2Qh/FleDcGqz/MK565l8z9s97HQIX\n4T4nl1hrH7bW/gEXQLwHqGH8fe+x1vbtIV+p6r6u4H3+MHCttfZ71tqngsba3wE+kLHPA9ban1tr\n/2GtjePahL0zSIvFBX0twNLg/KPAYO5nOgiiXgS8MXgf7wdeBxweLE+Z7LrtGxx7s7X2meB9fQO6\nB0gFUhub8vFfjAcq+fwE+Isxph5XDH0g8Itg3VrczSzTHbji+7nSiqsSuSe1wFrrG2PuJPsmltkQ\nsj94jOY53lrgmcy2I9ba+40xQ8HxHp9GmnqstbuC50PB46aM9UNAdfB8HXCcMeaTGeurmLpdxaZJ\nlq8DYpk9Way1f8FdHw/4DXCTMeYp3I3u2ox0ZroDFxQchSut+DvQhQt2/hws/32e/fJ5KuN56ly1\nU2yfe53yXaPJzOR9h4zPTOA+YG3Qu2sVcJ0xJrOqrgZX2pWSeexca4ENOSU5d+D+9x2IK73aW+uA\nw4wxmaVPUbI/m7lp+1/gfGPMibg8HBMsD+/hXGuBIZvR889a2xF8hjK/D5Ndt5/gSr+eMsasx30G\nv5NT2idSERTYlI/O4Nd2Xtbau4O2MGfj/gneaq3tDlYP59klTP5/pvkaXU3nc5LvHPnOM5pnm3zB\nw4TjBUGBx55vAinxPMuSk2wbwf0Czw0UhvJsmzJZnvPlEXDBHnCOMeZo4JXB3zuNMa+w1v4xZ9tR\nY8xfcdVRx+BKOrpwVX1n4Eo68pX05JPvBjZV0Jabh9S20/l8zOR9h4lpC+Oq9lLHfT2uvVim/ozn\nk12HydaFcx73VgRXxXZ1zvKxfOcPqoz+hCuhuQFX1diLu657srffLw/SQdA63OfmbOAdwLuMMcda\na7dP4/wiZUPFkJXletyN8tXB85THcO0dMp0I2DzHGCWjO6sxphlXl5+St7W5tbYfaM9znhMmOc+e\nPAbsG7RJSTkK92t9b443nfOtDqqhngyCyHfjqsRgknxP4gmgLmi0CoAx5gJjzL3GmIONMV+y1t5n\nrb3cWns0rl3Fqyc51h+A03HVT3/H3QQPx1U7/m6SfQrZIyB148zs8rxfvg1n4Iic18cDjwalPjuA\ntoxr8hSundbx0zz2Y8CRxpjMEqITccHXU/l3mVTu+/oYsF/OZ+Y8Ji8JXYdrW/Mya+2nrWu83xKs\nmypwTJ2r1hiTfq+MMctwbYT2+H0wxpwNvMNa+wdr7XuBg3C9Ck/d074i5UYlNuWjKfhHlisWBBXg\nipvvZLwHUcrngZ8YYx7C3Rj/BferLbNuPmU98HpjzJm4YvTPkP0rcABYZozZN6irz/R54BPGmK3A\nQ7hfhQfhGtfO1J9wY9P8yBjzYVy1yTeA26y1D+3F8fbki8D3jTGP4nqknIcLbF4QrJ8q31mstY8Y\nY/4EXGOMeW+Q9k/j2kH1Am83xvTjxp3ZF9fw9IeTHO4WXOPUQeDBoHrvaVxgc8ok+0w7rXvhYVwp\n1seNMVfhPkMvxvV62lsnGGP+E/gxrh3VgbgqG3DX5XJjTDuu1Oa9wKuAT+Y7UB4/Crb9rjHmv3Bt\ni76M6222yxhTM4N0prp+Hx58xq8E7jLGfAj3fXs+8J+4EqZ8duFKrs43xvwM13j4K8G6VOA1AOyf\nE9BjrX3cGHMjcK0x5p24a3AlrifUreRvzJ/JAz4flOr+A9c+qwrXoBljTCuuqmtg8kOIlAeV2JSP\nz+FKRHL/vpTawFr7IK53ye9s0O07WP5rXEPej+G6MV8AvNpa+9c85/kBLkD6Ga631d1k37R+gStq\nfzj3ny/un/RVweMGXMPFM+xejJpqrU3ietAM4tpE3IQLul4102NN83w/xw1s9zHczfsi4DVBo2WY\nOt/5vB5X2nA77qb3G+DyoOroVbgebA/j3u+v4xqd5kvXU7jeO7cH1VjgAq9eJrZNSZlpWqctCKIv\nwV2bR3AB8kx6TOXzfVxp3AZcD6mXZ1SxfQnXm+sruM/u8bgSj2mVtlhrB3E90tpwbXe+j6sGesdM\nExlU7V6Lu2ZvsdauxwXtF+Ou5cdwjabzDppprd0anPe9wKO4HnGX4aoXjwo2+xauuuiWPId4Ey64\n+z3uM9CD+37l9sjLd+6bcF3cP48r4Xk/cEHQgBncd+tDezqOSDnQODYiIiJSMVRiIyIiIhVjj21s\njDErcG0DLrLjI1i24aoc2oHd1trLCplIERERkemYssTGGNOE6wKbO8bG24GrrbXvA1YHw3mLiIiI\nFNWUJTZBQ8H3B3OYZGpjfOCp1Giqm3K2IZFI+pXShicUCpFMTjUUR/lQXkpPpeQDlJdSVSl5KYd8\nRCJhTVVRRHvb3XszbkTQJ4CVuIndJvB9n97e2F6eorQ0N9cpLyWoUvJSKfkA5aVUVUpeyiEfra2N\ne95ICmZGgU0wnsjNwLeBrxpjzgWeTLW9ERERESmmaQU21tqLg6dfyFj8ujlPjYiIiMgsqLu3iIiI\nVAwFNiIiIlIxFNiIiIhIxVBgIyIiIhVDgY2IiIhUDAU2IiIiUjEU2IiIiEjFUGAjIiIiFUOBjYiI\niFQMBTYiIiJSMRTYiIiISMVQYCMiIiIVQ4GNiIiIVAwFNiIiIlIxFNiIiIhIxVBgIyIiIhVDgY2I\niIhUDAU2IiIiUjEU2IiIiEjFUGAjIiIiFUOBjYiIiFQMBTYiIiJSMRTYiIiISMVQYCMiIiIVQ4GN\niIiIVAwFNiIiIlIxFNiIiIhIxVBgIyIiIhVDgY2IiIhUDAU2IiIiUjEU2IiIiEjFUGAjIiIiFUOB\njYiIiFQMBTYiIiJSMRTYiIiISMVQYCMiIiIVQ4GNiIiIVAwFNiIiIlIxFNiIiIhIxVBgIyIiIhUj\nMtkKY0wbcBXQDuy21l4WLF8BfBLoA1YDb7TWDs9DWkVERESmNFWJzduBq6217wNWG2PWBMsNcBaw\nGKhSUCMiIiKlYtISG6AN2BQ83wqsCF53AGdbazcaY64yxpxqrf17vgOEQh7NzXVzmNziCYdDyksJ\nqpS8VEo+QHkpVZWSl0rJhxTOVIHNZmAV8ASwEtgWLL8UuCF43g40TnaAZNKntzc2B8ksvubmOuWl\nBFVKXiolH6C8lKpKyUs55KO1ddLbosyDqQKbbwNfNcacCzwJnGeMuRm4BrjCGPMEUAN8ofDJFBER\nEdmzSQMba20X8LpJVp9dmOSIiIiI7D119xYREZGKocBGREREKoYCGxEREakYCmxERESkYiiwERER\nkYqhwEZEREQqhgIbERERqRgKbERERKRiKLARERGRiqHARkRERCqGAhsRERGpGApsREREpGIosBER\nEZGKocBGREREKoYCGxEREakYCmxERESkYiiwERERkYoRKXYC5LnB932GRxOEQx6RSIiQ5xU7SSIi\nUoEU2MicGxqJs61rkK2dA2zb6R63dg4wOBxPbxMOeUTCISJhF+hEw6HgdYhoxMt4HqyLuG3T20Xc\nY1NDNWNjcbdtel32duPH9ohGQnmO7REOqfBSRKQSKLCRvZZIJtnRM5QOXLYGQUxX33B6m6a6KCta\nGzjp0DaWNNeS9H3iiSTxeJKxxPjzeCLJWCJJPOEH69yy4ZE4u1PbJZKMpdbFg20TSRJJf9Z58Txy\ngiEvK4CKTjMISwdQma8jkwRdGedKeB6DsdGsIMxTqZaIyIwpsJE98n2fXQOjbOscYGvnIFt2DrCt\nc4Dt3THiiSQAVZEQyxfXc/DqZla21rNySQMrWxtoqq8qaNqam+vo7hnMCI78nEApSTzuZzzPWJ7w\nXaAUT44HTsH22fsH2yaSjI0lGRqO5w3CXODlp9+T2YqEc4MmbzzIiuQ8BkHXpCVXucfILLnK3Tfr\n2OMBnqoPRaQcKLCRLMOjcbZ1pqqPBtm6M7sayQOWNNeysrWBIw9czMrWBlYtaaB1YS2hUHFufCHP\noyoapioaLsr5c/m+ny5NGssqkRoPunKDqKqaCH19w+Pb5QZhcZ+xRCIdOGUee2Q0kROETQzw/NkX\nahEJh6ipClMdDVNTHaYmGqa6KkxNVcQtq3J/C5tq8BNJaqrd8uqqMLVV7tFtF0kfp1ifGRGpXAps\nnqOyq5EG2dY5wJad2dVIjXVRVgbVSKlSmOUt9VRXlUYAUao8zyMace15aqe5T3NzHb29sYKlKZHM\nX3I1llGllx2ETSy5Gh1zQdTwWILh0bh7Ppqgd/cII2MJRkbjDAfLpls9WBUNZQdIVS5gqgkCoZpo\nhJrqcDpAqskTSKX2rYmGqYqGVIUn8hynwOY5Ipn0eXzLLtbbnWzu2M2zOwbSVSbRSIgVGdVIK4Jq\npAUFrkaS+RMOhQhXQTWFD0qbm+vo7BoIgpzxAGh4LBE8j48HSCMJRsYSWduOjCXo7BtLPx8OAqbp\nlDp54EqG0gFSJCMgygyQxkuN8gVIqX3qG6rxfV/BkkgZUWBTwZK+z1Pb+rj30Z3847Gd9A2OUlcd\n4bADFnPovotY2drAyiUNLCliNZJUpkg4RENtiIba6Jwcz/dde6i8AVLwlxkEjZcsBc9H4uwaGMla\nNzKamNa5Q56XEfiEMwKi8aBpPECaWOqUL2CKhNULT6RQyiqweXzLLsYSSVbNQ6PUcuX7Pps6dnPv\noztY/9hOevpHqKkKc9SBrRy3dgmH7LuI1sUNBa32EJlrXkY7qqY5OmbS9xlNB0n5gyMvEqK3b3hi\nSdNogoGhMbr6ht3y4Dij8ek1HI+Evez2RnsZIGVWy+nHiYhTNoHNWDzBF35yf7ruPtWNeGVrw3j7\nj8X1VJdIA9L55Ps+W3YOcO+jO1n/2A46dw1TFQ1x5AGLOe5FSzlsv0VEI8+990VkKq4kJkJNVYQF\nk2wz07ZPiWSSkdHkxNKj0QTDY+MlTNklSuNtk/oGR9nZm7HvWIJ4YprtlSKhjAApklXF1lAb5dwX\nHkhDGbePSyST/NN28sS2fmJDY/j4JJM+vu/+ByZTj0kfHxe4+smM5elHt09qX/d6/HnSJ9gve33W\n8X2C8+cc063gN1eeU+y36zmtbAKbHb1DJJI+rzltf+prIukGr3c93MHA0BiQ3WNnRWt9SfTYKaRt\nXYOsf3QH9z66k46eGNFIiMP3a+HcF+zPEfsvViNfkXkWDoWoqwlRVxMBqufkmPFEMicYCtoi5Wm7\nNF7qNL5soG+Mxzb3ctdDHfzL8/fjRc9bVVZd94dH4/z9gXZuXb+Frr5hFi+ooToaxvNccOp5HqGQ\nK9XLWuYRrEu99oJ1pB/H15OzLtg/5BHCwwuRtX/6nOTZvwLvNeWmbAKbjm73q+moAxfT1lKfXp47\nxkpqsLiNT3VPGGMls3RnPsZYKYQdPTHufXQH9z62k22dg4RDHoft18IrTl7DkQcspra6bC6piEzD\nXLRXGhga4+d/e5rr//Ik9z3RxZvPXkvrwun22SuOXQMj/PmfW/nrfduIjcQ5dL9FXPyygznpyBXs\n2jVU7ORJCSubu2B7T4xwyJvwZfQ8j+bGapobqzl0v5b08tzuzFt3DvDYs73c/mB7eptyqc7q2jXE\n+sd2cu+jO9m8Yzchz2PdmmZe8rxVHH1QK/U1c9NAU0QqU0NtlA9ccBSH7tPM9/9g+Y/v3Mv5ZxzA\nC45YXnI9vrZ1DfKHe5/l7oc78H04Yd1SXnrcalYuaQAoufRK6SmbwKajO0brwtpp9yYIh1wpzfLF\n9Ry3dnz50Eic7cE8RqkpAO58qD3vAHSZ1VkNjdUMj8bTY36kRq7NfJ7ImRYgHowd4h6DbZJu7JBE\n0vXySAQj2qb2TeQcc3A4ztbOATwPDl7dzBvONBxzUCuNdeVX2iQixXXswUs4cNVCrvv9Y3z/Fst9\nj3fyppetpblxbqrN9pbv+9hnd3HLvc/ywFPd1FaHefGxq3jRsauKnjYpP+UT2PQM0tZSN+vj1FZH\n2H/FAvZfMd5cMLM6a0sQ8GzrHGDjU13Tbrg3XSHPyxgqPxjKPhTMSRTyCAdD3IfDriFgQ12UFxy5\nnGNNKwsa9AUXkdlZUF/Fe849jDsf6uDHf3qcT1xzDxe95CBOWLd03ktDEskk/3isk1vufZbNHbtZ\n1FTNBS88gFOPWK5qddlrZfHJ8X2f9u4YB69uLsjxJ6vOiieS7OgdYlvnACMJn7HR+HhAkjFZYTiY\nYyccTH6YCk4yt0k9V8MyESk2z/M4+bA21u7TzHd/9yhX//YR7nu8k399qaFpHkqDh0bi3P5AO39c\nv4Xu/mFWL23gba9cx7Fmicb4kVkri8Cmb3CU4dEEy+agxGYmImE3Iu+KxfUFH/JeRGS+LWqq4YPn\nH8lt92/jhr8+ySeuuYc3nnkwRx/UWpDz9e52DYJvu981CD5svxYuOetgDt6nWW1nZM6URWDTHvSI\naltUv4ctRURkJjzP4/SjV7Ju30V85+ZH+dovH+SkQ5dx4YsOpG6OOiZs6xzgD/du4a6HOwA44ZCg\nQXBrw5wcXyRTWQQ2HT0usJnvEhsRkeeKpc11fPTCo/nj+i388m9P8ejmXi45ay2H7Ltor47n+z6P\nbe7llnu38ODT3dRWR3jpcas545iVahAsBVUWgU179yCNddE5m3dGREQmCoU8zjx+NYftt4hrbnqU\nK2/YwOlHreC80/enpmp6t4t4Isk/7E7+cM8WNu/YTUtTNReccSCnHt6mBsEyL8riU9bRE2PZIpXW\niIjMhxWtDXz8Dcdw812buenOTTz8TA+XnL2Wg1YtnHSfoZE4f9+4nVv/sYXu/hH2WdrI2195CMce\n3Eo4pAbBMn/KI7DpjrFuTWF6RImIyESRcIhzTtmXIw5o4ZqbHuVzP7qPlx63mlc/f9+sued6d4/w\np39u4bb7tzM0Eufw/Vu45Ox1HLx6oRoES1FMGtgYY9qAq4B2YLe19rJgeT3wP8AuYBnwPmttT6ES\nODqWoLtvmGVqOCwiMu/WLGvikxcfy6/+/gx/uOdZHni6m7e8fC2RcCgYIXgHACcesoyXHreKFWoQ\nLEU2VYnN24GrrbW3GmO+b4xZY63dBLwV2A3UAY8UMqgBN/mljxoOi4gUSzQS5rWnH8CRByzmOzc/\nwqev+we+D3XVEc483jUIXqgBRKVETBXYtAGbgudbgRXBawM8Zq29yhjzRWPMC6y1/5fvAKGQR3Pz\n7AKSR57dBYDZd9GsjzUb4XCoqOefS8pL6amUfIDyUqrmIi/HN9dxuFnCb//+NHW1Uc44dtW8Nwiu\npGsihTHVJ3IzsAp4AlgJbAuWbwf6g+c7gUnLHZNJf9aD2j25pZdwyKPKo6gD5FXSAH3KS+mplHyA\n8lKq5jIvLz5mJQDDsVGGY6NzcszpKodr0traWOwkPKdNFdh8G/iqMeZc4EngPGPMzcA3gW8aY47C\nzRl5ZSET2NEdY0lzrVrVi4iIyB5NGthYa7uA102y+tzCJGei9p4YbS1qOCwiIiJ7VtLFIL7v09ET\nm5NZvUVERKTylXRgs2tglJHRhAbnExERkWkp6cCmvXsQUFdvERERmZ4SD2xSs3orsBEREZE9K+nA\npqMnRlN9FXU1mvxSRERE9qy0A5vuQbWvERERkWkr6UkwO3piHLpfS7GTISIiMide8cEb1wGzndW5\n97dXnvPIXKSnEpVsYDMymqC7f0Tta0REpCK84oM3rgQenqNjrfrtledszV1ujLkcWA10AvsAG4Fj\ngC/jJq9+DxAHksAHgMuApmCf64Fu4OPAbcCxwMXW2r7g2K3AF4AuYC3wBmAp8BFgFNgSnOfLwADQ\nClwOfNxae7Ex5uIgmWuAw4F7gD4gFeytt9Z+xRjzaaAe2B/4D+C/rbVnGWNOAs601v7HVO9NyQY2\nO4Ihs9UjSkREKsFvrzxn6ys+eOMhzE2JzYSgJsMvrbU3GWPWW2svMMa8BHghcAjQDozggob9gT8D\nUeBk4OXAdcBGa+1/GWOuAI4A/hYcNwF8F1gYLD8UNzH2h6y1HcaY44DzgT9Za683xhwaHDufH1hr\nfxXs8zAuuDnfGPM7oNVa+w5jzIpg2weD7d6FC8amVLKBTapH1DKNOiwiIhVinqqQBoLHweAxCYRx\n7Wq/a6190BhzIbAbuAb4DHAvcEDO/qNkt8U9HXgx8D3cPJIeUB0cH1wJUTXj80m24CbPrsp43R08\n7w0eP4sr4bkLuCDneFXBPv8DXA08a63duafMl3BgM0gkHGJxU02xkyIiIlIJPgN82hizBReU/BwX\n3LwMqMOVxEylCzc59muBg4HFwOeALxtjBoGncfNJft0Yc2pwzHcCDcaYr+Oqpn6Xc8wdwIuAGiBs\nrX3YGDNsjPkK0IYrDWo3xiSBq6aTSc/3/elst1fi8YS/t7OwfvPGh9jWNcin33z8HKdq75TDjLLT\npbyUnkrJBygvpapS8lIO+WhtbfSKnYZKEgQ5A9baj01n+5Itsenojqmrt4iIyHOctfa9M9m+JMex\nSfo+Hb2a/FJERERmpiQDm97+EUbHkiqxERERkRkpycCmoyeYI0o9okRERGQGSjKwSc/qrRIbERGR\nimKM2aeQxy/JxsMdPTEWNFRRW12SyRMREdkrr73h0jmZUuGn53+jLKdUCIKaTwEXZyw7DTjNWnv5\nXJyjJCOH9u6YplIQEZGK8tobLp2zKRVee8Olq356/jfyTalwEvBm3OB8C6y1bzTGnA+chhvs7jrg\nGbKnQXgWwFp7rTHmWtw0CN/AjVL8W9zAfKO4Afy+ADxA9rQJ1+ICk48bYz4G3G2t/UuQnlfhRjSO\nAH/ETd9wlDHmKOBC3Fg3y4AH5+J9gRINbDp6YhxxwOJiJ0NERGTO/PT8b2x97Q2XzsmUCvmCmsAO\nXKCxHDd4HsBbrLUvNsbUAAcBHyN7GoR1eY5TC3wQN4BfN24AvShuML21ZE+b8DTwAWPMQuAka+0V\nGcdZhWuSGe0rAAAgAElEQVT2ciNwH25gwBNwQVGDtfZSY8xZwHF78T7kVXKBzdBInN7dmvxSREQq\nzzxUIb0PeAw3v1Nq2oLUlAZxXKlL7jQIcVzJCbhSHQCstb3GmNW4STE/C9yPm8Ayd9qETtzUDDfg\nSoQy3Qn8HngB8L+4kYgJzp8ayDC+NxmdTMkFNpr8UkREZK9tBk7ElcJUGWMWA9caY67GzZj9fSZO\ng/AT4LvGmHVkBDaBGC4AeTmu2qkfN3dT5rQJ78CVyHwK+GnO/vviJsbcCtwK9OCmY1gN9Bhjvoqb\nmsHO1RtQcoFNRzD5pUpsREREZsZae2Wexd8L/jJdlPP6jJzXpwXH6wJekueY6f2NMUuAbwGXW2uz\nSl+stT/HzUmV6YTg8a95jjtrJRfYtHfHiEZCLFqgyS9FRERKXTDj9quLnY6UkhvHpqMnxtLmOkKe\n5hATERGRmSm5wKa9W3NEiYiIyN4pqcAm6fvs6NWs3iIiIrJ3Siqw6ekbZiyeVImNiIiI7JWSajzc\n3qOu3iIiUrnuOOfcOZlS4eQbf5F3PBxjjAGuxvVSugE3wnDEWvupWZ6zbJRUYJPq6q2qKBERqTR3\nnHPunE2pcMc55646+cZf5Bt9+H3B40PAe3GD8/mZGxhjFgFfAzqARdbai/NMuzCAG5+mF+gNpkt4\nBtdF+8vA24EEsAj4iLV2+1zkay6UVGDT3hOjubGamqqSSpaIiMisnXzjL7becc65czKlwiRBDbgB\n8nZaazcCG1MTTOZsUw+0AbcBG4wxHhOnXfgi8Epr7bAx5upg8L5+a+0lxph3BPs/AYwAp+JKh0pC\nSUUQHd2DKq0REZGKNVkV0jwbwZXstOECmIuZOO1CiOySHo/xKRpCwG3W2q8aY07Gle6UjJIKbNp7\nYhxzUGuxkyEiIlLJfODTwFPBXwcTp124ArjGGLMT2Gatfdg13wHgB7gpGPYHVgJvnef0T8nzfX/P\nW+2leDzh9wZzP+3J0Eicd335b1z4ogN50bGrCpamvdXcXMd081LqlJfSUyn5AOWlVFVKXsohH62t\njRphtohKprt3R9Ajqq2lvsgpERERkXJVMoFNe/cgoB5RIiIisvdKJrDp6IlRFQ3R3FRd7KSIiIhI\nmSqZwKa9O8YyTX4pIiIis1AygU1HT0wjDouIiMisTNrd2xjTBlwFtAO7rbWX5az/PBCz1l4+20Qk\nkz47eobU1VtERERmZaoSm7cDV1tr3wesNsasSa0wxrwNiM5VIrr6h4knkuoRJSIiIrMy1QB9bcCm\n4PlWYAWwyRjzMqAOuJGJwzRnCYU8mpv3XL309I7dAJh9F01r+2IIh0Mlm7aZUl5KT6XkA5SXUlUp\neamUfEjhTBXYbAZW4eaCWAlsC5ZfAnQDLwSWG2N+FcxJMUEy6U9rIKUnNrtRmmvDoZIdeKkcBoWa\nLuWl9FRKPkB5KVWVkpdyyEdra2Oxk/CcNlVg823gq8aYc4EngfOMMTdba88DSE2sNVlQMxMdPTFa\nmqqprgrP9lAiIiLyHDZpYGOt7QJeN8X623Azg85ae3dMA/OJiIjIrJVEd2/X1VsNh0VERGR2ih7Y\nxIbH6B8cVYmNiIiIzFrRA5v29OSXCmxERERkdooe2HR0a1ZvERERmRvFD2x6YlRXhVnYUFXspIiI\niEiZK3pgk+oR5WnySxEREZmlEghsBmlTw2ERERGZA0UNbBLJJDt7hzSrt4iIiMyJogY2XbuGSSR9\ndfUWERGROVHUwGa8q7d6RImIiMjsFTWw6eiO4QFLm2uLmQwRERGpEMUNbHoGaVlQQ1VUk1+KiIjI\n7BW3Kqo7pobDIiIiMmeKH9io4bCIiIjMkaIFNgNDYwwMjanhsIiIiMyZogU2qTmiVGIjIiIic6Vo\ngU17zyCgWb1FRERk7hS1xKamKsyCek1+KSIiInOjeIFNT4y2Fk1+KSIiInOneFVR6hElIiIic6wo\ngU08kaRz1xDL1CNKRERE5lBRApvOXUMkkj5tKrERERGROVSUwCbd1Vs9okRERGQOFSew6YnheZr8\nUkREROZWUQKb9u4YixfUEI1o8ksRERGZO0UrsdFUCiIiIjLXilRiM6iu3iIiIjLn5j2w2R0bZXA4\nrobDIiIiMufmPbBpD3pEqau3iIiIzLV5D2w6elJdvdXGRkRERObW/Ac23THqqiM01UXn+9QiIiJS\n4YpQFTXIMk1+KSIiIgVQlKoota8RERGRQpjXwMZNfjmsHlEiIiJSEPMa2OzsHSLp+yxbpIbDIiIi\nMvfmNbBp1+SXIiIiUkDzGth09AwS8jyWLNTklyIiIjL35jew6Y7RurCGaKQoMzmIiIhIhZvfqqie\nmOaIEhERkYKZt8DG9306ujWrt4iIiBROZLIVxpg24CqgHdhtrb0sWH4Q8Klg+RLgLdba4T2dqD82\nRmxEk1+KiIhI4UxVYvN24Gpr7fuA1caYNcHyJcDHrLX/BuwCzHRO1NE9CKCqKBERESmYSUtsgDZg\nU/B8K7AC2GStvR3AGPNqwLPWbpzsAKGQR3OzC2T6H+8E4OD9WljQUD3rhM+3cDiUzku5U15KT6Xk\nA5SXUlUpeamUfEjhTBXYbAZWAU8AK4FtAMaYCHAF8LS19l1THTyZ9OntdWPXPLVlF/U1ERKjcXp7\nE3OR9nnV3FyXzku5U15KT6XkA5SXUlUpeSmHfLS2NhY7Cc9pU1VFfRt4qzHm68CTwHnGmHXAJ4Az\ngVOMMT80xhwynRN19MQ0+aWIiIgU1KQlNtbaLuB1eVZ9MvibkY7uGAetWjjT3URERESmbV66e4/F\nE3T2DdGmHlEiIiJSQPMS2OzoHcL31SNKRERECmteApsOTX4pIiIi82BeApv2nhjhkEerJr8UERGR\nAipoYBN79ln8ZDKY/LKWSFiTX4qIiEjhTDWOzaw9+IEPUXvgQcQaj2HZkrZCnkpE9oLv+ySHh0n0\n95Ho7yeefuwn0deXfp3o78eLRIgsbiXaupjo4laii4PH1lbCdZoDTkRKQ0EDmwM+9G88c/V3ecmT\nN9B1yEkkxw4mFK0q5CllCiPbtrLj+9eyZXSEmkMOo+HIo6jZb3+8kErSKonv+ySHhjIClT4XqGQG\nLallfX34Y2MTjuFFo4QXLCDS1ES4aQHVK1eSHB1lrKuLgfvvI7FrV9b2obq6dLDTv7KNZMNCFwQF\ny0JV+t6LyPwoaGDTcuIJ7GpZzR//+3856qE72Hz5Eyx9/RupW7uukKeVHL7vs+vPf6Lr5zcQbmqi\nYd817PrLn+i95XeEGxqpP/xw6o84kvpDDiVUo3ZQpcj3fZKxWEaQklm60heUrvQHpSt9+PH4hGN4\nVVVEmhYQTgcrqzOCl6asdaGamikH00yOjRLv7masq5Oxzi732NXJWFcXXY9b4gMDWduHFyzILuXJ\nKO2JNC/CC4fn/D0TkeemggY2ADtiPn9ccgLHnHcW3i0/Z+uVn6fpxJNpfe0FhBs17HShxft20fG9\n7xB76EEajzueJa9/A4tXtNLd0UvssUcZ3LiBgQc20H/nHXiRCLXmYOqPOJKGI44k2rK42MmvaL7v\nkxwcJN7fT/+2EXZv25kdrGSUriT6+/MHK9XV6YAk2ryImn3WZAUpkaYF6eAlVFMzZ2kPRauoWtZG\n1bKJVczNzXV0besi3t3FaGcn8YygZ2TLswzcfx/+6GjGwUJEF7UQyQx6WsdLe8JNCzRiuYhMW8ED\nm44e19W77chDqH/eofT+8Ra6f3sjAw9upPW8C2g66WT90yqQgQ33s+Pa7+LHx1j25rfReMKJ6fc6\nVF1NQxDALEkmGXn2WQY23s/gxg10/viHdP74h1StWEnDEUdSf8SR1Oy7n6qspsFPJknGYjmlKNml\nKZnLSEycNy1UU0M4FawsWkTNmn1dqUpG1VAqaAlVl+aEsuG6OsJ1q6letXrCOt/3SfT3p4OdVGlP\nvKuLIfso/Xf1ZL0vXlUV0ZbFRBcvDqq3xkt7ootbCddpGAkRGef5vl+wg8fjCf/rP9vA3Q/v4Cvv\nOzW9fHTnTnb+8DpijzxM7cFrWfr6N1K1bFnB0jEXymHitZTkyAidP/0Jff93GzUHHEjbm99GtLU1\nvX5PeRnr6WHwgQ0MbtxA7NFH8ONxwo1N1B9+hKuyWnfInP76n435uC5+MklicCDdiNYFLbntV4LX\nu3fnD1ZqaydU90QyHheuXMIQ1YSbmko2WJmu2V4TP5EgvquXsc5U4LMzo7qri0Rfbvue+oxSnlSp\nT2sQCC2eVbu+cvre70ml5KUc8tHa2qhf60VU8MDmE9+6k5HRBP/v9cdkrfN9n9333k3n9T8mOTTE\norNfwaKXnY0XKXgh0l4phy8TwPCmTbRf/U3GOnfS8spXufc0p/3CTPKSHBkh9sjDDGzcwOADG9K9\nY2oPXpsuzYkuailEVqZlb6+Ln0ySGBjILkHp68sOXlJVQbv7IZmccIxQXV1OtU92aUq4aQGRBU2E\nG5v22Hi2XD5f01HovCRHR4l3d2WV9ox1daUDoWRsMGv78IKFGaU8OYHPHtr36LqUnnLIhwKb4ip8\nVVR3jLX7NE9Y7nkeTcefSP0hh9H5i5/SfeOv2H3vPSx9w8XUHnhQoZNVcfxkkt5bfkfXjb8iuqiF\nVR/9OLX77T/r44aqq2k46mgajjoaP5lkeNMmBh9wVVY7f/QD+NEPqF61ypXkHH4UNWvWFK3Kyk8m\nSezePaG6J6tUpc+VsiR274Y8QX2orj7dmDa6dBm1B5msoGU8eGlUD78iCVVVUdW2nKq25XnXJ2Kx\ndLAT7+pMt/MZeXYzA/f/M7t9TzhMdNEiootbM9r4jJf8+BpUVKTsFLTEZnBo1H/dZb/nvNMP4Mzj\nJ9a1Z4o9btn5/WsZ7Win6dTn0/qa8wnXl87YGKX8K2Gsu5uO73ybocctTac8nyUXXDhlVdFc5WWs\nu5vBBzYwsHEDQ4896qqsFiyg/vAjaDjiKOrWrpt1tYqfSJDYvTtrPJV0aUpfH97QAMPdvW7dwCTB\nSn09kQULMgKTzCqhJiJNC9PBS7FKDEv58zVTpZwX176nb7y0p3O8nU+8q4uxnu6s0rlQdTWRlpas\nYCeznU85te8p5esyE+WQD5XYFFdB/4u3dw3iM73JL+sOMqz+5H/Se8vv6Ln5twxu2EDrBRfSeNzx\nalw8hf5772bnD64DL0Tbpe+i8Zjnzdu5oy0tLDz9DBaefgbJ4SEGH36YwY0bGNxwP/1//xteNErd\n2nVBac6RRJtdyZ0fj5MY2D2hCijfWCuJgYG8wUq4oZFwUxM1Lc1Ur1hJeG1TELzktGNpbCzZ6k2Z\nf57nEVmwkMiChdTuf8CE9X4iQby3Jx3shAd2sXvLdsa6uhjevIlEX1/W9qH6+pxu7OPVXZGWFpXq\niRRBQf/jb+t0Y1m0TXPyy1A0SssrzqHxecex4wfX0XH1N+m/83aWvv6NWY1fxRW37/zxD9h9913U\nrV3H0kvemg4ciiFUU0vjMcfSeMyxrsrqmaddV/KNG1zgxXVEF7eSGB4imTPGCQCeR7ihIShVWUD1\nylVEmg4Zb7OyoCm9LtzYmG4XUQ6/3qR8eOFwOkCBtRM+X6lBCl0JT2dWyU/s0UdIxrI/i+GFCycG\nPkFD50jzIvU0FCmAwgY2OwcIhzwWL5xZD5qqZW2s/NBH6L/zdjp/ej2bPvlxWl5xDs0vfql+fQND\nTzxO+zXfItHXx+LzznfvSwn9g/RCIWr3P4Da/Q9g8b+8hrGuTgY2bmBk8+YgeMkeXyXc1ES4oVGD\ntEnJC1VVUb18OdXLJ2vfM5jRkHm8nc/Ipk0M/PMf2aM8p9v3LBkv5ckIgMKNjSqtFtkLBS+xWdJc\nS3gvbrqe57Hg5FOpP/wIOn96PV2/+Bn999zNgpNPoWr5Cqra2twvnufQF9+Px+m+6UZ6br6JqrY2\nVnz8P/KOE1JqootbaT7jxcVOhkjBhevqCa+up2b1PhPW+b5Poq8vuydX8Dj4yMPEe3qy2vd41dXp\n8XsmVHMtbiVcq4bNIvkUOLAZpK1ldg2AI41NtL35bSw46RR2Xv9jOm/4SXpdqKaG6LI2qpcvT/eS\nqGpbTrS1taRKMObC6I4ddFzzLYafeZqFLzyDxa85X/PviJQRz/OILFxIZOFCag84cMJ6Px4n3tub\nN/CZtH1P65KJ01Sk2/dE5ytrIiWl4CU2ZtXKOTlW3dp1rPnUZ0gODzPa0c7o9u2MtG9ntH07Q089\nSf9dd6YbmXqRiAt42tpcsBMEPtElS8vuy+77Pv23/42d1/+YUFU1y9/7ARoOP6LYyRKROeZFIm6Q\nwUnaEyZHRhjrHg924qkeXTt2EHv4IZJDQxkHc0FUunprUQuDddUMDQ7jJ5OQSOAnk/jJBCSyH/1E\nEpKJ7O0SCch9BAiF8EIeeCH3Y9LzgmXuuRcKQc5zL9hmT+s9LwSZxw6eD9ZVMzzkuuxP6NXr+0Cw\nzM9cRvb2mfuln/sT9vGztsl3XD/ncD7g0fr+d+a9hjI/ChrYDI8mpt1weLpCNTXUrNmXmjX7Zi1P\njo0ytmMHo+3tjGzfFjxud/PSpObYCYWIti6hqq2N6qA6q6rNPRZjtNfU8PuJWIxkbDD7cTBGcijG\n8OZNxB5+iPrDj2DpxW8m0tQ07+kUkeILVVdTvXwF1ctX5F2fGBycWNrT2cXIM88wcN8/2QUuiAiF\n8EJhCLtHLxyCzMdQyLV3y3n0olEIh8eDEQDfdwFQMjnxeSoo8n0XDAXLJzxPPyYh6YMfrE/6wbLs\nY/emoohUM4Tc5gieB3jZqzK3GV+Y2iy9zMtclnqStU32Mi/PsgnpkXlX8Ja4K1rnZyyaULSK6pWr\nqF65isypNf1EgrGuTka3u9IdV8rT7qYKGBlJbxdpaUlXZVWnq7XaCDc0THne5OhoVjCSGBoMHmMk\nB1PBSoxEbJBkKnAJgpasX1j58lRbS7ihgSUXvYEFp53+nGpPJCIzE66vJ1xfT80+a/Kur5QehJWS\nDymcggY2V1x6EssWFHdOIS8cpmrpMqqWLoOjjk4v95NJ4r29jLZvd0FPx3ZGtm+n/87bSQ6OD8ke\nbmqiqm053YsXMdQ3QDIIUFKBSr4Zl8d3DhOuqyNUV0+4vo5QbR1VbcsJ1dYRrq8nVFtHqL6OcF29\nG54//VhHqLZWvYRERERmqKCBzdo1i0o2svZCIaItLURbWqg/9LD0ct/3SezezWhQnTXa7h6Htm7D\nr65xgc6yZS5YyQlawvUuMEmt86qqVMoiIiIyjzQoTA7P84gEw+vXHbw2vVzFnyIiIqWvsvpEi4iI\nyHOaAhsRERGpGApsREREpGIosBEREZGKocBGREREKoYCGxEREakY6u49B5J+Mvjzx5/jhgpP+Ing\n0S33/SQJP4lPxrbpfRPukfHlvu8TDoWJeBEioQjRkHuMhMJEQ9H080goSsQLa9wcEZmRRDLBUHyY\nWHyIofjQhOdDY0PE4sPE4jH3Ou5eR7wwdZFa6qJ11EdrqYvUURetTS+ri9RSHzzWReuoCVfr/5PM\ni7IKbHaPDjAwNkg8mSCejJPw3WM8GSfuJ0gk426dHycRbBMPtkkkE4xNsjzux/McM5H1PEmSRDKR\nEbiMByalJOIFQU4onBMIBc+9CDXVVSTjEA6FCXshwl6YUPCYuSzruRcmFBp/Hs54HvJCE7bNd5yQ\nF3JTxvl+auo4kunJ5lKTyY2v84NJ6bKWZK73oSFRze7dw8Hkdtnr/JxjjK+b4pzB88x16Vd51vk5\n63LTPVV60mnwfWo6o8RyJvbLOn5wjIyj5JxzfPusY0x6TkimcuNnHCNIz8RzTl91dYSRkTgzu4VN\nf+uZ3RtncNw8y1J5mWUi3BxEs0jHbLf28fHDcXYNDmQFL7H4EKOJ0Sn3rQnXUBupoS5aS22khvpo\nPYtrW0gkE8TiQ/QO72LrwHZiYy7wmezTEvJCLshJBT5BEJQZ/KRe1+Ysj4Zmfqvygx+aCT/1v9w9\nuv/57jHf8mRqWTJOwk+m7ycJP5m+R0w4hp8gmd4/yb+f9vYZp1fmTskHNrGxITZ0Psj6jvt5YtfT\nM/wX60RCkeCGHyEcPKZu/OnXnisBqQnVZG8fChP1ItTVVjM6kiDsuVlnQ55HyAsHj+6mHSKU9drz\nQu6Gjpd+PmFfQuP7Zx7LC+HhXqcCrLFUEJfxN+YniCfH0oHZWMbzuB9nLBEPArfM/ROMxEdJpL70\nwZc29TzpJ/Kuk/nnpSfzG5+kz3MLyJi2D89LrfGCe276VTDXX/Ykf6H0tt6E82Qt8bKPMx3hUIhE\nMslMv6p7892e8R65s0HvgRfySCZnepaZbT/zXOeZ1XoaGqprqfKqqYvW0lK7iNpILXWRGmojtdnP\nozXUBctqIzWEvOm3WEj6SYaDoGlwLEZsbIhYPMbgmAuiYmOx9ONgPMau4b7069Hk2KTHjYai44FO\ndQ2jY/E8QUmcZDI7ENmbz9SehLwQkfQPttR9JJQuWQ+HNBVOsXl78wWZrng84e/NaL1jyTiPdD/G\nvR3381D3o8STcVY2LOfYpUeyrH5J+sOTG4C46prxgCUaigQBxuyLPytp5OG9zUvST6Z/oUwWFLkq\nt8QkAVMy64YcynNTneqGnHXTDfZtaqxhYGCErC29VACQeYzxACD3nJnHzAwfMgOEzOOO7+dNXJ5x\nzon75Z6D9PvR3FzPrt5YdgBTpsX2+q6UplLPy1gynhEIxRhKBUcZAdHgWAw/lCQZ9zMCiRDh4H9/\nqoQ4kllynApAsoKRnG2D4CR170hvk9oueJ360bknra2N5fnlrRAlU2KT9JM8tWsT63fcx307H2Qo\nPsSimmZeuOpUnrf0KJY3LCt2Ep/zQl6IUDhElGixk5LW3FxHb7h0/1lPV+qfqMhzVTQUYUF1Iwuq\nG6fcrtQDNCm+ogc22wc6WL/jftZ33E/vyC7qIrUcs+RwnrfsaPZbsM+MikFFRETkua0ogU3v8C7+\nsWMD63fcz7aBdiKhCIctXsd5S8/hkBZDZC8aiomIiIjMWwSR2wgY4KDm/Xn9wedx5JJDqY3UzldS\nREREpEJNGtgYY9qAq4B2YLe19rKpluczlhhjY+dDExoBv+qAszh26ZEsrF4wt7kRERGR57SpSmze\nDlxtrb3VGPN9Y8waa+2mKZZP8J7ff4L+kQEW1TRzxqrn87xlR9FWv3TucyEiIiLC1IFNG7ApeL4V\nWBG8nmz5BNe86gsV1eWttXXq1vrlRHkpPZWSD1BeSlWl5KVS8iGFMVWXo83AquD5SmDbHpaLiIiI\nFNWkA/QZYxYDXwV6gB3AEHAzsDNzubX2P+cnqSIiIiJTK+jIwyIiIiLzSaPfiYiISMVQYCMiIiIV\no2AD9M1kvJtSZoxZBfwO+CeQsNa+uchJ2ivGmBXAT4GLgBOAM3HX/5vW2tuLmbaZysnLvwJHA33A\n/1lrv1fMtE2XMeYg4FO478cS4CbK9JrkyYulDK8JgDFmf+AKXF52AM9QhtclTz4ilOk1STHGfB6I\nAY9RhtdE5k8hRx6e9ng3Je50oANIAncVOS17xRjTBHwY2BUsehfwfKAW+AXwsiIlbcby5OX5wONA\nI3B3sdK1F5YAH7PWPmOM+RrwDWARZXhNmJiXFwKPUH7XBGAhcBnwJPAX3P/Icvyu5OYjTnl+TwAw\nxrwN0rPvlu3/L5kfhayKyjfeTTm6F3iTtfYS4Mzgl1BZsdb2W2vfD3QGi+LWWt9aGwPqipi0GcuT\nl/8G3oMLdv6naAmbIWvt7UEg8GrAAzaU8TXJzcunKcNrAmCt/SfQDfwaeIIy/a7kyUdZfk8AjDEv\nw733NwaLyvKayPwpZGBTKePdHA3UBM97KIEZ0edA3BjjGWPqcN34y5UHnGStTeKK2MuGMSYSFK0v\ntda+izK+Jpl5Ad5NmV4TAGPMkUDEWnsOrpSjLK9LTj4SlPE1AS4BDgb+DXg5QDleE5k/hbxJfxv4\nqjHmXODJMq2GAtde4HPGmK3AVmutLXaC5sD/At/F/dr5RJHTMhs+kDDGfBuoBi4vbnJm5BO4dgIP\nGGNOAW6jfK9JOi/AKcCjZXpNwP1P/IYxZgsQxv0fK8frkpmPEOCV6zWx1p4HYIw5DTgN2Eh5XhOZ\nJxrHRkRERCqGunuLiIhIxVBgIyIiIhVDgY2IiIhUDAU2IkVmjFlToOPuU4jjioiUskrouixSdowx\nq3GD2O0PPAzUz/Hx3wMcA1ycOpe1tmEuzyEiUooU2IgUgbX2WaAhKK0pxCBjLbnnKsA5RERKjrp7\nixRBENA8A2zBDWQ5CBwE9ANfAl4JjABfs9Z+IdhnE3ArcC5uJNnrgK8DxwOtwD24+bMOA36DG8Dw\nFtyIs89Ya73gOG8EPh7scy/wLmvtk8E4IZ8FHgTOx00l8gFr7c3GmIbgfKcHafwZ8O/WWv0DEZGS\nojY2IsX1fABrbYO1djsuqGnFBTkvAC4xxpyfsf1C3HQl3wQ+jxs2f3/cqL8+8G5r7S24CRB/ZK19\nRebJjDFnAVfiAqClwB3ATcaYqmCT43EDoLUA3we+HCy/FDeC7VLgROC1wKlz8xaIiMwdBTYiJcIY\n4+FmK//3YE6sTbjA4uKMzX5prR2x1u4GPgp8CDdC7irclB/L9nCaC4FvW2vXW2tHgc/gJkY8Llg/\naK39mrV2DLgeOCBY3odrs/NaXEnSvtbav80mvyIihaDARqR0tOLmJVtvjNlljNkFfBFYnrFNR8bz\nVcCfgWeBr+CCGm8a59icehHMH7QFN58bQFfGtvGM410DXA38P2An8GtjTOu0cyYiMk8U2IiUjm5g\nDDjIWrvQWrsQ2Bc4K2ObzDYtPwKusdYusdaeAayfxjm2Aulu4MaYELAaF6xMZS1wvbX2UFzVVyOu\nxEhEpKQosBEprhEAY0yTtTYB3AB81hhTb4xpBn4BfGySfRfiGh1jjDkJeD0QzThuU559fgS81Rhz\nTA2KrFoAAACySURBVNCu5jJc25k79pDOC4FvGWMagU5cANYzvSyKiMwfBTYixdWB67m03RhzCPDu\nYPlTwJPAJuDDk+x7KfDfxpg+4CrcjMdrg3U3AycZY7LawVhr/wJ8BPgJroTodOCl1tqRPaTzClwg\nswnYDrTjGjqLiJQUdfcWERGRiqESGxEREakYCmxERESkYiiwERERkYqhwEZEREQqhgIbERERqRgK\nbERERKRiKLARERGRivH/AYauCIVY6GRoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11418a8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average accuracy: 64.05%\n",
      "Standard deviation: 4.87%\n",
      "F1 score: 21.08%\n"
     ]
    }
   ],
   "source": [
    "# Results on test set\n",
    "m_ts, s_ts, f1_ts = model.assess_performance(Kx_ts, y_label_ts, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Correlation thresholding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t8.47%\t0.00%\t0.156\n",
      "\t\t0\t3\t1.6e+00\t9.8e-02\t0.0e+00\t4.9e+05\t\t82.37%\t1.55%\t0.055\n",
      "\t\t1\t5\t3.6e+02\t8.5e-02\t0.0e+00\t5.0e+05\t\t85.08%\t2.25%\t0.120\n",
      "\t\t1\t16\t6.2e+02\t9.1e-02\t0.0e+00\t5.0e+05\t\t84.58%\t1.93%\t0.165\n",
      "\t\t1\t19\t9.3e+02\t8.5e-02\t0.0e+00\t5.0e+05\t\t85.76%\t1.55%\t0.067\n",
      "\t\t2\t21\t1.2e+03\t9.3e-02\t0.0e+00\t5.0e+05\t\t84.41%\t1.83%\t0.164\n",
      "\t\t2\t25\t1.5e+03\t8.8e-02\t0.0e+00\t5.0e+05\t\t84.75%\t2.63%\t0.182\n",
      "\t\t2\t28\t1.8e+03\t9.6e-02\t0.0e+00\t5.0e+05\t\t86.44%\t3.03%\t0.259\n",
      "\t\t2\t39\t2.1e+03\t8.5e-02\t0.0e+00\t5.0e+05\t\t87.12%\t1.55%\t0.156\n",
      "\t\t3\t42\t2.4e+03\t9.4e-02\t0.0e+00\t5.0e+05\t\t87.12%\t1.55%\t0.156\n",
      "\t\t3\t45\t2.7e+03\t8.7e-02\t0.0e+00\t5.0e+05\t\t85.42%\t1.55%\t0.140\n",
      "\t\t3\t45\t2.7e+03\t8.7e-02\t0.0e+00\t5.0e+05\t\t85.42%\t1.55%\t0.140\n",
      "Descent done.\n"
     ]
    }
   ],
   "source": [
    "model = model_MMCRF(parameters)\n",
    "model.learn(Kx_tr, y_label_tr, G_correl_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average accuracy: 85.42%\n",
      "Standard deviation: 1.55%\n",
      "F1 score: 14.00%\n"
     ]
    }
   ],
   "source": [
    "# Results on training set\n",
    "m_tr, s_tr, f1_tr = model.assess_performance(Kx_tr, y_label_tr, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average accuracy: 77.29%\n",
      "Standard deviation: 5.54%\n",
      "F1 score: 10.53%\n"
     ]
    }
   ],
   "source": [
    "# Results on test set\n",
    "m_ts, s_ts, f1_ts = model.assess_performance(Kx_ts, y_label_ts, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test of stratified k-fold to estimate the risk of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "parameters = params(C = 100,\n",
    "                    epsilon = 1,\n",
    "                    tolerance = 1e-8, # 10e-8\n",
    "                    max_iter = 5,\n",
    "                    max_CGD_iter = 5, #100\n",
    "                    max_LBP_iter = 5,\n",
    "                    profiling = True,\n",
    "                    profile_tm_interval = True,\n",
    "                    verbosity = 1,\n",
    "                    debugging = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stratified kfold with maximum-weight spanning tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fold 1\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t37.32%\t38.83%\t0.544\n",
      "\t\t5\t12\t1.3e+01\t4.1e-01\t0.0e+00\t1.7e+10\t\t52.26%\t9.03%\t0.382\n",
      "\t\t12\t31\t1.1e+03\t4.0e-01\t0.0e+00\t1.7e+10\t\t53.20%\t12.07%\t0.358\n",
      "\t\t18\t53\t1.7e+03\t4.5e-01\t0.0e+00\t1.7e+10\t\t55.09%\t17.71%\t0.309\n",
      "\t\t23\t86\t2.6e+03\t2.3e-01\t0.0e+00\t1.7e+10\t\t53.53%\t14.31%\t0.330\n",
      "\t\t28\t124\t2.8e+03\t3.2e-01\t0.0e+00\t1.7e+10\t\t54.45%\t12.85%\t0.370\n",
      "\t\t28\t124\t2.8e+03\t3.2e-01\t0.0e+00\t1.7e+10\t\t54.45%\t12.85%\t0.370\n",
      "Descent done.\n",
      "Fold 2\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t36.63%\t38.54%\t0.536\n",
      "\t\t4\t13\t1.1e+01\t4.1e-01\t0.0e+00\t1.7e+10\t\t54.57%\t13.14%\t0.347\n",
      "\t\t8\t33\t8.2e+02\t3.8e-01\t0.0e+00\t1.7e+10\t\t54.56%\t13.39%\t0.349\n",
      "\t\t12\t57\t1.5e+03\t2.1e-01\t0.0e+00\t1.7e+10\t\t53.50%\t15.57%\t0.309\n",
      "\t\t16\t86\t1.8e+03\t2.3e-01\t0.0e+00\t1.7e+10\t\t55.57%\t16.15%\t0.331\n",
      "\t\t21\t112\t2.1e+03\t2.8e-01\t0.0e+00\t1.7e+10\t\t53.92%\t10.86%\t0.373\n",
      "\t\t21\t112\t2.1e+03\t2.8e-01\t0.0e+00\t1.7e+10\t\t53.92%\t10.86%\t0.373\n",
      "Descent done.\n",
      "Fold 3\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t36.03%\t38.33%\t0.530\n",
      "\t\t4\t19\t1.2e+01\t4.6e-01\t0.0e+00\t1.7e+10\t\t50.76%\t4.34%\t0.418\n",
      "\t\t8\t60\t1.2e+03\t3.0e-01\t0.0e+00\t1.7e+10\t\t55.50%\t15.50%\t0.324\n",
      "\t\t12\t92\t1.5e+03\t2.9e-01\t0.0e+00\t1.7e+10\t\t54.52%\t15.01%\t0.319\n",
      "\t\t16\t110\t2.2e+03\t2.1e-01\t0.0e+00\t1.7e+10\t\t55.06%\t14.41%\t0.333\n",
      "\t\t20\t136\t2.4e+03\t2.9e-01\t0.0e+00\t1.7e+10\t\t54.24%\t14.07%\t0.324\n",
      "\t\t20\t136\t2.4e+03\t2.9e-01\t0.0e+00\t1.7e+10\t\t54.24%\t14.07%\t0.324\n",
      "Descent done.\n",
      "Fold 4\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t35.79%\t38.23%\t0.527\n",
      "\t\t4\t9\t1.2e+01\t4.2e-01\t0.0e+00\t1.7e+10\t\t56.13%\t14.24%\t0.351\n",
      "\t\t8\t39\t9.0e+02\t2.8e-01\t0.0e+00\t1.7e+10\t\t54.63%\t13.87%\t0.334\n",
      "\t\t12\t82\t1.3e+03\t1.2e-01\t0.0e+00\t1.7e+10\t\t54.04%\t14.07%\t0.325\n",
      "\t\t17\t114\t1.4e+03\t3.5e-01\t0.0e+00\t1.7e+10\t\t53.82%\t12.09%\t0.345\n",
      "\t\t21\t152\t1.9e+03\t3.0e-01\t0.0e+00\t1.7e+10\t\t54.34%\t12.79%\t0.348\n",
      "\t\t21\t152\t1.9e+03\t3.0e-01\t0.0e+00\t1.7e+10\t\t54.34%\t12.79%\t0.348\n",
      "Descent done.\n",
      "Fold 5\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t36.23%\t38.40%\t0.532\n",
      "\t\t4\t6\t6.9e+00\t4.0e-01\t0.0e+00\t1.7e+10\t\t57.46%\t15.57%\t0.365\n",
      "\t\t8\t19\t5.2e+02\t4.8e-01\t0.0e+00\t1.7e+10\t\t58.51%\t19.82%\t0.334\n",
      "\t\t13\t50\t1.1e+03\t1.1e-01\t0.0e+00\t1.7e+10\t\t58.94%\t19.73%\t0.340\n",
      "\t\t17\t71\t1.3e+03\t3.5e-01\t0.0e+00\t1.8e+10\t\t59.34%\t20.06%\t0.343\n",
      "\t\t21\t89\t1.7e+03\t1.6e-01\t0.0e+00\t1.8e+10\t\t57.58%\t17.21%\t0.356\n",
      "\t\t22\t89\t1.7e+03\t1.6e-01\t0.0e+00\t1.8e+10\t\t57.58%\t17.21%\t0.356\n",
      "Descent done.\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "skf = StratifiedKFold(5)\n",
    "listOfAccuracies_tr = []\n",
    "listOfAccuracies_ts = []\n",
    "listOfF1_tr = []\n",
    "listOfF1_ts = []\n",
    "\n",
    "#Little trick to stratify according to equivalence classes based on the number of cells against which the molecules\n",
    "#are active.\n",
    "\n",
    "listOfPercent = np.percentile(np.sum(y_label, axis = 1),list((np.arange(10))*10))\n",
    "bins = np.digitize(np.sum(y_label, axis = 1), listOfPercent)\n",
    "\n",
    "k = 1\n",
    "\n",
    "for train, test in skf.split(Kx, bins):\n",
    "    print('Fold %d\\n----------------------------------------------------------------------------------------------' %k)\n",
    "    \"\"\"\n",
    "    We run a MMCRF model for each fold\n",
    "    \"\"\"\n",
    "    Kx_tr, y_label_tr = Kx[np.ix_(train, train)], y_label[train,:]\n",
    "    Kx_ts, y_label_ts = Kx[np.ix_(train, test)], y_label[test,:]\n",
    "    \n",
    "    model = model_MMCRF(parameters)\n",
    "    model.learn(Kx_tr, y_label_tr, G_max_spanning_tree)\n",
    "    \n",
    "    \"\"\"\n",
    "    We compute the accuracy of the model on both the train and test datasets\n",
    "    \"\"\"\n",
    "    \n",
    "    m_tr, s_tr, f1_tr = model.assess_performance(Kx_tr, y_label_tr, False)\n",
    "    m_ts, s_ts, f1_ts = model.assess_performance(Kx_ts, y_label_ts, False)\n",
    "    \n",
    "    if k==1:\n",
    "        Y_pred = model.Ypred\n",
    "        Y_true = y_label[test,:]\n",
    "    else:\n",
    "        Y_pred = np.vstack((Y_pred, model.Ypred))\n",
    "        Y_true = np.vstack((Y_true, y_label[test,:]))\n",
    "    \n",
    "    listOfAccuracies_tr.append(m_tr)\n",
    "    listOfAccuracies_ts.append(m_ts)\n",
    "    listOfF1_tr.append(f1_tr)\n",
    "    listOfF1_ts.append(f1_ts)\n",
    "    \n",
    "    k +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training set:\n",
      "-------------------------\n",
      "Average accuracy: 54.91%\n",
      "Average F1 score: 35.43%\n",
      "\n",
      "Validation set:\n",
      "-------------------------\n",
      "Average accuracy: 54.76%\n",
      "Average F1 score: 35.24%\n"
     ]
    }
   ],
   "source": [
    "print('Training set:\\n-------------------------')\n",
    "print('Average accuracy: %.2f%%' %(np.mean(listOfAccuracies_tr)*100))\n",
    "print('Average F1 score: %.2f%%' %(np.mean(listOfF1_tr)*100))\n",
    "print('\\nValidation set:\\n-------------------------')\n",
    "print('Average accuracy: %.2f%%' %(np.mean(listOfAccuracies_ts)*100))\n",
    "print('Average F1 score: %.2f%%' %(np.mean(listOfF1_ts)*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With middle Active dataset (+ maximum-weight spanning tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fold 1\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t46.69%\t20.19%\t0.637\n",
      "\t\t1\t44\t7.9e+01\t9.4e-01\t0.0e+00\t8.5e+08\t\t48.55%\t6.50%\t0.441\n",
      "\t\t2\t65\t5.7e+03\t1.8e+00\t0.0e+00\t8.8e+08\t\t53.15%\t6.02%\t0.527\n",
      "\t\t2\t97\t3.1e+04\t4.5e+00\t0.0e+00\t9.3e+08\t\t50.25%\t6.21%\t0.431\n",
      "\t\t3\t119\t9.3e+04\t7.1e+00\t0.0e+00\t1.0e+09\t\t52.64%\t5.93%\t0.490\n",
      "\t\t4\t135\t2.5e+05\t1.1e+00\t0.0e+00\t9.8e+08\t\t48.02%\t6.77%\t0.433\n",
      "\t\t4\t135\t2.5e+05\t1.1e+00\t0.0e+00\t9.8e+08\t\t48.02%\t6.77%\t0.433\n",
      "Descent done.\n",
      "Fold 2\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t46.29%\t20.12%\t0.633\n",
      "\t\t1\t54\t1.0e+02\t9.7e-01\t0.0e+00\t8.2e+08\t\t50.24%\t7.14%\t0.437\n",
      "\t\t2\t88\t8.2e+03\t2.6e+00\t0.0e+00\t8.9e+08\t\t48.41%\t5.18%\t0.454\n",
      "\t\t2\t126\t3.8e+04\t5.8e-01\t0.0e+00\t8.7e+08\t\t50.85%\t7.39%\t0.422\n",
      "\t\t3\t156\t4.5e+04\t6.2e+00\t0.0e+00\t1.1e+09\t\t50.77%\t6.62%\t0.458\n",
      "\t\t4\t188\t4.9e+05\t1.0e+01\t0.0e+00\t1.2e+09\t\t52.31%\t7.26%\t0.434\n",
      "\t\t4\t188\t4.9e+05\t1.0e+01\t0.0e+00\t1.2e+09\t\t52.31%\t7.26%\t0.434\n",
      "Descent done.\n",
      "Fold 3\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t45.88%\t20.04%\t0.629\n",
      "\t\t1\t31\t3.4e+01\t4.8e-01\t0.0e+00\t8.2e+08\t\t53.30%\t9.93%\t0.402\n",
      "\t\t2\t93\t2.5e+03\t4.5e-01\t0.0e+00\t8.3e+08\t\t51.48%\t6.24%\t0.464\n",
      "\t\t3\t157\t4.4e+03\t4.1e-01\t0.0e+00\t8.4e+08\t\t51.70%\t6.39%\t0.446\n",
      "\t\t4\t187\t5.8e+03\t4.4e-01\t0.0e+00\t8.4e+08\t\t52.00%\t8.05%\t0.412\n",
      "\t\t6\t222\t6.6e+03\t3.5e-01\t0.0e+00\t8.4e+08\t\t51.34%\t7.94%\t0.425\n",
      "\t\t6\t222\t6.6e+03\t3.5e-01\t0.0e+00\t8.4e+08\t\t51.34%\t7.94%\t0.425\n",
      "Descent done.\n",
      "Fold 4\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t45.84%\t19.96%\t0.629\n",
      "\t\t1\t30\t1.6e+02\t1.5e+00\t0.0e+00\t9.0e+08\t\t52.78%\t5.88%\t0.468\n",
      "\t\t2\t71\t1.1e+04\t5.6e-01\t0.0e+00\t8.7e+08\t\t50.20%\t7.66%\t0.394\n",
      "\t\t4\t97\t1.3e+04\t1.4e+00\t0.0e+00\t8.8e+08\t\t51.81%\t8.03%\t0.427\n",
      "\t\t6\t115\t2.1e+04\t3.4e-01\t0.0e+00\t8.8e+08\t\t50.36%\t6.84%\t0.439\n",
      "\t\t8\t174\t2.2e+04\t1.9e+00\t0.0e+00\t8.7e+08\t\t51.33%\t6.54%\t0.441\n",
      "\t\t8\t174\t2.2e+04\t1.9e+00\t0.0e+00\t8.7e+08\t\t51.33%\t6.54%\t0.441\n",
      "Descent done.\n",
      "Fold 5\n",
      "----------------------------------------------------------------------------------------------\n",
      "alg: M3LBP\ttm:\titer:\tobj:\tmu:max\tmu:min\tdgap:\ttrain\tacc:\tstd:\tf1:\n",
      "\t\t0\t0\t0.0e+00\t0.0e+00\t0.0e+00\tinf\t\t45.94%\t20.02%\t0.630\n",
      "\t\t1\t36\t3.6e+01\t8.1e-01\t0.0e+00\t8.9e+08\t\t49.89%\t5.55%\t0.472\n",
      "\t\t3\t76\t3.5e+03\t3.9e-01\t0.0e+00\t8.8e+08\t\t53.23%\t9.31%\t0.429\n",
      "\t\t4\t139\t5.2e+03\t2.8e-01\t0.0e+00\t8.8e+08\t\t54.14%\t8.49%\t0.462\n",
      "\t\t6\t181\t6.0e+03\t6.2e-01\t0.0e+00\t8.8e+08\t\t50.93%\t7.14%\t0.454\n",
      "\t\t8\t200\t7.5e+03\t8.3e-01\t0.0e+00\t8.8e+08\t\t53.42%\t7.92%\t0.424\n",
      "\t\t8\t200\t7.5e+03\t8.3e-01\t0.0e+00\t8.8e+08\t\t53.42%\t7.92%\t0.424\n",
      "Descent done.\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "skf = StratifiedKFold(5)\n",
    "listOfAccuracies_tr = []\n",
    "listOfAccuracies_ts = []\n",
    "listOfF1_tr = []\n",
    "listOfF1_ts = []\n",
    "\n",
    "#Little trick to stratify according to equivalence classes based on the number of cells against which the molecules\n",
    "#are active.\n",
    "\n",
    "listOfPercent = np.percentile(np.sum(y_label_temp, axis = 1),list((np.arange(10))*10))\n",
    "bins = np.digitize(np.sum(y_label_temp, axis = 1), listOfPercent)\n",
    "\n",
    "k = 1\n",
    "\n",
    "for train, test in skf.split(Kx_temp, bins):\n",
    "    print('Fold %d\\n----------------------------------------------------------------------------------------------' %k)\n",
    "    \"\"\"\n",
    "    We run a MMCRF model for each fold\n",
    "    \"\"\"\n",
    "    Kx_tr, y_label_tr = Kx_temp[np.ix_(train, train)], y_label_temp[train,:]\n",
    "    Kx_ts, y_label_ts = Kx_temp[np.ix_(train, test)], y_label_temp[test,:]\n",
    "    \n",
    "    model = model_MMCRF(parameters)\n",
    "    model.learn(Kx_tr, y_label_tr, G_max_spanning_tree)\n",
    "    \n",
    "    \"\"\"\n",
    "    We compute the accuracy of the model on both the train and test datasets\n",
    "    \"\"\"\n",
    "    \n",
    "    m_tr, s_tr, f1_tr = model.assess_performance(Kx_tr, y_label_tr, False)\n",
    "    m_ts, s_ts, f1_ts = model.assess_performance(Kx_ts, y_label_ts, False)\n",
    "    \n",
    "    if k==1:\n",
    "        Y_pred = model.Ypred\n",
    "        Y_true = y_label_temp[test,:]\n",
    "    else:\n",
    "        Y_pred = np.vstack((Y_pred, model.Ypred))\n",
    "        Y_true = np.vstack((Y_true, y_label_temp[test,:]))\n",
    "    \n",
    "    listOfAccuracies_tr.append(m_tr)\n",
    "    listOfAccuracies_ts.append(m_ts)\n",
    "    listOfF1_tr.append(f1_tr)\n",
    "    listOfF1_ts.append(f1_ts)\n",
    "    \n",
    "    k +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training set:\n",
      "-------------------------\n",
      "Average accuracy: 51.28%\n",
      "Average F1 score: 43.13%\n",
      "\n",
      "Validation set:\n",
      "-------------------------\n",
      "Average accuracy: 51.51%\n",
      "Average F1 score: 43.32%\n"
     ]
    }
   ],
   "source": [
    "print('Training set:\\n-------------------------')\n",
    "print('Average accuracy: %.2f%%' %(np.mean(listOfAccuracies_tr)*100))\n",
    "print('Average F1 score: %.2f%%' %(np.mean(listOfF1_tr)*100))\n",
    "print('\\nValidation set:\\n-------------------------')\n",
    "print('Average accuracy: %.2f%%' %(np.mean(listOfAccuracies_ts)*100))\n",
    "print('Average F1 score: %.2f%%' %(np.mean(listOfF1_ts)*100))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}