evaluating_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating Approximated Ratios\n",
    "\n",
    "In order to evaluate the quality of the approximated ratios we can use the fact that \n",
    "$$\\frac{p(x|\\theta_1)}{p(x|\\theta_0)}\\cdot p(x|\\theta_0) = p(x|\\theta_1)$$\n",
    "The algorithms works as follow:\n",
    " 1. Compute the approximated density ratios for $p(x|\\theta_1)$ and $p(x|\\theta_0)$.\n",
    " 2. Compute the values of the ratio $p(s(x)|\\theta_1)/p(s(x)|\\theta_0)$ for data $\\in \\theta_0$.\n",
    " 3. Resample the data $\\in \\theta_0$ using the ratio values as weight.\n",
    " 4. Train a classifier to separate between data $\\in \\theta_1$ and the resampled dataset. Is expected that the classifier will works badly if the original data and the resampled data follows the same distribution. This only happens if the density ratios are being approximated correctly.\n",
    " 5. Evaluate the ROC curves of the classifier in the data. AUC = 0.5 means that the density ratios were approximated correcly. Better AUC show that the approximation is bad.\n",
    " \n",
    "In what follows we will show how to compute and evaluate the approximated density ratios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ROOT\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn import cross_validation\n",
    "from xgboost_wrapper import XGBoostClassifier\n",
    "from sklearn.externals import joblib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from mlp import make_predictions, train_mlp\n",
    "from utils import printFrame,makePlotName,makeSigBkg,saveFig, makeROC\n",
    "import pdb\n",
    "from xgboost_wrapper import XGBoostClassifier\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import os\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting some global variables (features, classifier used to approximate the ratios)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "workspace_file = 'workspace_evaluating.root'\n",
    "#vars_g = ['x']\n",
    "vars_g = ['x0','x1','x2','x3','x4','x5','x6','x7','x8','x9']\n",
    "dir='evaluating'\n",
    "model_g='xgboost'\n",
    "epochs_training = 40\n",
    "clf = XGBoostClassifier(num_class=2, nthread=4, silent=1,\n",
    "          num_boost_round=epochs_training, eta=0.5, max_depth=3)\n",
    "verbose_printing=False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the models from where the data will be sampled. In this case 2 multivariate normal distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating initial distributions\n"
     ]
    }
   ],
   "source": [
    "mu_l = []\n",
    "cov_l = []\n",
    "\n",
    "mu_l.append([5.,5.,4.,3.,5.,5.,4.5,2.5,4.,3.5])\n",
    "#mu_g.append([7.,8.,7.,6.,7.,8.,6.5,5.5,7.,6.5])\n",
    "mu_l.append([4.,4.5,3.6,4.,6.,4.5,4.2,2.2,4.1,3.3])\n",
    "mu_l.append([1.,0.5,0.3,0.5,0.6,0.4,0.1,0.2,0.1,0.3])\n",
    "\n",
    "cov_l.append([[3.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,2.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,14.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,6.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,17.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,10.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,5.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,1.3,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,1.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,0.,9.3]])\n",
    "cov_l.append([[3.5,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,3.5,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,9.5,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,7.2,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,4.5,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,4.5,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,8.2,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,9.5,3.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,3.5,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,0.,4.5]])\n",
    "cov_l.append([[13.,0.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,12.,0.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,14.,0.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,6.,0.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,1.,0.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,10.,0.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,15.,0.,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,6.3,0.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,11.,0.],\n",
    "         [0.,0.,0.,0.,0.,0.,0.,0.,0.,1.3]])\n",
    "# Statistical model\n",
    "w = ROOT.RooWorkspace('w')\n",
    "\n",
    "print 'Generating initial distributions'\n",
    "cov_m = []\n",
    "mu_m = []\n",
    "mu_str = []\n",
    "cov_root = []\n",
    "vec = []\n",
    "argus = ROOT.RooArgList()\n",
    "#features\n",
    "for i,var in enumerate(vars_g):\n",
    "    w.factory('{0}[{1},{2}]'.format(var,-25,30))\n",
    "    argus.add(w.var(var))\n",
    "\n",
    "for glob in range(2):\n",
    "    # generate covariance matrix\n",
    "    cov_m.append(np.matrix(cov_l[glob]))\n",
    "    cov_root.append(ROOT.TMatrixDSym(len(vars_g)))\n",
    "    for i,var1 in enumerate(vars_g):\n",
    "      for j,var2 in enumerate(vars_g):\n",
    "        cov_root[-1][i][j] = cov_m[-1][i,j]\n",
    "    getattr(w,'import')(cov_root[-1],'cov{0}'.format(glob))\n",
    "    # generate mu vector\n",
    "    mu_m.append(np.array(mu_l[glob]))\n",
    "    vec.append(ROOT.TVectorD(len(vars_g)))\n",
    "    for i, mu in enumerate(mu_m[-1]):\n",
    "      vec[-1][i] = mu\n",
    "    mu_str.append(','.join([str(mu) for mu in mu_m[-1]]))\n",
    "    # multivariate gaussian\n",
    "    gaussian = ROOT.RooMultiVarGaussian('f{0}'.format(glob),\n",
    "          'f{0}'.format(glob),argus,vec[-1],cov_root[-1])\n",
    "    getattr(w,'import')(gaussian)\n",
    "  # Check Model\n",
    "w.Print()\n",
    "\n",
    "\n",
    "if verbose_printing == True:\n",
    "    printFrame(w,['x0','x1','x7','x8'],[w.pdf('f0'),w.pdf('f1')],'distributions',['f0','f1']\n",
    "    ,dir=dir,model_g=model_g,range=[-15,20],title='Distributions',x_text='x0',y_text='p(x)',print_pdf=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sampling data from the distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Making Data\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD9CAYAAABTJWtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UU3eeN/DPpUltuzUgtASaSxuPXKSxUVKdmPNMPZtu\nBNGZBjz0ATk7ENs4pwe2swz1zINyTlvsnhX67HFdx2l8embSbaTPMbB2BWZXMyhjZlt3hacMbFmj\n43UWNAkhnUrFtKUCcp8/7L3FiPyG8OP9Oud7mny53+v3m2refL/3FyMIAgEAABARRUW6AwAAMH8g\nFAAAQIJQAAAACUIBAAAkCAUAAJAgFAAAQDKhULh9+/YDOp2u7YUXXvg1EVFvb29senr66ZSUlMsZ\nGRmNN27ciBG3rays3MtxHJ+amnqpsbExQ6xvbW1dr9VqOziO40tKSg6J9bdu3VqWl5dXw3EcbzAY\nzl+9evWpmRwgAABM3IRC4dChQyUajcbDMIxARFRVVbUnPT399OXLl1NMJlNTVVXVHiIij8ejqamp\nyfN4PBqXy5VZXFxsEwSBISIqKio6YrfbrTzPczzPcy6XK5OIyG63W+Pi4q7zPM+VlpYeLCsre3u2\nBgsAAOMQBGHM4vV6WZPJdOa3v/3t8z/84Q9/LQgCrV69+lJPT49SEAQKBAIJq1evviQIAu3fv39v\nVVVVmdh2y5Ytrv/4j/8wdHd3J6ampl4U648dO7bjlVde+T/iNufPn98oCAINDg7KHnvssT+N1ycU\nFBQUlNkpsvFCo7S09ODf/d3f/ezmzZsKsS4YDCqVSmWQiEipVAaDwaCSiKi7u/sJg8FwXtyOZVmf\n3+9XyeXyQZZlfWK9SqXy+/1+FRGR3+9XJSUleYmIZDLZUHR0dF9vb29sbGxsr7i9OEMBAIDJEb5d\nrZmoMZeP/uVf/uWH8fHxn+l0urb77ZhhGGEuvrQjnZ7zpbz55psR78N8Kfgs8Fngsxi7TMWYM4V/\n//d//x8NDQ3mkydPbvvmm28eunnzpqKgoKBaqVQGe3p6EhISEnoCgUBifHz8Z0R3ZgBerzdJbO/z\n+ViWZX0qlcrv8/nY8HqxzbVr15584oknuoeGhmR9fX3RI2cJAAAwd8acKezfv7/c6/UmdXZ2rnQ6\nnTv+4i/+4rfV1dUFZrO5weFwWIiIHA6HJTs7u46IyGw2Nzidzh0DAwMPdnZ2ruR5ntPr9S0JCQk9\nCoXiZnNz80ZBEJjq6uqCrKyserGNuK/jx4+/aDKZmmZ70AAAMLpxjymMJC4T7dmzpyo3N7fWbrdb\n1Wp1V21tbS4RkUaj8eTm5tZqNBqPTCYbstlsxWIbm81WvHPnzvf7+/sf3rZt28nMzEwXEZHVarUX\nFBRUcxzHx8XFXXc6nTtmepCLidFojHQX5g18Ft/BZ/EdfBbTw0x13WkuMQwjLIR+AgDMJwzDkDCT\nB5oBAGBpQSgAAIAEoQAQYYoYBTEMc1dRxCjGbwgwC3BMASDCGIYhqgirrKApn2cOIMIxBQAAmBaE\nAsB8FEVYUoKImNR1CgAwPYoYBYX6QuNvOEz3LCmFKibQDmCaEAoAcyjUFxr1+AHAfIHlI4AFSqGI\nvXeJSREb6W7BAodQAFigQqEviEi4q4RCIYQETAuWjwAWlSG6ExB3hEKTOhsRAKEAsJAwDL7kYXYh\nFAAWlJEXtCEgYObhmAIAAEgQCgAAIEEoACwUUUR3lozEQkRRuMoZZhaOKQAsFKNc5Uy4yhlmGGYK\nAAAgQSgAzKLwZyUAzHdjhsI333zz0MaNG5vT0tLaNRqNZ+/evZVERBUVFRUsy/p0Ol2bTqdrO3Xq\n1FaxTWVl5V6O4/jU1NRLjY2NGWJ9a2vreq1W28FxHF9SUnJIrL9169ayvLy8Go7jeIPBcP7q1atP\nzcZAASJButeRWADmuTFD4aGHHvrm7Nmzz7e3t6d9+umna8+ePfv8xx9//BzDMMJrr732921tbbq2\ntjbd1q1bTxEReTweTU1NTZ7H49G4XK7M4uJim/iAh6KioiN2u93K8zzH8zzncrkyiYjsdrs1Li7u\nOs/zXGlp6cGysrK3Z3/YAEuFDLe9gEkZd/nokUce+ZqIaGBg4MHbt28/sGLFii+IaNSn+dTX12fl\n5+cfk8vlg2q1uis5OflKc3PzxkAgkBgKhZbr9foWIqLCwsKjdXV12UREDQ0NZovF4iAiysnJ+bCp\nqck0kwMEWNrE216I90b6IsL9gflu3LOPhoeHo5599tnf//GPf1xVVFR0ZM2aNReOHz/+4uHDh39y\n9OjRwg0bNnxy4MCB3TExMTe6u7ufMBgM58W2LMv6/H6/Si6XD7Is6xPrVSqV3+/3q4iI/H6/Kikp\nyUtEJJPJhqKjo/t6e3tjY2Nje0f2o6KiQnptNBrJaDROe/AAAIuJ2+0mt9s9rX2MGwpRUVHD7e3t\naX19fdFbtmz5jdvtNhYVFR1544033iIiev311/9m9+7dB+x2u3VaPRnHyFAAAIB7hf/CvG/fvknv\nY8JnH0VHR/f94Ac/+NdPPvlkQ3x8/GcMwwgMwwi7du36VUtLi57ozgzA6/UmiW18Ph/LsqxPpVL5\nfT4fG14vtrl27dqTRERDQ0Oyvr6+6PBZAgAAzI0xQ+Hzzz9/7MaNGzFERP39/Q+fPn06XafTtfX0\n9CSI25w4cWK7VqvtICIym80NTqdzx8DAwIOdnZ0reZ7n9Hp9S0JCQo9CobjZ3Ny8URAEprq6uiAr\nK6tebONwOCxERMePH3/RZDI1zd5wAQBgLGMuHwUCgUSLxeIYHh6OGh4ejiooKKg2mUxNhYWFR9vb\n29MYhhFWrlzZ+e67775CRKTRaDy5ubm1Go3GI5PJhmw2WzHDMAIRkc1mK965c+f7/f39D2/btu1k\nZmami4jIarXaCwoKqjmO4+Pi4q47nc4dsz9sAAAYDSMIwvhbRRjDMMJC6CdAOIZh7r4+oYJGf0bz\ndOruuZ322O/xb2npYBhm1DNFx4IrmgEAQIJQAAAACUIBAAAkCAUAAJAgFAAWMjx4B2YYHrIDsJDh\nwTswwzBTAAAACUIBAAAkCAWAGRL+lLX5+aQ12T19xDMWYCQcUwCYIdJT1kYKfx9x4vMVvhMKzcfw\ngkjBTAEAACQIBQAAkCAUAABAglAAAAAJQgEAACQIBQAAkCAUAABAglAAAADJmKHwzTffPLRx48bm\ntLS0do1G49m7d28lEVFvb29senr66ZSUlMsZGRmNN27ciBHbVFZW7uU4jk9NTb3U2NiYIda3trau\n12q1HRzH8SUlJYfE+lu3bi3Ly8ur4TiONxgM569evfrUbAwUAADGN2YoPPTQQ9+cPXv2+fb29rRP\nP/107dmzZ5//+OOPn6uqqtqTnp5++vLlyykmk6mpqqpqDxGRx+PR1NTU5Hk8Ho3L5cosLi62ic8H\nLSoqOmK32608z3M8z3MulyuTiMhut1vj4uKu8zzPlZaWHiwrK3t79ocNAACjGXf56JFHHvmaiGhg\nYODB27dvP7BixYovGhoazBaLxUFEZLFYHHV1ddlERPX19Vn5+fnH5HL5oFqt7kpOTr7S3Ny8MRAI\nJIZCoeV6vb6FiKiwsPCo2GbkvnJycj5samoyzdZgAQBgbOPe+2h4eDjq2Wef/f0f//jHVUVFRUfW\nrFlzIRgMKpVKZZCISKlUBoPBoJKIqLu7+wmDwXBebMuyrM/v96vkcvkgy7I+sV6lUvn9fr+KiMjv\n96uSkpK8REQymWwoOjq6r7e3NzY2NrZ3ZD8qKiqk10ajkYxG47QGDgCw2LjdbnK73dPax7ihEBUV\nNdze3p7W19cXvWXLlt+cPXv2+ZE/ZxhGYBhGuF/7mTIyFAAA4F7hvzDv27dv0vuY8NlH0dHRfT/4\nwQ/+tbW1db1SqQz29PQkEBEFAoHE+Pj4z4juzAC8Xm+S2Mbn87Esy/pUKpXf5/Ox4fVim2vXrj1J\nRDQ0NCTr6+uLDp8lAADA3BgzFD7//PPHxDOL+vv7Hz59+nS6TqdrM5vNDQ6Hw0JE5HA4LNnZ2XVE\nRGazucHpdO4YGBh4sLOzcyXP85xer29JSEjoUSgUN5ubmzcKgsBUV1cXZGVl1YttxH0dP378RZPJ\n1DS7QwZY5O56bvO3r/HcZpigMZePAoFAosVicQwPD0cNDw9HFRQUVJtMpiadTteWm5tba7fbrWq1\nuqu2tjaXiEij0Xhyc3NrNRqNRyaTDdlstmJxaclmsxXv3Lnz/f7+/oe3bdt2MjMz00VEZLVa7QUF\nBdUcx/FxcXHXnU7njtkfNsAihuc2wzQwgjDrhwOmjWEYYSH0E5Y2hmFGf8hOxRjv57KOBLozewj/\nt8QQ/n0tTgzDkHhZwEThimYAAJAgFAAAQIJQAAAACUIBYMmTEcMwUlEoYiPdIYigcS9eA4DFbohG\nHnwOhSZ1XBIWGcwUAABAglAAAAAJQgFgChQxirvW4RkGSy6wOOCYAsAUhPpC97lADGBhw0wBAAAk\nCAUAAJAgFAAAQIJQAAAACUIBAAAkCAUAAJAgFAAAQIJQAAAACUIBAAAkCAUAAJCMGQperzfp+eef\nP7tmzZoLzzzzzH/9/Oc//2siooqKigqWZX06na5Np9O1nTp1aqvYprKyci/HcXxqauqlxsbGDLG+\ntbV1vVar7eA4ji8pKTkk1t+6dWtZXl5eDcdxvMFgOH/16tWnZmOgAAAwvjFDQS6XDx48eLD0woUL\na86fP2945513/urixYtPMwwjvPbaa3/f1tama2tr023duvUUEZHH49HU1NTkeTwejcvlyiwuLraJ\nD40uKio6YrfbrTzPczzPcy6XK5OIyG63W+Pi4q7zPM+VlpYeLCsre3v2hw0AAKMZMxQSEhJ60tLS\n2omIHn300S+ffvrpi36/X0VEJH7Zj1RfX5+Vn59/TC6XD6rV6q7k5OQrzc3NGwOBQGIoFFqu1+tb\niIgKCwuP1tXVZRMRNTQ0mC0Wi4OIKCcn58OmpibTTA8SAAAmZsJ3Se3q6lK3tbXpDAbD+XPnzn3/\n8OHDPzl69Gjhhg0bPjlw4MDumJiYG93d3U8YDIbzYhuWZX1+v18ll8sHWZb1ifUqlcovhovf71cl\nJSV5iYhkMtlQdHR0X29vb2xsbGzvyD+/oqJCem00GsloNE550AAAi5Hb7Sa32z2tfUwoFL788stH\nX3zxxeOHDh0qefTRR78sKio68sYbb7xFRPT666//ze7duw/Y7XbrtHoyjpGhAAAA9wr/hXnfvn2T\n3se4Zx8NDg7Kc3JyPvzRj370QXZ2dh0RUXx8/GcMwwgMwwi7du36VUtLi57ozgzA6/UmiW19Ph/L\nsqxPpVL5fT4fG14vtrl27dqTRERDQ0Oyvr6+6PBZAgBMUxQRkbjiyxBFKSLYGZjPxgwFQRAYq9Vq\n12g0np/+9Kf/INYHAoFE8fWJEye2a7XaDiIis9nc4HQ6dwwMDDzY2dm5kud5Tq/XtyQkJPQoFIqb\nzc3NGwVBYKqrqwuysrLqxTYOh8NCRHT8+PEXTSZT0+wMFWAJG6Y7DwESy3BojI1l9zxVTqGInf0+\nwrww5vLRuXPnvv/BBx/8aO3atZ/qdLo2IqL9+/eXHzt2LL+9vT2NYRhh5cqVne++++4rREQajcaT\nm5tbq9FoPDKZbMhmsxUzDCMQEdlstuKdO3e+39/f//C2bdtOZmZmuoiIrFarvaCgoJrjOD4uLu66\n0+ncMduDBoCxDBGRcFdNKITHjS4VY4bCc8899/Hw8PA9swnxFNTRlJeX7y8vL98fXr9+/frWjo4O\nbXj9smXLbtXW1uZOtMMAADB7cEUzAABIEAoAACBBKAAAgAShADAORYzinrNxABarCV/RDLBUhfpC\nd07jHCn8PcAigZkCAABIEAoAACBBKAAAgAShAAAAEoQCAABIEAoAACBBKAAAgAShAAAAEoQCAABI\nEAoAACBBKAAAgAShAAAAEoQCAABIEAoAACAZMxS8Xm/S888/f3bNmjUXnnnmmf/6+c9//tdERL29\nvbHp6emnU1JSLmdkZDTeuHEjRmxTWVm5l+M4PjU19VJjY2OGWN/a2rpeq9V2cBzHl5SUHBLrb926\ntSwvL6+G4zjeYDCcv3r16lOzMVAAABjfmKEgl8sHDx48WHrhwoU158+fN7zzzjt/dfHixaerqqr2\npKenn758+XKKyWRqqqqq2kNE5PF4NDU1NXkej0fjcrkyi4uLbYIgMERERUVFR+x2u5XneY7nec7l\ncmUSEdntdmtcXNx1nue50tLSg2VlZW/P/rABlrgoIiLm2/Lt6yhF5PoD88aYoZCQkNCTlpbWTkT0\n6KOPfvn0009f9Pv9qoaGBrPFYnEQEVksFkddXV02EVF9fX1Wfn7+MblcPqhWq7uSk5OvNDc3bwwE\nAomhUGi5Xq9vISIqLCw8KrYZua+cnJwPm5qaTLM5YAAgomG686CgkWU4NEYD2V1PnlMoYme7hxAh\nE37yWldXl7qtrU23cePG5mAwqFQqlUEiIqVSGQwGg0oiou7u7icMBsN5sQ3Lsj6/36+Sy+WDLMv6\nxHqVSuX3+/0qIiK/369KSkryEhHJZLKh6Ojovt7e3tjY2NjekX9+RUWF9NpoNJLRaJzSgAFgKoaI\nSJDehUJ4JOl85Ha7ye12T2sfEwqFL7/88tGcnJwPDx06VLJ8+fK7fp1gGEZgGEa4X9uZMjIUAADg\nXuG/MO/bt2/S+xj37KPBwUF5Tk7OhwUFBdXZ2dl1RHdmBz09PQlERIFAIDE+Pv4zojszAK/XmyS2\n9fl8LMuyPpVK5ff5fGx4vdjm2rVrTxIRDQ0Nyfr6+qLDZwkAADA3xgwFQRAYq9Vq12g0np/+9Kf/\nINabzeYGh8NhISJyOBwWMSzMZnOD0+ncMTAw8GBnZ+dKnuc5vV7fkpCQ0KNQKG42NzdvFASBqa6u\nLsjKyqoP39fx48dfNJlMTbM3XAAAGMuYy0fnzp37/gcffPCjtWvXfqrT6dqI7pxyumfPnqrc3Nxa\nu91uVavVXbW1tblERBqNxpObm1ur0Wg8MplsyGazFYtLSzabrXjnzp3v9/f3P7xt27aTmZmZLiIi\nq9VqLygoqOY4jo+Li7vudDp3zPagAQBgdGOGwnPPPffx8PDwqLOJM2fObB6tvry8fH95efn+8Pr1\n69e3dnR0aMPrly1bdksMFQAAiCxc0QwwgiJGcdeplwyDs2xgaZnwKakAS0GoL3TnnP2Rwt8DLGKY\nKQAAgAShAAAAEoQCAABIEAoAACBBKAAAgAShAAAAEoQCAABIEAoAACBBKAAAgAShAAAAEoQCAABI\nEAoAACBBKAAAgAShAAAAEoQCAEyB7K5nTigUsZHuEMwQPE8BAO6IIqJh8aFC3/43ajnR8M1RNh4i\nIkF6FwrhYUSLBUIBAO4YplEeMBSKQEcgksZcPnr55ZffUyqVQa1W2yHWVVRUVLAs69PpdG06na7t\n1KlTW8WfVVZW7uU4jk9NTb3U2NiYIda3trau12q1HRzH8SUlJYfE+lu3bi3Ly8ur4TiONxgM569e\nvfrUTA8QAAAmbsxQeOmll/7R5XJljqxjGEZ47bXX/r6trU3X1tam27p16ykiIo/Ho6mpqcnzeDwa\nl8uVWVxcbBMEgSEiKioqOmK32608z3M8z3PiPu12uzUuLu46z/NcaWnpwbKysrdna6AAADC+MUNh\n06ZNH61YseKL8Hrxy36k+vr6rPz8/GNyuXxQrVZ3JScnX2lubt4YCAQSQ6HQcr1e30JEVFhYeLSu\nri6biKihocFssVgcREQ5OTkfNjU1mWZmWAAAMBVTOqZw+PDhnxw9erRww4YNnxw4cGB3TEzMje7u\n7icMBsN5cRuWZX1+v18ll8sHWZb1ifUqlcrv9/tVRER+v1+VlJTkJSKSyWRD0dHRfb29vbGxsbG9\n4X9mRUWF9NpoNJLRaJxK1wEkihgFhfqwZg6Lh9vtJrfbPa19TDoUioqKjrzxxhtvERG9/vrrf7N7\n9+4DdrvdOq1eTMDIUACYCaG+0CgHViPQEYAZEv4L8759+ya9j0lfpxAfH/8ZwzACwzDCrl27ftXS\n0qInujMD8Hq9SeJ2Pp+PZVnWp1Kp/D6fjw2vF9tcu3btSSKioaEhWV9fX/RoswQAAJgbkw6FQCCQ\nKL4+ceLEdvHMJLPZ3OB0OncMDAw82NnZuZLneU6v17ckJCT0KBSKm83NzRsFQWCqq6sLsrKy6sU2\nDofDQkR0/PjxF00mU9NMDQwAACZvzOWj/Pz8Y7/73e/+/PPPP38sKSnJu2/fvjfdbrexvb09jWEY\nYeXKlZ3vvvvuK0REGo3Gk5ubW6vRaDwymWzIZrMVMwwjEBHZbLbinTt3vt/f3//wtm3bTmZmZrqI\niKxWq72goKCa4zg+Li7uutPp3DH7QwYAgPsZMxSOHTuWH1738ssvv3e/7cvLy/eXl5fvD69fv359\na0dHhza8ftmyZbdqa2tzJ9pZAACYXbj3EQAASBAKAAAgQSgAAIAEoQAAABKEAgAASBAKAAAgQSgA\nAIAEoQAAABKEAgAASBAKADADZMQwzF1FoYiNdKdgCvCMZgC4vygiGhafqfXtf6OW33me812GiEi4\nqyYUuudZXLAAIBQA4P6GaZRnTuDBRIsZlo8AAECCUAAAAAlCAQAAJAgFAACQIBQAAECCUIAlQRGj\nuOc8egC4F05JhSUh1Bca5dTKCHQEYJ4bc6bw8ssvv6dUKoNarbZDrOvt7Y1NT08/nZKScjkjI6Px\nxo0bMeLPKisr93Icx6empl5qbGzMEOtbW1vXa7XaDo7j+JKSkkNi/a1bt5bl5eXVcBzHGwyG81ev\nXn1qpgcIAAATN2YovPTSS//ocrkyR9ZVVVXtSU9PP3358uUUk8nUVFVVtYeIyOPxaGpqavI8Ho/G\n5XJlFhcX2wRBYIiIioqKjtjtdivP8xzP85y4T7vdbo2Li7vO8zxXWlp6sKys7O3ZGigAAIxvzFDY\ntGnTRytWrPhiZF1DQ4PZYrE4iIgsFoujrq4um4iovr4+Kz8//5hcLh9Uq9VdycnJV5qbmzcGAoHE\nUCi0XK/XtxARFRYWHhXbjNxXTk7Oh01NTabZGCQAAEzMpI8pBINBpVKpDBIRKZXKYDAYVBIRdXd3\nP2EwGM6L27Es6/P7/Sq5XD7IsqxPrFepVH6/368iIvL7/aqkpCQvEZFMJhuKjo7u6+3tjY2Nje0N\n/3MrKiqk10ajkYxG42S7DgCwqLndbnK73dPax7QONDMMIzAMI4y/5fSNDAUAALhX+C/M+/btm/Q+\nJn1KqlKpDPb09CQQEQUCgcT4+PjPiO7MALxeb5K4nc/nY1mW9alUKr/P52PD68U2165de5KIaGho\nSNbX1xc92iwBAADmxqRDwWw2NzgcDgsRkcPhsGRnZ9eJ9U6nc8fAwMCDnZ2dK3me5/R6fUtCQkKP\nQqG42dzcvFEQBKa6urogKyurPnxfx48ff9FkMjXN5OAAAGByxlw+ys/PP/a73/3uzz///PPHkpKS\nvG+99dYbe/bsqcrNza212+1WtVrdVVtbm0tEpNFoPLm5ubUajcYjk8mGbDZbsbi0ZLPZinfu3Pl+\nf3//w9u2bTuZmZnpIiKyWq32goKCao7j+Li4uOtOp3PH7A8ZAOaG7K6LBJcvX0E3b2IhYL4bMxSO\nHTuWP1r9mTNnNo9WX15evr+8vHx/eP369etbOzo6tOH1y5YtuyWGCgAsNnc/eAcP3VkYcJsLAACQ\nIBQAYHKkbw3mu4JvkkUD9z4CgMkZ9RGdc98NmB3IdwAAkCAUYNHBbbIBpg7LR7Do4DbZAFOHmQIA\nAEgQCgAAIEEoAACABKEAAAAShAIAAEgQCgAAIEEoAACABKEAAHNEdtcFhQpFbKQ7BKPAxWsAMEdw\nK+2FADMFAJi+KKK77pqKO6cuWJgpAMD04c6piwayHAAAJAgFAACQTDkU1Gp119q1az/V6XRter2+\nhYiot7c3Nj09/XRKSsrljIyMxhs3bsSI21dWVu7lOI5PTU291NjYmCHWt7a2rtdqtR0cx/ElJSWH\npjccWIrCb5UNAFM35VBgGEZwu93GtrY2XUtLi56IqKqqak96evrpy5cvp5hMpqaqqqo9REQej0dT\nU1OT5/F4NC6XK7O4uNgmCAJDRFRUVHTEbrdbeZ7neJ7nXC5X5swMDZYK6VbZYgGAKZvW8pH4xS5q\naGgwWywWBxGRxWJx1NXVZRMR1dfXZ+Xn5x+Ty+WDarW6Kzk5+Upzc/PGQCCQGAqFloszjcLCwqNi\nGwAAmHtTPvuIYRhh8+bNZx544IHbr7zyyrs//vGPfxkMBpVKpTJIRKRUKoPBYFBJRNTd3f2EwWA4\nL7ZlWdbn9/tVcrl8kGVZn1ivUqn8fr9fNdqfV1FRIb02Go1kNBqn2nUAgEXJ7XaT2+2e1j6mHArn\nzp37fmJiYuBPf/rT4+np6adTU1Mvjfw5wzACwzDC/dpP1shQAACAe4X/wrxv375J72PKy0eJiYkB\nIqLHH3/8T9u3bz/R0tKiVyqVwZ6engQiokAgkBgfH/8Z0Z0ZgNfrTRLb+nw+lmVZn0ql8vt8PnZk\nvUql8k+1TwCwkMjueZY2bn0ReVMKha+//vqRUCi0nIjoq6+++rPGxsYMrVbbYTabGxwOh4WIyOFw\nWLKzs+uIiMxmc4PT6dwxMDDwYGdn50qe5zm9Xt+SkJDQo1AobjY3N28UBIGprq4uENsAwGIn3vbi\nuxIKfRHZLsHUlo+CwaBy+/btJ4iIhoaGZH/5l3/5fzMyMho3bNjwSW5ubq3dbreq1equ2traXCIi\njUbjyc3NrdVoNB6ZTDZks9mKxaUlm81WvHPnzvf7+/sf3rZt28nMzEzXzA0PACImioiGw04RjlIQ\nDd+MSHdgYqYUCitXruxsb29PC6+PjY3tPXPmzObR2pSXl+8vLy/fH16/fv361o6ODu1U+gEA89io\nt74IRaAjMBm4ohkAACQIBVhQwq9exhXMADMLd0mFBUW6enmk8PcAMGWYKQAAgAShAAAAEoQCAMwj\neI5zpOHLk5E3AAAGRklEQVSYAgDMnfBrF+65bgHPcY40hAIAzJ3waxdw3cK8g+UjmLdw+inA3MNM\nAeYtnH4KMPcwUwAAAAlCAQDmMdxee65h+QgAImfcO6nefTYSEc5Imm0IBZgXFDGKO8cQYGnBnVTn\nHYQCzAs4qAwwP+CYAgAsMLjqeTZhpgAA88s9xxkYoqjl9z3OgGMMMwuhAHMOxw9gTDjOEFFYPlpg\n3G53pLswbdLxg5FlKrpmpDuLQ1ekOxBJdy8nPfLI8kh3aEGbF6HgcrkyU1NTL3Ecx7/99ttlke7P\nfLbQQmFWb1XRNXO7WvC6It2BWRZFRMR8W759HaX49rW4nHSn9Pd/jWsbpiHiy0e3b99+4NVXX/3F\nmTNnNqtUKv/3vve9/2c2mxuefvrpi5HuG0zOfZeFKsZ5DzCe0ZaU3grRXSFB9O2xh2G699oG+V2/\nkCxfvoJu3uydnb4ucBEPhZaWFn1ycvIVtVrdRUS0Y8cOZ319fRZCYX4Z9Qv/ASK6HbZhxTjvAWbK\nfYOC6LuwoG8PXN+9WeirL+6etY6yzZINDkEQIlr+6Z/+6cVdu3b9UnxfXV39o1dfffXwyG1o5NwQ\nBQUFBWXCZbLfyRGfKTAMI4y3jSAIOOcMAGAORPxAs0ql8nu93iTxvdfrTWJZ1hfJPgEALFURD4UN\nGzZ8wvM819XVpR4YGHiwpqYmz2w2N0S6XwAAS1HEl49kMtnQL37xi1e3bNnym9u3bz9gtVrtOMgM\nABAhkT7QPFapra39nxqN5kJUVNTt1tbWZ0f+bP/+/XuTk5P51atXX/rNb36TEem+zmV58803K1Qq\nlS8tLa0tLS2t7dSpU5mR7tNcl1OnTmWuXr36UnJyMl9VVVUW6f5Esjz11FNdWq3207S0tLbvfe97\nLZHuz1yWl1566b34+PjgM8880yHWXb9+PXbz5s2nOY67nJ6e3vjFF1/ERLqfkfospvJdEfGBjFUu\nXryY+oc//CHFaDSeHRkKFy5c0Kxbt659YGBA3tnZqV61atWV27dvR0W6v3NVKioq3jxw4MBrke5H\npMrQ0NADq1atutLZ2akeGBiQr1u3rt3j8Twd6X5FqqjV6s7r16/HRrofkSj/9m//tun3v/+9buQX\n4c9+9rP//fbbb/8vQRCoqqqqrKysrCrS/YzUZzGV74qIH1MYS2pq6qWUlJTL4fX19fVZ+fn5x+Ry\n+aBare5KTk6+0tLSoo9EHyNFWMJnZI28tkUulw+K17ZEul+RtFT/PmzatOmjFStWfDGyrqGhwWyx\nWBxERBaLxVFXV5cdmd7NrdE+C6LJ/92Y16FwP93d3U+MPEOJZVmf3+9XRbJPc+3w4cM/Wbdu3X9a\nrVb7jRs3YiLdn7nk9/tVSUlJXvH9Uvz/PxLDMMLmzZvPbNiw4ZNf/vKXP450fyItGAwqlUplkIhI\nqVQGg8GgMtJ9iqTJfldEPBTS09NPa7XajvDy61//+oXJ7Gci1zssJPf7XBoaGsxFRUVHOjs7V7a3\nt6clJiYGdu/efSDS/Z1Li+3/9XSdO3fu+21tbbpTp05tfeedd/7qo48+2hTpPs0XDMMIS/nvy1S+\nKyJ+9tHp06fTJ9sm/NoGn8/HqlQq/8z2LLIm+rns2rXrVy+88MKvZ7s/8wmubblbYmJigIjo8ccf\n/9P27dtPtLS06Ddt2vRRpPsVKUqlMtjT05OQkJDQEwgEEuPj4z+LdJ8iZeTYJ/pdEfGZwkSNXBcz\nm80NTqdzx8DAwIOdnZ0reZ7n9Hp9SyT7N5cCgUCi+PrEiRPbtVptRyT7M9dwbct3vv7660dCodBy\nIqKvvvrqzxobGzOW2t+HcGazucHhcFiIiBwOhyU7O7su0n2KlCl9V0T6iPlY5Z//+Z+3syzrfeih\nh/qVSmVPZmbmKfFnf/u3f1u+atWqK6tXr77kcrm2RLqvc1kKCgqOarXaT9euXfufWVlZdT09PcpI\n92muy8mTJ7empKT8YdWqVVf279+/N9L9iVT57//+75Xr1q1rX7duXfuaNWv+a6l9Fjt27DiWmJjY\nLZfLB1iW9b733nsvXb9+PdZkMp1Zaqekhn8Wdrv95al8VzCCsGSX2wAAIMyCWT4CAIDZh1AAAAAJ\nQgEAACQIBQAAkCAUAABAglAAAADJ/wcx/WRww8r7zQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x155e8950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Start generating data\n",
    "num_train = 500000\n",
    "num_test = 5000\n",
    "no_train = False\n",
    "\n",
    "print 'Making Data'\n",
    "# Start generating data\n",
    "\n",
    "   \n",
    "# make data from root pdf\n",
    "def makeDataFi(x, pdf, num):\n",
    "    traindata = np.zeros((num,len(vars_g))) \n",
    "    data = pdf.generate(x,num)\n",
    "    traindata[:] = [[data.get(i).getRealValue(var) for var in vars_g]\n",
    "        for i in range(num)]\n",
    "    return traindata\n",
    "  \n",
    "  # features\n",
    "vars = ROOT.TList()\n",
    "for var in vars_g:\n",
    "    vars.Add(w.var(var))\n",
    "x = ROOT.RooArgSet(vars)\n",
    "\n",
    "\n",
    "train_data = np.zeros((num_train*2,len(vars_g) + 1))\n",
    "test_data = np.zeros((num_test*2,len(vars_g) + 1))\n",
    "\n",
    "if not no_train:\n",
    "    #traindata[:num_train,0] =  makeDataFi(x,w.pdf('f0'), num_train).reshape(num_train)\n",
    "\n",
    "    #traindata[num_train:,0] = makeDataFi(x,w.pdf('f1'), num_train).reshape(num_train)\n",
    "    data_f0 =  makeDataFi(x,w.pdf('f0'), num_train)\n",
    "    data_f1 = makeDataFi(x,w.pdf('f1'), num_train)\n",
    "\n",
    "    train_data[:num_train,:len(vars_g)] =  data_f0\n",
    "\n",
    "    train_data[num_train:,:len(vars_g)] = data_f1\n",
    " \n",
    "    train_data[:num_train,-1] = np.ones(num_train)\n",
    "\n",
    "#testdata[:num_test,0] = makeDataFi(x, w.pdf('f0'), num_test).reshape(num_test)\n",
    "#testdata[num_test:,0] = makeDataFi(x, w.pdf('f1'), num_test).reshape(num_test)\n",
    "test_data[:num_test,:len(vars_g)] =  makeDataFi(x,w.pdf('f0'), num_test)\n",
    "\n",
    "test_data[num_test:,:len(vars_g)] = makeDataFi(x,w.pdf('f1'), num_test)\n",
    "\n",
    "test_data[num_test:,-1] = np.ones(num_test)\n",
    "\n",
    "plt.hist(data_f0[:,0],bins=50)\n",
    "plt.hist(data_f1[:,0],bins=50)\n",
    "plt.plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tran classifiers on the sampled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 2.53032553995 (+/- 0.0135329382044)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "model_file = 'train'\n",
    "traindata = train_data[:,:-1]\n",
    "targetdata = train_data[:,-1]\n",
    "if model_g == 'mlp':\n",
    "    train_mlp((traindata, targetdata), save_file='{0}/{1}_F0_F1.pkl'.format(dir,model_file),n_epochs=epochs_training)\n",
    "else:\n",
    "    rng = np.random.RandomState(1234)\n",
    "    indices = rng.permutation(traindata.shape[0])\n",
    "    traindata_random = traindata[indices]\n",
    "    targetdata_random = targetdata[indices]\n",
    "    scores = cross_validation.cross_val_score(clf, traindata_random, targetdata_random)\n",
    "    print \"Accuracy: {0} (+/- {1})\".format(scores.mean(), scores.std() * 2)\n",
    "    clf.fit(traindata_random,targetdata_random)\n",
    "    #clf.plot_importance_matrix(vars_names)\n",
    "    joblib.dump(clf, '{0}/{1}_F0_F1.pkl'.format(dir,model_file))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Prediction function \n",
    "def predict(filename, traindata,model_g='mlp', sig=1):\n",
    "  sfilename,k,j = filename.split('/')[-1].split('_')\n",
    "  sfilename = '/'.join(filename.split('/')[:-1]) + '/' +  sfilename\n",
    "  j = j.split('.')[0]\n",
    "  sig = 1\n",
    "  if k <> 'F0':\n",
    "    k = int(k)\n",
    "    j = int(j)\n",
    "    sig = 1 if k < j else 0\n",
    "    filename = '{0}_{1}_{2}.pkl'.format(sfilename,min(k,j),max(k,j))\n",
    "  if model_g == 'mlp':\n",
    "    return make_predictions(dataset=traindata, model_file=filename)[:,sig]\n",
    "  else:\n",
    "    clf = joblib.load(filename)\n",
    "    if clf.__class__.__name__ == 'NuSVR':\n",
    "      output = clf.predict(traindata)\n",
    "      return np.clip(output,0.,1.)\n",
    "    else:\n",
    "      return clf.predict_proba(traindata)[:,sig]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute score distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating Score Histograms\n",
      "Estimating F0 F1\n"
     ]
    }
   ],
   "source": [
    "\n",
    "bins = 80\n",
    "low = 0.\n",
    "high = 1.  \n",
    "c1_g=''\n",
    "  \n",
    "print 'Generating Score Histograms'\n",
    "\n",
    "w.factory('score[{0},{1}]'.format(low,high))\n",
    "s = w.var('score')\n",
    "  \n",
    "def saveHisto(w,outputs,s,bins,low,high,k='F0',j='F1'):\n",
    "    \n",
    "    print 'Estimating {0} {1}'.format(k,j)\n",
    "    for l,name in enumerate(['sig','bkg']):\n",
    "      data = ROOT.RooDataSet('{0}data_{1}_{2}'.format(name,k,j),\"data\",\n",
    "          ROOT.RooArgSet(s))\n",
    "      hist = ROOT.TH1F('{0}hist_{1}_{2}'.format(name,k,j),'hist',bins,low,high)\n",
    "      values = outputs[l]\n",
    "      #values = values[self.findOutliers(values)]\n",
    "      for val in values:\n",
    "        hist.Fill(val)\n",
    "        s.setVal(val)\n",
    "        data.add(ROOT.RooArgSet(s))\n",
    "      norm = 1./hist.Integral()\n",
    "      hist.Scale(norm) \n",
    "        \n",
    "      s.setBins(bins)\n",
    "      datahist = ROOT.RooDataHist('{0}datahist_{1}_{2}'.format(name,k,j),'hist',\n",
    "            ROOT.RooArgList(s),hist)\n",
    "      histpdf = ROOT.RooHistFunc('{0}histpdf_{1}_{2}'.format(name,k,j),'hist',\n",
    "            ROOT.RooArgSet(s), datahist, 1)\n",
    "\n",
    "      getattr(w,'import')(hist)\n",
    "      getattr(w,'import')(data)\n",
    "      getattr(w,'import')(datahist) # work around for morph = w.import(morph)\n",
    "      getattr(w,'import')(histpdf) # work around for morph = w.import(morph)\n",
    "      score_str = 'score'\n",
    "      # Calculate the density of the classifier output using kernel density \n",
    "      #w.factory('KeysPdf::{0}dist_{1}_{2}({3},{0}data_{1}_{2},RooKeysPdf::NoMirror,2)'.format(name,k,j,score_str))\n",
    "\n",
    "\n",
    "numtrain = traindata.shape[0]       \n",
    "\n",
    "outputs = [predict('{0}/{1}_F0_F1.pkl'.format(dir,model_file),traindata[targetdata==1],model_g=model_g),\n",
    "            predict('{0}/{1}_F0_F1.pkl'.format(dir,model_file),traindata[targetdata==0],model_g=model_g)]\n",
    "\n",
    "saveHisto(w,outputs,s, bins, low, high)\n",
    "\n",
    "if verbose_printing == True:\n",
    "    can = printFrame(w,['score'],[w.function('sighistpdf_F0_F1'),w.function('bkghistpdf_F0_F1')], makePlotName('full','all',type='hist',dir=dir,c1_g=c1_g,model_g=model_g),['signal','bkg'],\n",
    "  dir=dir,model_g=model_g,y_text='score(x)',print_pdf=True,title='Pairwise score distributions',use_rootnotes=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def findOutliers(x):\n",
    "  q5, q95 = np.percentile(x, [5,95])  \n",
    "  iqr = 2.*(q95 - q5)\n",
    "  outliers = (x <= q95 + iqr) & (x >= q5 - iqr)\n",
    "  return outliers\n",
    "def singleRatio(f0,f1):\n",
    "  ratio = f1 / f0\n",
    "  ratio[np.abs(ratio) == np.inf] = 0 \n",
    "  ratio[np.isnan(ratio)] = 0\n",
    "  return ratio\n",
    "\n",
    "def evalDist(x,f0,val):\n",
    "  iter = x.createIterator()\n",
    "  v = iter.Next()\n",
    "  i = 0\n",
    "  while v:\n",
    "    v.setVal(val[i])\n",
    "    v = iter.Next()\n",
    "    i = i+1\n",
    "  return f0.getVal(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute approximated density ratios and do reweighting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating ratios\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEgCAYAAACJomVNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcFPX+P/D3LjcRUK5y2UVXBAVNkeSiKboqglghap7w\nmGJhkl/JPJWplbpoWXa6WNLxkHfS1EwFLUVCw9AUBAVMRUFchYVVwBs3gZ2d3x/+lrPCgtz2Bq/n\n4zGPx87nM5+Z94y47535zHyGw7IsAQAAdDautgMAAICuCQkGAADUAgkGAADUAgkGAADUAgkGAADU\nAgkGAADUAgkG1GL37t2zg4KCjrdm2R07dszz9/dPVVcsbV2/QCAQnzx5cgIR0bp16z588803N3dW\nLBYWFhVisVhARDRv3rwdK1euXNtZ6164cOGmTz755OPOWl9rHDp0aJqzs3OhhYVFRXZ2tqcmtw26\nDwkGGnz22WcrpkyZclS5zM3NLU9V2c8///yPltY1e/bs3cePHw/qjLiEQmHK1q1bIzpjXa3B4XAa\nHg778MMP123evPnNZ7VpbYwVFRUWAoFArNiO8rbaQlXS3LRp08KPP/74k/asr73ef//9L//zn//8\nX0VFhYWnp2d2R9bF5XLlBQUFLp0V27MEBwcfs7CwqFBMJiYmtcOGDcvR1Pa7AyQYaDBu3LhTf/31\n1wssy3KIiEpKShxlMplhVlbWcLlczlWU3bhxY8DYsWP/1FRc7f0S1qRnxSiTyQxVlSuOtT5iWZZz\n+/btvoMHD77Smetsb1vF32hrHTt2LLiiosJCMb3wwgt//eMf//i5vduHppBgoIG3t3dGfX29UVZW\n1nAiotTUVP/x48f/MXDgwOvKZQMGDLjh4OAgffjwYe+IiIitTk5OxXw+v2jlypVrFf/JG//CTkpK\nChw0aNA1S0vLB4sWLfp+3Lhxpxr/4l+6dOm/ra2t77m4uBQkJiZOJiL66KOPPk1NTfWPioqKsbCw\nqFi8ePF3RES5ubnukyZN+t3Gxqbc3d09d//+/TMV6ykvL7cJCQk53Lt374d+fn5pN27cGNDSfv/4\n449z+vXrd8vW1rZs3bp1HyrXiUQi0Zw5c34kInr8+HGP1157bZetrW2ZlZXVfV9f3/S7d+/2aS5G\nLpcr/89//vN/bm5ueYMGDbqmKFP+lV5WVmYbGBiY1KtXr0dCoTDl9u3bfYmIxGKxgMvlypW/NBVn\nSbm5ue5vvfXWf8+ePTvKwsKiwtra+h5R00tumzdvftPNzS3PxsamfOrUqQklJSWOijoulyuPjY2N\nHDhw4HUrK6v7UVFRMYq6/Px813Hjxp2ytLR8YGdnVxoWFra38TGrra01sbCwqGAYxsDT0zPbzc0t\nj4iouLjYacaMGQf69Olz18XFpWDjxo1vK9qkp6f7jho16qyVldV9Jyen4rfffntjfX29ERGR4geL\np6dntoWFRcXPP//8D1VnacrHb968eTsWLly4acqUKUfNzc0rU1JShC1tvyVisViQmprqP3fu3LjW\nLA+txLIsJkwN0/jx409+8803S1iWpUWLFsVs27bt9Y8++ugT5bKIiIgtLMtSaGjoobfeemtTdXW1\n6d27d+18fX3TYmNjF7AsS9u3b583ZsyYVJZlqbS01LZXr14PDx06FMowDPfbb79dbGRkVLd169Y3\nFMsaGRnVbdmyJUIul3M2bdr0lpOTk0QRk1Ao/EOxLMuyVFlZacbn8wt37NgRzjAM9+LFi8NtbW1L\nr1y54sGyLL366qt7X3311b3V1dWmf//99xAej1fk7+//p6r9vXz58mBzc/OK1NTUMbW1tcbvvvvu\nV4aGhvUnTpyYwLIsiUSi1XPmzIljWZb++9//Rr788suHa2pqesjlcs6FCxe8Hj16ZKEqRpZlicPh\nyAMDA4/fv3/f8vHjxyaKshs3briwLEvh4eE7LCwsHim2/c4772xQHLObN28KOByOnGEYrqrjsGPH\njnDFsopp3rx521euXLmGZVk6ceLEBFtb29KLFy8Or62tNX777be/Gzt27Cnl2F5++eXDDx8+7HX7\n9m1nOzu7u8ePHw9kWZbCwsL2rFu3bgXLslRbW2t85syZF5r7e1HeH4ZhuM8//3zm2rVrP66vrzcs\nKCjo7+LickOx3szMzOfT0tJ8GYbhisXifh4eHlc2bNjwjqp1Nf4bUrVMeHj4jt69ez/466+/RrEs\nS9XV1aYtbb+lKTo6etX48eNPavv/X1ebcAYDTxk3btypP//8cywR0enTp8eMHTv2T39//1RFWWpq\nqv+4ceNO3blzx/7YsWPB33zzzb9MTU1r7OzsSpcsWbJh7969YY3XefTo0SnPPffc36GhofFcLle+\nePHi7xwcHKTKy/Tr1+9WRETEVg6Hw86dOzeupKTE8e7du30U9azSpZNff/31pf79+98MDw/fyeVy\n5cOHD8+aPn36wf37989kGMbg4MGD09esWbPK1NS0ZsiQIZfDw8N3ss1cevnll19eefnll4+MGTPm\ntLGxcd3atWtXcrlcufJ2FW2NjY3rysvLbfLy8tw4HA7r5eV10cLCokJVjAorVqz4zNLS8oGJiUmt\nqu2/9NJLvyq2/emnn3509uzZURKJhNf8v1Dz21K2e/fu2REREVuHDx+eZWxsXPfZZ5+tOHv27CjF\nGRIR0fLlyz/v1avXI2dn58Lx48f/oThLNTY2rhOLxQKJRMIzNjaue+GFF/56VjxEROfPn/cpKyuz\n/fjjjz8xNDSU9e/f/+b8+fO3KP4mnn/++Qu+vr7pXC5X3q9fv1sLFiz44dSpU+Nas+7mhIaGxo8a\nNeosEVFOTs6wlrbfkri4uLnz5s3b0ZFYoCmV14Wh+xo7duyf33///aL79+9blZaW2g0YMOCGnZ1d\naXh4+M779+9bXb58ecjYsWP/vHXrVr/6+nojR0fHEkVbuVzO7du37+3G6ywuLnbi8/lFymWN55UT\nTs+ePauJiCorK8379Olzl+jpPo5bt271S0tL87OysrqvKJPJZIZz586NKysrs5XJZIbOzs6FijpV\nMSmUlJQ4KsfSs2fPahsbm3JVy86ZM+fHwsJC57CwsL0PHjywfO2113Z9+umnHxkaGsoax6igHEdj\nHA6HVd62mZlZlbW19b3i4mInOzu70ubatUZJSYmjt7d3hvK6bWxsyiUSCU9xPBof84qKCgsioi++\n+OKDlStXrvX19U23srK6/9577331+uuvb3/WNm/dutWvuLjYSfnfhWEYA8Xlr+vXrw989913v87M\nzBxRXV3dUyaTGSrH2FYcDofl8XiS1m6/OadPnx5z584d+1deeeWX9sYCqiHBwFNGjhx57uHDh703\nb9785ujRo88QEfXq1euRk5NT8Q8//LDAycmpuF+/freMjY3rTExMasvLy22Uf/Gr4uTkVHzkyJGX\nFfMsy3KKior4rY2p8Rd33759b48bN+5UUlJSYONlGYYxMDQ0lN2+fbuvot9D+Vd7Y46OjiVXr171\nUMxXV1f3LC8vt1G1rKGhoWzVqlVrVq1atebWrVv9pkyZcnTQoEHX3njjjW3NdfK31PnPsiynsLDQ\nWTFfWVlpfu/ePWsnJ6diY2PjOkU85ubmlUREUqnUoTXrJXpyzBW3QxMRVVVVmZWXl9sofyE3x97e\n/s4PP/ywgIjozJkzowMCApLHjRt3ysXFpaCldn379r3dv3//m9evXx+oqn7hwoWbRowYkblv375X\nzczMqjZs2LDkwIEDM5pbn5mZWVV1dXVPxbzy/isoH4dnbb85O3fuDJ8xY8YBxQ8b6Dy4RAZPMTU1\nrfH29s74+uuv31X+5TdmzJjTX3/99bvjxo07RfTkizkwMDDp3Xff/bqiosJCLpdzb9y4MUBxKU3Z\nlClTjl66dGloQkLCVJlMZvj9998vUvVl0Rx7e/s7yh31L7300q/Xr18fuGvXrtfq6+uN6uvrjc6f\nP++Tm5vrbmBgwEyfPv2gSCQS1dTUmF65cmXwzp07w5v7Qn7llVd++fXXX186c+bM6Lq6OuNVq1at\nae5upJSUFOGlS5eGMgxjYGFhUWFkZFRvYGDAqIqxtY4ePTpFse2VK1euHTVq1Fkejyexs7Mr5fF4\nkh9//HEOwzAG27Zte0N5/fb29neKior4ik5yoqcv582aNWvP9u3bX8/Ozvasra01+fDDD9eNHDny\nXHNnc8qX3Pbv3z9T8QPA0tLyAYfDYZ/1I4KIyNfXN93CwqLiiy+++KCmpsaUYRiDv//++7mMjAxv\noicJ1MLCoqJnz57Vubm57ps2bVqo3L7xMfT09My+fPnykOzsbM/Hjx/3EIlEouZibs32VampqTHd\nv3//TFweUw8kGGhi3Lhxp0pLS+3GjBlzWlHm7++fWlZWZqucdOLi4ubW1dUZDx48+Iq1tfW9mTNn\n7lckDuVnPGxtbcv2798/84MPPvjC1ta27OrVqx7e3t4Zin4JVc+DKM+/88473/7yyy+vWFtb31uy\nZMkGc3PzyqSkpMC9e/eG8Xg8iaOjY8mKFSs+q6urMyYiiomJiaqsrDR3cHCQvvHGG9veeOONbc3t\n6+DBg698//33i/75z3/+5OTkVGxtbX1P+bKWcmxSqdRh5syZ+3v37v1w8ODBV4RCYYriDrPGMTa3\nPeX94nA47OzZs3dHR0evtrGxKb948aLXrl27XlPUb968+c1///vfS21tbcuuXLkyWHFGSUQ0ceLE\nE0OGDLns4OAgVb6MqFj/xIkTT6xdu3bljBkzDjg5ORXfvHmzv3JfhKrjrSjLyMjwHjly5DkLC4uK\nqVOnJnz33XeLFc/utLQ/XC5X/uuvv76UlZU13MXFpcDOzq50wYIFPzx69KgXEdGXX375/k8//fTP\nXr16PVqwYMEPYWFhe5Xbi0QiUXh4+E4rK6v7v/zyyysDBw68vmrVqjUBAQHJgwYNuubv75/a+Pi1\nZfuqxMfHh1pZWd0XCoUpzS0D7cdhWZ1/xAC6GLlcznV2di786aef/qk4IwKArkcnz2ASExMnu7u7\n57q5ueWtX79+WeP6hISEqZ6entleXl4XR4wYkakY1oPoyTAfw4YNy/Hy8rro6+ubrtnIoTlJSUmB\nDx48sKytrTVRPGsycuTIc9qOCwDUSNv3STeeZDKZwYABA/Jv3rwpqKurM/L09MxSPN+gmCorK80U\nn3NycoYOGDAgXzEvEAhulpeXW2t7PzA9PYlEotU2NjZlFhYWj0aOHHk2PT3dR9sxYcKESb2Tzp3B\npKen+7q6uuYLBAKxkZFRfVhY2N6EhISpysuYmZlVKT5XVlaa29ralinXs3o8/EZXtXr16uiysjLb\nR48e9Tp79uwoHx+f89qOCQDUS+duU5ZIJDzlTlY+n1+Ulpbm13i5+Pj40BUrVnxWUlLiqHy7KofD\nYQMCApINDAyYyMjI2MYj4erDuFYAALqorT/edS7BtDYBhIaGxoeGhsanpqb6z5kz58dr164NInpy\n376jo2NJaWmp3aRJk353d3fPbTyeEcsix3QWkUhEje4ehQ7Q9eO5fvlyqpFKm5SbOjjQss8/10JE\nzdP1Y6lvOJy2XxjSuQTD4/Ekyg+fFRYWOjd+6luZv79/qkwmMywvL7exsbEpVzxZbmdnVzpt2rRD\n6enpvup81whAV7N8+XqSSmtU1uVn/kWnZ0xoUi4Si9UcFegjnUsw3t7eGXl5eW5isVjg5ORUvG/f\nvlf37NkzS3mZGzduDHBxcSngcDjshQsXnicisrGxKa+uru6peAiuqqrKLCkpKXD16tXR2tkTAP0k\nldaQQCBSWZd1+rTKcgBVdC7BGBoaymJiYqKCgoKOMwxjEBERsdXDw+NqbGxsJBFRZGRk7IEDB2bE\nxcXNNTIyqjc3N69UPEAmlUodpk+ffpDoydhUs2fP3h0YGJikzf3p6oRCobZD6FJwPDsPjqX2dbsH\nLTkcDtvd9hmgLebNEzV7BvOfb71oSr9+TcrzmXt0+pLG3kEHWsDhcPS/kx8AdJeMMSWB5ZIm5VlF\nGn1TM+gJJBgAeMrfmZkkzopXWVdX+1jD0YA+Q4IBgKfIampIyLdUWXdejsvL0HpIMADwlOKqfIqX\nqn6mpc64WMPRgD5DggHQQ8098NhRpg4OJDOqI8uRql/Xw96Vdfo2oetCggHQQzVSKYkEgk5fb2sf\nmMwtO00nbm4hI64JzfD4mPKKztOgQYPIxsaG4uLiyNXVtdNjA/2jc4NdAoDuyyw+QuMFb9CCEbH0\na943ZNvbma5du0aLFi2iyMhIbYcHOgIJBgDa5Hj+93T70SVKLoilndnvkrQyjyzN7YmIKCwsjC5c\nuEDl5eVajhJ0ARIMALRJkOsicjIfRMGub1Ogy1tkYWzbMBCigYEBOTk5UWFh4TPWAt0BEgwAtAtL\nRCzhtmVoHhIMALQLh4h6m/ShirqyhldgMAxDxcXF5Ozs3HJj6BaQYACg3cyMrcjBzJUeVN4hIqI9\ne/bQ888/TzY2NlqODHQBblMG0EOmDg5qeQeLqYPq519a8uLAf9H2nLdp0KBBZG1tTXFxcZ0eF+gn\njKYMAE+x7utEwyYEqqw7ffgIfTzsQJPy+KJPKCs/Wd2hgRZhNGUA0Iry8nISzZunsk4XX6cMmoEE\nAwAdxmGYZkcWwOuUuy908gMAgFrgDAYAWk3OVqscaVnKFNG8+CyVbfKZeyRSc1ygm5BgAKD1jOUq\nR1qW3+WofNMlEd522Z0hwQDoofXrl1NNTfuH68/MyqSaupom5XUyU6qWPepIaAANkGAA9FBNjZRE\nIkG728cnZpGlO79J+YbPHxB7Wf7M9mW3b9PNCxeIa2BAPczNiXn4gKJPTaD/895Gdmbtjwu6FiQY\nAGiz4uvXSeDlRXb9+tHDu3fp3sNH1Ku2h7bDAh2DBAMAbZKfnk6P7tyhmocPqTg3lzyDgojDxQ2p\n0JRO/lUkJiZOdnd3z3Vzc8tbv379ssb1CQkJUz09PbO9vLwujhgxIvPkyZMTWtsWADrG1deXzG1t\nydXPjzyDgojoyd1l1cxDOlG2heKlnz81FVflazli0BadO4NhGMYgKioqJjk5OYDH40l8fHzOh4SE\nHPbw8LiqWCYgICB56tSpCUREly5dGjpt2rRD+fn5rq1pCwBqYCwnrsyYLDxtyczS8qmq2ydztBQU\naJvOJZj09HRfV1fXfIFAICYiCgsL25uQkDBVOUmYmZlVKT5XVlaa29ralrW2LRGRSCRq+CwUCkko\nFKpxjwAA9E9KSgqlpKR0aB06l2AkEgnP2dm54XV4fD6/KC0tza/xcvHx8aErVqz4rKSkxDEpKSmw\nLW2VEwxAd1RcfJek9KBJ+e3bj4mRP/suMuj6Gv/4jo6ObvM6dC7BcDicVg11HBoaGh8aGhqfmprq\nP2fOnB9zc3Pd1R0bgK4wNXUgkUjc7vYHf31MPfl9mpQzBubEshltWld+WhqxZXVUSzLKSUoiIxMT\n8p46td2xQdehcwmGx+NJCgsLG16HV1hY6Mzn84uaW97f3z9VJpMZ3rt3z5rP5xe1pS2Avlq2rGOj\nE+/6dThZDw9tUt6DiOj8sxPM8P/fuU9E5OrnR8X3Cmhs8KwOxQRdj87dRebt7Z2Rl5fnJhaLBXV1\ndcb79u17NSQk5LDyMjdu3BigeC/BhQsXnicisrGxKW9NWwAA0AydO4MxNDSUxcTERAUFBR1nGMYg\nIiJiq4eHx9XY2NhIIqLIyMjYAwcOzIiLi5trZGRUb25uXrl3796wltpqd48AALonvNESoBuy7e9M\n5n6jVdbd/uNQs5e7/jy2R2Vdc+VERDknk+je7eL2Bws6AW+0BIBWkTMMCXqoHtrltoZjga5L5/pg\nAACga0CCAQAAtcAlMgA9tHz5epJKm77PpbWqyi0p93TTO/iNTes6EhbAU5BgAPSQVFpDAoGo3e25\nhhbUw3xAk/LHlTGtaq/8PhiPsWOJrZBR2oED9LiykrynTm0yHllnWb16NQ0ZMoT+8Y9/tLicSCSi\nqqoq+ve//92k7tSpU1RXV0eTJk1q8/ZTUlJo6dKldP78eZX133//PVVWVtKyZcsoOzubrl+/TjNn\nzmzzdp5l6dKl5Ovrq5Z1dyYkGABoM+X3wRARkQmXhgdMpqzERLVut7XDlXA4zd/s9Mcff1BVVVW7\nEkxL6urq6JtvvqHs7GwiIrp48SL99ttvzSYBhmHIwMCgXdt69913KTAwUOcTDPpgAKBNFO+DKcjI\noOzjx4mIiGPMJRMzs2e25fP5VFZWRkREU6ZMoZdeeomIiO7evUvOzk8G4airq6OlS5eSn58fDR8+\nnObOnUtVVU/Gt503bx59//33RET08OFDmjFjBnl4eFBAQADNnTuXli5d2rAtiURCL774Inl4eNBL\nL71ENTU1dOnSJYqNjaW4uDjy8vKiL774goiIjh49SmPGjCFvb2964YUXKC0trWE9H3/8Mbm5uZGv\nry/99ttvze7bkSNHyNfXl8zMzKi8vJxWr15NycnJ5OXlRUuWLCEiIi6XS9HR0eTr60vR0dEUHR39\nVMwikahhvqXj4OjoSLa2tvTXX38985hrExIMALSJqvfBtNb48ePpxIkTJJPJ6ObNmyQWi0kmk9GJ\nEydowoQnr3X64osvyNLSktLS0igrK4scHR3ps88+I6InZyaKs5M1a9aQjY0NXb16lfbv30+nT59u\nqGNZljIyMmjPnj109epVqq+vp927d9PQoUPprbfeovDwcLp48SJ98MEHdOPGDfrkk0/o2LFjlJGR\nQZs3b264BHfkyBE6cuQIZWdn07lz5+jatWvNnh2dPHmSRo0aRURENjY2tGbNGgoICKCLFy/Shg0b\nGpbr2bMnpaen05o1a5qsQ3n/WjoOREQvvPACnThxok3HX9NwiQwANGbixImUnJxMPB6v4cs4LS2N\nkpOTGxLM4cOHqaKign755RciIqqtraXhw4c3WVdKSgrFxDzpM7KysqLQ0KfHVps8eTL16tWLiIj8\n/Pzoxo0bRPQk+Sg/bH38+HG6ceMGjR07tqGMYRi6e/cu/fHHHxQWFkY9e/YkIqKIiAj65JNPVO6b\nWCymiRMnNsw390B3eHh4c4fnKc86Dnw+n9LT01u1Lm1BggEAjZkwYQKtWbOG+Hw+TZw4kViWpeTk\nZDp58uRT/SubNm1q1XualL/ElT9zOBwyMTFpmDcwMKCamubvups8eTLt3LmzSfn/f3pd5TaeFU9z\nzM3NGz4bGhqSXOn1CI1jfNZx0PVRSXCJDAA0pm/fvsTlcmnnzp0UEBBAEydOpB07dpCxsTHx+Xwi\nIgoJCaGvvvqKHj9+TEREFRUVlJub27AOxZeqUCikuLg4IiJ68OABHT58uMkyquZ79+5NDx8+bJgP\nDAykxMREunLlSkOZ4i6xCRMm0M8//0zV1dXEMAxt37692X0TCAQkkUia3Y4qrq6ulJmZSSzLUkVF\nBf36668Ndc86DkVFReTi4tLi+rUNZzAAesjBwZTEYlG728tlP9PjSqcm5camdUQVbV8fWyGjc/v3\nU93jx03eCVNdW0vz5v0vVgMDSyovL6RlyzaRg4Mp9ezZ86nLU8uXLyeRSEQ+Pj7E5XKJw+GQSCQi\nd/cnr3xS9FGsWrWKXn/9dfLw8CBHR0fy9vam3r17Nyyj3FeiPD9t2rSGTv5Zs2bRBx98QLt27aKI\niAiqqamhuro6GjNmDPn4+NCLL75IZ8+eJU9PT7KysiKhUEglJSUqj8H48ePp4MGDtHjxYiJ6cjnw\nyy+/pOHDh5NQKKQNGzY06b+ZPn067du3jzw8PKhv377k4+PT6uNw9uxZWrt2bdv/sTQIg10CdEPW\nfZ1o2IRAlXUtDVzZnsEuz/52lFYsuquyTiwW0Y4dotYF3YhMJiOGYcjExIQePXpE/v7+9M033zT0\n5WhabW0tDRkyhLKzs8msFXfUdURJSQkFBATQ5cuX1bodZe0Z7BKXyABAL927d49Gjx5NXl5e5Ofn\nRzNnztRaciEiMjExoX/9618Nt1Gr09dff60Xr37HJTIAUCsZU0nxKfNU1jFV+UQkatd6+/TpQxkZ\nbXu9s7otWrRII9tRNUKBLkKCAQD1MpaTpVCgsqro1yzNxgIahUtkAACgFkgwAACgFrhEBqCHlouW\nk/SBtN3tq5gHlPv36SblxhzTjoQF8BQkGAA9JH0gJUGooN3tuTcNqEd/8yblj09Xtqo9hutX/3D9\nrd3Xtpg3bx75+PjQokWLNDLkPxIMALQZhutvqrOH62/tvraF8gOnmhjyH30wANAmHRmu/+zXX1N9\ndTUREV3avZtK8guIqHsO179mzRr6+++/aezYsTRixAgaMmQIffvttw3rU95XkUhEs2bNarI/zzpe\nEomEJk6cSEOGDKEXX3yRysrKGobN0cSQ/zqZYBITEye7u7vnurm55a1fv35Z4/rdu3fP9vT0zB42\nbFjO6NGjz+Tk5AxT1AkEAvGwYcNyvLy8Lvr6+ur2UKMAeqitw/WzLEvi3FwS5+aSkbU1XT99mm5e\nuUKP7tyh2uqabjtcf3R0NAkEAkpOTqbMzExKS0ujH374ga5du9ZkX4mIMjMzm+zPs47X4sWLSSgU\n0uXLlykmJoZOnTr11DrVPeS/zl0iYxjGICoqKiY5OTmAx+NJfHx8zoeEhBz28PC4qljGxcWl4M8/\n/xzbu3fvh4mJiZMXLFjww7lz50YSPRkKJiUlRWhtbX1Pe3sBAMoEPXoQEVEPPp8e3r1L1r17U529\nPd0rkXTr4fqrqqrorbfeopycHOJyuVRcXEzZ2dk0aNCgJu2a25+WjpfyMerfv/9T8RGpf8h/nUsw\n6enpvq6urvkCgUBMRBQWFrY3ISFhqnKCGTVq1FnFZz8/v7SioiK+8jraOl4OAGiGpYMD3crOJpOe\nPcnS0ZEe3Svr1sP1f/jhh+Tk5ERxcXHE5XIpKCioYfTkxnE03h/l5Vo6Xs+KX51jM+rcJTKJRMJz\ndnYuVMzz+fwiiUTCa275rVu3RkyZMuWoYp7D4bABAQHJ3t7eGZs3b35TVRuRSNQwpaSkdGr8ANC8\nHubmRBzXDS0qAAAa0UlEQVQO3blxg6wcHcmoh3G3Hq7/4cOHxOfzicvl0t9//02pqakq41a1P4qy\nlo7XhAkTGmK+efMmnTx58qn1tDTkf0pKylPfle2hc2cwHA6n1en0jz/+GL9t27Y3zpw5M1pRdubM\nmdGOjo4lpaWldpMmTfrd3d0919/f/6l/NX0YJA6gJQ6WDiSOF7e7vTyDoceSprckt/c5mJaG62/M\nytGRHt29S8ampsQ1MOjWw/V//PHHNGfOHNq6dSsNHDiQxo0b91S9YvmW9qel4/Xtt9/S3Llz6aef\nfqL+/fvT+PHjn1p/S0P+C4XCp86K2nNXm84N13/u3LmRIpFIlJiYOJmI6LPPPlvB5XLly5YtW6+8\nXE5OzrDp06cfTExMnOzq6pqval3R0dGrzc3NK997772vFGUYrh9As8P1t1SXczKJ7t0ubmXUT+vO\nw/V3hrYO+d8lhuv39vbOyMvLcxOLxYK6ujrjffv2vRoSEnJYeZnbt2/3nT59+sFdu3a9ppxcqqur\ne1ZUVFgQEVVVVZklJSUFDh069JKm9wEA1K87D9ffGTQx5L/OXSIzNDSUxcTERAUFBR1nGMYgIiJi\nq4eHx9XY2NhIIqLIyMjYNWvWrLp//77VwoULNxERGRkZ1aenp/tKpVKH6dOnHyQikslkhrNnz94d\nGBiYpM39AQD16M7D9XcGTQz5r3MJhogoODj4WHBw8DHlssjIyFjF5y1btszfsmXL/MbtXFxcCrKy\nsprezwgAABqnc5fIAACga0CCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAA\ntUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtUCCAQAAtdDJ4foBoHOsX76caqTSJuV1tbVa\niAa6GyQYgC6sRiolkUDQpPzfcrnmg4FuB5fIAABALXAGA9CFJV/OJLE4q0l5HadGC9FAd4MEA9CF\nVbI1JBDym5SzuawWooHuBpfIAABALZBgAABALZBgAABALdAHA9CFlZRXUHxu0+dg5OiCAQ3QyTOY\nxMTEye7u7rlubm5569evX9a4fvfu3bM9PT2zhw0bljN69OgzOTk5w1rbFqA7qWcMyLLHyCYTgCbo\nXIJhGMYgKioqJjExcfKVK1cG79mzZ9bVq1c9lJdxcXEp+PPPP8fm5OQMW7ly5doFCxb80Nq2AACg\nGTqXYNLT031dXV3zBQKB2MjIqD4sLGxvQkLCVOVlRo0adbZ3794PiYj8/PzSioqK+K1tCwAAmqFz\nfTASiYTn7OxcqJjn8/lFaWlpfs0tv3Xr1ogpU6YcbUtbkUjU8FkoFJJQKOyk6AEAuoaUlBRKSUnp\n0Dp0LsFwOJxWdz/+8ccf47dt2/bGmTNnRrelrXKCAQCAphr/+I6Ojm7zOnQuwfB4PElhYaGzYr6w\nsNCZz+cXNV4uJydn2Jtvvrk5MTFxspWV1f22tAUAAPXTuT4Yb2/vjLy8PDexWCyoq6sz3rdv36sh\nISGHlZe5fft23+nTpx/ctWvXa66urvltaQsAAJqhc2cwhoaGspiYmKigoKDjDMMYREREbPXw8Lga\nGxsbSUQUGRkZu2bNmlX379+3Wrhw4SYiIiMjo/r09HTf5tpqd48AALonnUswRETBwcHHgoODjymX\nRUZGxio+b9myZf6WLVvmt7YtAABonk4mGADoHqofP6J5S+aprHOwdKDPRZ9rNiDoVEgwAKA1rJGc\nBKEClXXieLFGY4HOp/ZO/sePH/eora01Ufd2AABAt3R6gpHL5dyDBw9Onzlz5n4ejyfp37//zX79\n+t3i8XiSV1555ZdDhw5NY1mW09nbBQAA3dLpCUYoFKZkZmaOeP/9978sKChwKSkpcZRKpQ4FBQUu\n77///pfnz5/3GTdu3KnO3i4AAOiWTu+D+f333yeZmJjUNi43MTGpHTly5LmRI0eewyUzAICur9MT\njCK5JCcnBwQEBCQr1+3cuTM8PDx8p6oEBADdDyNnKCU+XmVdZSaj4Wigs6mtkz86Onr1woULN1VV\nVZlJpVKHl19++cjhw4dD1LU9ANBDLJHQ0lLlJKup0XZ00EFqSzCnTp0a5+LiUuDp6Znt7++fOmvW\nrD0HDhyYoa7tAQCAblHbczD379+3On/+vM+AAQNuFBUV8W/fvt2XZVlOW0ZLBoCuTc7KVb7SmejJ\n655Bv6ntDGbUqFFng4KCjh8/fjzo/PnzPhKJhDd69Ogz6toeAOgjrspXOlv2GEn1jIG2g4MOUtsZ\nzO+//z6pX79+t4iIevbsWb1x48a3T506NU5d2wMAAN3S6WcwN27cGEBEpEguyhTPvyiWAQCArqvT\nz2A+/PDDdVVVVWYhISGHvb29MxwdHUvkcjlXKpU6ZGRkeB8+fDjEwsKiYu/evWGdvW0AANAdnZ5g\n9u3b92p+fr7r3r17wz766KNPb9261Y/oyRnNmDFjTm/cuPFtFxeXgs7eLgAA6Ba19MG4urrmv/fe\ne1+ZmprWpKam+nO5XPmYMWNOL1y4cJOpqSlubgcA6AbU1sk/d+7cuF69ej165513vmVZlvPTTz/9\nc+7cuXH79++fqa5tAgCA7lBbgrl8+fKQK1euDFbMT5gw4eTgwYOvqGt7AACgW9T2HMzzzz9/4ezZ\ns6MU8+fOnRs5YsSITHVtDwAAdIvazmAyMjK8R48efcbZ2bmQw+Gwt2/f7jto0KBrQ4cOvcThcNic\nnJxh6to2AABon9oSTGJi4mR1rRsAAHSf2hKMQCAQq2vdAACg+9TWBwMAAN2bTiaYxMTEye7u7rlu\nbm5569evX9a4Pjc3133UqFFne/To8firr756T7lOIBCIhw0bluPl5XXR19c3XXNRAwCAMrVdImsv\nhmEMoqKiYpKTkwN4PJ7Ex8fnfEhIyGEPD4+rimVsbGzKN27c+HZ8fHxo4/YcDodNSUkRWltb39Ns\n5AAAoEznEkx6erqvq6trvqIPJywsbG9CQsJU5QRjZ2dXamdnV/rbb7+9qGodLMtyNBQugNYtX76e\npFLVA2TU1NZrOBqA/9G5BCORSHjOzs6Fink+n1+Ulpbm19r2HA6HDQgISDYwMGAiIyNj33zzzc2N\nlxGJRA2fhUIhCYXCDkYNoD2/nThCBmauKuvq2CoNRwNdRUpKCqWkpHRoHTqXYDr6xsszZ86MdnR0\nLCktLbWbNGnS7+7u7rn+/v6pyssoJxgAfVfDVhJfKFBdmSvXaCzQdTT+8R0dHd3mdehcJz+Px5MU\nFhY6K+YLCwud+Xx+UWvbOzo6lhA9uYw2bdq0Q+np6b7qiBMAAFqmc2cw3t7eGXl5eW5isVjg5ORU\nvG/fvlf37NkzS9WyjftaqqurezIMY2BhYVFRVVVllpSUFLh69eq2p10APfKgvJxkubkq61i2QxcE\nADpE5xKMoaGhLCYmJiooKOg4wzAGERERWz08PK7GxsZGEhFFRkbGSqVSBx8fn/OPHj3qxeVy5d9+\n++07V65cGXz37t0+06dPP0hEJJPJDGfPnr07MDAwSbt7BKBecoYhQY8eKutuazgWAGU6l2CIiIKD\ng48FBwcfUy6LjIyMVXx2cHCQKl9GUzA3N6/MysoarokYAQCgZTrXBwMAAF0DEgwAAKgFEgwAAKgF\nEgwAAKgFEgwAAKgFEgwAAKgFEgwAAKiFTj4HAwDawxKRuLKsTXUYLwBUQYIB6IZaSiJERAK+6q+G\nW83UYcQAUAUJBqAb4lDLSQSgMyDBAHRhmrqk9azLatA9IcEAdGHNnamo4yyluTMiXD7rvpBgAPRc\nbX2tTp894Oym+0KCAdB3HLnO96e05+ym+vEjmrdkXpNyB0sH+lz0eecEBmqFBAMAOok1kpMgVNCk\nXBwv1ngs0D540BIAANQCZzAAoJNkcobi41OalDOZlZoPBtoFCQYAtKalGwDkcqIHUocm5ZVFZ9Qc\nFXQWJBgAPbB+/XKqqZGqrJPLGQ1H07laugFA0KNHk/IcRr/3tztBggHQAzU1UhKJBCrr1v9Xs7EA\ntBY6+QEAQC2QYAAAQC1wiQwAdFJzNwDU1tdqPhhoF51MMImJiZOXLFmygWEYg/nz529ZtmzZeuX6\n3Nxc99dff337xYsXvT799NOP3nvvva9a2xZAlzXXmX/pUiYRCTQej7apugGglCPXQiTQHjqXYBiG\nMYiKiopJTk4O4PF4Eh8fn/MhISGHPTw8riqWsbGxKd+4cePb8fHxoW1tC6DLmuvMDw09rflgADpI\n5/pg0tPTfV1dXfMFAoHYyMioPiwsbG9CQsJU5WXs7OxKvb29M4yMjOrb2hYAADRD585gJBIJz9nZ\nuVAxz+fzi9LS0vw6s61IJGr4LBQKSSgUdihmAICuJiUlhVJSUjq0Dp1LMBwOp90jeLe2rXKCAQCA\nphr/+I6Ojm7zOnQuwfB4PElhYaGzYr6wsNCZz+cXqbstgKa09FR+d+3MbwuGYUgkmqeyztTUgZYt\nw1D+ukLnEoy3t3dGXl6em1gsFjg5ORXv27fv1T179sxStSzLspz2tgXQlpaeykdn/rNxONTs8ROJ\nxBqNBVqmcwnG0NBQFhMTExUUFHScYRiDiIiIrR4eHldjY2MjiYgiIyNjpVKpg4+Pz/lHjx714nK5\n8m+//fadK1euDDY3N69U1Vbb+wQA0B3pXIIhIgoODj4WHBx8TLksMjIyVvHZwcFBqnwp7FltAfTd\n5YJymrckXmVdPYsHD0E36WSCAYCn1RJDglBL1ZX78OAh6CYkGAA9UF4uo107r6mskz/dFQmgM5Bg\nAPSAnDEkfm/Vj4PdYHFFGHQTEgyAGuBWZO3IzMzELcw6BAkGQA1wK7J2GBjU4BZmHaJzY5EBAEDX\ngAQDAABqgQQDAABqgQQDAABqgQQDAABqgbvIADoArzjWH83dwozbl9UHCQagA/CKY/3R3C3MuH1Z\nfXCJDAAA1AIJBgAA1AKXyABAr9TLHjf76oLLBeUajgZaggQDAPrFiG321QUnTxZoOBhoCRIMwDNg\n4ErdwhKRODdXZV1dbdtfvoYBMtUHCQbgGTBwpW5hWaKsDNVfXRUVbe9WxgCZ6oMEAwB6hWW5ZGnu\nrbJOLj+p4WigJbiLDAAA1AIJBgAA1AKXyAAIHfkA6oAEA0DoyAfVcIdZx+hkgklMTJy8ZMmSDQzD\nGMyfP3/LsmXL1jdeZvHixd8dO3YsuGfPntU7duyY5+XldZGISCAQiHv16vXIwMCAMTIyqk9PT/fV\n/B4AQFeAO8w6RucSDMMwBlFRUTHJyckBPB5P4uPjcz4kJOSwh4fHVcUyR48enZKfn++al5fnlpaW\n5rdw4cJN586dG0lExOFw2JSUFKG1tfU97e0FAADoXCd/enq6r6ura75AIBAbGRnVh4WF7U1ISJiq\nvMzhw4dDwsPDdxIR+fn5pT148MDyzp079op6lmU5mo4bAACepnNnMBKJhOfs7FyomOfz+UVpaWl+\nz1pGIpHw7O3t73A4HDYgICDZwMCAiYyMjH3zzTc3N96GSCRq+CwUCkkoFKplXwAA9FVKSgqlpKR0\naB06l2A4HA7bmuWaO0s5ffr0GCcnp+LS0lK7SZMm/e7u7p7r7++fqryMcoIBAGiPrv4Cs8Y/vqOj\no9u8Dp1LMDweT1JYWOismC8sLHTm8/lFLS1TVFTE5/F4EiIiJyenYiIiOzu70mnTph1KT0/3bZxg\noHvSh1uRMy8X0eTp3zUpr2cfayEaaAleYPZsOpdgvL29M/Ly8tzEYrHAycmpeN++fa/u2bNnlvIy\nISEhh2NiYqLCwsL2njt3bqSlpeUDe3v7O9XV1T0ZhjGwsLCoqKqqMktKSgpcvXp129MudEn6cCty\nHYehkSG9mpT/frZCC9EAdIzOJRhDQ0NZTExMVFBQ0HGGYQwiIiK2enh4XI2NjY0kIoqMjIydMmXK\n0aNHj05xdXXNNzMzq9q+ffvrRERSqdRh+vTpB4mIZDKZ4ezZs3cHBgYmaXN/AAC6K51LMEREwcHB\nx4KDg48pl0VGRsYqz8fExEQ1bufi4lKQlZU1XN3xAahLZaUhxZ+WNSnHjZGgj3QywQB0V3K5kcqR\ngln2mIqlAXQbEgx0Oc115utKRz50bRhe5n+QYKDLaa4zX1c68qFrw/Ay/6NzT/IDAEDXgAQDAABq\ngQQDAABqgT4Y0Ev68FQ+aB7LspQSH6+y7n55uYajaaq73QCABAN6SR+eygfNkzFEO/7bQ2VdQYGJ\nhqNpqrvdAIAEAwBdBsPUU1aZ6q+1e49rNRwNIMEAaNjlgnKat0T1ZZx6Fl+CHWLMkuUYW5VV8iJG\nw8EAEgzorK7az1JLDAlCLVVX7pNrNhjQGV2xfwYJBnQW+lmgO+mK/TNIMAAaVl4uo107r6msk2NQ\nS+hCkGBA67rb2GFyxpD4vf1U1t3AoJaggr6+PRMJBrQOY4cBtExf356JJ/kBAEAtcAYDGtFV7wgD\n0CZdv/MMCQY0AneEgbY1N4yMLgwh0166fucZEgx0GpylgC5rbhgZXRhCpqtCgoFOg7OU/8HT+rqn\nuWFkuuoQMrpw+QwJBkANistr6fTDOyrr5BxWw9EAETU7jExXHUJGFy6fIcFAm3W351baA8+6gC7T\n1HM1SDCg0rP6Uw4cmEFERCkpYhIKBUTU/S6DqUUdxiLrLA+KysiSr3rgy+6uubObl18+0Oz/+/bQ\nyQSTmJg4ecmSJRsYhjGYP3/+lmXLlq1vvMzixYu/O3bsWHDPnj2rd+zYMc/Ly+tia9vCE61NIo0p\nJxLlBAOdoB6XzzpLaxNMbd1j4nlsVFn3qKxr9s80p6XLatHRbV+fziUYhmEMoqKiYpKTkwN4PJ7E\nx8fnfEhIyGEPD4+rimWOHj06JT8/3zUvL88tLS3Nb+HChZvOnTs3sjVtu6rmkkV2di55erqrbNPa\nJALNO3FGqvLLqZbpXl9M+ow14pBbkL/KurS9JzUcTdeicwkmPT3d19XVNV8gEIiJiMLCwvYmJCRM\nVU4Shw8fDgkPD99JROTn55f24MEDS6lU6nDz5s3+z2qrK9qTENqTLEJDT+POrg5qLokQEVVRPY0I\nCmhSXrwH/SxdQUtnN3fvPKTJ079rUn75crW6w9IfLMvq1LR///5X5s+fv1kx/+OPP74WFRW1UXmZ\nl1566ciZM2deUMxPnDgxOSMjY8Qvv/wy41ltiYjFhAkTJkxtn9r6fa5zZzAcTuvu4WTZ9o1r3t52\nAADQNjqXYHg8nqSwsNBZMV9YWOjM5/OLWlqmqKiIz+fzi+rr642e1RYAADRD50ZT9vb2zsjLy3MT\ni8WCuro643379r0aEhJyWHmZkJCQw3FxcXOJiM6dOzfS0tLygb29/Z3WtAUAAM3QuTMYQ0NDWUxM\nTFRQUNBxhmEMIiIitnp4eFyNjY2NJCKKjIyMnTJlytGjR49OcXV1zTczM6vavn376y211e4eAQB0\nU9ru1NfU9PPPP88cPHjwZS6Xy2RmZj6vXLdu3boVrq6ueYMGDco9fvx4oLZj1bdp9erVIh6PVzR8\n+PCLw4cPv3js2LHJ2o5JH6djx45NHjRoUK6rq2ve559/vkzb8ejz1K9fP/HQoUNzhg8fftHHxydd\n2/Ho2/T6669v69Onz53nnnvukqKsvLzcOiAg4Hc3N7frkyZNSrp//77ls9aj9R3R1HT16lX3a9eu\nDRQKhX8oJ5jLly8P9vT0zKqrqzO6efOmYMCAAfkMw3C1Ha8+TSKRaPVXX331rrbj0OdJJpMZDBgw\nIP/mzZuCuro6I09Pz6wrV654aDsufZ0EAsHN8vJya23Hoa/Tn3/+6X/hwgUv5QSzdOnSL9avX/8B\ny7L0+eefL1u2bNnnz1qPzvXBqIu7u3vuwIEDrzcuT0hImDpr1qw9RkZG9QKBQOzq6pqfnp7uq40Y\n9RmLu/M6RPn5LyMjo3rFM1zajkuf4W+y/fz9/VOtrKzuK5cpP38YHh6+Mz4+PvRZ6+k2CaY5xcXF\nTsp3mvH5/CKJRMLTZkz6aOPGjW97enpmR0REbH3w4IGltuPRNxKJhOfs7FyomMffYcdwOBw2ICAg\n2dvbO2Pz5s1vajueruDOnTv29vb2d4iI7O3t79y5c8f+WW10rpO/IyZNmvS7VCp1aFy+bt26D19+\n+eUjrV1Pa5/F6U6aO7affvrpRwsXLty0atWqNUREK1euXPvee+99tXXr1gjNR6m/8DfXuc6cOTPa\n0dGxpLS01G7SpEm/u7u75/r7+6dqO66ugsPhsK35m+1SCeb333+f1NY2qp6p4fF4ks6NTP+19tjO\nnz9/S1uSOTzRmue/oPUcHR1LiIjs7OxKp02bdig9Pd0XCaZj7O3t70ilUgcHBwdpSUmJY58+fe4+\nq023vESmfG02JCTk8N69e8Pq6uqMb9682T8vL8/N19c3XZvx6ZuSkhJHxedDhw5NGzp06CVtxqOP\n8AxX56muru5ZUVFhQURUVVVllpSUFIi/yY4LCQk5vHPnznAiop07d4aHhoaqfmWrMm3fraCp6eDB\ng9P4fH5hjx49auzt7aWTJ08+pqj79NNPPxwwYED+oEGDchMTE4O0Hau+TXPmzIkbOnRozrBhw7Kn\nTp0aL5VK7bUdkz5OR48eDR44cOC1AQMG5K9bt26FtuPR16mgoKC/p6dnlqenZ9aQIUP+xrFs+xQW\nFrbH0dGx2MjIqI7P5xdu27bt9fLycuuJEycmt+U2ZQ7L4tIvAAB0vm55iQwAANQPCQYAANQCCQYA\nANQCCQYAANQCCQZAy86fP+/j6emZXVtba1JVVWX23HPP/X3lypXB2o4LoKNwFxmADli5cuXax48f\n96ipqTF1dnYuXLZs2XptxwTQUUgwADqgvr7eyNvbO8PU1LTm7NmzozB0DHQFuEQGoAPKyspsq6qq\nzCorK81rampMtR0PQGfAGQyADggJCTn8z3/+86eCggKXkpISx40bN76t7ZgAOqpLDXYJoI/i4uLm\nmpiY1IaFhe2Vy+XcF1544a+UlBShUChM0XZsAB2BMxgAAFAL9MEAAIBaIMEAAIBaIMEAAIBaIMEA\nAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBa\nIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBaIMEAAIBa/D8TJqTXimw4KgAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x382a8950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "print 'Calculating ratios'\n",
    "true_dist = True\n",
    "npoints = 50\n",
    "\n",
    "score = ROOT.RooArgSet(w.var('score'))\n",
    "getRatio = singleRatio\n",
    "\n",
    "# NN trained on complete model\n",
    "F0pdf = w.function('bkghistpdf_F0_F1')\n",
    "F1pdf = w.function('sighistpdf_F0_F1')\n",
    "\n",
    "\n",
    "# Transforming f1 into f0\n",
    "data_f1 = train_data[train_data[:,-1] == 0.]\n",
    "data_f0 = train_data[train_data[:,-1] == 1.]\n",
    "xarray = data_f1[:,:-1]\n",
    "testtarget = data_f1[:,-1]\n",
    "\n",
    "outputs = predict('{0}/{1}_F0_F1.pkl'.format(dir,model_file),xarray,model_g=model_g)\n",
    "\n",
    "F1fulldist = np.array([evalDist(score,F1pdf,[xs]) for xs in outputs])\n",
    "F0fulldist = np.array([evalDist(score,F0pdf,[xs]) for xs in outputs])\n",
    "\n",
    "completeRatio = getRatio(F0fulldist,F1fulldist)\n",
    "\n",
    "if len(vars_g) == 1:\n",
    "    F1dist = np.array([evalDist(x,w.pdf('f1'),[xs]) for xs in xarray])\n",
    "    F0dist = np.array([evalDist(x,w.pdf('f0'),[xs]) for xs in xarray])\n",
    "else:\n",
    "    F1dist = np.array([evalDist(x,w.pdf('f1'),xs) for xs in xarray])\n",
    "    F0dist = np.array([evalDist(x,w.pdf('f0'),xs) for xs in xarray])\n",
    "\n",
    "trueRatio = getRatio(F1dist, F0dist)\n",
    "\n",
    "#trueIndexes = findOutliers(trueRatio)\n",
    "#completeIndexes = findOutliers(completeRatio)\n",
    "#indexes = np.logical_and(trueIndexes,completeIndexes)\n",
    "#trueRatio = trueRatio[indexes]\n",
    "#completeRatio = completeRatio[indexes]\n",
    "#data_f1_red = data_f1[indexes]\n",
    "data_f1_red = data_f1\n",
    "\n",
    "bins = 50\n",
    "feature = 7\n",
    "x_range = (-10.,10.)\n",
    "# Transfering distributions\n",
    "# Doing histogram manipulation\n",
    "fig,ax = plt.subplots()\n",
    "colors = ['b-','r-','k-']\n",
    "colors_rgb = ['blue','red','black']\n",
    "    \n",
    "hist,bins = np.histogram(data_f1[:,feature],bins=bins, range=x_range,density=True)\n",
    "\n",
    "\n",
    "hist_transfered,bins_1 = np.histogram(data_f1_red[:,feature],weights=trueRatio,bins=bins, range=x_range,density=True)\n",
    "hist_transfered_clf,bins_2 = np.histogram(data_f1_red[:,feature],bins=bins,weights=completeRatio, range=x_range,density=True)\n",
    "hist0,bins0 = np.histogram(data_f0[:,feature], bins=bins, range=x_range,density=True)\n",
    "\n",
    "widths = np.diff(bins)\n",
    "\n",
    "ax.bar(bins[:-1], hist0,widths,label='f0',alpha=0.5,color='red')\n",
    "ax.bar(bins[:-1], hist,widths,label='f1',alpha=0.5,color='yellow')\n",
    "ax.bar(bins[:-1], hist_transfered,widths,label='f1 weighted (true)',\n",
    "    alpha=0.5,color='blue')\n",
    "ax.bar(bins[:-1], hist_transfered_clf,widths,label='f1 weighted (trained)',\n",
    "        alpha=0.5,color='green')\n",
    "\n",
    "ax.legend(frameon=False,fontsize=11)\n",
    "ax.set_xlabel('x') \n",
    "ax.set_ylabel('p(x)') \n",
    "if len(vars_g) > 1:\n",
    "    ax.set_title('Weighted distributions feature {0}'.format(feature))\n",
    "else:\n",
    "    ax.set_title('Weighted distributions')\n",
    "file_plot =  makePlotName('all','transf',type='hist_v{0}'.format(feature),model_g=model_g) \n",
    "fig.savefig('{0}/plots/{1}/{2}.png'.format(dir,model_g,file_plot))\n",
    "\n",
    "#saveFig(xarray, [true_transfer, data_f0[:,0]], makePlotName('all','transf',type='hist'),title='Transfered distribution',labels=['Transfer True', 'Truth'],hist=True, print_pdf=True,dir=dir)\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Resampling data from f1 using the ratios as weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x4152af90>"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD9CAYAAABTJWtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1Yk/e9P/B3FIquCgqTBBJqPBKk6VAsFbm2ekyLUMZZ\nIx0OdB2EFlt/8vNo7dl+UNZtuK0au2vtbFe6zcWWsnNEpi2wTinCxHZdhR7UyoqUWHlKSNIOBFF5\n5v794ZIhhCclhIf367ru60q+9wOfO7V55/5+7weRIAggIiICgDnOLoCIiKYOhgIREdkwFIiIyIah\nQERENgwFIiKyYSgQEZHNiKHw2WefrVi9evU56+Th4dH2yiuv7GxpafGMiIg4GRAQUBMZGVnU2tq6\nyLrOvn37nlMoFPrAwMDqoqKiSGt7RUVFSFBQUKVCodDv2rXrgLW9q6vLLT4+/ohCodCHhYWdqa+v\nX+qYXSUiolEJgjCmqa+vb45EIjE1NDT4/eAHP3hx//79/08QBGi12tTU1FStIAj49NNPlatWrTrf\n3d3tWltbK1++fPml/v5+kSAIWLNmTXlZWVmoIAj45je/efzEiRNRgiDgtddeS9m+fXumIAjIycmJ\nj4+PzxlrTZw4ceLEaWKnMS/43nvvRT744IMfCIKAFStWVJvNZrEgCDCZTJIVK1ZUC4KAvXv3PqfV\nalOt6zzyyCOFH330UVhTU5NPYGDgRWv74cOHN2/btu031mXOnDmzVhAE9PT0uHz1q1/90tkfCidO\nnDjN1sllrEcUOTk5m7ds2XIYACwWi1gsFlsAQCwWWywWixgAmpqafMPCws5Y15HJZAaj0Sh1dXXt\nkclkBmu7VCo1Go1GKQAYjUapn59fIwC4uLj0enh4tLW0tHh6enq2WJcXiUS87JqI6DYIgiAaz/Jj\nCoXu7u67/vSnPz26f//+1MHzRCKRMBlf2oLAXJgoGRkZyMjIcHYZM8JU/yzT0vbDbO6wO08imQ+t\ndsj/0k411T/P6UYkGlceABhjKJw4ceKbISEhFUuWLPkSuHl0YDabJRKJxGwymXy8vb2/AG4eATQ2\nNvpZ1zMYDDKZTGaQSqVGg8EgG9xuXaehoeEeX1/fpt7eXpe2tjaPgUcJRHT7zOYOyOUZdufV1dlv\np9ltTKekHj58eIu16wgA1Gp1QVZWlgYAsrKyNDExMXnW9pycnM3d3d131dbWLtPr9YrQ0NByiURi\ndnd3v1pWVrZWEARRdnZ2wsaNG/MHb+vo0aObwsPDSyZ+N4mIaCxGPVK4fv363cXFxRsOHjz4lLUt\nLS1NGxcXl6vT6ZLlcnldbm5uHAAolcqquLi4XKVSWeXi4tKbmZmZYu1ayszMTElKSnqzo6NjfnR0\n9PGoqKhCAEhOTtYlJCRkKxQKvZeXV3NOTs5mR+0s3aRSqZxdwozBz3Ji8fN0PtF06KsXiUTCdKiT\naKpJSsoYsfvozTftz6OZQSQSOWagmYhmJk9PT1y5csXZZdAdWrx4MVpaJmYolqFANAMMd5ZRRUUl\n5PLh17ty5QrP7JsBbucso+EwFIhmgL/9+U/wnOs/pP264YITqqHpjKFANAPc1XENz8jkQ9p3Xv7L\n5BdD0xrvkkpERDYMBSKasj777DMEBwfD3d0dv/71r51dzqzA7iMiusVIt8aYCOO5vcaLL76I8PBw\nnD9/HgCQmpoKnU4HANi6dSu0Wq3D6pytGApEdIuRbo0xEcZze436+np8/etfBwD89re/RX5+Pi5c\nuDl4HhERgWXLlmHbtm2OKHPWYvcREU1JDz/8MEpLS7Fjxw4sXLgQP/3pT/H9738fvr6+8PX1xfe/\n/328+eabzi5zxmEoENGU9Je//AXr1q3Da6+9hvb2dly/fh2rVq2yzV+5ciU+/fRTJ1Y4MzEUiGha\nuHbtGjw8PGzv3d3dce3aNSdWNDMxFIhoWliwYAGuXr1qe9/W1oYFCxY4saKZiaFARNPCfffdZzsL\nCQA++eQTfO1rX3NiRTMTQ4GIpoXExES89NJLaGpqgtFoxEsvvYSkpCRnlzXj8JRUolmqoqLCbrtE\nMt+hT2WTSObf1nrbtm3D5cuXERQUBAB46qmn8PTTT09kaQSGAtGs1dEx1277VHpu86lTp255v3//\nfuzfv99J1cwO7D4iIiIbhgIREdmw+4hoBuvsbMX5vCS78/qbeeEXDcVQIJom9qelocNstjuvvdkE\nyIa23y304ZlFcrvr8FkLZA9DgWia6DCbkTHMszWL+nomtxiasTimQERENgwFIiKyGTUUWltbF23a\ntOnovffee1GpVFaVlZWtbWlp8YyIiDgZEBBQExkZWdTa2rrIuvy+ffueUygU+sDAwOqioqJIa3tF\nRUVIUFBQpUKh0O/ateuAtb2rq8stPj7+iEKh0IeFhZ2pr69fOvG7SUREYzFqKOzatetAdHT08YsX\nL9574cKFlYGBgdVarTYtIiLiZE1NTUB4eHiJVqtNA4CqqirlkSNH4quqqpSFhYVRKSkpmYIgiABg\n+/btr+t0umS9Xq/Q6/WKwsLCKADQ6XTJXl5ezXq9XrF79+6XU1NTeWUK0STo7Gx1dgmjGvg4zldf\nfdXZ5cwKIw40t7W1eXzwwQfrsrKyNADg4uLS6+Hh0VZQUKA+ffr0egDQaDRZKpWqVKvVpuXn52/c\nsmXLYVdX1x65XF7n7+9/qaysbO3SpUvr29vbF4aGhpYDQGJi4lt5eXkxUVFRhQUFBeo9e/b8BABi\nY2OP7dixw+6DWDMyMmyvVSoVVCrVhHwARLPV3UKf3fa0jDSYW+2f5TQRJIsk0GaM7TGaAx/HWVpa\nioceegjnzp3D4sWLUVtb67Aap6vS0lKUlpbe0TZGDIXa2tplS5Ys+fKJJ55445NPPlkVEhJS8atf\n/eoZi8UiFovFFgAQi8UWi8UiBoCmpibfsLCwM9b1ZTKZwWg0Sl1dXXtkMpnB2i6VSo1Go1EKAEaj\nUern59cI/Ct0WlpaPD09PVsG1jIwFIjIccytZshj5A7bfl1e3ZiXHfg4zrvvvhtbt27FjRs3sHfv\nXgdVN70N/sG8Z8+ecW9jxO6j3t5el7Nnz96fkpKSefbs2fvvvvvu69auIiuRSCSIRCJh3H+ZiGgE\nAx/H6e7ujsWLF+Pxxx/HsmXLnF3ajDZiKMhkMoNMJjOsWbPmYwDYtGnT0bNnz94vkUjMZrNZAgAm\nk8nH29v7C+DmEUBjY6OfdX2DwSCTyWQGqVRqNBgMssHt1nUaGhruAW6GUFtbm8fgowQimn0GPo7z\n6tWr8Pf3d3ZJs8KIoSCRSMx+fn6NNTU1AQBQXFy84b777vv00Ucf/ZN1nCErK0sTExOTBwBqtbog\nJydnc3d39121tbXL9Hq9IjQ0tFwikZjd3d2vlpWVrRUEQZSdnZ2wcePGfOs61m0dPXp0U3h4eIlj\nd5mIphNBYEfEZBr1iuZXX331Px9//PH/7u7uvmv58uWfv/HGG0/09fXNjYuLy9XpdMlyubwuNzc3\nDgCUSmVVXFxcrlKprHJxcenNzMxMsXYtZWZmpiQlJb3Z0dExPzo6+nhUVFQhACQnJ+sSEhKyFQqF\n3svLqzknJ2ezY3eZiKYTkUjk7BJmlVFDYdWqVZ98/PHHawa3FxcXb7C3fHp6+t709PQho0AhISEV\nlZWVQYPb3dzcuqyhQkTDq6j4O/LO19md19nVPbnF0IzFex8RTRMdHb1YJFPZnSf0fzy5xUwia/eR\nIAjo6upCT0+P7bVIJMJdd93l5ApnFoYCEd1CskgyrtNGb2f742HtPjp9+jQefvhhW9v8+fOhUqnw\nl7/wbq8TiaFARLcY64Vlk2Hg4zhVKhX6+/udWM3swBviERGRDUOBiIhsGApERGTDUCAiIhuGAhER\n2TAUiIjIhqFAREQ2DAUiIrJhKBDRjLZ9+3b8/Oc/H9OySUlJ+NGPfuSwWkbb/vPPP48lS5bA19fX\nYTWMhlc0E9Et9qelocPsuMdxzpdIkKqdvKumX3/99TEvKxKJRrwr65w5c3Dp0iX827/9223VMtL2\nGxoa8NJLL6GxsRFeXl4AgKeffhrvv/8+9Ho9Dh06BI1Gc1t/dzwYCkR0iw6zGRlyucO2n1FX57Bt\nT4Y7fb7DcOs3NDTAy8vLFggAEBwcjM2bNyM1NXXSbiHO7iMimpLeeOMNqNVq23uFQoG4uH/dZd/P\nzw8XLlwAAFRXVyMiIgJeXl4IDAzEH//4R9tyg7tsXnzxRfj6+kImk+H3v/895syZg8uXL9vmt7S0\n4Fvf+hbc3d0RFhZmm/fv//7vAIBVq1Zh4cKFtr/x7rvvIjg4GIsXL8Y3vvENVFZW2rZ17tw53H//\n/XB3d8fmzZvR2dlpd1+Li4sRGRmJpqYmLFy4EE8++SQAICUlBQ8//DDmzZt3ex/ibWAoENGUpFKp\n8MEHHwAAmpqa0NPTgzNnzgAALl++jOvXr2PlypW4fv06IiIi8L3vfQ9ffvklcnJykJKSgosXLwK4\ntcumsLAQL7/8MkpKSqDX61FaWnrL3xQEATk5OcjIyMCVK1fg7++PH/7whwCA999/HwBw4cIFtLe3\n4zvf+Q7OnTuH5ORkHDx4EC0tLdi2bRvUajV6enrQ3d2NmJgYaDQaXLlyBd/5zndw7Ngxu7/4N2zY\ngBMnTsDX1xft7e04dOiQQz7TsWAoENGUtGzZMixcuBDnzp3D+++/j0ceeQS+vr747LPPcPr0adsv\n93fffRfLli2DRqPBnDlzEBwcjG9/+9u3HC1Y5ebm4sknn8S9996L+fPnY8+ePbfMF4lE+Pa3v40H\nHngAc+fOxeOPP47z588PW+Pvfvc7bNu2DWvWrIFIJEJiYiLc3Nzw0Ucf4cyZM+jt7cWuXbswd+5c\nxMbGYs2aIc8rs5kqjx3lmAIRTVnr169HaWkpLl26hPXr12PRokU4ffo0PvroI6xfvx4AUF9fj7Ky\nMixevNi2Xm9vLxITE4dsz2QyITQ01PZeJpMNWUYsFttez58/H9euXRu2vvr6erz11lt49dVXbW09\nPT0wmUwQBAFSqfSW5ZcuXTplvvyHwyMFIpqy1q9fj1OnTuGDDz6ASqWyhcTp06dtoXDPPfdg/fr1\nuHLlim1qb2/Ha6+9NmR7Pj4+aGxstL0f+Pp23HPPPfjhD394y9++du0a4uPj4ePjA6PReMvy9fX1\nU/6Z0wwFIpqyrKHQ2dkJX19fPPjggygsLERLSwtWr14NAPjWt76Fmpoa/OEPf0BPTw96enrw8ccf\no7q6GsDNbhnrr/O4uDi88cYbqK6uxo0bN/Czn/3slr832q94sViMzz//3Pb+qaeewm9+8xuUl5dD\nEARcv34df/7zn3Ht2jV8/etfh4uLC1555RX09PTg7bffxscfj++xqT09Pejs7ER/fz+6u7vR2dnp\n8CMNdh8R0S3mSyQOPW10vmTsj+NUKBRYuHAh1q1bBwBwd3fH8uXL4e3tbfvFvWDBAhQVFeHZZ5/F\ns88+i/7+fgQHB+Oll14CcOtAc1RUFHbu3ImHHnoIc+fOxfPPP4/s7Gy4ubkNWdZq4PuMjAxoNBp0\ndHTg4MGD2LRpEw4ePIgdO3ZAr9dj/vz5WLduHdavXw9XV1e8/fbbeOqpp/D8888jOjoasbGxI+7v\n4L8dERGB999/HyKRCB999BGefvpplJaW2sZTHEE01fu3AEAkEgnToU4iR9rgH4znZTF25+3+8Bd4\n+Rs/GHO7dd753htTvo/bkS5evIigoCB0d3djzpzp23EiEons/nf8Z/u4+qtG/RTkcnndypUrL6xe\nvfpcaGhoOQC0tLR4RkREnAwICKiJjIwsam1tXWRdft++fc8pFAp9YGBgdVFRUaS1vaKiIiQoKKhS\noVDod+3adcDa3tXV5RYfH39EoVDow8LCztTX1y8dzw4QEY3HO++8g66uLly5cgWpqalQq9XTOhAm\n2qifhEgkEkpLS1Xnzp1bXV5eHgoAWq02LSIi4mRNTU1AeHh4iVarTQOAqqoq5ZEjR+KrqqqUhYWF\nUSkpKZnWlNq+ffvrOp0uWa/XK/R6vaKwsDAKAHQ6XbKXl1ezXq9X7N69++XU1NT9jtxhIprdfve7\n30EsFsPf3x+urq7jug3GbDCmeBx8+FFQUKDWaDRZAKDRaLLy8vJiACA/P3/jli1bDru6uvbI5fI6\nf3//S2VlZWtNJpNPe3v7QuuRRmJi4lvWdQZuKzY29lhJSUn4RO4gEdFAJ06cQGtrK5qbm3Hs2LFb\nTkGlMQw0i0QiYcOGDcVz587t27Zt22+feuqpgxaLRSwWiy0AIBaLLRaLRQwATU1NvmFhYWes68pk\nMoPRaJS6urr2yGQyg7VdKpUajUajFACMRqPUz8+vEQBcXFx6PTw82lpaWjw9PT1bBtaRkZFhe61S\nqaBSqe5ox4mIZprS0tIhV2mP16ih8OGHH37Dx8fH9OWXXy6JiIg4GRgYWD1wvkgkEkQikcNHqgaG\nAhERDTX4B/PgK7bHYtTuIx8fHxMALFmy5MvHHnvsnfLy8lCxWGwxm80SADCZTD7e3t5fADePABob\nG/2s6xoMBplMJjNIpVKjwWCQDW63rtPQ0HAPAPT29rq0tbV5DD5KICKiyTFiKNy4ceMr7e3tCwHg\n+vXrdxcVFUUGBQVVqtXqgqysLA0AZGVlaWJiYvIAQK1WF+Tk5Gzu7u6+q7a2dpler1eEhoaWSyQS\ns7u7+9WysrK1giCIsrOzEzZu3JhvXce6raNHj24KDw8vcewuE5GVu9vdtnPzOU3faeAtPu7UiN1H\nFotF/Nhjj70D3PwV//jjj/93ZGRk0QMPPPC/cXFxuTqdLlkul9fl5ubGAYBSqayKi4vLVSqVVS4u\nLr2ZmZkp1q6lzMzMlKSkpDc7OjrmR0dHH4+KiioEgOTkZF1CQkK2QqHQe3l5Nefk5GyesL0johHl\nh33fbvtI1zf83JCH4kvD3ySOprcRQ2HZsmW158+fDx7c7unp2VJcXLzB3jrp6el709PT9w5uDwkJ\nqaisrAwa3O7m5tZlDRUiInIuXrFBREQ2DAUiIrJhKBARkQ1DgYiIbBgKRERkw1AgIiIbhgIREdkw\nFIiIyIahQERENgwFIiKyYSgQEZENQ4GIiGwYCkREZMNQICIim1Efx0lEkyctbT/M5g6785qbrwAy\nu7OIJgxDgWgKMZs7IJdn2J13tu/Q5BZDsxK7j4iIyIZHCkRTyOWKYrSer7M7T+hqndxiaFZiKBBN\nIXd1XMMzMrndebv7+ya3GJqV2H1EREQ2DAUiIrJhKBARkQ1DgYiIbMYUCn19fXNXr1597tFHH/0T\nALS0tHhGREScDAgIqImMjCxqbW1dZF123759zykUCn1gYGB1UVFRpLW9oqIiJCgoqFKhUOh37dp1\nwNre1dXlFh8ff0ShUOjDwsLO1NfXL53IHSQiorEbUygcOHBgl1KprBKJRAIAaLXatIiIiJM1NTUB\n4eHhJVqtNg0AqqqqlEeOHImvqqpSFhYWRqWkpGQKgiACgO3bt7+u0+mS9Xq9Qq/XKwoLC6MAQKfT\nJXt5eTXr9XrF7t27X05NTd3vqJ0lIqKRjRoKBoNBdvz48eitW7f+3voFX1BQoNZoNFkAoNFosvLy\n8mIAID8/f+OWLVsOu7q69sjl8jp/f/9LZWVla00mk097e/vC0NDQcgBITEx8y7rOwG3FxsYeKykp\nCXfUzhIR0chGvU5h9+7dL//iF7/4wdWrV92tbRaLRSwWiy0AIBaLLRaLRQwATU1NvmFhYWesy8lk\nMoPRaJS6urr2yGQyg7VdKpUajUajFACMRqPUz8+vEQBcXFx6PTw82lpaWjw9PT1bBtaRkZFhe61S\nqaBSqW5zl4mIZqbS0lKUlpbe0TZGDIV33333W97e3l+sXr36XGlpqcreMiKRSLB2KznSwFAgIqKh\nBv9g3rNnz7i3MWIo/O1vf/t6QUGB+vjx49GdnZ3zrl696p6QkJAtFostZrNZIpFIzCaTycfb2/sL\n4OYRQGNjo591fYPBIJPJZAapVGo0GAyywe3WdRoaGu7x9fVt6u3tdWlra/MYfJRARESTY8Qxhb17\n96Y3Njb61dbWLsvJydn88MMP/yU7OztBrVYXZGVlaQAgKytLExMTkwcAarW6ICcnZ3N3d/ddtbW1\ny/R6vSI0NLRcIpGY3d3dr5aVla0VBEGUnZ2dsHHjxnzrOtZtHT16dFN4eHiJo3eaiIjsG9e9j6zd\nRGlpadq4uLhcnU6XLJfL63Jzc+MAQKlUVsXFxeUqlcoqFxeX3szMzBTrOpmZmSlJSUlvdnR0zI+O\njj4eFRVVCADJycm6hISEbIVCoffy8mrOycnZPNE7SUREYzPmUFi/fv3p9evXnwYAT0/PluLi4g32\nlktPT9+bnp6+d3B7SEhIRWVlZdDgdjc3ty5rqBARkXPximYiIrLhrbOJaFyam5uQkZQ0pH2+RIJU\nrXbyC6IJxVAgonGZ39eLDLl8SHtGXd2k10ITj91HRERkw1AgIiIbhgIREdlwTIGIxqWzqwt5eaVD\n2iv6rk1+MTThGApENC5CP7BokWpIe4chb/KLoQnH7iMiIrJhKBARkQ1DgYiIbBgKRERkw1AgIiIb\nhgIREdkwFIiIyIahQERENgwFIiKyYSgQEZENQ4GIiGwYCkREZMNQICIiG4YCERHZMBSIiMhmxFDo\n7Oyct3bt2rLg4ODzSqWy6rnnntsHAC0tLZ4REREnAwICaiIjI4taW1sXWdfZt2/fcwqFQh8YGFhd\nVFQUaW2vqKgICQoKqlQoFPpdu3YdsLZ3dXW5xcfHH1EoFPqwsLAz9fX1Sx2xo0RENLoRQ2HevHmd\np06deuj8+fPBFy5cWHnq1KmH/vrXvz6o1WrTIiIiTtbU1ASEh4eXaLXaNACoqqpSHjlyJL6qqkpZ\nWFgYlZKSkikIgggAtm/f/rpOp0vW6/UKvV6vKCwsjAIAnU6X7OXl1azX6xW7d+9+OTU1db/jd5uI\niOwZtfvoK1/5yg0A6O7uvquvr2/u4sWLrxQUFKg1Gk0WAGg0mqy8vLwYAMjPz9+4ZcuWw66urj1y\nubzO39//UllZ2VqTyeTT3t6+MDQ0tBwAEhMT37KuM3BbsbGxx0pKSsIdtbNERDSyUR/H2d/fP+f+\n++8/+/nnny/fvn376/fdd9+nFotFLBaLLQAgFostFotFDABNTU2+YWFhZ6zrymQyg9FolLq6uvbI\nZDKDtV0qlRqNRqMUAIxGo9TPz68RAFxcXHo9PDzaWlpaPD09PVsG1pGRkWF7rVKpoFKp7mjHiZxp\nf1oaOszmIe3tzSZA5oSCaEYoLS1FaWnpHW1j1FCYM2dO//nz54Pb2to8HnnkkfdOnTr10MD5IpFI\nEIlEwh1VMQYDQ4Fouuswm5Ehlw9pL+rrmfxiaMYY/IN5z549497GmM8+8vDwaPuP//iPP1dUVISI\nxWKL2WyWAIDJZPLx9vb+Arh5BNDY2OhnXcdgMMhkMplBKpUaDQaDbHC7dZ2GhoZ7AKC3t9elra3N\nY/BRAhERTY4RQ+Ef//jHV61nFnV0dMw/efJkxOrVq8+p1eqCrKwsDQBkZWVpYmJi8gBArVYX5OTk\nbO7u7r6rtrZ2mV6vV4SGhpZLJBKzu7v71bKysrWCIIiys7MTNm7cmG9dx7qto0ePbgoPDy9x7C4T\nOV9Fxd+Rl1c6ZOrs6nZ2aTTLjdh9ZDKZfDQaTVZ/f/+c/v7+OQkJCdnh4eElq1evPhcXF5er0+mS\n5XJ5XW5ubhwAKJXKqri4uFylUlnl4uLSm5mZmWLtWsrMzExJSkp6s6OjY350dPTxqKioQgBITk7W\nJSQkZCsUCr2Xl1dzTk7OZsfvNpFzdXT0YpFMNaRd6P948oshGmDEUAgKCqo8e/bs/YPbPT09W4qL\nizfYWyc9PX1venr63sHtISEhFZWVlUGD293c3LqsoUJERM7FK5qJiMiGoUBERDYMBSIismEoEBGR\nDUOBiIhsGApERGTDUCAiIptR731ERNNXn9CN6rq8YecRDcZQIJoBhvvyF9CHefJFdtYAYHD4fSxp\nGmIoEE0TI/3qH/bL3zC0iWgkDAWiKeS2vviB2/ryF4R+di3REAwFoilFmNAv/hGJwK4lGoJnHxER\nkQ2PFIicoLnVhOpeewPD/U6oZqjb6Vpqbm5CRlKS3XnzJRKkarUTVR45EEOByAn60DO1B4Zvo2tp\nfl+v3UeMAkBGXd3E1EUOx+4jIiKyYSgQEZENu4+IaEJ0dnUhL6/U7ryKvmuTWwzdNoYCkYOkZaTB\n3Gq2O6+zu2OSq5k4ww1C9/Z1Y9Eild11Ogz2B61p6mEoEDmIudUMeYzc7rwPD02Ns4xuy3CD0Ly2\nYUbgmAIREdkwFIiIyIahQERENiOOKTQ2NvolJia+9cUXX3iLRCLh6aef/t3OnTtfaWlp8YyPjz9S\nX1+/VC6X1+Xm5sYtWrSoFQD27dv33KFDh56cO3du3yuvvLIzMjKyCAAqKipCkpKS3uzs7JwXHR19\n/MCBA7sAoKuryy0xMfGts2fP3u/l5dV85MiR+KVLl9Y7fteJ7txIg8kV5yuGHVOYiUa6Crr5mmmS\nq6HbNeKRgqura8/LL7+8+9NPP73vzJkzYa+99tr/vXjx4r1arTYtIiLiZE1NTUB4eHiJVqtNA4Cq\nqirlkSNH4quqqpSFhYVRKSkpmYIgiABg+/btr+t0umS9Xq/Q6/WKwsLCKADQ6XTJXl5ezXq9XrF7\n9+6XU1NT9zt+t4kmhnUw2d7UMY3PMLot/xyAtjf1ocfZ1dEYjRgKEonEHBwcfB4AFixYcO3ee++9\naDQapQUFBWqNRpMFABqNJisvLy8GAPLz8zdu2bLlsKura49cLq/z9/e/VFZWttZkMvm0t7cvDA0N\nLQeAxMTEt6zrDNxWbGzssZKSknBH7jAREQ1vzKek1tXVyc+dO7d67dq1ZRaLRSwWiy0AIBaLLRaL\nRQwATU1NvmFhYWes68hkMoPRaJS6urr2yGQy211dpFKp0Wg0SgHAaDRK/fz8GgHAxcWl18PDo62l\npcXT09NgKqUKAAANoUlEQVSzZeDfz8jIsL1WqVRQqVS3tcNERDNVaWkpSktL72gbYwqFa9euLYiN\njT124MCBXQsXLmwfOE8kEgkikcjhJygPDAUiIhpq8A/mPXv2jHsbo4ZCT0+Pa2xs7LGEhITsmJiY\nPODm0YHZbJZIJBKzyWTy8fb2/gK4eQTQ2NjoZ13XYDDIZDKZQSqVGg0Gg2xwu3WdhoaGe3x9fZt6\ne3td2traPAYfJRA523ADyrNtMJlmvhFDQRAEUXJysk6pVFY988wzv7K2q9XqgqysLE1qaur+rKws\njTUs1Gp1wXe/+93/efbZZ18yGo1SvV6vCA0NLReJRIK7u/vVsrKytaGhoeXZ2dkJO3fufGXgtsLC\nws4cPXp0U3h4eIljd5lo/Ia7Ovmv5X+d/GKmoc6uG0h6JsnuPMkiCbQZfNbCVDFiKHz44Yff+MMf\n/vC9lStXXli9evU54OYpp2lpadq4uLhcnU6XbD0lFQCUSmVVXFxcrlKprHJxcenNzMxMsXYtZWZm\npiQlJb3Z0dExPzo6+nhUVFQhACQnJ+sSEhKyFQqF3svLqzknJ2ezo3eaiCaXIBKGPaKqy6ub1Fpo\nZCOGwoMPPvjX/v5+u2coFRcXb7DXnp6evjc9PX3v4PaQkJCKysrKoMHtbm5uXdZQIZpJ+uubcV5r\n/7x9l67eSa6GaGx4QzwiB1nQ04dnJPafXrZzGt8Pj2Y2hgKRg3R2daO6us7uPEHgHUVpamIoEDmI\nIAiYN08+zNyGySyFaMwYCkT/xPsYOUdFRYXdM5N4VpJzMBSI/mmkh+Lw1FPH6ejrsPu586wk5+Ct\ns4mIyIahQERENgwFIiKyYSgQEZENQ4GIiGx49hERTUnDnaoK8HRVR2Io0KzCaxGmj+FOVQV4uqoj\nMRRoVuG1CEQjYygQkcP1dfYOe8fY/vrmSa6GRsJQICKHWyBg+DvGnr08ydXQSHj2ERER2fBIgWYk\nPlN5ahEEYdjbiHd2dY97ezwzyXEYCjQj8ZnKU89wtxEXBNO4t8UzkxyH3UdERGTDUCAiIhuGAhER\n2TAUiIjIhgPNNG3xlhVkDx/veWdGPFJ48sknD4nFYktQUFClta2lpcUzIiLiZEBAQE1kZGRRa2ur\n7YqUffv2PadQKPSBgYHVRUVFkdb2ioqKkKCgoEqFQqHftWvXAWt7V1eXW3x8/BGFQqEPCws7U19f\nv3Sid5BmLusZRvamju4OZ5dHTmI9M2nwNNwPCLrViKHwxBNPvFFYWBg1sE2r1aZFREScrKmpCQgP\nDy/RarVpAFBVVaU8cuRIfFVVlbKwsDAqJSUlUxAEEQBs3779dZ1Ol6zX6xV6vV5h3aZOp0v28vJq\n1uv1it27d7+cmpq631E7SkREoxsxFNatW/fB4sWLrwxsKygoUGs0miwA0Gg0WXl5eTEAkJ+fv3HL\nli2HXV1de+RyeZ2/v/+lsrKytSaTyae9vX1haGhoOQAkJia+ZV1n4LZiY2OPlZSUhDtiJ4mIaGzG\nPaZgsVjEYrHYAgBisdhisVjEANDU1OQbFhZ2xrqcTCYzGI1Gqaura49MJjNY26VSqdFoNEoBwGg0\nSv38/BoBwMXFpdfDw6OtpaXF09PTs2Xw383IyLC9VqlUUKlU4y2diGhGKy0tRWlp6R1t444GmkUi\nkSASiYQ7qmCMBoYCERENNfgH8549e8a9jXGHglgstpjNZolEIjGbTCYfb2/vL4CbRwCNjY1+1uUM\nBoNMJpMZpFKp0WAwyAa3W9dpaGi4x9fXt6m3t9elra3Nw95RAs1eU/0Mo+vVRhSk/MbuvLkdPZNc\nDY2E90sam3GHglqtLsjKytKkpqbuz8rK0sTExORZ27/73e/+z7PPPvuS0WiU6vV6RWhoaLlIJBLc\n3d2vlpWVrQ0NDS3Pzs5O2Llz5ysDtxUWFnbm6NGjm8LDw0smegdpepvqD8VZ0NuPnd6Bduf9suaL\nSa6GRsL7JY3NiKGwZcuWw6dPn17/j3/846t+fn6NP/3pT3+clpamjYuLy9XpdMlyubwuNzc3DgCU\nSmVVXFxcrlKprHJxcenNzMxMsXYtZWZmpiQlJb3Z0dExPzo6+nhUVFQhACQnJ+sSEhKyFQqF3svL\nqzknJ2ez43eZiIiGM2IoHD58eIu99uLi4g322tPT0/emp6fvHdweEhJSUVlZGTS43c3NrcsaKkRE\n5Hy8zQUREdnwNhfkdFN9MJlmPg5C/wtDgZxuqg8m08zHQeh/YfcRERHZMBSIiMiGoUBERDYcU6BJ\nM9yAMgeTZ7f+vj7k5ZUOaW9ubp38YuyYbYPQDAWaNMMNKHMweXYTBGDRItWQ9st9Fya/GDtm2yA0\nu4+IiMiGRwpEY9Bf34zz2rwh7S5dvU6oZmYRuvvRmTf0s/1K8w0nVEMMBZpQM/VCtAU9fXhGsmhI\n+85+JxQzwywA8H8WDf1sX7g8KXflvyMz8XnQDAWaULwQjWaT4cYbpvNYA0OBaAw6u7pRXV03pF0Q\npv6vWaLxYCjQuM3ULqKRCIKAefPkduY0THYpRA7FUKBxYxcR0cim87UNDAUiogk2na9tYCjQsHgF\nMtHEm+pHEQwFGhavQCZnGu72F8DUuQXG7ZjqRxEMhVluNg4a0/TQ39WPeaX2v/zdvrg2ydXMHgyF\nWY6Dxv8y3FXLAK9cdobhLmoDpseFbbdjKlwMx1Ag+ieXax2IarX/y/Tzvpn5JURTy1S4GI6hMAuw\ni2hshr8WAeD1CORMkzk4zVCYQUY6Wyg2I9b2vu58HeTBcgCzr4tool3nFc0Tqrq1FYHDdBnNZiMN\nTh/70bFhf/TdjikRCoWFhVHPPPPMr/r6+uZu3br196mpqfudXdNUNdqv/oFf/laDv/gHhgLdmevM\nhAk11lDo7ehBQcpv7M67bm6b6LKmtJECAwfGvz2nh0JfX9/cHTt2/Lq4uHiDVCo1rlmz5mO1Wl1w\n7733XnR2bY400pd79d+rEfi1QLvzhvviB/irfyyuVxuH/TKZ29EzydXQ7VogADu97f8/8oLho0mu\nZmZxeiiUl5eH+vv7X5LL5XUAsHnz5pz8/PyNUy0UbvdLfLh5o32584ygOzPcl79Hexd2rrT/3+qX\nNV84uiyaBMMdRbQZmmfkdQ8TThAEp05//OMfN23duvWg9X12dvb3duzY8erAZQAInDhx4sRp/NN4\nv5OdfqQgEomE0ZYRBEE0GbUQEc12Tn9Gs1QqNTY2NvpZ3zc2NvrJZDKDM2siIpqtnB4KDzzwwP/q\n9XpFXV2dvLu7+64jR47Eq9XqAmfXRUQ0Gzm9+8jFxaX317/+9Y5HHnnkvb6+vrnJycm6qTbITEQ0\nazh7oHmkKTc39ztKpfLTOXPm9FVUVNw/cN7evXuf8/f3169YsaL6vffei3R2rdNt+slPfpIhlUoN\nwcHB54KDg8+dOHEiytk1TbfpxIkTUStWrKj29/fXa7XaVGfXM92npUuX1gUFBV0IDg4+t2bNmnJn\n1zPdpieeeOKQt7e35Wtf+1qlta25udlzw4YNJxUKRU1ERETRlStXFo22HafvyEjTxYsXAz/77LMA\nlUp1amAofPrpp8pVq1ad7+7udq2trZUvX778Ul9f3xxn1zudpoyMjJ/88pe/fNbZdUzXqbe3d+7y\n5csv1dbWyru7u11XrVp1vqqq6l5n1zWdJ7lcXtvc3Ozp7Dqm6/T++++vO3v27OqBofCDH/zgxf37\n9/8/QRCg1WpTU1NTtaNtx+ljCiMJDAysDggIqBncnp+fv3HLli2HXV1de+RyeZ2/v/+l8vLyUGfU\nOJ0JPKvrtg28vsbV1bXHen2Ns+ua7vhv8vatW7fug8WLF18Z2FZQUKDWaDRZAKDRaLLy8vJiRtvO\nlA6F4TQ1NfkOPENJJpMZjEaj1Jk1TUevvvrqf65ateqT5ORkXWtrK284Mw5Go1Hq5+fXaH3Pf4N3\nTiQSCRs2bCh+4IEH/vfgwYNPObuemcBisYjFYrEFAMRiscVisYhHW8fpA80REREnzWazZHD73r17\n0x999NE/jXU7Y7neYbYZ7rN94YUXfrh9+/bXf/zjH/8UAH70ox/97L/+679+qdPpkie/yumJ/94m\n3ocffvgNHx8f05dffrkkIiLiZGBgYPW6des+cHZdM4VIJBLG8u/W6aFw8uTJiPGuM/jaBoPBIJNK\npcaJrWz6G+tnu3Xr1t+PJ4CJ19c4go+PjwkAlixZ8uVjjz32Tnl5eShD4c6IxWKL2WyWSCQSs8lk\n8vH29h71Xi7TpvtoYF+jWq0uyMnJ2dzd3X1XbW3tMr1erwgNDS13Zn3Tjclk8rG+fueddx4LCgqq\ndGY90w2vr5lYN27c+Ep7e/tCALh+/frdRUVFkfw3eefUanVBVlaWBgCysrI0MTEx9h8tOJCzR8xH\nmt5+++3HZDJZ47x58zrEYrE5KirqhHXeCy+8kL58+fJLK1asqC4sLHzE2bVOtykhIeGtoKCgCytX\nrvxk48aNeWazWezsmqbbdPz48W8GBAR8tnz58kt79+59ztn1TOfp8uXLy1atWnV+1apV5++7776/\n8/Mc/7R58+bDPj4+Ta6urt0ymazx0KFDTzQ3N3uGh4cXj+eUVJEgsGuUiIhumjbdR0RE5HgMBSIi\nsmEoEBGRDUOBiIhsGApERGTDUCAiIpv/DwQoNKH2U/IlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4152ae50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Weighting with approximated densities\n",
    "ratio_norm = completeRatio / completeRatio.sum()\n",
    "#ratio_norm = trueRatio/trueRatio.sum()\n",
    "data_f1_idx = np.random.choice(range(data_f1.shape[0]),data_f1.shape[0],p=ratio_norm)\n",
    "data_f1_re_clf = data_f1[data_f1_idx]\n",
    "traindata_weighted_clf = traindata.copy()\n",
    "traindata_weighted_clf[targetdata==0] = data_f1_re_clf[:,:-1]\n",
    "# Weighting with true densities\n",
    "ratio_norm = trueRatio / trueRatio.sum()\n",
    "data_f1_idx = np.random.choice(range(data_f1.shape[0]),data_f1.shape[0],p=ratio_norm)\n",
    "data_f1_re_true = data_f1[data_f1_idx]\n",
    "traindata_weighted_true = traindata.copy()\n",
    "traindata_weighted_true[targetdata==0] = data_f1_re_true[:,:-1]\n",
    "#traindata[targetdata == 0] = data_f0[:,:-1]\n",
    "plt.hist(data_f0[:,feature],label='f0',alpha=0.5,bins=bins,range=x_range)\n",
    "plt.hist(data_f1[:,feature],label='f1',alpha=0.5,bins=bins,range=x_range)\n",
    "plt.hist(data_f1_re_clf[:,feature],label='weighted f1',alpha=0.5,bins=bins,range=x_range)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the approximation by training a classifier to classify between the weighted data and the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 1.47303094775 (+/- 0.000782237553734)\n",
      "Trained\n"
     ]
    }
   ],
   "source": [
    "xgb = XGBoostClassifier(num_class=2, nthread=4, silent=1,\n",
    "          num_boost_round=50, eta=0.1, max_depth=3)\n",
    "rng = np.random.RandomState(1234)\n",
    "indices = rng.permutation(traindata.shape[0])\n",
    "traindata = traindata[indices]\n",
    "targetdata = targetdata[indices]\n",
    "traindata_weighted_clf = traindata_weighted_clf[indices]\n",
    "traindata_weighted_true = traindata_weighted_true[indices]\n",
    "scores = cross_validation.cross_val_score(xgb, traindata_weighted_clf, targetdata)\n",
    "print \"Accuracy: {0} (+/- {1})\".format(scores.mean(), scores.std() * 2)\n",
    "xgb.fit(traindata_weighted_clf,targetdata)\n",
    "predictions_weighted_clf = xgb.predict_proba(traindata_weighted_clf)[:,1]\n",
    "xgb.fit(traindata_weighted_true,targetdata)\n",
    "predictions_weighted_true = xgb.predict_proba(traindata_weighted_clf)[:,1]\n",
    "xgb.fit(traindata,targetdata)\n",
    "predictions = xgb.predict_proba(traindata)[:,1]\n",
    "print 'Trained'\n",
    "# Fit also on original data and another on weighted data with true ratios\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEVCAYAAADgh5I1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdczfsfB/BXE6W0iJZIIiWFbLKuvW+Ea147I3tfXNfI\nvLKyLy4/4qKQSiqjUkZDC6UllZY67c45n98fn5ublOqM5uf5eHwfnfE9n/M+0fd9PluCEAKGYRiG\nkazpABiGYZjagSUEhmEYBgBLCAzDMMy/WEJgGIZhALCEwDAMw/yLJQSGYRgGAEsIDFNrzJ49+6+t\nW7furOk4mIaLJQSm2ujq6sbIycnlKigocFq2bJk0Y8aMy1lZWYolz/Hx8ek9aNAgD0VFxSwlJaUv\nY8eOdQoPD+9Y8pysrCxFGxubP1u3bh2roKDAadeuXeTKlSsPp6WlqVbvJxItCQkJIiEhIbKJQenp\n6SrNmzdP6dev39OSjwcGBnbp2rXrK3l5+Zxu3bq9DAoKMhHVezJ1G0sITLWRkJAg9+7dG83hcBSC\ngoJM3rx5Y/zHH39sKX7e19e317Bhw1wnTJhwOzExsVV0dHQbExOToD59+nhHR0e3AYDCwkLZwYMH\nPwoPD+/o6uo6jMPhKPj6+vZSU1NL9ff3N6+5T1f7rF+/3tbQ0DCsZJIpLCyUHTdunOPMmTMvffny\nRWnWrFkXx40b51hUVCRTk7EytQQhhB3sqJZDV1c3+tGjR4OK769du3bfyJEj7xff79u371Nra+tj\npV83YsQI55kzZ14khODMmTPz1NXVk3JycuQq+74SEhJ8e3v7hfr6+u+UlJQyynqPso5z587N7dix\nY5iysnL6sGHDXGJjY3VKlmlnZ7esbdu2UWpqailr167dx+fzJQgh4PP5Ejt37tzSunXrmBYtWiTP\nnDnzYmZmpmLxa58+fdq3V69ePkpKShna2tpxFy9enEkIwezZsy9YW1sfGzVq1D0FBYWsHj16PI+K\nimpb/DobG5vDLVq0SFZUVMw0NjYODgkJ6VRe7N7e3r179erlc+HChdl9+/Z9Wvy4q6vrT5qamh9L\nnqujoxPr4uIyrKb/f7Cj5o8aD4AdDefQ1dWNdnd3H0wIQXx8vJaxsXHwjh07fiOEICcnR05KSorr\n5eU1oPTrLly4MLtVq1afCCGYMmXKtdmzZ1+oyvtKSEjwx4wZ45SZmakYFxen3bx5888VXQDv3Lkz\nrl27du8jIiIMeDye5B9//LG5d+/e3iXLHDRo0KOMjAyluLg47fbt2789e/bsr4TQRNKuXbv30dHR\nutnZ2fITJ078Z8aMGZcIIYiJiWmtoKCQde3atSlcLlcqLS1NJTAw0IQQglmzZv2lqqqa+uLFi25c\nLldq+vTpf1tZWf2PEAIXF5dhXbt2fVmcWCIiIgwSExNblhU7l8uVMjMze/X69WvT0gnh0KFDK0eM\nGOFc8vwxY8Y4HTx4cFVN//9gR80frMmIqTaEEInx48ffUVRUzNLR0YnT09OL2rJlyx8Abe/m8/mS\nrVq1Siz9upYtWyalpqaqAUBaWppqWedUZMOGDXsVFRWztLW14wcOHOgZGBjY5Ufn29vbL9q4ceMe\nAwODt5KSkvyNGzfuCQwM7BIfH69dfM769ettlZSUvmhra8fb2Nj8+b///W8qAFy5cmX66tWrD+rq\n6sbIy8vn7NmzZ+O1a9eseDye1NWrV6cNHTr04ZQpU65LSUnxVFRU0k1MTIIA2qQ2ceLEW926dXsp\nJSXFmz59+pXiOGVkZIo4HI5CeHh4Rz6fL2lgYPC2ZcuWSWXFbmdnt7xnz57PTU1NA0o/l52d3bRZ\ns2aZJR9TVFTM4nA4ClX9nTL1D0sITLWRkJAgjo6O47KyshS9vLwsPDw8Br18+bIbACgrK2dISkry\nExMTW5V+XWJiYqvmzZunAICamlrqp0+fNKr63iUvnnJycrnZ2dlNf3R+bGxs6xUrVhxRVlbOUFZW\nzlBVVU0DgISEBM3ic7S1teOLb+vo6MQVx5WYmNiqdevWsSWf43K50snJyeofP37Uatu27Yfy3ldd\nXT25+HaTJk3yiuMcNGiQx9KlS49ZW1sfV1dXT164cOGpsi7inz590jh69Oiykn0zJSkoKHBKd+Rn\nZmY2U1RUzPrR74NpGFhCYGpE//79nyxbtuzo+vXrbQFAXl4+p1evXr4ODg6TS5/r4OAwefDgwY8A\nYMiQIe6urq7DcnNz5cQZn46OTtzp06cXZGRkKBcfOTk58j179nxefE5cXJxOyduampoJAKChofEp\nJiZGt+Rz0tLS3JYtWyZpa2vHR0VF6QkS07Jly46+fPmyW1hYmOG7d+/a79+/f23pc/z9/c0TExNb\nGRoahrVq1SrRxsbmT39/f3MNDY1PhBAJQ0PDsODg4M4lXxMcHNy5U6dOoYLExNQzNd1mxY6Gc5Tu\nVE5JSVGTk5PLef78eQ9CCJ49e9ZHXl4+287ObllWVpZCenq68ubNm/9QVlZOj4yM1COEoKCgQLZ7\n9+7+w4cPf1Dcvp+amqq6a9euTc7OziPKel8JCQl+yc7ZWbNm/bVly5adP4r19u3b442MjN6EhoYa\nEkLw5cuXZg4ODpYlyxwyZMjD4j6EDh06hJ85c2YeIQRnz579VV9f/110dLQuh8NpOmnSpJvFfQix\nsbE6CgoKWQ4ODpZFRUXSqampqiX7EErG5enpaaGlpRVPCMGLFy+6PX/+vEdhYaFMdna2/PDhwx9s\n3759W+m4CwoKZJOTk1sUH0eOHFneo0eP58nJyS0IISgsLJRp3bp1zJEjR5bn5+c3OnLkyHJdXd3o\noqIi6Zr+/8GOmj9YDYGpMWpqaqmzZs26aGtrux4A+vTp4+3q6jrs1q1bEzU0ND7p6urGBAUFmTx7\n9qyvnp5eFADIysoWuru7D+nQoUPE0KFDHzZr1iyzR48efunp6Solv72XVHpsf2XG+48fP/7O+vXr\nba2srK41a9Ys09jY+I2rq+uwkueMGzfOsWvXrq9MTU0DRo8efW/u3LnnAWDu3LnnZ8yYcbl///5P\n2rZt+0FOTi736NGjywBa83B2dh558ODB1aqqqmmmpqYBxd/Yy4qr+H5WVpbiggULTquoqKTr6urG\nqKmppa5du3Z/6bhlZWULW7Ro8bn4aNasWWbxYwDti7hz5874S5cuzVRWVs64dOnSzDt37oyXlpbm\n/uj3wTQMEoSwDXIYpqokJSX5kZGR7X7UH8AwdQ2rITAMwzAAWEJgGrBFixbZKygocEofS5YsOVHR\na0W5xATD1BasyYhhGIYBAEjXdADlYd/AGIZhBEMIkRDkdbW6yaimh2DVlmPbtm01HkNtOdjvgv0u\n2O/ix4cwanVCYBiGYaoPSwgMwzAMAJYQ6gQLC4uaDqHWYL+L/7DfxX/Y70I0au0oIwkJCVJbY2MY\nhqmtJCQkQOpjpzLDMAxTfUSeEObOnXteXV092djY+E155yxfvtxOX1//vYmJSVBAQICpqGNgGIZh\nqk7kCWHOnDkXXFxchpf3vLOz88jIyMh279+/1z99+vSCxYsXnxR1DAzDMEzViTwh9OvX76mysnJG\nec87OTmNnTVr1kUA6NGjh9+XL1+UkpOT1UUdB8MwDFM11T5TOSEhQbPkTlNaWlofP378qFVyp6hi\n27dv/3rbwsKCjSRgmLqAEIDH+/7g8//7yefT80rfLuux8p4vfq/iwSdl/RTmsZI/K/u5q/lcQgju\n+vvA7UUQsvOLkFdQVPlyy1AjS1eU7gEvb5mKkgmBYRo0QoD8fCAnB8jKAgoKgNxceuTl0SM//9uj\noIAeJW8XFACFhUBR0fdHYSHA5X57FBV9/xiP99/P0rd5PBqrpCQgJfX9ISn5/SEh8e3P8m6Xfkzi\n38tI6dulfwrzWMmflSHmcwmAzOxCfExPQVphEgolM9CEL4VxfBlIQBISRAoOlS/1O9WeEDQ1NRNK\nblT+8eNHreKtBxmmXikqAr58ATIy6E8Oh17Miw8O578jO5sexfdzcuj9nBx65OYC0tJA03+3gpaT\nA5o3B5o0obebNAEaN6ZHkyZAo0b/HY0bA82affuYjMz3h6ws/Skt/e0hI0Mv5sU/paTo4yV/lr7o\nV+Vix5SLxwOePwdevAAee+fhYe5+5AzfBh3OZPTTWIkZPUehl1ErKCqW+H0L8buv9oQwduxYp2PH\nji21srK69vz5855KSkpfymouYphapagISEkBPn+mR0oKPVJTgfR0ICaGXrSzsuj99HR6v1kzQEWF\n/mzWDFBQABQVv/2poQHIy9OLffFj8vL/HU2b0ou+dK1di5IRIS4XuH8f8PIC7twB5OQJNIfcQkCX\n1Rig1R0nxsWgtVJrsby3yCemTZ069X+PHz8ekJqaqqaurp68Y8eObUVFRTIAsHDhwlMAsHTp0mMu\nLi7D5eXlcy5cuDDHzMzs9XeBsYlpTHXIyQE+fvzvSEgAkpLokZxML/5JSfTbuqoqoK5Ov5mXPFRU\nADU1etHW0KD3lZXpxV2STfVhKkYI4O8PXL8OODgAbdoAw4cD7fu9wcnoFUjJTYHdcDsMbDMQAMDn\n83H+/HmMGjUKrVq1+qYsYSamsZnKTP3F4wGJiUBsLBARQb/Rh4XRn8UJoKAA0NL679DQAFq1ohd+\ndXWgRQv6U0WFXdwZkeLxgNevgRs3aBJo3BiYMoUeLdukY5vXNlwPuY5tA7ZhYbeFkJakNcSIiAgs\nWLAAhYWF+Pvvv9GuXbtvyhUmIbA6KFO3FRUBHz4Ab98C794BkZH0/ocPQFQU0LIl0Lo1vdBraADm\n5vQxfX1AW5t+k2ft3Uw1io4GTp8GLl6krYgTJwJOToCxMcAnPJx5fQbbjm/DpI6TEGYdBjU5NQBA\nQUEB9u7di6NHj2L79u1YvHgxpKSkRBobSwhM3fDlC73oR0R8e0RH02/27dvTw9gYGDeO1rlbt6Yd\nrAxTw/LzAUdH4MIF2jQ0ezbg7g4YGv53ztPYp1jushyKjRTh9osbTFqafH2usLAQ3bt3R5s2bRAQ\nEABtbe3v30QEWJMRU7twufTC//o1EBxMjzdv6MibDh3oYWDw3+127Whdm2FqGUKAly9pErh+HTAz\nA+bMASZM+PZ7SnxmPNa5r4N3nDf2D92PyZ0mQ6KMWuubN29gZGRU5nMlsT4Epm4qKqJt+q9e0QRQ\nnAQ0NABTU8DEBOjcmX7r19FhTTtMnZCSAvz1Fz3y82ltYNYs+l+4pHxuPg74HMDh54dh3d0a6/us\nh7ysvNDvz/oQmLohIQF49gzw9QV8fIDQUEBXl351MjMDLC2BLl1owyrD1DFJScD+/bRGMHYsYG8P\n9O37/fcYQgjuRNzBarfVMG1lipfzX6KNcpuvz3/58gVKSkrVHD3FEgIjHoTQpp+nT4GHD+nsmrQ0\noFMnWmc+eJAmAXnhvxExTE0qmQh++YW2cGpqln1u6OdQrHBZgaTsJJwZcwaD2w7++hyfz8fJkyex\nfft2vHz5Eq1bi2euwY+whMCITkIC8OgR7S1zd6czW/v3B4YMAbZsAYyM2NBNpt5ISgL27aNNQzNm\nACEhtLWzLBl5Gdj+eDuuvrmK3/r/hsXdF38dRgoAoaGhmD9/PiQlJfH48eMaSQYASwiMMHJy6IX/\n/n1aE/j8GRg0CBg6FNi+HWjbtqYjZBiRq0oi4PF5OBdwDr95/obxHcYjbEkYmss3//p8QUEBdu3a\nhZMnT2Lnzp1YsGABJGvwSxNLCEzVxMcD9+7RgdPPntERP5MmAQsX0vZ/EY+LZpjaomQimDnzx4kA\nALzjvLHswTLIy8rjwfQHMG31/V5gBQUFSEpKQmBgIDTLa2eqRmyUEfNjfD4d/XP3Lk0C8fHAiBHA\nTz8BY8YANdT5xTDVJTGRJoKLF2kiWLfux4ngY9ZHrHdfjyexT7BvyD5YGVlVOFRUlNgoI0a08vNp\nX4CTE60NKCjQi7+dHdCrF1tkjWkQSieC0FC6qkl58rn5OOR7CAd9D2Jxt8U4Pfq0SIaRVif2l81Q\nfD7w+DFw6RJw+zYd/z9uHODpSWcAM0wDkZgI2NrSP4VZsypOBIQQOL11wiq3Veis3hkv5r9AW+Vv\n+8/i4uKwf/9+HDhwAI0aNRLzJxAcSwgN3adPwLlzdMycvDydSrlnD13vh2EakI8fgQMHKp8IACA8\nJRwrXFbgY9ZH2I+yx1C9od88z+PxcOzYMezcuRM2NjbV2nQkCJYQGiJCAA8P4ORJ+nPKFLrkopkZ\nmw3MNDi+vsDWrXTC/OzZlUsEX/K/YMfjHfg7+G9s6bcFS7ovgYyUzDfnBAcHY/78+WjSpAm8vb1h\nYGAgvg8hIiwhNCTp6bRB1N6e7pq1eDGtGSgo1HRkDFOt4uOBq1dpF1lsLLBtG/DgAZ068yM8Pg8X\nAi9gq+dWjGk/BqFLQtFCvsV35wUHB2PIkCHYvXs35s6dW6NDSauCjTKq74p33rC3p9svjRpFE0Hv\n3qw2wDQohNAusUuX6KA5S0v65/DTT/T7UUV84n2w/MFyNJJuBLvhduiq0fUH70WQnp4OVVVVEX6C\nymGL2zHf43KBf/6hc+rT04FFi2j/QPPmFb+WYeoRPp/+KezaRQfQzZtHm4bU1Cr3+k+cT1jvvh6e\n0Z6wHWKLacbTanVfABt2yvwnPx/4+286Xq55c+C334DRo9mSEUyDw+fTSvH27bQG8McfwMiRlf9T\nKOAW4PDzwzjgcwALui5AxNIINJVt+s05hBC8e/euTvQPVAZLCPVFbi5w/Ph/i8adPg0MGMCahZgG\nycMDWLuWNhPt3k2bhir7p0AIwb1397DSdSU6tegEv3l+0FPR++68mJgYLF68GF++fIG3t3ed6Sf4\nkbr/CRq6oiI6Wqh9e8DPj64s6uwMWFiwZMA0OCEhtBYwbx6dUfzqFa0gV/ZPISI1AiOujMA693U4\nPvI4HK0cv0sGXC4Xhw4dQrdu3dC/f388efKkXiQDgNUQ6i5CaH143Tq6VaSjI9C1/E4uhqnPPn2i\nraNOTsCmTXRuZVXmf2XmZ+L3J7/jYuBFbO63GUvNl343jBQAwsLCMHPmTDRr1gy+vr7Q19cX4aeo\neSwh1EWBgcDKlUBqKnDiBF1dlGEaIA6Hjps4fhz49Vfg3buqLa/FJ3z8FfgXNntsxij9UQhdEgr1\npurlni8lJYWlS5di1qxZtbpjWVBslFFdkpRE9xW4dw/YsYP+BbB1hZgGiMsFzp6lfwZDhtAO46pu\nIfD843Msf7Ac0pLSsBthh24a3cQTbDVjo4zqu/x84PBh2mE8dy7diYxtM8k0QITQZqH16+muZPfv\n0zEUVZHIScSGRxvg/sEdewfvxfTO0yEpUT/6AITFEkJtd+sWsHo1/V/v5wfofT/agWEaAn9/YM0a\nICODfj8aPrxq4yYKuAU44ncE+7z3YZ7ZPERYR0Ch0fez9AkhuHz5Mry9vXHq1CkRfoLajyWE2io5\nmU4mi4igy0tYWNR0RAxTIz58oB3FT58Cv/9OJ5VVdR+m++/uw8bVBh3UOsD3V1/oq5bdGRwVFYWF\nCxciPT0dZ86cET74OobVk2qjO3fo7mMdOtAOZJYMmAYoLY2OnejeHejUiXYY//pr1ZLB29S3GHll\nJFa5rYLdcDvcnXq3zGRQVFQEW1tb9OjRA8OHD4e/vz+6NsBRe6yGUJtkZgIrVtCtKW/eBPr0qemI\nGKba5ecDR4/SyfaWlkBYGKBe/sCfMmUVZGHnk524EHABG/tuxB2rO5CVki33/GPHjuHRo0fw9/dH\n2wa8FzgbZVRbeHrStYaGD6eLsjdtWvFrGKYe4fOB//0P2LyZVpD37qWV5CqVQfi4FHQJmx5twvB2\nw7F78G60bFrx3h5cLhdSUlL1YigpG2VUl+Xl0b+A69fpOLoRI2o6IoapVnl5NBEcP05HUV++DPTr\nV/Vy/BP8sezBMkhAAnes7sBc07zSr5Vmw7cBsIRQs0JD6eY0hoZAcDBQA0vlMkxNCggApk4FdHWB\nDRuAn3+u+oorSdlJ2PhoI1wjXbFn8B7MMJlR7jDS5ORkJCQkwKyqY1UbCNapXFPOn6edxatX09oB\nSwZMA5KbS1ch/eknuluZiwvtL6hKMijkFeKAzwEYnTBCc7nmiFgagVldZpWZDAghOH/+PIyNjeHl\n5SWyz1HfsBpCdcvPB5YtA7y9gSdPgI4dazoihqk2eXl0r6Z9+2iz0OvXgLZ21ct58P4BbFxt0E6l\nHXx+9UF71fblnvvu3TssXLgQ2dnZcHNzQ5cuXYT4BPUbqyFUp9hYuiR1ejqdZcOSAdNAFBYCdnZ0\nXuXTp7RG4OBQ9WQQkRqBUVdHYYXLChz66RDuT7v/w2Rgb2+P3r17Y9y4cXj+/DlLBhUQS0JwcXEZ\n3qFDhwh9ff33tra260s/n5qaqjZ8+HCXLl26BBoZGYX89ddfs8URR63i4gKYmwOTJ9MhpWwUEdMA\nFC/K26kT/RO4f59OvjcxqVo5GXkZsHGxQd/zfTFIdxBCloRgVPtRFb7OxMQEr169go2NDaSqOput\nISKEiPTgcrlSenp6kdHR0bqFhYUyJiYmgWFhYR1LnrNt27btGzZs2EMIQUpKipqKikpaUVGRdMlz\naGj1AJ9PyP79hGhoEPL0aU1HwzDVgs8nxN2dkJ49CTE2JsTVVbByinhF5JjfMdJifwuy8O5Ckpyd\nLNpA66F/r50CXb9F3ofg7+9v3q5du0hdXd0YALCysrrm6Og4rmPHjuHF57Rq1SoxODi4MwBkZWUp\nqqqqpklLS3NFHUuNKyigy08EBQHPnwvWWMowdYy3N12UNyGBrkY6ZYpgO7g+jHqIla4r0UK+BR7O\neIjO6p1/eD6fz683G9XUFJEnhISEBE1tbe344vtaWlof/fz8epQ8Z/78+WcGDRrkoaGh8YnD4Sg4\nODhMLqus7du3f71tYWEBi7q0hENGBjBuHB099PQpIC9f0xExjFi9ekUTQXg4sG0bMGOGYKuzv097\nj9VuqxGaEooDQw9gfIfxP5wwlpiYiGXLlqF3795YtWqVEJ+gbvLy8hLZyCmRJwQJCYkKpxfv3r17\nU5cuXQK9vLwsoqKi9IYOHfowKCjIREFBgVPyvJIJoU5JTgYGD6Yb1xw8yDa4Z+o1Hx+6SY2/P51j\n6egIyJa/SkS5vuR/wc4nO3Ex8CLW9VmHG5Y30Ei6/G3P+Hw+zp49i82bN2PBggVYvHixEJ+i7ir9\nZXnHjh0ClyXyK5WmpmZCfHz817aR+Ph4bS0trY8lz/Hx8eltaWl5AwD09PSi2rRpE/327VsDUcdS\nIz5/BgYNojNsDh1iyYCpt169AsaPB6ys6CY1kZHAkiVVTwY8Pg+nXp5Ch2MdkJmfiZAlIVjXZ90P\nk0FERAQsLCxw/vx5eHh4YNeuXWjSpImQn4gR+dWqW7duL9+/f68fExOjW1hYKHv9+vUpY8eOdSp5\nTocOHSLc3d2HAEBycrL627dvDdq2bftB1LFUu5QUWjP4+Wc666YerIvCMKVFR9PZxWPG0O8+ERGA\ntTUgyPXYI9oDZqfNcDXkKh5Mf4CzY89Wau2hvXv3wtLSEt7e3jA2NhbgUzBlErQ3+keHs7PziPbt\n27/V09OL3L1790ZCCOzt7Rfa29svJP+OLBo9evTdzp07BxkZGb25cuXKtNJloK6NMkpJocMptmyh\nQywYpp5JSyNk1SpCVFQI2b6dEA5H8LIi0yLJ+Gvjie6fuuRG6A3CZ38zIgMhRhmx1U5FIS2N1gxG\njgR27WI1A6Zeyc8Hjh0DbG2BSZNo5bdlxV/iy5RVkIVdT3fh7OuzWN1rNVb1WoXG0o1FGm9DJ8xq\np6yBW1jp6bQBdfhwlgyYeoXPB65coUtQP31KV1qxtxcsGfD4PJx9fRYGxwzwOeczQhaHYFO/TRUm\ng1u3biEyMlLAT8BUFVvLSBgZGXQk0ZAhwJ49LBkw9YaHB7B2LR02eukS0L+/4GU9iX0CGxcbyMnI\n4e7Uu+im0a3C1yQkJGDp0qWIiIjAlStXBH9zpkpYDUFQWVl0qUYLC7pSF0sGTD0QEkJbPufPB9av\np/MpBU0G0RnRsLxhiRm3Z2B9n/V4OudphcmAz+fjxIkT6NKlC0xMTBAYGMiWqq5GrIYgiIICOt6u\nWze6uxlLBkwdl5BAJ5M5OdG5BLdvA43KH/X5Q5wCDvY824NTr07BpocNLo2/hCYyFQ9BIoRg6NCh\nKCgowOPHj2FoaChYAIzAWEKoKkLo1yclJdrTxpIBU4dxOLSCe+IEMG8e3cheSUmwsoq3r9zssRmD\n2wxG0KIgaClqVfr1EhISOHjwIDp37syWoKghLCFU1e+/07n5Xl4AWz2RqaOKioAzZ4CdO2nLZ0AA\noKMjeHnP4p7BxsUGMlIyuDX5Fnpo9aj4RWVgy1PXLJYQqsLBAbhwAfDzY2sTMXVS8XLUGzbQBODs\nDJiaCl5e7JdYrHdfD+94b+wdvBdTjaeWu31lSRwOB02bNq0Xm9rXJ6xeVlmvXtHpmHfuAOrqNR0N\nw1SZry/dpWzbNrpZjZub4MkgpzAHWz23wuy0GQzUDBBhHYHpnadXmAwIIXBwcICBgQECAwMFe3NG\nbFgNoTISE4EJE4BTpwBWpWXqmMhIYONGOmJo5066CqmgrZ18wseV4CvY+Ggj+rfuj4CFAdBpVrm2\npri4OFhbWyM6Oho3b96EqTBVE0YsKl1DyM3NlRNnILVWURFdm2jePGDixJqOhmEqLTUVWL4c6NkT\nMDMD3r4FZs8WPBk8//gcvc71gp2/HRwsHXB10tVKJQMejwc7OzuYmZmhR48eeP36NXr37i1YEIxY\nVZgQfHx8ehsaGoYZGBi8BYDAwMAuS5YsOSH+0GqJtWsBFRW60DvD1AF5eXSeZIcO9H54OK0hyAn4\nlS4+Mx7Tb03HJIdJsO5uDb95fuitXfkLelFREQICAuDt7Y0tW7ZAVpC1sZnqUdFiR927d/ePjY3V\n6dKlS0DxY4aGhqGCLp5U2QO1YXG727cJ0dUlJD29piNhmApxuYRcuECIlhYhP/9MyLt3wpWXU5hD\ntntuJyqdpR8bAAAgAElEQVS2KmTzo82EUyDEanZMtYG4t9DU0dGJK3m/Xm53WdqHD8CCBcDdu4Cy\nck1HwzA/FBoK/PorbQ5ycAB69RK8LEII/hfyP2xw34Be2r3wasEr6CrpiixWpvaqsMlIR0cnztvb\nuw8AFBYWyh44cGBNyf2R66WiImDaNFrP7iHYeGqGqQ5cLrB6NTBwIPDLL3QROmGSwYuEF+hzvg8O\n+h7ElYlXcP3n65VOBmlpaVi9ejWys7MFD4CpURUmhJMnTy4+fvy4dUJCgqampmZCQECA6fHjx62r\nI7gas20b3QvZxqamI2GYcnl60s7iN29oDWHpUsE36EvISsCsO7Mw7to4zDObB/95/ujXul+lXksI\nwdWrV2FkZAQut/43HtRrFbUpPXv2rE9lHhP1gZrqQ3j0iBANDUKSk2vm/RmmArGxhEydSoi2NiHX\nrxNSWCh4WbmFuWTn451ExVaFbHDfQLLys6r0+ujoaDJ8+HBibGxM/Pz8BA+EERkI0YdQ4feJpUuX\nHqvMY/UChwPMnQucOwe0aFHT0TDMN7Kzgd9+o5PJ9PXp6KHJkwEZmaqXRQiBQ6gDOh7viICkALyY\n/wJ7Bu+BQiOFSpcRHx+Pbt26oX///nj16hXMzc2rHghTq5Tbqezr69vLx8end0pKSvNDhw6tIv/u\nwMPhcBT4fH79nOG8bt1/m90wTC3B5wN//w1s2kRXWw8MBLS1BS/v1adXsHG1AaeAg7/G/wULXQuB\nytHW1kZoaCjU2cz9eqPchFBYWCjL4XAUeDyeFIfD+fq1QVFRMevmzZs/V0941cjTE7h3jy4IzzC1\nhL8/nVxGCHDzJp1kJqjPOZ+x8dFG3H93HzsH7sRc07mQkhRugUaWDOqXCvdUjomJ0dXV1Y2pnnD+\nU617KufkAJ070wVeRo2qnvdkmB9ITKSD3Nzc6CSzGTME7zDm8rk48eIEdj7ZiRmdZ2DbgG1o1rhZ\nlcr48OED2rZtK1gATLUSZk/lCuchyMnJ5a5Zs+ZAWFiYYV5eXpN/35B4eHgMEuQNa6VNm4A+fVgy\nYGpcQQFw5Ajdo+DXX4GICEBRUfDyHsc8xrIHy9Bcvjm8ZnmhU4tOVXp9SkoKVq1aBT8/P7x58waN\nBN01h6kTKvzOMX369CsdOnSI+PDhQ9vt27dv19XVjenWrdvL6giuWjx7Bty4Afz5Z01HwjRghNA5\nkEZGdC6Bry9gayt4MkjISsC0f6Zhxu0Z2Np/K9xnuFcpGRBCcOnSJRgZGUFdXR0BAQEsGTQEFQ1D\nMjU1fU0IgbGxcXDxY127dn0p6LCmyh6ojmGnhYWEGBsT4uAg/vdimHKEhxMybBghBgaEPHggXFkF\n3AKy9+leomqrSjY/2kyyC7KrXEZsbCwZMmQIMTU1JS9fvhQuIKbaQZxLV8jKyhYCQMuWLZPu3bs3\nWkND41NGRkb9WMvhzz+BVq3oaqYMU82+fKEb8F2+TFstly4VbAhpMddIVyx3WY52Ku3wfN5ztFNp\nJ1A5UlJSGDFiBJYvXw5pabZCfoNSUcZwcnIak5GRoRQcHGw8YMAAL1NT09eOjo5jBc1AlT0g7hpC\nUhIhqqqEvH8v3vdhmFK4XELOnCGkZUtC5s0Tfg7kh/QPZPy18aTtkbbk7tu7ogmSqbMgRA2hwlFG\nZfH39zc3Nzf3F3VyKknso4yWLAEaNQIOHxbfezBMKd7edBhpkya087hrV8HLyivKwz7vfbDzt8PK\nniuxpvcaNJZuLLpgmTpJLKOM+Hy+5O3btydERUXpGRkZhYwcOdL55cuX3TZt2rT78+fPLQIDA+vu\n1mFv39IlISMiajoSpoH4+JHOe3z6lHYWT50KCLqdMCEETm+dYONqg24a3aq0a1lJzs7OuHbtGi5e\nvMj2NmYA/GAewrx5885GR0e3MTc393/8+PGAVq1aJUZERHTYtWvX5nHjxjlKSEiIdZKAWGsIY8YA\n/fvTzW8YRozy8oCDB2lFdMkSurm9vLzg5b1NfYsVLisQmxmLoyOOYkjbIVUuIzk5GStWrMCLFy9g\nb2+PoUOHCh4QU+sIU0Moty2pU6dOITweT5IQgry8vMbNmjX7kpqaqipo21RVD4irD+HhQ0LatiUk\nP1885TMMIYTPJ+Sff+j+ShMnEvLhg3DlcQo4ZP3D9UTVVpUc8D5ACrgFAsTEJ2fPniXNmzcn69ev\nJzk5OcIFxdRKEMcoIxkZmSJJSUk+ADRu3Di/TZs20aqqqmkCZZ3agscD1qyhdXY2ppoRkzdv6Mrp\nnz/TdRIHCTGFkxCC66HXsfbhWgzUHYg3i9+glUIrgcq6evUq7O3t4ebmhi5d6m6LLyM+5TYZNWnS\nJK9du3aRxfejoqL09PT0ogDanBMcHNxZrIGJo8no/Hngr7+Ax48Fb8BlmHKkp9PVSB0c6JYaCxcC\nwozafJP8BsseLMOX/C84NvIY+ur0FSo+LpcLCQkJSEkJt34RU7uJpVM5PDy8o+Ah1UI8HrB7N3Dh\nAksGjEhxucDp08COHXRKS3g43V9JUFkFWfjN8zdcfXMV2y22Y2HXhUIvQgeAzSlgKlTu/5CaWNBO\nrG7cANTVgb7CfctimJK8vOgwUlVV4OFDukaiMDyjPTHHcQ4Gtx2MMOswqMmpVbkMDoeDiIgIdO/e\nXbhgmAZHLPsauLi4DO/QoUOEvr7+e1tb2/VlnePl5WVhamoaYGRkFGJhYeEljji+4vOBXbvodFBW\nO2BEIDYWsLQE5syhzUMeHsIlg/CUcIz931j8cvsXnBh1AufGnhMoGTg5OaFTp064ceOG4MEwDZeg\nvdHlHVwuV0pPTy8yOjpat7CwUMbExCQwLCysY8lzMjIylAwNDUPj4+O1CCFISUlRK10ORDnK6J9/\nCOnalQ79YBghZGcTsnUrneT++++E5OYKVx6ngEPWPVxH1PapkcO+h0lGXoZA5Xz69IlMmjSJ6Ovr\nEw8PD+GCYuo0iHMLTQDIzc2Ve/v2rUFlzvX39zdv165dpK6uboyMjEyRlZXVNUdHx3Elz7l69eq0\nSZMm/aOlpfURANTU1FKrmMcqjxBg507a28dqB4yACAGuXQM6dgTevwcCAoCtW+mMY8HKI7gRegMd\nj3dEIicRIYtDYNPTBkqNlapc1s2bN9G5c2cYGBggKCgIAwcOFCwopsGrsJfJyclp7Nq1a/cXFBQ0\niomJ0Q0ICDDdtm3bDicnp7FlnZ+QkKCpra0dX3xfS0vro5+fX4+S57x//16/qKhIZuDAgZ4cDkdh\nxYoVR2bMmHG5dFnbt2//etvCwgIWFhZV+Gj/cnOjvX5jxlT9tQwDevFfsYJuuX3lCtCvn3DlvU19\ni2UPliExOxFXJl5B/9b9hSqvTZs28PDwgLGxsXCBMXWSl5cXvLy8RFNYRVUIU1PT1xkZGUpdunQJ\nKH6sU6dOIeWdf/PmzUnz5s07U3z/8uXLvyxduvRoyXOsra2P9erVyyc3N7dJamqqqr6+/rt3797p\nlzwHomoyGjiQkIsXRVMW06B8/kzIggWEqKsTcvo0XZROGNkF2WSj+8avzUOF3ELRBMowJUCcy1/L\nyMgUKSkpfSn5WPGEtbJoamomxMfHf90CPD4+Xru4aaiYtrZ2vJqaWmqTJk3ymjRpkte/f/8nQUFB\nJvr6+u+rnNF+5M0bul6RlZVIi2Xqt6Ii4PhxOg7hl1/ofyGlqrfkfEUIwZ2IO7BxtUFfnb4IXhQs\n8OQyQghbd4gRmwr7EDp16hR65cqV6VwuV/r9+/f6y5YtO9q7d2+f8s7v1q3by/fv3+vHxMToFhYW\nyl6/fn3K2LFjnUqeM27cOMdnz5715fF4Urm5uXJ+fn49DA0Nw0Txgb5x5AhgbQ3Iyoq8aKZ+cnMD\nTEyABw+AJ0/oGkTCJIP3ae8x8upIbPHcgovjL+LKxCsCJYPMzEwsXrwYmzdvFjwYhqlIRVWI7Oxs\n+Y0bN+7u2rXry65du77ctGnTrry8vMY/eo2zs/OI9u3bv9XT04vcvXv3RkII7O3tF9rb2y8sPmf/\n/v1rDA0NQ42MjN4cOXJkeekyIGyTUUYGIUpKhCQmClcO0yBERhIydiwhenqEODoKPyAtpzCHbPXY\nSlRtVcl+7/1CNQ/9888/RFNTkyxYsICkp6cLFxhT70Gc+yG8fv3azMzM7HV1JKeShF66Yvdu4N07\nulQFw5SDw6H/Vc6coYvf2tgIv8yV01snrHBZgR6aPXDgpwPQUtQSqJyEhAQsXboU4eHhOH36NPr3\nF67zmWkYhFm6osKEYGFh4ZWUlNTS0tLyxpQpU64bGRmFCBRlVQMTJiFwuUDbtoCTE8AW8WLKwOcD\nf/8NbNwIDBkC7NkDaGgIV+aHjA9Y4bIC79Pe49jIYwItTV3S8uXLoaKigo0bN7IN7plKE2tCAIDE\nxMRWDg4Okx0cHCZnZWUpTp482WHr1q07BXnDSgcmTEK4eZP2Hzx9KtqgmHrB358uN0EI/W/Ss6dw\n5eVz82H7zBZH/Y9ibe+1WNlrJWSlhO+3IqwDmRGA2BNCsTdv3hjb2tquv379+pSioiIhtgOvmFAJ\nYfhwOjzkl19EGxRTpyUl0RqBqyutEcyYAUgKuXiL83tnLH+wHF1adsGhYYcE2rmMYURJLKudFgsL\nCzN0cHCYfPPmzZ9VVVXTpkyZcv3QoUOrBHmzapGQQL8C3r5d05EwtURBAWBnR7fB+PVXOoxUUVG4\nMmO+xMDGxQZhKWE4PvI4hrUbJnBZjx8/hrKyMjoLuzIewwipwoQwd+7c81ZWVtdcXV2HaWpqJlRH\nUEK5fBmYNEnwNQWYeiU0lO5frKMD+PoC+vrClVfALcABnwM4/PwwVvZcies/X0cjacHa9zMyMrBu\n3Tq4uLjg0qVLwgXGMCJQYUJ4/vy5kC2s1YjPp8NFrl2r6UiYGlZUBPz5J7BvH7B3LzB3rvBLWblG\numLZg2Xo1KITXi54CV0lXYHKIYTgxo0bsLGxwcSJExEaGgpFYassDCMC5SYES0vLGzdu3LA0NjZ+\nU/q56tgxTSAeHrQtoFu3mo6EqSGEAM7OdKdULS3Az48OOBNGXGYcVrmuQkBSAOyG22FU+1FClTdz\n5kwEBATg5s2b6N27t3DBMYwIldup/OnTJw0NDY1PsbGxrUt3UEhISJDWrVvHijUwQTqVJ08GBg4E\nFi8WT1BMrfb2LbB0KRAfDxw8CIwcKVytoJBXiEO+h3DA5wCWmS/D+r7r0Vi6sdBxvnr1CsbGxpBl\nM+gZMRCmU7nCmWvr1q2zrcxjoj5Q1ZnKGRmEKCoS8uVL1V7H1HlcLiGHDtE9Cg4fJqRQBGvGPYx6\nSAyOGpBRV0aRqPQo4QtkmGoCce6H4Obm9lPpx5ydnUcKlH3EydkZGDAAaNaspiNhqtH79/Sf/fZt\n2jxkYwPICDEg+mPWR0y5OQXz787H/qH7cW/aPbRVFqzNKS8vD3x+uetAMkytU25COHny5GJjY+M3\nb9++NTA2Nn5TfOjq6sZ07tw5uDqDrJQ7d4Bx4yo+j6kX+Hw6lLRXL7qVpZcXoKcneHlFvCIc8DmA\nLvZdYKBqgNAloRhjIPgeGo8ePYKxsTHc3d0FD4phqlm5fQiZmZnNMjIylDds2LDX1tZ2Pfm3TUpB\nQYGjqqqaJvbAqtKHUFAAqKvTtYtatBBvYEyNi4yk8wl4PODCBeGHknpGe2Lpg6XQVtTG0RFHoa8q\neIFpaWlYvXo1PD09cfz4cYwePVq44BimioTpQyi3hiAhIUF0dXVjjh8/bq2goMBRVFTMUlRUzJKQ\nkCDp6ekqgocrBp6egJERSwb1XFoaMHMmXWpi3Djg8WPhksEnzidM+2caZjvOxh8D/8CD6Q8ETgaE\nEFy9ehVGRkZo1qwZQkJCWDJg6pxyh51OnTr1f/fv3x/VtWvXVxISEt99VY+Ojm4j3tCq4PZt1lxU\nz7m60lrBlCm030BZWfCyinhFOOZ/DLue7sLCbgtxZswZyMvKCxUfn8+Hm5sbHB0dYW5uLlRZDFNT\nqrSWUXWqdJMRIXTAuacn0L69+ANjqtWHD8Dq1UBQEGBvD/z03RCHqnka+xTWztZQb6qOYyOOwUDN\nQDSBMkwtIZYmo2Le3t59srOzmwLA5cuXZ6xatepQbGxsa0HeTCyCg+kyFSwZ1Cs5OcCWLYC5OdC9\nOxAWJlwySMpOwszbMzHt1jRs7b8Vbr+4sWTAMKVUmBAWLVpkLycnlxsUFGRy6NChVW3btv0wc+bM\n2rPwirs7MHRoTUfBiNCTJ0DnzrR2EBgIbNoENBZwPhiXz4Wdnx2MTxqjlUIrhFuHw7KTpcDLSufm\n5mLr1q1ITU0VLCCGqcUqTAjS0tJcSUlJ/p07d8ZbW1sfX7p06TEOh6NQHcFVirs7MHhwTUfBiEBq\nKl1zaPp0OtP46lXaGigon3gfdDvdDXci7uDJ7CewHWKLprJNBS7Pzc0NRkZG+PDhg+BBMUwtVuHi\ndgoKCpzdu3dv+vvvv395+vRpPx6PJyXuvRAqraAAePYMuHKlpiNhhMDn051ON26kK5OGhQEKQnzl\n+JzzGRvcN8Atyg0HfjqAKZ2mCLXRTEpKClatWoVnz57hxIkTGDFihODBMUwtVmEN4fr161MaNWpU\ncP78+bktW7ZMSkhI0Fy7du3+6giuQs+e0eGmKrVrFCxTeQ8e0Mllp07R23/+KXgy4PF5OPHiBIxO\nGEG5iTLCrMNgZWQlVDLIzMyEiYkJWrRogZCQEJYMmHqtUqOMkpKSWr548aK7hIQEMTc392/RosVn\nsQdWmVFGGzfSdQp+/13c4TAilptLVyR1dqZLVP/8s3C7l/l99MMS5yVoKtsUx0ceh1ELI5HF+vHj\nR2gJ03bFMNVIrKOMHBwcJvfo0cPvxo0blg4ODpPNzc39b9y4YSnIm4mcuzvdIZ2pU1xdaadxVhYd\nTjp5suDJIDU3FfPvzseE6xOwsudKeM3yEmkyAMCSAdNgVFhD6Ny5c7C7u/uQ4lpBSkpK88GDBz8S\n934IFdYQOBygZUsgPR1oJNiOVUz1Skigi8+9fg0cOwYI0/rC4/NwLuActnpuxVSjqdhhsQPNGgu3\nsGFcXBx0dNieyEzdJtYaAiFEonnz5inF91VVVdMEfTOR8vGhG+GwZFDrcbnA4cOAiQnQsSMQEiJc\nMniR8AK9zvXCpaBLcPvFDX8O/1OoZJCdnY2VK1eiZ8+eyMjIEDwwhqnjKhxlNHz4cJdhw4a5Tps2\n7SohROL69etTRowY8aA6gvuhJ0+Afv1qOgqmAj4+dL+i5s0Bb2/AQIi5YOl56dj0aBMc3zpi7+C9\nmGEyA5ISQnQ8AHB2dsaSJUswYMAABAcHQ1mYNTEYpo6rVKfyrVu3Jj579qwvAPTr1+/phAkTbos9\nsIqajLp2BWxtWR9CLZWWBmzYANy/T+cUWFkJvnsZn/BxIeACNntsxs+GP2PnwJ1QbiLchTs1NRVL\nly7FixcvYG9vj6FsciNTTwjTZFRuDeHdu3ft165duz8yMrJd586dg/fv379WS0vro+BhilBBAV2y\nonv3mo6EKaXknIIpU4DwcOH2LHqd+BrWztYAAOfpzjBrZSaSOCUlJWFgYIDz589DTk5OJGUyTF1X\nbg2hb9++z2bNmnWxX79+T+/evTvG19e3161btyZWW2A/qiH4+QGLFgEBAdUVDlMJb97Q5qHCQuDk\nSVqJE1RGXga2em7FjbAb2D1oN+aYzhG6eYhhGgKx1BCys7Obzp8//wwAdOjQIcLU1LT2XH39/IAe\nPWo6CuZf2dnAjh20ZvD778CCBYCUlGBl8Qkfl4IuYYP7BkzoOAHh1uFQacImHjJMdSg3IeTn5zd+\n/fq1GUBHGuXl5TV5/fq1GSFEQkJCgpiZmb2uvjBLef6cLWhXCxBCt6KwsQEsLOjoIXV1wcsLTg7G\nkvtLUMArwN2pd9FdU/gmwefPn+PYsWP466+/IC1d4RgKhmnQym0ysrCw8Cq5MU5xIii+7+npOVCs\ngf2oyahdO8DJCTA0FGcIzA98+AAsWwZERwMnTtCEIKisgixs89qGK8FXsHPgTswzmwcpSQGrGMVl\nZmVh8+bN+Oeff3D48GFMnjxZqCUsGKauEEuTkZeXl4XAEYlTZiaQlCTc+EVGYAUFwP79dF7BmjW0\nhiArK1hZhBBcC7mGNQ/XYHi74QhdEorm8s2FjtHJyQnW1tb46aefEBISAhW21hXDVErdq0MHBdEF\n7QRtpGYE9ugRYG1N9yJ69QrQ1RW8rPCUcFg7WyM9Lx03LG+gt3ZvkcTo7u6ONWvW4NKlSxg4UKyV\nWIapd+peQggMBExNazqKBiU9HVixAnj6FLCzA8aOFbys7MJs7HyyE+cDzmNr/61Y0n0JpCVF999w\n8ODBCA4ORmNBd9RhmAZMLOP4XFxchnfo0CFCX1//va2t7fryznvx4kV3aWlpbpWGswYF0TUQmGpx\n7x5gbEw3tQ8NFTwZEELwT9g/MDxuiISsBAQvCsbyHstFmgwA2n7KkgHDCIgQ8sODx+NJXrp0acaO\nHTt+I4QgNjZWx8/Pz7y887lcrpSenl5kdHS0bmFhoYyJiUlgWFhYx7LOGzhwoMeoUaPu3bx5c1Lp\n52loZTAzI8TXt+znGJHJyCBk1ixC2rQhxNNTuLLepb4jwy4PI4bHDYlntJCF/Ss/P5/4+/uLpCyG\nqU/+vXZWeG0v66iwhrBkyZITvr6+va5evToNAJo2bZq9ZMmSE+Wd7+/vb96uXbtIXV3dGBkZmSIr\nK6trjo6O40qfd/To0WU///zzzZIL51WoqIhOfe0s1oVWG7wHD2itQF6eTggXdARRXlEefvP8Db3O\n9cLgNoMRuDAQFroCFlbCs2fPYGpqiiNHjghdFsMw/6mwvu7n59cjICDAtHhimoqKSvqPttBMSEjQ\n1NbWji++r6Wl9dHPz69H6XMcHR3HeXh4DCreeKessrZv3/71toWFBSxUVYHWrQG21IBYZGYCq1fT\nzuOLF4FBgwQv6967e1j+YDm6anRF4KJAaCkKv6dAZmYmNmzYACcnJxw5cgSTJk0SukyGqeu8vLzg\n5eUlkrIqTAiysrKFPB7v65CelJSU5pKSkvzyzi/v4l6SjY3Nn3v37t3w71wDCVLOmNmSCQEAcOkS\n0KVLRcUzAnj4EJg3jy5LHRws+DaWMV9isMJlBcJTwmE/2h4/6f0kkvg8PDwwc+ZMjBo1CqGhoVBS\nUhJJuQxT11lYWMCiRDV+x44dApdVYUJYtmzZ0QkTJtz+/Plzi02bNu2+efPmz3/88ceW8s7X1NRM\niI+P1y6+Hx8fr116UbxXr151tbKyugYAqampag8ePBghIyNTNHbsWKcfBhMYyBKCiHE4wNq1tJno\nzBngJwGv33lFedjvsx9H/I5gZc+VcPjZAY2kRbdXhbq6Oq5evYr+/fuLrEyGYb5VqeWvw8PDOz56\n9GgwAAwePPhRx44dw8s7l8vlShsYGLx99OjRYA0NjU/m5ub+//vf/6aW95o5c+ZcGDNmzN2JEyfe\n+iawsmYqDxsGLF8OjBpViY/GVMTLiy5G17s3cOhQ1VclTchKgEe0BzxiPOAa6Yre2r1x8KeDaK3U\nWizxMgxTMbHMVC4WFxenIy8vnzNmzJi7/74ZiYuL09HR0Ykrs0Bpae6xY8eWDhs2zJXH40n9+uuv\n5zp27Bh+6tSphQCwcOHCU4IECgAICwM6dRL45QyVkQGsXEkTwt69dK+CqkjkJGK122q4RrlioO5A\nDG4zGBv7bkR71fYiiY8QwpaZYJgaUGENwcjIKKS4XyA/P79xdHR0GwMDg7ehoaFivTJ/V0PIzAQ0\nNenO7ILuyM7g7l1aKxg/niaDpk0r/1oen4cTL07g9ye/Y57ZPGzptwXysvIiiy0/Px+7du1Camoq\nTp48KbJyGaYhEWsNISQkxKjk/devX5sdP37cWpA3E0pgIKCnx5KBgIpnG/v4AH//XfWhpP4J/lh0\nbxEUGyni8ezHMGwu2oUFHz9+jAULFsDIyAh2dnYiLZthmMqp8jRRMzOz16WHkVaLDx8Aff1qf9v6\n4M4dYMkSwNKSjiCSr8KX+oy8DGzy2IQ7EXewb8g+/NL5F5E252RkZGDdunVwcXHB0aNHMX78eJGV\nzTBM1VSYEA4ePLi6+Dafz5d8/fq1maamZoJ4wyrDu3dsQloVpabSJapfvgSuXwf69av8awkhuBx8\nGesersOEjhMQtiRM6H2My3L48GE0atQIoaGhUFRUFHn5DMNUXoUJITs7+2srs7S0NHf06NH3Jk2a\n9I94wypDVBQwYUK1v21ddfMmTQbTptHln6oyly8sJQxL7i8Bp5ADp6lOMNc0F1ucO3bsYB3IDFNL\n/DAh8Hg8qaysLMWStYQa8+ED0LZtTUdR66Wl0U7joCDgn3/okNLKyinMwc4nO3Eu4By2DdiGxd0W\nC71RTUVYMmCY2qPcHloulystJSXF8/b27iNoj7VIxcYKtwB/PUcI4OBA1yDS0qJ98FVJBo4RjjA8\nYYi4zDgELwrGUvOlIk0GwcHB8Pf3F1l5DMOIXrk1BHNzc//Xr1+bdenSJXDcuHGOlpaWN+Tk5HIB\nOiS09EQyscrJoTu5t2hRbW9Zl8TF0U7jmBjaVFSVRBDzJQbLHyzH27S3ODf2HIa0HSLS2PLy8rBz\n506cPXsWJ06cgLm5+JqfGIYRTrkJobhWkJ+f31hVVTXNw8Pjm6XOqjUhxMTQRe1Y88I3eDzg2DFg\n50460ezWrcpvZ1nIK8Qh30M44HMANj1tcMPyhkiXmgCAR48eYeHChejatSuCg4PRsmVLkZbPMIxo\nlZsQUlJSmh86dGiVsbHxm+oMqEwxMay5qJSgIGD+fDqE1MeHbmtZWV4xXlhyfwnaKLeB/3x/tFUW\nfd/MunXrcP36dRw/fhyjR48WefkMw4heuQmBx+NJcTgcAde8FLHYWFpDYJCbC/z+O3D+PJ1pPGdO\n5Z8391oAACAASURBVCtOydnJWPNwDR7HPMafw//EhA4TxNapa2Vlha1bt0JB0GVTGYapduUmhJYt\nWyZt27ZN8HVURYklBAB0iepFi4AePYA3bwB19cq9jsfn4fSr0/jN6zfM7jIbYdZhaCpbhTUrBGBm\nZibW8hmGET3RbmgrLrGxwJgxNR1FjUlNBVatAp48AU6epHsWVNarT6+w+P5iyErJwmOmB4zVjUUa\nG5fLBSEEMjLl7pnEMEwdUe6wU3d3d9EONxFGXFyDrCEQAly+DBgZAc2bAyEhlU8GmfmZWPZgGUZe\nHYnF3RbjyZwnIk8GAQEB6NmzJ65duybSchmGqRnl1hBUVVXTqjOQH/r4kQ6ub0CiomjzUFoacP8+\n0LVr5V5HCMH/Qv6HNW5rMLr9aIQtCYOqnKpIY8vNzcW2bdtw6dIl2Nra4pdffhFp+QzD1IxKbZBT\nE74uf83nA02a0GWvG4l2WGR1UFFRQUZGRk2HwTBMPaOsrIz09PTvHhfr8tc1LiODjq2sg8kAoKt5\n1takyzBM3SWOEYK1f3OBpCSATWhiGIYRu9qfEBITKz++kmEYhhFY7U8Inz4BGho1HQXDMEy9V/sT\nAmsyYhiGqRa1PyF8/swSAlNnPH36FB06dKjUuV5eXtDW1hZbLOIu/0cKCgrQqVMnJCcn18j71yXJ\nyckwNDREYWFhTYdSBxJCdDSdlcWIha6uLuTk5KCgoICWLVtixowZyMrK+uYcHx8fDBo0CIqKilBS\nUsLYsWMRHh7+zTlZWVmwsbFB69atoaCggHbt2mHlypVIS6s901mqQ79+/RARESGSsmbPno2tW7eK\npKzqdvr0aQwYMADqdbj/r6CgAHPnzkWzZs3QqlUrHD58+Ifnp6SkYNq0aVBSUoKKiso383N+VJa6\nujoGDhyI06dPi+2zVFbtTwgcDtBUvOvuNGQSEhK4d+8eOBwOgoKC8ObNG/zxxx9fn/f19cWwYcMw\nYcIEJCYmIjo6GiYmJujTpw+io6MBAIWFhRg8eDDCw8Ph6uoKDocDX19fqKmpiXVTHC6XK7ayGeGc\nOnUKM2bMEOi1teXfdfv27YiKikJcXBw8PT2xb98+uLq6lnv+xIkToaGhgfj4eKSkpGDNmjWVLmv6\n9Ok4deqUWD9PpRBCauVBQyOEmJoS8uIFqWsyMghZtIiQr5+jltLV1SWPHj36en/t2rVk5MiRX+/3\n7duXWFtbf/e6ESNGkJkzZxJCCDlz5gxRV1cnOTk5lX7fkJAQMmTIEKKiokLU1dXJnj17CCGEzJo1\ni2zZsuXreZ6enkRLS+vr/datWxNbW1tibGxMGjVqRGxtbcnPP//8TdnLly8ny5cvJ4QQ8uXLFzJ3\n7lzSqlUroqmpSbZs2UJ4PN538eTl5ZHGjRuTtLQ0Qgghf/zxB5GWliYcDocQQsiWLVuIjY0NIYSQ\n/Px8snr1aqKjo0PU1dXJokWLSF5eXpnxvnr1inTp0oUoKCgQS0tLMnny5K+fr/jcgwcPkhYtWpBW\nrVqRCxcuEEIIOXXqFJGRkSGysrKkadOmZOzYsYQQQhISEsjEiRNJ8+bNSZs2bYidnd3X98rNzSWz\nZs0iysrKxNDQkOzbt++bWEpbvnw50dbWJoqKiqRr167k6dOnX5/btm0bmTRpEpkyZQpRUFAgZmZm\nJCgo6Jt/hz179hBDQ0OirKxM5syZQ/Lz8wkhhMTGxpImTZp883u+d+8e6dKlC1FUVCTa2tpk+/bt\nX5+Ljo4mEhIS5Ny5c0RHR4cMGDCAEELIuXPnSMeOHYmysjIZNmwYiY2NrVTsoqKhoUEePnz49f5v\nv/1GrKysyjzX1dWV6Orqlvl/qzJlFRUVETk5ORIXF1fp+Mq7tvz7uEDX3dpfQ0hOrlPDTgmhexl3\n6kRv1wXk30A/fvwIFxcX9OjRAwBdosLX1xeWlpbfvWby5Ml4+PAhAMDd3R0jRoyAnJxcpd6Pw+Fg\nyJAhGDlyJBITExEZGYnBgwcDoDWWiibcXLt2DQ8ePEBmZiasrKzg7OyM7OxsAACPx8ONGzcwffp0\nALTZRVZWFlFRUQgICICbmxvOnj37XZmNGzeGubk5vLy8AACPHz+Grq4unj179vW+hYUFAGDDhg2I\njIxEUFAQIiMjkZCQgN9///27MgsLCzFhwgTMnTsXGRkZmDp1Ku7cufPN50tKSkJWVhY+ffqEc+fO\nwdraGpmZmViwYAGmT5+O9evXg8PhwNHREXw+H2PGjIGpqSk+ffqER48e4c8//4SbmxsAYMeOHYiO\njsaHDx/g6uqKixcv/vB3aW5ujqCgIGRkZGDatGmwtLT8ph3byckJkydP/vr8+PHjwePxvj5/9epV\nuLm5ISoqCu/evftas3zz5g3atm0LScn/Li9NmzbF33//jczMTNy/fx8nT56Eo6PjN/E8efIEERER\ncHFxgaOjI/bs2YPbt28jNTUV/fr1w9SpUysde0l79+6FsrJymYeKikqZr8nIyEBiYiJMTEy+Pta5\nc2eEhoaWef7z589hYGCAWbNmQU1NDebm5njy5Emly5KWlka7du0QGBhYZvnVRtBMIu4DACF8PiEy\nMoT8++2rtouLI2TsWEI6diSk+AsLKlFDoKlDuENQrVu3Jk2bNiUKCgpEQkKCjB8//uu3nPj4eCIh\nIUHevn373esePHhAZGRkCCGEDBkyhGzcuLHS73n16lViZmZW5nOzZ8/+YQ1BV1f367foYn379iWX\nLl0ihBDi5uZG9PT0CCGEJCUlkUaNGn399l783gMHDizzvbdu3UqWL19OuFwuadmyJbGzsyMbNmwg\neXl5pEmTJiQ9PZ3w+XwiLy9PoqKivr7Ox8eHtGnT5rt4Hz9+TDQ1Nb+LdevWrV/PLf1NukWLFsTP\nz6/M38Xz58+Jjo7ON+Xt3r2bzJkzhxBCSNu2bYmrq+vX506f/n975x3W5PXF8W9EUGQGgYoMURAF\nIYADwYkLcA+qUveoq7V1/Nq6Z63VulocVetW6qxbBK2KIKJWQUChCooIOEBW2GGc3x9pYgIJRgwE\n9H6e5z6SvPe973lP4nty1/nurLSHUB4ul0tRUVFEJOwhuLm5iY+VlZWRiYkJ3bhxg4iEn8OOHTvE\nx/39/cV+P3ToELm6ulZ6rVmzZtGcOXOI6G0PISEhQXzcy8uLdu/eLX5dWlpa6S9oSduVwfPnz4nD\n4VBRUZH4vUuXLpGlpaXM+lOmTCEOh0N79uyhkpISOnLkCOnr61N6errCbXXu3JkOHjyosI3yni34\naHsIhYWAmhrQsKGqLakUkZRl27bCJHQREUCXLoqfr4yQUFU4HA7OnDkDPp+PoKAgXL16FXfv3gUg\nzJVSr149vHz5ssJ5L1++hNF/k/2GhoZ48eKFwtdMSkpCixZVV2krv3Jm1KhROHz4MADhr1ZR7yAx\nMRHFxcUwMTER/yKcPn060tLSZLbbvXt3BAUFITw8HA4ODujduzeuX7+O27dvw9raGlwuF2lpacjP\nz0e7du3Ebfbt2xdv3ryp0N6LFy9gampaqe2NGzeW+iXdqFEjcW+nPImJiXjx4oXUL9yff/4Zqamp\n4utJtm9hYSGzHRHr16+HnZ0d9PX1weVykZ2dLXUfZhIJJTkcDszMzKQ+5/LXEh3jcrnIycmRutbt\n27fRo0cPGBsbQ19fHzt27Kiw4ECyvcTERMyaNUt8n40bCxM0pqSkKGT7h6L937yl5AKL7OxsuYJP\nmpqaaN68OSZOnAg1NTWMHDkS5ubmCA0NVbitnJwc6OvrK+0eqkLtDggZGQCXq2orKiU6GujcGTh2\nTKhXsHRpnU27hG7duuGbb77BvHnzAABaWlpwc3PDsWPHKtQ9duyYeJind+/eCAwMRH5+vkLXsbCw\nwNOnT2Ue09LSkmrn1atXFeqUHwb5/PPPERQUhJSUFJw+fRqjRo0CIHzANGjQAOnp6cjMzERmZiay\ns7MRHS1bFdbNzQ2PHj3CqVOn4O7uDltbWzx//hz+/v7i4SJDQ0NoamoiJiZG3GZWVlaFlVkAYGJi\nIn6AiXj+/Ll8x7zjPi0sLNC8eXPxdTMzM8Hn83H+/Hnx9STbr+xaISEhWLduHY4fP46srCxkZmZC\nT09PPHwICAO3iLKyMiQnJ6OpxCbR8tcSHePxeEhISEBZWZn4uGjIKTk5GVlZWZg+fbrU8fL3a2Fh\ngZ07d0rda15eHlxdXRWyXZLVq1dDR0dHZtHV1ZV5DpfLhYmJidQQTmRkJOzt7WXWlxwOKn8/irRV\nUlKC+Ph4me3UKFXtWlR3AUAUGUnUpo3CXaiaJD+faMECIiMjop07ieTMJdW5SeW0tDRq1KgR3bp1\ni4iIbty4QVpaWuTr60t8Pp8yMjJo0aJFxOVyKT4+noiIioqKqEOHDuTl5UX//vsvlZaW0ps3b+in\nn34if3//CtfMyckhExMT+vXXX6mwsJD4fL54mOSPP/6g1q1bU0ZGBr18+ZI6duxYYchI0l4Rffv2\npd69e1cYiho8eDDNmjWL+Hw+lZaWUnx8PF2/fl2uPzp16kS6urrioZHhw4eTrq4unThxQlxn1qxZ\nNGLECEpNTSUiouTkZPFQjeSQUVFREVlYWNDmzZupuLiYTp8+TRoaGlJDRuWHdCTvb/78+TRq1Cjx\nsdLSUmrbti2tXbuW8vPzqaSkhKKjo+mf/xZdzJs3j7p3706ZmZmUlJREDg4OcoeMLly4QE2bNqVX\nr15RUVERrVixgtTU1MTXXrZsGamrq9PJkyepuLiYNmzYQM2bN6eSkhIiEg418ng8Sk5OpvT0dOrc\nuTMtWrRI3D6Px6ObN2+KXxsbG9P+/fuJiOj27dtkbGxMY8eOJaK3Q0aSQ2enTp0ie3t7evjwIREJ\nFwccO3ZMIduVxfz588X+jImJoSZNmkgNyUmSkZFBXC6X9u/fTyUlJXT8+HFq3LixeJHCu9oKDQ0l\nOzu797JP3rMFH+2QUVZWrewhXL0K8HhCzQKR2H292u1JhTE0NMT48eOxdu1aAEDnzp0RGBiIkydP\nomnTprC0tERkZCRu3LgBKysrAICGhgb+/vtvtG7dGn369IGenh46duyIjIwMuLq6VriGtrY2Ll++\njHPnzsHExAQ2NjbiydyxY8fC0dERlpaW8PLygo+Pj0JZHUeNGoUrV66IewciDhw4AIFAADs7OxgY\nGGD48OEyex0iunfvjpKSEri4uIhf5+bmolu3buI6a9euhbW1NVxdXaGnp4c+ffrg8ePH4uMiezU0\nNHDy5Ens3r0bXC4Xfn5+GDBgADQ0NCrUlcXkyZMRExMDLpeLYcOGoV69ejh//jzu37+PFi1awMjI\nCFOnThX3TpYtW4ZmzZqhefPm8PLywrhx4+S27+XlBS8vL9jY2MDS0hKamppSQ0wcDgeDBw/G0aNH\nYWBgAD8/P5w8eRJqamri46NGjYKHhwesrKzQsmVLLF68WHz+tGnTcPDgQfHrbdu2YenSpdDV1cWP\nP/6IkSNHStlT3s4hQ4Zg3rx58PHxgZ6eHhwcHMTLNN9lu7JYsWIFrKys0KxZM/To0QPz5s2Dh4eH\n+LiOjg5CQ0MBCHsBZ8+exfr166Gvr49ffvkFZ86cEU9av6stPz8/zJgxQ+n38N5UNZJUdwFAdPo0\n0cCB7xU1q5OiIqIffiBq2pTo3DnFzkEt7yEwahYXFxfat2+fqs14J8uXL6cxY8bIPS6vpyaiqKiI\n7Ozs6NWrV9Vh3kfF69evydbWVmrSWRHkPVvwAT2E2q2HkJUFqHiSRcTjx8CoUYCJCXD//rs3T+cU\n5WD1jdU1Yxyj1hIcHAwbGxsYGhrCz88PDx48gJeXl6rNeif0ISsVIOwdyVuiyZDG2NgYMTExqjYD\nQG2fVM7KAvT0VGoCEbBrl3DieNIk4OzZyoNBGZXhQOQBtN7aGin8FPkVGZ8Ejx49gpOTE7hcLjZt\n2oQTJ07UiXQOiuwHYXx8VIuEZkBAgNfs2bN/LS0tVfvyyy93zZs3b63kcT8/v9G//PLLD0TE0dHR\nyfn9999n8Hi8KCnDOByiFSuA4mLgxx+VbqMiZGQAU6cCcXHAn38KN5tVRmJWIiadnQR+ER+b+26G\nq5mrSM6uZgxmMBifDPKeLR8ioan0HkJpaanazJkztwQEBHjFxMTYHT58+IvY2FhbyTotWrR4Ghwc\n3C0qKoq3ZMmSH6dOnSo7q5MKewjXrgGOjoC5OXD7duXBoIzKsP3udrT/oz36tOiDsMlhcDWrOJnK\nYDAYtRmlzyHcuXPHxdraOt7S0vIZAPj4+Bw5c+bMYFtbW3F6TDc3tzDR3x07drydnJxsJqMpgM8H\nbG1lHqouiouFewn27wf27AHeNdz7OP0xppybAkGpAEHjg9DG+B3dCAaDwailKD0gpKSkmJqbm4t3\ntJiZmSXfvn27o7z6u3fvntyvXz9/WceWh4UJewkpKXB3dxdvDqou4uKEE8effSacODY2ll+3uLQY\nG8I2YP3N9VjSbQlmusyEWj21arWPwWAwyhMUFCRetv2hKD0gcDgchQfMr1271mPPnj2TQkNDO8s6\nvrxZM2DCBGDAAKXZJws+HzhwAFixAli+HPjqK0DefNrTzKfYELYBV55egYWeBf6Z8g+ac5tXq30M\nBoMhj/I/llesWFHltpQeEExNTVOSkpLESUmSkpLMzczMksvXi4qK4k2ZMuWPgIAALy6Xmymzsdxc\nQE7uEGVx4wYwZoxwo1lQkPy5AkGpAOtC12HTrU342uVr7Bm8B25mbmwlBoPB+Hio6gYGeaW4uLh+\nixYtniQkJFgWFRVpODo63o+JibGVrJOYmGhhZWUVHxYW5iqvHQDVqoVQWkq0bBlRkybv3mSWVZBF\n/fz6kcdBD3qW+ey9rgO2Me2TIjg4mFq1aqVQXVmpK5RJdbdfFVJTU6l169Zi7QSGfCIjI6lTp05y\nj8t7tqA2pa6oX79+yZYtW2Z6enoG2tnZxYwcOfKora1t7I4dO6bt2LFjGgCsXLlyaWZmJnfGjBm/\nOzs7R7i4uMiW1aqmHsKLF0C3bsIeQUSE/BGpnKIc/BT8E6w3W8NS3xLnvziPZvrNlG6PKmESmsql\nLklo1qtXT26SwepizZo1mDhxIhrU1QyQADIyMjB06FBoa2vD0tJSnGlXFvv27YOamppUQj2RTgIg\nHO7R1NQUH7OVWETD4/Ggr68vTl5YI1Q1klR3AUBkYkKUnKxALFWcBw+ImjUjWrKE6L88XRXILcql\ntTfWkvE6Yxr912h69KaiHoCioJb3ECRTELx69YocHR3p+++/Fx+/efMmaWtrk6+vL+Xm5lJGRgYt\nXryYuFwuPX36lIiEaQrat29PHh4eFBsbS0TCX4KrVq2SmdxOWRQXF1db2zXBu37Bl1ePU3b7HA5H\nnKBQFsr2b2FhIRkaGlJKSkqVzq8tn7ePjw/5+PhQXl4e3bhxg/T09MRJ+Mqzd+9e6tq1q9y23N3d\npXQfyuPn50cDBgyQeUzeswUf0ENQ+YNfrmEAkY4OUVaWXGe9D2VlRCdOCLOTytOgyBfk06awTdRk\nfRMafmw4PUyV/SG/D3UpIBAxCc1PRUKza9euxOFwSEtLi7S1tenYsWN07do1MjU1pbVr11KTJk1o\n7NixtG/fPurSpYvUuRwORywQVJkvynP9+nWytraWem/Pnj1ka2tLOjo61KJFCynRnfL2jBs3jsrK\nyujnn38mKysraty4MY0YMYIyMjLE53z++efUpEkT0tPTo27dusl9UFeV3Nxc0tDQoLi4OPF748aN\no/nz58usv3fv3gr+k8Td3Z127dol93hycjJpamqSQCCocKw6AkLtTl2RlwdoaX1wM0+eAL17C1cQ\nHT8unESWJLswG5tvb0bLzS0R9CwIAaMDcGz4MdgZ2X3wtesCwu8Qk9D8lCQ0RcMWUVFRyMnJEX/G\nr1+/RmZmJp4/f46dO3eKvxvyUNQXgFBas1WrVlLvffbZZ7hw4QL4fD727t2LOXPmICIiQnxc0p4d\nO3bA19cXZ8+eRXBwMF6+fAkul4uvv/5aXL9///6Ij49HWloa2rZtK/4eyOKrr76SK63p5OQk85zH\njx+L5S5FODo6ys3bxOFwEBERASMjI7Rq1QqrVq2SkiEFgAULFsDIyAhdunTB9evXpY6ZmppCXV0d\njx49knsfSqWqkaS6CwCi+vXlRk5F2bWLqHFjonXriMr3OPMEebQyaCVx13DJ+6g33U25+8HXKw8U\nkdBcjg8uVYVJaL7lU5PQlPylL7JHQ0NDKuumrF+4ovPe5YvyrFq1Sq5IvYghQ4bQb7/9JtceW1tb\nqR7tixcvSF1dXWavLzMzkzgcDvH5/Eqv+T4EBwdTkyZNpN7buXMnubu7y6z/9OlTevbsGRERRUdH\nk52dnbg3TCTUhsjNzSWBQED79+8nHR0dKX8SEZmamlKISJNXAnnPFny02U4/UDpz927h3oLQUEDy\nh0kZlcEvyg8Lry6Em5kb7k69ixbcqks6fii0THW5jkQSmj179kRwcDAGDhyIu3fvwsXFRUpC08bG\nRuq82iihOXbsWLkSmiLKysrk5s7v3r075s6dKyWhOXnyZCkJzdTUVLGEpggiqqD+BVSvhKaI0tJS\nsVbD+0poysLIyEhKr6EyJOVERcjzBQAYGBhUkNa8ePEiVqxYgbi4OJSVlSE/Px88Hk+uPc+ePcPQ\noUOlfFa/fn28fv0axsbGWLRoEU6cOIG0tDTUq1cPHA4Hb968kSt9+b5oa2tXWHRRmbRm8+Zv9yjZ\n29tj6dKlWLduHebPnw8AYt0NABg3bhwOHz4Mf39/zJw5U/x+TUpr1u4hI03NKp1WVATMnAn8/DPg\n7y8dDJ5lPUPHXR2x9Z+tOOJ9BMeGH1NpMKhNMAnNT0dCU9FrVvZ5vI8vAOGqGUkhoaKiInh7e+OH\nH35AamoqMjMz0a9fP6lhKlk+CAgIkPJBfn4+TExM8Oeff+Ls2bO4cuUKsrOzkZCQIDniUIHp06fL\nldZ0cHCQeY6NjY1Y7lJEZdKaspBnjyxSUlIgEAgqDLVVG1XtWlR3AUBUrnusCE+fErVrRzR0qPR8\ndGFxIa0OXk0Gaw1oTcgaKisre++2qwLq2KQyk9D8NCQ0iYiaNGlCly5dEr+WZc+jR4+oQYMGdP/+\nfSooKKBp06ZJDTVV5ovyFBUVkZGRkXiVEZ/PJzU1Nbp+/TqVlZWRv78/NWrUqFL/bNq0idzd3Skx\nMZGIhKvZzpw5Q0RE27ZtIycnJ+Lz+ZSbm0szZsyoMCymDHx8fOiLL76gvLw8CgkJIT09PYqJiZFZ\n19/fXywSFBsbS/b29rRy5UoiEi54CAgIoIKCAiouLqZDhw6RlpaW1IS1n58f9e/fX2bb8p4t+Ggn\nld9zyOj0aaBjR2D0aOCvv94mSr0YdxH2v9vjVsot/DPlH8zrMo/tMJYDk9D8NCQ0AWD58uUYP348\nuFwuTpw4IXNC38bGBkuXLkXv3r3RqlUrdO3aVarOu3whiYaGBiZMmIBDhw4BEEpQ+vr6YsSIETAw\nMMDhw4cxePBgqXPK2zNr1iwMGjQIHh4e0NXVhZubG+7cEW5jGjduHJo1awZTU1PY29vDza16Mgls\n27YNBQUFMDY2xpgxY7B9+3bx/oHnz59DR0cHycnC5AxXr16Fo6MjtLW10b9/f3h7e2PhwoUAgOLi\nYixZsgTGxsYwMjLC1q1bcebMGakJaz8/P0yfPl3p9yCXqkaS6i4AiBwcZEZAWfz8s3B/wX8/bImI\nKD0/nQYdHkTWvtZ04fEFhdtSJqjlPQRGzVJXJDSri7S0NLZTWUE+ip3KSkXB3YxHjwJr1gA3bwp7\nCADw99O/0X5ne5jpmuHBjAfo17JfNRoq5NmzZ+JligwGIFze+erVK5SUlGD//v11RkKzujA0NERs\nbGyd3qlcU/B4PISGhtboNWt3QFBgyGjPHuB//wOCg4GmTYGSshL8cPkHTDozCVv7bcXWflvRoH71\nfvlKSkqwYcMGtG/fHg8ePKjWazHqFnVVQpPxaVItEprKgMPhEPXuDfy3+UkW+/cD8+cD168DNjZA\nniAPg48MFi4rHeYHEx0Tuecqi/DwcEyZMgX6+vrYsWOH1PgfIF/mjsFgMD6EOiGhqVQq6VZu2AAs\nWwZcuiQMBgXFBei+rzuMtIxweezlGgkGGzduRN++ffHtt9/i77//rhAMGAwGoy5RJwPCL78AO3cC\nISGAgwNwJ+UOTDeaoo1xG/w57M8aUy7r2bMnoqOjMX78eLZqicFg1Hlq905ldfUKb/36qzAYXL8O\nmJoCV55ewayAWZjrNheLuy2uUfPk5TthMBiMukidCgj/+x9w7hxw5QrQgPsGw4/PQPjLcKztvRbe\ntt7VZgYRoaSkBOoyAhSDwWB8LNTuIaP6b+PVkyfCSeSgEAGOJv+C1ltaw0LPAjFfxeBzu8+rbcjm\nyZMn8PDwgK+vb7W0z2AwGLWFOhEQSkqA774DJk0rwOJb03Eh7gJuTr6JDR4bqm1JaXFxMX755Rd0\n7NgRnp6emDVrVrVch1G7mDhxIgwMDGTusJZFYGAghg4dWs1WfRxs2bJFnNSNUUup6o626i4AiGbM\noLIyoqlTiVp4XiDzDc1oxPERlJqbWvkWvw/kn3/+IUdHR+rTp88H50FBLd+p3KxZM9LU1CRtbW36\n7LPPaMyYMZSdnS1VJzQ0lHr06EE6Ojqkp6dHAwcOrJC7JTs7m2bNmkUWFhakra1NVlZWNHv2bHrz\n5k1N3s4HERwcTGZmZpSfn09ERAKBgLy9vcnS0pI4HA4FBQVVOKddu3biPEx1lYiICGrbti01atSI\n2rVrR/fv35dbNzk5mQYNGkQGBgZkZmZG27dvFx979OgRDRo0iIyMjMjAwIA8PT2lUqcXFhaSmZmZ\nOO8R48OQ92zBR7tTuX59HPZ/Dr/88cjv/SV2DPwdRz8/CiMto2q97L59+/Ddd98hMDDwg9I00rQr\noAAAFSpJREFU1wU4HA7Onz+PnJwcREZGIjo6GqtWrRIfDwsLg6enJ4YOHYqXL18iISEBjo6O6Ny5\nMxISEgAIhWB69eqF2NhYBAYGIicnB2FhYTA0NBTnmakOSkpKlNpeYmIiLC0toSmRZbdbt244dOgQ\nmjRpUmFY8p9//gGfz5dKYfw+lBdKUQUCgQCDBw/GuHHjkJWVhfHjx2Pw4MEoLi6WWX/MmDGwsrJC\namoqLly4gIULF4rzUGVnZ2PIkCF4/PgxXr9+DRcXF6ncRA0aNEDfvn1x4MCBmrg1RlWoaiSp7gKA\n8mdOJ84yNXJbO4H4hcoTuahJUMt7CExCU8iuXbuoYcOGpKamRtra2rR8+XKp42ZmZhWypK5YsYKm\nTJlS4drm5uakq6tL7dq1kxI2WbZsGXl7e9OYMWNIV1eXdu/eXal98fHx1KNHD2rcuDEZGhrS6NGj\nKUtJkrIiAgMDK4j4WFhYUEBAQIW6OTk5xOFwKC0tTfze1KlTaezYsTLbTk9PJw6HIyVx6efnJ1eg\niPF+yHu24GPtIcSkPoFa/EBc+movdBooR+CCUREiJqE5efJkbN++HW5ubsjJycGyZcveeR8PHjyo\nkKfexcUFkZGRyMzMxKhRozB8+HAIBALx8bNnz2L48OHIzs7GqFGj3mnfokWL8PLlS8TGxiIpKQnL\nly+Xaw+Px5MrCSkpuCLJw4cPpQRpAPmSkKLviehfQCg4JC9dS3BwMExMTKQEfVq3bo3IyEi598BQ\nLbU6INxKuYtWJSOhra38tokIe/bsQUxMjPIbf184nA8vVYSIMGTIEOjq6sLCwgJWVlZYvFi4nyMj\nIwNlZWVSimMimjRpgjdv3gAA0tPTZdaRx/nz59G0aVPMmTMHGhoa0NbWRocOHaRskgeHw8G3334L\nU1NTNGjQABYWFmjbti1OnToFQJhuuFGjRnBxccHr169x8eJFbNq0CZqamjAyMsLs2bNx5MgRub54\nH7KysiooZY0ePVqsNDd37lwUFRVJ6eF26tQJgwYNAiAcYqnMPisrK/Tq1Qvq6uowNDTEnDlzKmju\nShIVFSUlHCNZtmzZIvOc3Nxc6InyxP+Hrq5uBWUzQJiuunPnzvjxxx9RVFSE8PBwnDx5EgUFBRXq\nJicnY+bMmdi4cWOFNrKzs+XeA0O11OqAUCAoxbavRiq93cePH6Nnz574/fffld52lSD68FJFRBKa\nfD4fQUFBuHr1Ku7evQsAUhKa5amNEpoA5Epoin4pT58+HWlpaVW+tiRcLreCOtj69ethZ2cHfX19\ncLlcZGdniwMnAJiZmYn/fpd9r1+/ho+PD8zMzKCnp4exY8ciPT1dKbaL0NHRkSkJqaurK7O+n58f\nEhISYG5ujq+//hpjxoypIBOalpYGDw8PfP311xg5Uvr/b05OToUAxKg91OqAUFqki06dlLe/QCAQ\n4KeffkKnTp0wePBg3Lp1C3Z2dkprv67zqUtovi/lJSFDQkKwbt06HD9+HFlZWcjMzISenp5Uz0PS\n9nfZt3DhQqipqeHBgwfIzs7GwYMH5eoVA0CbNm3kSkJ+9dVXcs+JioqSei8qKgpt2rSRWd/CwgLn\nzp1DamoqwsLCkJaWJh5iBIDMzEx4eHhgyJAhWLBgQYXzY2Nj2Q7/WkytDggN6+lJ7k37IIgI7u7u\nCA0Nxb179zB79myoqdVMzqO6xOzZs3Hnzh3cvn0bALBmzRrs378fmzdvRk5ODjIzM7F48WLcvn1b\nPM4+duxYmJubw9vbG48ePUJZWRnS09OxevVqXLx4scI1BgwYgJcvX+K3335DUVERcnJyxKuRnJyc\n4O/vj8zMTLx69Qq//vrrO202MjKCu7s7JkyYgBYtWojH9U1MTODh4YG5c+ciJycHZWVlePLkCYKD\ngxX2R1FREQoLCyv8DQD9+vWTGsLJyclB/fr1YWhoCIFAgJUrV8rVF1bEvtzcXGhpaUFXVxcpKSlY\nt25dpbY+fPgQOTk5Msu2bdtknuPu7g41NTX4+vqiqKgIvr6+qFevHnr27Cmz/r///oucnBwIBAIc\nOnQIly9fxty5cwEAfD4fnp6e6NKlC1avXi3z/OvXr6Nv376V3gdDddTqgNCogZbS2uJwONi/fz8u\nXLiAZs2aKa3dj41PWUJT1oR2q1at0KhRI7x48QKenp7Q0tISi9c7OztDT09PHMy8vLzg5eUFGxsb\n8fJVCwuLStuvzL5ly5YhPDwcenp6GDhwILy9vZW+I19dXR2nT5/GgQMHwOVyceDAAZw+fRr1//sl\n5ufnJyUgHxgYCCsrKxgYGGDnzp0IDAxE48aNAQCnTp3C3bt3sXfvXnHPRFdXVywnWVhYiIsXL2L8\n+PFKvQeG8qjVegi7O7lhUuhNVZvyQTA9hI+by5cvY9u2beJJbYZ8tmzZguTkZKxZs0bVpnwUVIce\nQq0OCLu6dMbkkPeXpExLS4OhoWGtSEnNAgKDwagOPj2BnPd8oJeVlWHnzp2ws7Nja50ZDAbjPand\n6a85iserf//9F1OnToVAIMDVq1fh4OBQjYYxGAzGx0ed7yEIBAKsWLECXbp0wYgRIxAaGsqCAYPB\nYFSBWt5DePeyUA6Hg6ysLERERFTYsFQb4HK5tWIug8FgfFxIpgRRFrU7INR794NUXV0dmzZtqgFj\nqkZGRoaqTWAwGAyFqNVDRlSvVptXY4jW6DOYLyRhvngL84VyqJYnbkBAgFfr1q3/bdmyZdzatWvn\nyarz7bff+rZs2TLO0dExMiIiwllWHY7EpHJycjImT56MrKys6jC5VsO+7G9hvngL88VbmC+Ug9ID\nQmlpqdrMmTO3BAQEeMXExNgdPnz4i9jYWFvJOv7+/v3i4+Ot4+LiWu7cuXPqjBkzZGeZ49RDWVkZ\ntm7dCicnJ5iZmaFhw4bKNpnBYDAYqIY5hDt37rhYW1vHW1paPgMAHx+fI2fOnBlsa2sbK6pz9uzZ\nQePHj98PAB07drydlZWl//r1688+++yz15JtpeTlokuXLuBwOAgODmaJ6BgMBqMaUXpASElJMTU3\nN08SvTYzM0u+fft2x3fVSU5ONisfEJaFhoj/lpd98VNhxYoVqjah1sB88Rbmi7cwX3w4Sg8IHA5H\noTwN5bdWlz+vqluvGQwGg1E1lD6HYGpqmpKUlCTeEJCUlGRuZmaWXFmd5ORkM1NT0xRl28JgMBgM\nxVF6QGjfvv3duLi4ls+ePbMUCAQaR48eHTlo0KCzknUGDRp09sCBA+MA4NatW676+vpZ5YeLGAwG\ng1GzKH3IqH79+iVbtmyZ6enpGVhaWqo2efLk3ba2trE7duyYBgDTpk3b0a9fP39/f/9+1tbW8Vpa\nWnl79+6dqGw7GAwGg/GeEJFKy8WLF71atWr1r7W1ddyaNWvmyarzzTff+FpbW8fxeLzI8PBwZ1Xb\nrCpfHDp0aDSPx4t0cHCI6tSpU2hkZCRP1Tar8ntBRLhz504HNTW1kr/++muYqm1WpS+uXbvm7uTk\nFNGmTZsH3bt3D1K1zaryRVpamqGnp2eAo6Pj/TZt2jzYu3fvBFXbXB1l4sSJe4yNjV/b29tHy6tT\nleemSm+qpKREzcrKKj4hIcFSIBCoOzo63o+JibGVrHPhwoV+ffv29Sci3Lp1q2PHjh1vqfrDUJUv\nbt686ZaVlaVHJPyP8Sn7QlSvR48eV/v373/+xIkT3qq2W1W+yMzM1Lezs3uYlJRkRiR8KKrablX5\nYtmyZcvnz5//s8gPBgYG6cXFxfVVbbuyS3BwcNfw8HBneQGhqs9NleaGkNyzoK6uXizasyBZR96e\nBdVYXH0o4gs3N7cwPT29bEDoi+TkZDPVWFu9KOILANi8efM3n3/++QkjI6M0VdhZEyjiiz///HOU\nt7f3X6LFG4aGhm9UY231oogvTExMXvL5fF0A4PP5uo0bN06vX79+iWosrj66du0awuVyM+Udr+pz\nU6UBQdZ+hJSUFNN31fkYH4SK+EKS3bt3T+7Xr59/zVhXsyj6vThz5sxg0S53RZc71zUU8UVcXFzL\njIwMgx49elxr37793YMHD46teUurH0V8MWXKlD8ePnzYpmnTpi8cHR0jf/vtt1k1b6nqqepzU6XZ\nTpW1Z+Fj4H3u6dq1az327NkzKTQ0tHN12qQqFPHF7Nmzf12zZs18DodDRMQp/x35WFDEF8XFxerh\n4eFtr1y50is/P7+Rm5tbmKur662WLVvG1YSNNYUivli9evVCJyen+0FBQe5Pnjyx6tOnz+XIyEhH\nHR2dnJqwsTZRleemSgMC27PwFkV8AQBRUVG8KVOm/BEQEOBVWZexLqOIL+7du9fOx8fnCAC8efPG\n8OLFi33V1dWLyy9xruso4gtzc/MkQ0PDN5qamgWampoF3bp1C46MjHT82AKCIr64efNmp0WLFv0E\nAFZWVk+aN2+e8OjRo1bt27e/W9P2qpIqPzdVOTFSXFxcv0WLFk8SEhIsi4qKNN41qRwWFub6sU6k\nKuKLxMRECysrq/iwsDBXVdural9IlgkTJuz9WFcZKeKL2NjY1r169fq7pKRELS8vr5G9vX30w4cP\n7VRtuyp8MWfOnI3Lly9fRkR49erVZ6ampsnp6ekGqra9OkpCQoKlIpPK7/PcVPlN+fv797WxsXlk\nZWUVv3r16gVEhO3bt0/bvn37NFGdr7/+eouVlVU8j8eLvHfvXltV26wqX0yePHmXgYFBupOTU4ST\nk1NEhw4d7qjaZlV+L0TlYw4Iivpi3bp139nZ2T20t7eP/u23375Vtc2q8kVaWprhgAEDzvF4vEh7\ne/toPz+/Uaq2uTqKj4/PYRMTkxfq6uoCMzOzpN27d09SxnOTQ/TRDcczGAwGowowSTIGg8FgAGAB\ngcFgMBj/wQICg8FgMACwgMBgMBiM/2ABgaFy1NTUSp2dnSNE5fnz5xby6mpra+d+6PUmTJiwr0WL\nFk+dnZ0j2rVrd+/WrVuu79vGlClT/vj3339bA8LNUJLHOnfuHPqhNgJv/cLj8aKGDRt2Mjc3V7uy\n+pGRkY4XL17sq4xrMz5RVL18ihVWtLW1c6qjrrwiuUz10qVLfXg8XmRN2V/VdsePH79v/fr1/6us\n/t69eyfMnDlzs6o/T1bqbmE9BEatIy8vT6t3795/t2vX7h6Px4s6e/bsoPJ1Xr58adKtW7dgZ2fn\nCAcHh+gbN250AYBLly55dOrU6Wa7du3ujRgx4lheXp6WrGvQf9v6u3btGhIfH28NABs3bpzr4OAQ\n7eDgEC3KgZOXl6fVv3//C05OTvcdHByijx8/PhwA3N3dg+7du9du/vz5awoKCjSdnZ0jxo4dexB4\n24vx8fE54u/v3090zQkTJuw7efLksLKysnrff//9OhcXlzuOjo6RO3funPoun7i5uYU9efLEChAm\neevUqdPNtm3bhnfu3Dn08ePHNgKBQGPp0qUrjx49OtLZ2Tni+PHjw/Py8rQmTZq0p2PHjrfbtm0b\nLsuPDIYUqo5IrLCipqZWItpsN2zYsL9KSkrU+Hy+DpFwo5G1tXWcqK7oV/P69ev/99NPPy0kIpSW\nltbLycnRTktLM+zWrdv1/Px8TSLCmjVr5q1cuXJJ+etNmDBhryhd9rFjx4a7urqG3bt3r62Dg0NU\nfn6+Zm5urlabNm0eREREOJ04ccJ7ypQpO0XnZmdn6xIR3N3dr4k2+5TvIYhenzp1asj48eP3ERGK\nioo0zM3NnxcWFjbYsWPH1FWrVi0iIhQWFjZo3779PwkJCZbl7RS1U1JSojZs2LC/tm7d+hURgc/n\n65SUlKgRES5fvtzb29v7BBFh375947/55htf0fkLFixYfejQodFEwhTZNjY2j/Ly8hqp+vNmpfYW\nleYyYjAAQFNTsyAiIsJZ9Lq4uFh9wYIFP4eEhHStV69e2YsXL5qmpqYaGxsbp4rquLi43Jk0adKe\n4uJi9SFDhpx2dHSMDAoKco+JibHr1KnTTQAQCAQaor8lISLO999/v27VqlWLjY2NU3fv3j358uXL\nfYYNG3ZSU1OzAACGDRt2MiQkpKuXl1fAd999t37+/PlrBgwYcL5Lly43FL0vLy+vgFmzZv0mEAg0\nLl682Ld79+7XGzRoUHTp0iWP6OhohxMnTnwOCNM0x8fHW1taWj6TPF/U80hJSTG1tLR8Nn369O0A\nkJWVpT9u3LgD8fHx1hwOh0pKSuqL7oskEppdunTJ49y5cwPXr1//HQAUFRU1SEpKMm/VqtUjRe+B\n8WnBAgKj1uHn5zf6zZs3huHh4W3V1NRKmzdvnlBYWNhQsk7Xrl1DQkJCup4/f37AhAkT9s2dO3cj\nl8vN7NOnz+U///xzVGXtczgcWr9+/XfDhg07KXrv77//7i35MCUiDofDoZYtW8ZFREQ4X7hwof/i\nxYtX9erV68qSJUt+VOQ+GjZsWOju7h4UGBjoeezYsRFffPHFYdGxLVu2zOzTp8/lys4XBcqCggJN\nT0/PwDNnzgweOnToqSVLlvzYq1evK6dOnRqamJjYzN3dPUheGydPnhz2sSW5Y1QfbA6BUevg8/m6\nxsbGqWpqaqXXrl3rkZiY2Kx8nefPn1sYGRmlffnll7u+/PLLXREREc6urq63QkNDO4vG2vPy8rTi\n4uJayroGlUsN3LVr15DTp08PKSgo0MzLy9M6ffr0kK5du4a8fPnSpGHDhoWjR4/2++6779ZL9mRE\nqKurF4t+pZdn5MiRR/fs2TNJ1NsAAE9Pz8Bt27Z9JTrn8ePHNvn5+Y3k+UNTU7PA19f320WLFv1E\nRBw+n6/btGnTFwAgqUeuq6vLz8nJ0RG99vT0DPT19f1W9FqW7QyGFKoes2KFFR0dHb7k6zdv3jR2\nc3O76eDgEDVx4sQ9dnZ2DxMTEy0k6+7bt2+8vb19tLOzc3i3bt2uP3v2rBkR4erVqz06dOhwh8fj\nRfJ4vMhz584NKH89ecnwNm7cOMfe3j5aMkFcYGCgB4/HixQlExTNG0jOIcybN2+Nra1tzJgxYw6W\nv5/i4uL6BgYG6ZMmTdoteq+srIyzcOHCnxwcHKLs7e2je/bseUU0N1GZXwYOHHj2yJEjI8PCwlxt\nbGweOTs7hy9evPjH5s2bPyUiZGRkcDt06HDHyckp4tixY8MLCgoaTps2bbuDg0NUmzZtHgwcOPCs\nqj9rVmp3YcntGAwGgwGADRkxGAwG4z9YQGAwGAwGABYQGAwGg/EfLCAwGAwGAwALCAwGg8H4DxYQ\nGAwGgwEA+D8qs8NwAc2L4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x41546a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fpr_clf, tpr_clf, _  = roc_curve(targetdata.ravel(),predictions_weighted_clf.ravel())\n",
    "roc_auc_clf = auc(fpr_clf, tpr_clf)\n",
    "fpr_true, tpr_true, _  = roc_curve(targetdata.ravel(),predictions_weighted_true.ravel())\n",
    "roc_auc_true = auc(fpr_true, tpr_true)\n",
    "fpr, tpr, _  = roc_curve(targetdata.ravel(),predictions.ravel())\n",
    "roc_auc = auc(fpr, tpr)\n",
    "fig = plt.figure()\n",
    "plt.plot(fpr_clf, tpr_clf, label='ROC curve weighted app(area = %0.2f)' % roc_auc_clf)\n",
    "plt.plot(fpr_true, tpr_true, label='ROC curve weighted true (area = %0.2f)' % roc_auc_true)\n",
    "plt.plot(fpr, tpr, label='ROC curve f1(area = %0.2f)' % roc_auc)\n",
    "plt.plot([0, 1], [0, 1], 'k--')\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('{0} n_epochs {1}'.format('ROC',epochs_training))\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.savefig('{0}/plots/{1}/{2}_{3}.png'.format(dir,model_g,'roc_eval',epochs_training))\n",
    "\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}