{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imputation of missing values using Matrix Factorization\n",
    "\n",
    "This is a short example about imputation of missing values using the package [cmfrec](https://github.com/david-cortes/cmfrec) (low-rank matrix factorization), which among other things, can be used for imputation of tabular data using matrix factorization models.\n",
    "\n",
    "The basic idea is to produce an approximate factorization of a matrix $\\mathbf{X}$ with misisng data as the product of two lower-dimensional matrices, determined by finding two such matrices that minimize the squared error with respect to the non-missing entries:\n",
    "$$\n",
    "\\mathbf{X} \\approx \\mathbf{A} \\mathbf{B}^T\n",
    "$$\n",
    "\n",
    "Optionally adding other aspects such as centering and column means - e.g.:\n",
    "$$\n",
    "\\mathbf{X} \\approx \\mathbf{A} \\mathbf{B}^T + \\mu + \\mathbf{b}\n",
    "$$\n",
    "\n",
    "Missing values are then imputed according to the the values that such factorization (fitted to the non-missing values) would produce.\n",
    "\n",
    "Typically, values imputed in this way are not as good quality as those produce by iterative imputation methods, but for high-dimensional data, such a model is typically faster to fit and to produce imputations.\n",
    "\n",
    "Note that `cmfrec` is aimed primarily at recommender systems, and using it for imputation of tabular data benefits from changing the default parameters - particularly: do not add user (row) biases, and do not regularize the item biases (column intercepts). See the code below for an example of changing these parameters.\n",
    "\n",
    "For more information about the library and methods, see the package's webpage:\n",
    "https://github.com/david-cortes/cmfrec\n",
    "** *\n",
    "The example here is copy-pasted from SciKit-Learn's usage guide for their imputer, and just adds a few extra lines that do the job with this package.\n",
    "\n",
    "Original code was taken from this link:\n",
    "\n",
    "https://scikit-learn.org/stable/auto_examples/impute/plot_iterative_imputer_variants_comparison.html#sphx-glr-auto-examples-impute-plot-iterative-imputer-variants-comparison-py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAGqCAYAAAARNPnDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABV1ElEQVR4nO3deZhlVXn+/e8tICAIREFFRFoRFIPQDKIoYBsUQYjiSARHEgmO0YixHaJoTCTRKDj9FH0VJwQjiAgqaLQZZB4b0IAM7azIIIJMDTzvH3sVffpQY3f1rqb7+7muvjhnD2s9e+9TRd211t6VqkKSJEmSpL48YKYLkCRJkiStXAyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSNClJKsnj2uvPJPnXgXWvS/KHJLckeeg097tfkpOns81l7f5Ycx+S7Jzk8uW1/ySz2ud81Wnsc7HPQpKnJ/l5+1rZO8nDk5ya5OYk/z1d/WrmP29LIsm8JP8wTW0dkeSD09GWtCwYRCVpJZJk3yTntR+Cf5fke0l2mmo7VXVgVf1ba3M14KPAblW1dlVdP501V9XXqmq3Jdl3tB/ElkXYGLY0NU8kyYIkt7Vr+Pt2jGsvi76mW1WdVlWPX176b+fyWUvaXjv3d7YQeXOSS5N8KMm6A30OfxY+AHyyfa0cBxwAXAesU1VvW9JalrD+Vyc5fYJtpi0YTUWSOUl+PcV97v1lGSy7z9vA95ALhpav3z4PCybZzsFJvjrd9Un3FwZRSVpJJPln4FDgP4CHA48GPg08fymbfjiwBnDZEtSUJP6/aOr+tqrWBmYD2wDvnO4OlmVQX8H8V1U9GNgAeA3wVOAnSdYaY/tNWPxrZRPgp1VVU+3YazTj1kqy5cD7fYFrZqoY6f7G//lL0kqgjdB8AHhDVR1bVX+pqoVV9Z2qenvbZockZyb5Uxst/WSSB47R3hFJPphkc2Bk6tufkvyorX9aknOT3NT++7SBfecl+fckPwFuBR7bRhcObFMWb0zyqSRp2y82apPksCS/SvLnJOcn2Xlpz02SLyf5Y5JfJHnPSDgeHrEYHk1ttV3dRsOuSbLfGDWPd3yrJPnvJNe1Nt442RHbqvo9cBJdIB3p66lJzmjX8eIkcwbWPSaLpoH+sNXx1aFj+/skvwRGruX+SX7W6j4pySZteZJ8LMm17TrPH/mhPMlzk/y09fObJAe15YuNciXZon0e/pTksiTPG1h3RKvvxNbO2Uk2HeMafinJ29rrjdpxvL69f1ySG1q99/af5Ct0v4z5TrrR5X8ZaHK/JL9s1+TdE12Hdi1ur6pzgecBD6ULpYt9FpJcBTx2oM+vA68C/qW9f1aSBySZm+SqJNcn+UaShyzJNWrrRv3sJdkC+AywY+v7TxMd48j5S/Iv7br/Lt304ucmuaKd53cNbH9wkm8mObpdwwuSbD1U2+MG3o98X1kL+B7wyFbbLUkemXG+RyU5tTVzcdt+n2X1eRvwlXb9RrwS+PLQOXtkkmPSfX+5Jsmb2/LdgXcB+7R6Lx7YbZMkP2l1nJxk/YH2ntdq/1M7li0G1m3TzvHNSY6m+wXhyLr1k5zQ9rshyWnxl4CaYX4AJWnlsCPdDyXfGmebu4G3Auu37XcFXj9eo1V1BfDX7e16VfU37YfmE4GP0/1A/lHgxCx+7+gr6KYkPhj4RVu2F/BkYGvgpcBzxuj2XLrg9RDgSOB/kqwxxraT8QlgXbqA8Ay6HyZfM9FO7YfljwN7tBGxpwEXjbPLWMf3WmAPumPaFth7soUneVTb98r2fiO6c/9BuvNzEHBMkg3aLkcC59Bdl4PprsOwZwBbAM9JsjfdD8svpBvxOw34ettuN2AXYHNgPWAfYGRa9v8H/GM7L1vSAtNQ7asB3wFOBh4GvAn4WpLBqZQvA94P/FU7xn8f41ScAswZqP/q9l9ajacNjzhW1SuAX9JGl6vqvwZW7wQ8nu5r4L2DP+xPpKpuBn4A3OcXJFW16VCfLwO+RjequnZV/RB4M91n4BnAI4EbgU8NNTXZazTiPp+9qvoZcCBwZut7vUke4iPovpdsBLwX+BzwcmC7dszvTfLYge2fD/wPi75ej2vXfkxV9Re6z/VvW21rV9VvGed7VFXt0nbfum1/9GCb0/x5G/FV4O/S/TJpC7rvZ2cP9PmA1ufFdOdrV+AtSZ5TVd+nm51ydKt364F296X7HvQw4IF0X8ek+8Xf14G30F3r79L9UuOBLZAfRxeOH0J3zl800ObbgF+3/R5O95mZ8ii8NJ0MopK0cngocF1V3TXWBlV1flWdVVV3VdUC4LMs+mF+KvYEfl5VX2ltfR34P+BvB7Y5oqoua+sXtmWHVNWfquqXwI8ZGOUbqvOrVXV92/e/gdXpQsNYDmqjAH9qoz7zR1YkWYUuQL2zqm5ux/3fjB7QRnMPsGWSNavqd1U13vTksY7vpcBhVfXrqroROGQS/R6X5GbgV8C1wPva8pcD362q71bVPVX1A+A84LlJHk0XRt5bVXdW1enA8aO0fXAbMb8N+EfgQ1X1s/bZ+Q9gdhtxW0j3g/cTgLRtftfaWAg8Mck6VXVjVV0wSj9PBdZu5+XOqvoRcAJdGBhxbFWd0/r+GmN8JuiC6M7tB/9dgP8Cnt7WPaOtn4r3V9VtVXUxXYjYeqIdhvyWLgwsiX8E3t0+D3fQ/cLgxVl8hHyy12jEpL62Jmkh8O/t6/YoulB4WPv6uYxu2vFWA9ufX1XfbNt/lC7EPnVJOl7K71HT+Xkb8Wu6GSHPohsZ/fLQ+icDG1TVB1qfV9MF97+boN0vVtUV7fp+Y6COfYATq+oH7Xx+BFiT7pdgTwVWAw5ts12+SfdLuxELgQ2BTdr6+/xyRuqbQVSSVg7XA+tnnOmeSTZvU7d+n+TPdD/Qrj/W9uN4JItGOUf8gm5EYMSvRtnv9wOvb6X7oXG0Ot+WbhriTS1YrjtBnR+pqvVG/rH4D8nr0404DNY7XOuo2qjNPnSjSr9rU/qeMM4uYx3fI1n8fIx2bobt3UYb59AFwZHj3wR4yVDw3onuB9BHAjdU1a0T9DW4bBPgsIG2bgACbNR+kP8k3WjdH5IcnmSdtt+LgOcCv0hySpIdR+nnkcCvquqegWXD535Sn4mqugq4he4H9p3pAsZv22jXkgTRSfU7jo3oztWS2AT41sA5/xndSODDB7aZ1DUa2GZpj2fQ9VV1d3t9W/vvHwbW3zbU/r21tmv9a7prP2VL+T1q2j5vQ74MvJou0A4/eGgTuunFg1+P72Lxazma8b5X3Pu9qh3Lr9oxPBL4zVC4HPy+9mG6Ud6T091OMHfiQ5OWLYOoJK0czgRuZ/xpn/+PbuRys6pah+4HpixBX7+l+wFs0KOB3wy8X6LfxKe7H/QddKOIf9WC5U0sWZ3QPa10IYvXO1jrX4AHDax7xODOVXVSVT2bLuj9H91ox1T9DnjUwPuNJ7tjVZ0CHEE3MgLdD6VfGQzeVbVWVR3S+nlIksHjGa2vwWvzK7optoPtrVlVZ7T+P15V29FNz94ceHtbfm5VPZ9uauFxdKM6w34LbDx0n9rw52QqTgFeDDywqn7T3r+SbprlRWPsM+0jQumeYPwsuimyS+JXdNO9B8/5Gu2YRkz6Gk2gjxGxez9j7Vo/iu7aQxeyxvr6Gq22pfkeNd2ftxHH0M0Cubqqhn8B9yvgmqFr8+Cqem5bP9Xzv9j31iShO7+/ofv63qgtG/HokRdtxPptVfVYutkp/5xk1yn2L00rg6gkrQSq6ia6+7k+le7hIg9KslqSPZKM3Bv3YODPwC1tZO91S9jdd4HN0/2pmFWT7AM8kW6Uamk9GLgL+COwapL3AuuMv8vY2sjON4B/T/LgNp3xn1k0snERsEuSR6d74NO9T6dN9/cfn9fuFb2DbkTubqbuG8A/pXvIznp0QXsqDgWenWR2q/tvkzyn3be2RroHtjyq/ZB8HnBwu6dsRxafLj2azwDvTPLXcO+DnV7SXj85yVPavXd/oftFx92t7f2SrNumD/6Z0c/L2W2/f2mfxTmtnqOmePwjTgHeCIw8tGYe3X2Apw+M4A37A929wUstyepJtqML3jcCX1zCpj5D93kceSjUBknGe7L1mNdoEv4APCpjPJRsmmyX5IVtNsZb6L5WzmrrLgL2bZ/V3Vl8mu0fgIdm4E/hMPH3qPGu53R/3oB7Z0b8DTDan7g5B/hzknckWbMd55ZJnjxQ76xM/qFB3wD2TLJr+7p7G935PIPul413AW9u33dfCOwwsmOSvdI9uCss+ppcku9X0rQxiErSSqKqPkoXst5DF+R+RfeD+3Ftk4PoHpJxM93I3tH3bWVS/VxP93CUt9FNCf4XYK+qum4pyh9xEt3TNK+gm3Z2O5ObyjqeN9H9gHo1cDrdA1W+ANDusTya7r7S81k8TD+A7hh/Szcd8hlM8HCnMXyO7gEq84EL6YL8XUzyh8Sq+iPd9MB/rapf0T0c5l0susZvZ9H/7/eje8jL9XQPNDqa7gfZsdr+FvCfwFFtKuSldA+Rge4XAJ+jC12/aG2OjMy+AljQ9jmQ7t7V4bbvpHvC7B50I9OfBl5ZVf83meMexSl0QWUkiJ5ON9p26ph7wIeA97RpkwctYb//ku5+3RvorsP5wNNaQFkSh9Hdu3tya/cs4CljbTzBNZrIj+ju6fx9kun4+hzNt+mmsN9I97l4YS26L/yf6MLgn+g+m8eN7NQ+B18Hrm7X55FM/D3qYOBLbfuXDq5YBp+3wbbPa9PDh5ff3Y5vNt2fdbkO+Dzd7QTQPVAI4PoM/U3SMfq5nO5r6ROtrb+le/DVne34Xkg3TfhGunN+7MDumwE/pPuF2ZnAp6tq3lSOU5pu8T5lSZKWH0n2AD5TVcPTm5dFX0cD/1dV75twY2mKkhwMPK6q7vOLCElyRFSSpBnUpuw9t02n24juCbjj/ZmdpenryUk2Tfe3KnenGz09bln0JUnSeAyikiTNrND97cIb6abm/ozuft5l4RF0907eQvc3UF9XVRcuo74kSRqTU3MlSZIkSb1yRFSSJEmS1CuDqCRJkiSpV6vOdAGSYP31169Zs2bNdBmSJEnStDn//POvq6oNRltnEJWWA7NmzeK8886b6TIkSZKkaZPkF2Otc2quJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrg6gkSZIkqVerznQBkuCS39zErLknznQZkiRpHAsO2XOmS5BWGI6ISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXi1xEE1yS/vvrCT7Tl9JkORdQ+/PWMJ2ZiW5dHqqGrOP9ZK8fln2MdDXjkk+N7RssWNM8tokFyT5qyRHJPlNktXbuvWTLJhEPxOe7yQLkqw/yvKDkxw0qQOagiTzklye5OIk5yaZPd19SJIkSerHdIyIzgKmFESTrDLBJosF0ap62hRr6tN6wJSCaDpLcu53B74/TruvAN4E7FZVN7bFdwP7T6WTmTrfk/hc7FdVWwOfBj7cU5/TIsmqffQjSZIk3R9MRxA9BNg5yUVJ3ppklSQfbqNW85P8I0CSOUl+nORI4JK27Lgk5ye5LMkBbdkhwJqtva+1ZSOjr0cnee5Ix23E70Vj9Tkoyatbf99Jck2SNyb55yQXJjkryUPadvOSHJrkjCSXJtmhLV9spK+tm9WOf9NW74fburcP1PL+tmxWkp8l+TRwAbDxQFs7JDm2vX5+ktuSPDDJGkmuHjiMXYEfjnYRkrwUmEsXQq8bWHUo8NbRgtBodQ6d7wck+XS7Pick+W6SFw808aY2+npJkicMLN86yY+S/DzJa1tbadfo0rb9Pm35Yp+LJGslObGNfF46st2QM4GN2v5rJflCO44Lkzy/LX9Qkm+0Yzs6ydlJth85viQfSHI2sGOSlyc5p13Dz7bP0yrt8zVS71vbvm9O8tPW7lFt2UPaZ2t++yxt1ZYfnOTwJCcDXx7tukmSJEkro+kYpZkLHFRVewG0QHlTVT053ZTQn7QfxAF2ALasqmva+/2r6oYkawLnJjmmquYmeWNVzR6lr6OAfYDvJnkgXTB7HfD3Y/RZQ/tvCWwDrAFcCbyjqrZJ8jHglXShDWCtqnpakl2AL7T9xjv+LUfqTbIbsFk71gDHt3Z+CTweeE1VDY+gXtDqAtgZuBR4Mt31Obu1uz6wsKpuGqWGTYBPAttU1e+H1v0SOB14BfCdkYVj1VlVpw7s+0K6Ee8nAQ8DftbOx4jrqmrbdFOTDwL+oS3fCngqsBZwYZITgR2B2cDWwPp013ukr3s/F0leBPy2qvZsda47yvHuDhzXXr8b+FFV7Z9kPeCcJD+k+1zcWFVbJdkSuGhg/7WAS6vqvUm2AN4BPL2qFrZfFOwHXAZsVFVbtjrWa/vOBR5TVXcMLHs/cGFV7Z3kb+hC5+y2bjtgp6q6bZTjkCQBvz9y7kyXIE3KnLOmZUKW1It58+bNdAnjWhbTBXcDthoYOVuXLvDcCZwzEEIB3pzkBe31xm2768dp+3vAx1vY3B04tapua6FqtD6vGNr/x1V1M3BzkptYFMwuoQtPI74OUFWnJllnIHBMxm7t34Xt/dqtll8Cv6iqs4Z3qKq7klzZQtEOwEeBXYBVgNMG2j15eN/mj8ANwEuBj42y/j+A44ETJ1HnYBDdCfifqroH+H2SHw+1e2z77/l0oXXEt1vwuq3ts0Nr6+tVdTfwhySn0IXtP7P45+IS4CNJ/hM4oapOG2j3a0nWojsv2w4cx/OyaLR6DeDRrb/DAKrq0iTzB9q5Gzimvd6VLiyemwRgTeBaus/GY5N8op23kXM/v9VxHIvC8E7Ai1pfP0ry0IEAffxYIbT90uYAgFXW2WC0TSRJkqQV0rIIogHeVFUnLbYwmQP8Zej9s4Adq+rWJPPoQsSYqur2tt1z6EZGvz5Bn7OGmrhj4PU9A+/vYfFzMTySWsBdLD6VeaxaA3yoqj47Si1/GXWPzmnAHsBCuum3R9AFrpGAtQddQB3NrW396UmuraqvLVZ81ZVJLqILquPWOcqxjGfk/N3NxOdvvLbuPS9VdUWS7YDnAh9KcnJVfaCt3g+4mG469Kfowm+AF1XV5YsV3lLlGG5vgZi2/5eq6p3DGyXZmu6z9ga6c7c/sCfdLwmeB/xrkr8e49hGzsGY17yqDgcOB1h9w82Gz5kkrTQese8hM12CNCnzDtlzpkuQVhjTcY/ozcCDB96fBLwuyWoASTZvo1jD1qWbOnlru7/wqQPrFo7sP4qjgNfQTWEdCZ6T7XOyRu5f3Iluyu9NwALaKFySbYHHtG1HO/79k6zdtt0oycMm0eepwFuAM6vqj8BDgScAl7VQtRWLTy9dTNtnd+A/kjxnlE3+nUWhdrJ1ng68KN29og8H5kziOACen+7+1oe2fc5tx7dPu/dyA7owd87wjkkeCdxaVV8FPsKikc+R41wIvAd4ahtBPonuXtW0/UemOJ9OC95Jnkg3vXg0/wu8eOTY2/2em7Sp0A+oqmOAfwW2TfeAqY2r6sfAv9A9qGrtdmz7tf3n0E1Z/vMkz5UkSZK00pmOEdH5wF1JLqYbxTuM7r7CC1o4+COw9yj7fR84sE2ZvBwYnLJ6ODA/yQVVtd/QfiMPfjm+qu5syz4/yT4n68Z0f8JkHRY9cfYY4JVtZPFc2rTfqro+yU/S/QmV71XV21tAOrNlo1uAl9ONGo7nbODhLJoaOx+4tqoq3UN2LqyqcUfN2j2Wz6O7h/aFQ+suS3IBLdhV1clj1HntwG7H0E1dvbQd79nAaPeoDjuHbjrro4F/q6rfJvkW3X2iF9ONFv5LVf0+iz/kCLrA+OEk99CNDr9ulOO8Lcl/0wXrN9Ld2zu/XfsFwF50T9b9Uvt8XUh3Pu9Te1X9NMl7gJNb0FxINwJ6G/DFLHq68TvpRqi/2qbdBvhYVf0pycFt2/l0o9OvmsQ5kiRJklZamSDbrHTa1N+Dquq8ma5lRAtKV1bVUTPQ99pVdUsb3TyH7qE+ww9EWu6k+7Msq7Xp3JvSjXxuPvDLi+XK6htuVhu+6tCZLkOSJI1jgVNzpSlJcn5VbT/aOv+24f1AVX1wBrs/oT2s6YF0o5vLfQhtHgT8uE3XDvC65TWESpIkSSsbg+iQqpoz0zUsT+6v56M9HXnU375IkiRJmlnT8bAiSZIkSZImzSAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVz6sSFoOPGmjdTnPR8JLkiRpJeGIqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1atWZLkASXPKbm5g198SZLkOStBJZcMieM12CpJWYI6KSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlX0xpEk9zS/jsryb7T3Pa7ht6fsYTtzEpy6fRUNWYf6yV5/bLsY6CvHZN8bmjZrCS3Jblo4N8rx2ljiepNcnZr+5dJ/jjQ16wlOJSJ+rq7tX1pku8kWW+6+5AkSZLUj2U1IjoLmFIQTbLKBJssFkSr6mlTrKlP6wFTCnbpLMn12B34/ijLr6qq2QP/vjxOG+sxRr3jXZeqekpVzQbeCxw90NeCtu+qkzyGybittb0lcAPwhqVtcBKfuWkzzedCkiRJul9bVj8cHwJskeQi4EvAx9uyOcDqwKeq6rNJ5gDvA34HzAaemOQ4YGNgDeCwqjo8ySHAmq29y6pqvyS3VNXaSY4GvlRV3wVIcgTwHeC40focLDLJq4G9gVWALYH/Bh4IvAK4A3huVd2QZB5wEbADsA6wf1Wdk+Rg4Jaq+khr71Jgr9bvpq3eH1TV25O8HXhpq+VbVfW+NnL4PeDHwI6tll+0tnYA5lbVC5M8HzgKWJfulwc/rarHtsPYFfjoZC5Kkk2AH7a+bgBOAf4N2H+wXuDEyVyXMfo4GHgk3S8jrkvyT8BngEe3Td5SVT9JshbwCeBJdJ/Dg6vq20n+Gvgi3XV4APCiqvr5UDdnAlu1/jYFPgVsANwKvLaq/q8t/xrdtf0e8M/t8zJn8NiSPInRP5sbAkfTXe9VgdcBZwD/H7A9UMAXqupjSWa3Y3wQcBXd5+PG9rk5A3g6cDzd50vSCuz3R86d6RKkSZtz1odnugRpUubNmzfTJWgZWFZBdC5wUFXtBZDkAOCmqnpyktWBnyQ5uW27A7BlVV3T3u/fwt+awLlJjqmquUne2Ebfhh0F7AN8N8kD6YLZ64C/H6PPGtp/S2AbuoB1JfCOqtomyceAVwKHtu3WqqqnJdkF+ELbb7zj33Kk3iS7AZu1Yw1wfGvnl8DjgddU1fCI5AWtLoCdgUuBJ9Nds7Nbu+sDC6vqplFqGAmWI95UVacl+U+60HQ2XaA9OckVQ/XOYXLX5foxjn87YKequi3JkcDHqur0JI8GTgK2AN4N/Kiq9m/TbM9J8kPgQLqg+7V2PRcbtWyjmLvSBUKAw4EDq+rnSZ4CfBr4G+Cw1s7Xkxw4VN+9xzbOZ/OFwElV9e+tzwfRhfKN2qgsA9ODv9zO7ylJPkAXdN/S1q1XVc8Y7SS1vg8AWGWdDcY4lZIkSdKKp6/pgrsBWyV5cXu/Ll0wuxM4ZyDsALw5yQva643bdmMFHuhGuz7eQsTuwKktAI3V5xVD+/+4qm4Gbk5yE91oKsAltFG35usAVXVqknWmeI/ibu3fhe392q2WXwK/qKqzhneoqruSXJlkC7rg9FFgF7pgdtpAuycP79tcNVpwr6rPJ3kJXeC7z/oBS3Ndjq+q29rrZ9GNqI6sWyfJg1vtz0tyUFu+Bt2o6ZnAu5M8Cjh2YDR0ZER8FnA+8IMkawNPA/5noP3V239HRpgBjgQ+MsaxjfU5ORf4QpLVgOOq6qIkVwOPTfIJulHjk5OsSxc2T2n7fwn4n4G+jh7jHNFGlQ8HWH3DzYZ/QSLpfugR+x4y0yVIkzbvkD1nugRJK7G+gmjoRoxOWmxhN/L2l6H3zwJ2rKpb29TGNcZruKpub9s9h25k9OsT9DlrqIk7Bl7fM/D+HhY/P8NBoYC7WPw+27FqDfChUaYGz2Lg+EdxGrAHsJBuSu0RdEF0JLztwSSn5Q70+SDgUe3t2sDNY2y6NNdl8Jge0Pa7bXCDdMnxRVV1+dC+P0tyNrAncFKSf6iqH9HuEW3B7wS6e0SPAP40xkj5eAbrG/Vz0mrcpdXxlSQfrqovJ9ma7rP2Brqp1m+dQl+SJEmSWHYPK7oZePDA+5OA17XRJZJs3u4RHLYucGMLO08AnjqwbuHI/qM4CngN3RTWkUAx2T4na5/Wzk50UzlvAhYA27bl2wKPaduOdvz7txE8kmyU5GGT6PNUuimeZ1bVH4GHAk8ALmtBbiu6e1en4j/p7p18LzDytN3heoeNd10mcjLwxpE37X5K6M7Jm9pxkGSb9t/HAldX1cfp7qscHJWmnfc304Xx24Br2gjvyAOftm6bngW8qL3+u3HqG/Vz0u6nvbaqPkc3DXjbNhX6AVV1DPCvwLatnhuT7NzaewXdvbeSJEmSxrCsRkTnA3cluZhu1OowuimVF7Tg8UcWTZsc9H3gwCTzgcvpwsSIw4H5SS6oqv2G9juZ7j6946vqzrbs85Psc7JuTPcnY9ahe7gPwDHAK9uU0XNp036r6vokP0n38KLvtYcVbQGc2XLXLcDLgbsn6PNs4OF0gRS683ptVVWS7YELq2qsKZ3D94h+AbiY7j7Tp1fV3UlelOQ1VfXFwXrppp0OGu+6TOTNwKfavqu2YzmQ7iFJh9Jd09CF+r3oAv/LkywEfg98YLjBqrqwfbb+DtgP+H9J3gOsRvdLiYvpAvxXk7ytHc9o99HC2J+TOcDbWx230N0vvBHwxSx6uvE7239fBXymjTZfTfdLEUmSJEljyNg5RiPaVNSDquq8ma5lRAteV1bVUTNdy/KohcLbWmj/O+BlVfX8ma5rLKtvuFlt+KpDZ7oMSdJKZIH3iEpaxpKcX1Xbj7bOv214P1VVH5zpGpZz2wGfbKOcf2LRKLYkSZKkGWYQnYSqmjPTNWhqquo0YOsJN5QkSZLUu2X1sCJJkiRJkkZlEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrH1YkLQeetNG6nOdj9CVJkrSScERUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9WnWmC5AEl/zmJmbNPXGmy5Akabm04JA9Z7oESdPMEVFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktSrCYNoklvaf2cl2Xc6O0/yrqH3ZyxhO7OSXDo9VY3Zx3pJXr8s+xjoa8cknxtaNivJbUkuTPKzJOckedVS9PH5JE8cZ/0HkjxrCdp9TZKL2r87k1zSXh+ypLW2ducluTzJxUnOTTJ7adqTJEmSNHOmMiI6C5hSEE2yygSbLBZEq+ppU2m/Z+sBUwqi6SzJqPPuwPdHWX5VVW1TVVsAfwe8NclrlqB9quofquqn46x/b1X9cAna/WJVza6q2cBvgWe293NHtpnE52Is+1XV1sCngQ8vYRuLWYpaptrPqn30I0mSJN0fTCUkHQLs3Ea33ppklSQfbqNT85P8I0CSOUl+nORI4JK27Lgk5ye5LMkBbdkhwJqtva+1ZSOjr0cnee5Ix0mOSPKisfoclOTVrb/vJLkmyRuT/HMbSTwryUPadvOSHJrkjCSXJtmhLT84yUED7V2aZFY7/k1bvR9u694+UMv727JZbcTy08AFwMYDbe2Q5Nj2+vlthPOBSdZIcvXAYewKjBsCq+pq4J+BN7f21kryhVbPhUme35avkuQjbWRyfpI3DRz/9m39Ee04L0ny1oFz/uL2etfW5iWtj9Xb8gVJ3p/kgrbuCWPVm+SWNsp6NrBjkpe3Ud2Lknx2JBAm2S3Jma3N/0my9ijNnQlsNMFxPyjJN9oxH53k7CTbT7aWcc7Lm5P8tLV7VFv2kPaZm98+Y1u15QcnOTzJycCXx7uekiRJ0spkKqM0c4GDqmovgBYob6qqJ7dg8pP2AzfADsCWVXVNe79/Vd2QZE3g3CTHVNXcJG9sI2fDjgL2Ab6b5IF0wex1wN+P0WcN7b8lsA2wBnAl8I6q2ibJx4BXAoe27daqqqcl2QX4QttvvOPfcqTeJLsBm7VjDXB8a+eXwOOB11TV8AjqBa0ugJ2BS4En012Hs1u76wMLq+qmcWoZbG8k/L0b+FFV7Z9kPeCcJD9sx/sYYJuqumskiA+YDWxUVVu2/tcbXJlkDeAIYNequiLJl+muxaFtk+uqatt005YPAv5hjFrXAi6tqvcm2QJ4B/D0qlrYQvt+Sb4LvAd4VlX9Jck76ML2B4ba2h04boLjfh1wY1VtlWRL4KKp1AJcNsZ5mQs8pqruGFj2fuDCqto7yd/Qhc7Zbd12wE5VddsY50WSlju/P3LuxBtJPZpz1rRMhJKmxbx582a6hBXC0kwX3A3YamTUDFiXLpjdCZwzEEIB3pzkBe31xm2768dp+3vAx1vY3B04tapua+FvtD6vGNr/x1V1M3BzkpuA77TllwBbDWz3dYCqOjXJOsMhbAK7tX8Xtvdrt1p+Cfyiqs4a3qEFwStb+NkB+CiwC7AKcNpAuycP7zuGDNXzvCwazV0DeDTwLOAzVXVXq+GGoTauBh6b5BPAiaP0/XjgmqoaOcdfAt7AoiB6bPvv+cALx6n1buCY9npXuoB2bhKANYFrgacCT6T7BQPAA+lGP0d8LcladOdr2wmOeyfgsHbMlyaZP8VavjPGeZnf6jiORWF4J+BFra8fJXloknXbuuPHCqHtlzkHAKyyzgajbSJJkiStkJYmiAZ4U1WdtNjCZA7wl6H3zwJ2rKpbk8yjCwtjqqrb23bPoRsZ/foEfc4aauKOgdf3DLy/h8WPeXgktYC7WHzK8li1BvhQVX12lFr+MuoendOAPYCFdNNvj6ALViNBag+6gDoZ2wA/G6jnRVV1+VA94b7Hea+qujHJ1nTn+g3AS4H9B5uYoIaRc3s343+ebq+quwfa/FJVvXOo1r8FflBVLxujjf2Ai+mmSX+KLviOd9xLXEtrY7TzsifdLw+eB/xrkr9m9HM0cs7H/CxU1eHA4QCrb7jZmNdIkvr2iH2X6vly0rSbd8ieM12CpGk2lXtEbwYePPD+JOB1SVYDSLJ5G60ati7dFMlb2z2ETx1Yt3Bk/1EcBbyGbgrrSPCcbJ+TtU9rZye6Kb83AQtoo21JtqWb1gqjH//+I/cwJtkoycMm0eepwFuAM6vqj8BD6abXXtbC01YsPo10VC3wfgT4xEA9bxoJYElGpgCfDByY9rCc4am5bSrwA6rqGOBfWTTSOOL/gFlJHtfevwI4ZRLHOZ7/BV48cr7aPZabAGcBTx/pK919npsP7lhVC+mm7z61jSyPddyn04VH0j0d+ElTqWW085LuwVMbV9WPgX+he4DV2nTXdL+2/xy66cp/XorzI0mSJK3QpjIiOh+4K8nFdKN4h9E9SfeCFgL+COw9yn7fpwtC84HL6cLGiMOB+UkuqKr9hvYbecDL8VV1Z1v2+Un2OVk3pvuTMeuwaBTwGOCVSS4CzqVN+62q65P8JN2fifleVb29BaEzWwa6BXg53cjgeM4GHk4XXqA7r9dWVaV7mM6FVTXW6NimSS6kG6W9GfhEVX2xrfs3uumy89u5WQDsRXfONm/LFwKfAz450OZGwBez6Om+i40MttHp1wD/08LsucBnJjjGcVXVT5O8Bzi59bsQeENVnZXk1cDX27Rs6ELnFUP735bkv+lGkd84xnF/GvhS+9xdSHee73Pf7Vi1ALdx3/OyCvDVNu02wMeq6k9JDm7bzgduBZb4z+pIkiRJK4OMnXlWbG3q70FVdd5M1zKiBaIrq+qoma7l/i7dU3hXa0F6U7qRz80HfqmxXFl9w81qw1cdOtNlSJK0XFrg1FzpfinJ+VW1/Wjr/NuGy5Gq+uBM17ACeRDw4zaNO8DrltcQKkmSJK1sVtogWlVzZroGLTvtqcmj/vZFkiRJ0syaysOKJEmSJElaagZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvVppH1YkLU+etNG6nOej6SVJkrSScERUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknq16kwXIAku+c1NzJp74kyXIUlaQS04ZM+ZLkGSFuOIqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6tX9OogmuaX9d1aSfae57XcNvT9jCduZleTS6alqzD7WS/L6ZdnHQF87Jvnc0LJZSW5LclGSi5OckeTxPdXz3STrLeG+Byf5Tav7p0leNrDuA0meNco+c5KcsBQlS5IkSSu9+3UQHTALmFIQTbLKBJssFkSr6mlTrKlP6wFTCqLpLMn13x34/ijLr6qq2VW1NfAlhs7fslJVz62qPy1FEx+rqtnA84HPJlmttfveqvrhNJQoSZIkaciqM13ANDkE2CLJRXQh6ONt2RxgdeBTVfXZJHOA9wG/A2YDT0xyHLAxsAZwWFUdnuQQYM3W3mVVtV+SW6pq7SRHA1+qqu8CJDkC+A5w3Gh9DhaZ5NXA3sAqwJbAfwMPBF4B3AE8t6puSDIPuAjYAVgH2L+qzklyMHBLVX2ktXcpsFfrd9NW7w+q6u1J3g68tNXyrap6X5JZwPeAHwM7tlp+0draAZhbVS9M8nzgKGBdul9W/LSqHtsOY1fgoxNcj3WAG1u7s4CvAGu1dW+sqjOSfAX4ZlV9u233NeBo4MTRzmOSDdv6deg+t6+rqtOSLAC2r6rrRruWre1bgMPauboNeH5V/WGw4Kr6eZJbgb8Crm3X9YSq+maS3YFDgeuAC0b2SbIBcCTwUOBcupC+Xavl5cCb6a7v2cDrq+ruCc6bJK2wfn/k3JkuYaU256wPz3QJK7V58+bNdAnScmdFGRGdC5zWRuQ+Bvw9cFNVPRl4MvDaJI9p2+4AvLuqntje719V2wHbA29O8tCqmgvc1trbb6ivo4B9AJI8kC6YfXeCPgdtSTd6uwPw78CtVbUNcCbwyoHt1mqjsK8HvjCJ4x8ZkXx7kt2AzVofs4HtkuzStn088OWq2qaqfjHQxgXANu31zsCl7TieQhekSLI+sLCqbhqlhk3bFNergH9mUVi9Fnh2VW1Ld94+3pZ/HnhNa3dd4GmMfx73BU5qo5db0wX1Yfe5lm35WsBZbbT2VOC1wzsm2Rb4eVVdO7R8DeBzwN+28/KIgdXvA37Uju1bwKPbPlu0Y316q/duYPhzRJIDkpyX5Ly7bx3tlEqSJEkrphVlRHTYbsBWSV7c3q9LF8zuBM6pqmsGtn1zkhe01xu37a4fp+3vAR9PsjrdCNipVXVbC3+j9XnF0P4/rqqbgZuT3EQ3mgpwCbDVwHZfB6iqU5OsM8X7IHdr/y5s79dutfwS+EVVnTW8Q1XdleTKFqJ2oAuSu9CN3p420O7JY/R5VQtdJNkHOJzu/KwGfDLJbLpAtnnr75Qkn0ryMOCFwDGthrHO47nAF9rU2eOq6qJRahjrWt4JjNzXeT7w7IF93prktcBjW73DngBcU1U/b8f2VeCAtm4n4AXteL6f5Ma2fFdgO+DcJABr0gXyxbQR28MBVt9wsxqlb0laYTxi30NmuoSV2rxD9pzpEiRpMStqEA3wpqo6abGF3dTcvwy9fxawY1Xd2qbErjFew1V1e9vuOXSjXl+foM9ZQ03cMfD6noH397D49RgOJgXcxeKj2GPVGuBDo0wNnsXA8Y/iNGAPYCHwQ+AIuiB6UFu/BxNPywU4Hvhie/1W4A90o5gPAG4f2O4rdCOFfwfsP1D7fc5jq38XYE/gK0k+XFVfHlg3h7Gv5cKqGjmfd7P4ef5YVX0kyQuBLyfZtKoGa4T7Xot7ux1n+Zeq6p1jrJckSZJWaivK1NybgQcPvD8JeN3Ig2eSbJ5krVH2Wxe4sQWXJwBPHVi3cGT/URxFN61059bXVPqcrJHpvzvRTVW9CVgAbNuWbwuMTP0d7fj3T7J223ajNvI4kVOBtwBnVtUf6e59fAJwWbqhva0YfUrssJ2Aq9rrdYHfVdU9dPfCDj4k6ojWH1V12UDt9zmPSTYBrq2qzwH/H+08DBjvWk6oqo4FzgNeNbTq/4DHJNm0vX/ZwLrT6e7DpY3k/lVb/r/Ai0fOeZKHtPolSZIkseKMiM4H7kpyMV24OYzuSboXtAD1R7oH8wz7PnBgkvnA5cDglNXDgflJLhjlPtGTgS8Dx1fVnW3Z5yfZ52TdmO5PxqzDotHCY4BXtocSnUub9ltV1yf5SXt40ffafaJbAGe2qaG3AC+nGw0cz9nAw+kCKXTn9dqqqiTbAxcOjCwOG3lYUuimwv5DW/5p4JgkL6F7SNK9I7JV9YckP6N70NOIsc7jHODtSRa24xm8nxbGv5aT9QHgyAz8eZo2An4AcGKS6+jC55Zt9fuBr7epyKfQPQTr5vawovcAJ6d7MvFC4A20B0NJkiRJK7uMnSs0U9q00oOq6ryZrmVEC1ZXVtVR09jmg+jujd12jAcgLdfafcJ3t3tbdwT+38h9slO1+oab1YavOnQ6y5Mk6V4LvEdU0gxIcn5VbT/auhVlRFTLWFV9cDrbS/IsuqcBf/T+GEKbRwPfaKOedzLK03glSZIk3ZdBdDlUVXNmuoZlrap+SPtzJ/dX7Um620y4oSRJkqTFrCgPK5IkSZIk3U8YRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKhxVJy4EnbbQu5/lofUmSJK0kHBGVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF6tOtMFSIJLfnMTs+aeONNlSJK0zCw4ZM+ZLkHScsQRUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuD6CiSvDvJZUnmJ7koyVPa8s8neeI09XHLdGwzDXW8a1n30frZMMnJoyyvJF8ZeL9qkj8mOaG9f3V7f1H79+Ux2n9LkldOQ50bJDk7yYVJdl7KtvZK8v6lrUmSJEla0RhEhyTZEdgL2LaqtgKeBfwKoKr+oap+OpP1LQNTDqJJVlmCfnYHThpl+V+ALZOs2d4/G/jN0DZHV9Xs9u8+YTPJqsD+wJFLUNewXYH/q6ptquq0oX6metwnAs9L8qBpqEuSJElaYaw60wUshzYErquqOwCq6rqRFUnmAQdV1XlttPJTdEH1RrpA91/Ao4G3VNXxSV4NvABYHXgMcGRV3WeELMnbgZe27b5VVe8bWj8HeD/wB2A2cCxwCfBPwJrA3lV1VZIjgNuBvwYeDvxzVZ3Q6ti+qt7Y2jsB+AhdOFwzyUXAZVW1X5KXA28GHgicDby+qu5ux/tR4DnA24DTW1sPA75XVdsl2Rq4CNikqn6Z5CrgSVV1a+trrNHB7wF7At8EXgZ8HZjKaOTfABdU1V2tpscBnwE2AO4GXgJsPNE5BB5Mdw1HzsmOwB8HjzvJrDHOz+7AfwCr0H1+dq2qap+ZvYBvTOF4JK3kfn/k3JkuQZp2c8768EyXIE27efPmzXQJ91uOiN7XycDGSa5I8ukkzxhju7WAeVW1HXAz8EG60bwXAB8Y2G4HYD+68POSJNsPNpJkN2Cztt1sYLsku4zS39Z0oelJwCuAzatqB+DzwJsGtpsFPIMu2H0myRpjHWhVzQVuayON+yXZAtgHeHpVzaYLcfsNHO+lVfWUqjp9oI1rgTWSrEMXHs8Ddk6yCXBtVd3aRhIfP85o8lHA37Vat6ILeIP2GZia+5pR9n86cP7A+68Bn6qqrYGnAb9ry8c9h1V1EfBeFo3A3jZ43MD1o52fJBsAnwNe1Pp8yUAt5zFGqE5yQJLzkpx39603jXFqJEmSpBWPI6JDquqWJNvRhYdnAkcnmVtVRwxteifw/fb6EuCOqlqY5BK6MDjiB1V1PUCSY4Gd6MLJiN3avwvb+7XpgumpQ/2dW1W/a+1cRReYR/p+5sB236iqe4CfJ7kaeMJkj51uWup2wLlJoBspvLatuxs4Zoz9zqALg7vQjQruDgQYmdr6FO4bLu9VVfPbSOPLgO+OssnRI6O5Y9gQ+BlAkgcDG1XVt1rbt7flMPlzOGjwuMc6P08FTq2qa1qfNwzsfy3wyDGO+3DgcIDVN9ysxjk+SSuZR+x7yEyXIE27eYfsOdMlSFqOGERHUVV3A/OAeS1Yvgo4YmizhVU1Eh7uAUam8t7T7lm8t7nh5ofeB/hQVX12grLuGHh9z8D7e1j8Oo7W310sPvo91ihpgC9V1TtHWXd7Oy+jOY0uuG8CfBt4R+v3hLZ+DxaF9rEcTzddeA7w0Am2HXYbi44p42w32XM4aPC4Rz0/SZ7Hfc/7iDVafZIkSZIap+YOSfL4JJsNLJoN/GIpmnx2koe0h/HsDfxkaP1JwP5J1m79b9Tuu1xSL0nygCSbAo8FLgcWALPb8o3ppgGPWJhktfb6f4EXj/Tf6t5kEn2eCrwc+Hkbjb0BeC6LjnXX1vZ4vgB8oKoumUR/w34GPA6gqv4M/DrJ3gBJVp/GhwWNdX7OBJ6R5DEjywf22Ry4dJr6lyRJklYIjoje19rAJ5KsRzeSeCVwwFK0dzrwFbqgdGRVDU7LpapObvdmntmme95CF+quHW5oki4HTqF7WNGBVXV7kp8A19BNQb0UuGBg+8OB+UkuaPeJvgc4OckDgIXAG5ggiFfVglb7yHTi04FHVdWN7f7J21tAHK+NXwOHTfFYR3yP7hyPeAXw2SQfoDuGl4y61xRV1U9HOz9VdVaSA4Bj2/Jr6e4Xhm7K72gjzJIkSdJKK4tml2q6DT+ttof+jgBOqKpv9tHfZLSn8D6qqpbpDU9JvgX8S1X9fFn2MxVJHk73y4ddJ9p29Q03qw1fdeiyL0qSpBmywHtEpZVOkvOravvR1jkiqmWqqr7aU1dz6R5atNwEUbo/5fO2mS5CkiRJWt4YRJeh9qTdI3rs79V99bW8qarL6aYlLzeq6tyZrkGSJElaHvmwIkmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSeuXDiqTlwJM2WpfzfKy9JEmSVhKOiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrg6gkSZIkqVerznQBkuCS39zErLknznQZkiRpDAsO2XOmS5BWKI6ISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEe5Tk3UkuSzI/yUVJntKWfz7JE6epj1umY5tpqONdy7qP1s+GSU5exn3MS7J9e/3dJOsty/4kSZKkFd2qM13AyiLJjsBewLZVdUeS9YEHAlTVP8xoccvGu4D/mMoOSVapqrun2M/uwElT3GeJVdVz++pLkiRJWlE5ItqfDYHrquoOgKq6rqp+C/cZcbslyX8mOT/JD5Ps0NZfneR5bZtXJ/l2ku8nuTzJ+0brMMnbk5zbRmDfP8r6OUlOSfKNJFckOSTJfknOSXJJkk3bdkck+UyS09p2ew3U8cmB9k5obR4CrNlGfb/W1r28tXtRks8mWWXgeD+Q5Gxgx4G2Hpbk/PZ66ySV5NHt/VVJHtQ23R343tBxzUryf22k+dIkX0vyrCQ/SfLzJDu07dZK8oV2ji5M8vy2fM0kR7XzdjSw5kDbC9ovEUhyXLtOlyU5YGCbW5L8e5KLk5yV5OHjfTAkSZKklY0jov05GXhvkiuAHwJHV9Upo2y3FjCvqt6R5FvAB4FnA08EvgQc37bbAdgSuBU4N8mJVXXeSCNJdgM2a9sFOD7JLlV16lB/WwNbADcAVwOfr6odkvwT8CbgLW27WcAzgE2BHyd53FgHWlVzk7yxqma3WrYA9gGeXlULk3wa2A/4cjveS6vqvUNtXJtkjSTrADsD5wE7JzkduLaqbm1h9vFV9dNRyngc8BLgAOBcYF9gJ+B5dKO1ewPvBn5UVfu36bbnJPkh8I/ArVW1VZKtgAvGONT9q+qGJGvSXYNjqur6dkxnVdW7k/wX8Fq66yhJmka/P3LuTJeglcicsz480yVoJTJv3ryZLmGZM4j2pKpuSbIdXah6JnB0krlVdcTQpncC32+vLwHuaOHtErowOOIHLfSQ5Fi6kHXewPrd2r8L2/u16YLpcBA9t6p+19q5ii4wj/T9zIHtvlFV9wA/T3I18ITJHjuwK7AdXViDboTx2rbubuCYMfY7A3g6sAvdNN/d6UL1aW39U4Czx9j3mqq6pB3XZcD/VlUNncfdgOclOai9XwN4dOvv4wBVNT/J/DH6eHOSF7TXG9Od3+vpruEJbfn5dL9IuI82inoAwCrrbDBGF5IkSdKKxyDao3b/4zxgXgtErwKOGNpsYVVVe30PMDKV954kg9erhvYbfh/gQ1X12QnKumPg9T0D7+9h8c/HaP3dxeLTu9cYo48AX6qqd46y7vZx7gs9jS64bwJ8G3hH63ck5O3BotA+bDLHFeBFVXX5YsV2YXn4eBnaZg7wLGDHNjo7j0XHP3gN72aMr7OqOhw4HGD1DTcbtz9J0n09Yt9DZroErUTmHbLnTJcgrVC8R7QnSR6fZLOBRbOBXyxFk89O8pA2LXRv4CdD608C9k+ydut/oyQPW4r+XpLkAe2+0ccClwMLgNlt+cZ004BHLEyyWnv9v8CLR/pvdW8yiT5PBV4O/LyNxt4APJdFx7pra3tJnQS8KS15JtlmoN/92rItga1G2Xdd4MYWQp8APHUp6pAkSZJWKo6I9mdt4BPtXsS7gCtp0zKX0OnAV+juhTxy8P5QgKo6ud2beWbLWbfQhbprhxuapMuBU4CHAwdW1e1JfgJcQzeN91IWv5fycGB+kguqar8k7wFOTvIAYCHwBiYI4lW1oNU+Mp34dOBRVXVjkg3oRlP/vITHA/BvwKGtztAF672A/wd8sU3JvQg4Z5R9vw8c2La5HDhrKeqQJEmSVipZNINQ9xdJXg1sX1Vv7Km/I4ATquqbffQ3GUleThdKV4h5WatvuFlt+KpDZ7oMSZI0hgVOzZWmLMn5VbX9aOscEdX9UlV9daZrkCRJkrRkDKL3Q+1Ju0f02N+r++pLkiRJ0orPhxVJkiRJknplEJUkSZIk9cogKkmSJEnqlUFUkiRJktQrH1YkLQeetNG6nOdj4SVJkrSScERUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknq16kwXIAku+c1NzJp74kyXIUmSBiw4ZM+ZLkFaYTkiKkmSJEnqlUFUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQXQ4leXeSy5LMT3JRkqe05Z9P8sRp6uOW6dhmGup417Luo/WzYZKT++hLkiRJ0vhWnekCtLgkOwJ7AdtW1R1J1gceCFBV/zCjxS0b7wL+Yyo7JFmlqu6eYj+7AydNcR9JkiRJy4BBdPmzIXBdVd0BUFXXjaxIMg84qKrOa6OVnwKeBdxIF+j+C3g08JaqOj7Jq4EXAKsDjwGOrKr3D3eY5O3AS9t236qq9w2tnwO8H/gDMBs4FrgE+CdgTWDvqroqyRHA7cBfAw8H/rmqTmh1bF9Vb2ztnQB8hC4crpnkIuCyqtovycuBN9OF77OB11fV3e14Pwo8B3gbcHpr62HA96pquyRbAxcBm1TVL5NcBTypqm5tfS127ElmAd9vbT0VuBj4YtvuYcB+VXVOkrWATwBPovuaObiqvt32/wqwVmvyjVV1RjtfBwPXAVsC5wMvr6oaPveStDL6/ZFzZ7oEaVLmnPXhmS5BmtC8efNmuoQl4tTc5c/JwMZJrkjy6STPGGO7tYB5VbUdcDPwQeDZdMHzAwPb7QDsRxcgX5Jk+8FGkuwGbNa2mw1sl2SXUfrbmi54Pgl4BbB5Ve0AfB5408B2s4BnAHsCn0myxlgHWlVzgduqanYLoVsA+wBPr6rZwN2t9pHjvbSqnlJVpw+0cS2wRpJ1gJ2B84Cdk2wCXFtVtyZZBXh8Vf10lDIeBxwGbAU8AdgX2Ak4iC7cA7wb+FFVPRl4JvDhFk6vBZ5dVdu2uj8+0O42wFuAJwKPBZ4+3HGSA5Kcl+S8u2+9aazTJEmSJK1wHBFdzlTVLUm2owtVzwSOTjK3qo4Y2vROutE86EYn76iqhUkuoQuDI35QVdcDJDmWLmSdN7B+t/bvwvZ+bbpgeupQf+dW1e9aO1fRBeaRvp85sN03quoe4OdJrqYLd5O1K7AdcG4S6EZbr23r7gaOGWO/M+iC3i5003x3BwKc1tY/hW50dTTXVNUl7bguA/63qmroPO4GPC/JQe39GnQjz78FPplkdqtv84F2z6mqX7d2L2ptnT6wnqo6HDgcYPUNN3O0VNJK4xH7HjLTJUiTMu+QPWe6BGmFZRBdDrX7H+cB81ogehVwxNBmCwemet4DjEzlvSfJ4HUdDjjD7wN8qKo+O0FZdwy8vmfg/T0s/jkarb+7WHz0faxR0gBfqqp3jrLu9nHuCz2NLrhvAnwbeEfr94S2fg8WhfZhkzmuAC+qqssXKzY5mG668tZ0x3f7GO3ejV9rkiRJ0r2cmrucSfL4JJsNLJoN/GIpmnx2kockWRPYG/jJ0PqTgP2TrN3636jdd7mkXpLkAUk2pZuSejmwAJjdlm9MNw14xMIkq7XX/wu8eKT/Vvcmk+jzVODlwM/baOwNwHNZdKy7traX1EnAm9KGaZNs05avC/yu9fkKYJWl6EOSJElaaThKs/xZG/hEkvXoRhKvBA5YivZOp3ugzuPoHlY0OC2Xqjq53Zt5ZstZt9CFumuHG5qky4FT6B5WdGBV3Z7kJ8A1dNN4LwUuGNj+cGB+kgvafaLvAU5O8gBgIfAGJgjiVbWg1T4ynfh04FFVdWOSDehGU/+8hMcD8G/Aoa3O0AXrvYBPA8ckeQnwY+AvS9GHJEmStNKID/JccQ0/rbaH/o4ATqiqb/bR32S0p/A+qqqW6xuSVt9ws9rwVYfOdBmSJGnAAu8RlZZKkvOravvR1jkiqhVaVX11pmuQJEmStDiD6AqsPWn3iB77e3VffUmSJEm6//JhRZIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9cqHFUnLgSdttC7n+Yh4SZIkrSQcEZUkSZIk9cogKkmSJEnqlUFUkiRJktQrg6gkSZIkqVcGUUmSJElSrwyikiRJkqReGUQlSZIkSb0yiEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXq060wVIgkt+cxOz5p4402VIkpbSgkP2nOkSJOl+wRFRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBdAWV5FFJvp3k50muSnJYkgeOse0jk3xzEm1+N8l6S1jPwUkOWsJ9V0ty/ijLFyS5JMlF7d/TxmljTpIT2utXJ/nkKNu8Oskfk1zYzttJ47U5sN/eSZ441eOSJEmSVlYG0RVQkgDHAsdV1WbA5sDawL+Psu2qVfXbqnrxRO1W1XOr6k/TXe8k7AScMca6Z1bV7PZvrG2m4uiq2qadt0OAY5NsMcE+ewMGUUmSJGmSDKIrpr8Bbq+qLwJU1d3AW4H9kzyojfz9T5LvACcnmZXkUoC2/htJ5ic5OsnZSbZv6xYkWb9t/7Mkn0tyWZKTk6zZtnltknOTXJzkmCQPGqvIJKskuTqd9ZLck2SXtu60JI9rm+4OfG8yB55k3kC96ydZsATnD4Cq+jFwOHDAWMfWRkyfB3y4jcpuOpVzIEmSJK2MVp3pArRM/DWw2FTWqvpzkl8CI+FuR2CrqrohyayBTV8P3FhVWyXZErhojD42A15WVa9N8g3gRcBXgWOr6nMAST4I/D3widEaqKq7k1xBN5r4mFbzzknOBh5VVVe2TZ8JvH+MOn6c5G7gjqp6yhjbLI0LgH9sr+9zbFX1iSTHAydU1Tfbuj9N9hxIWnq/P3LuTJcg3WvOWR+e6RKkxcybN2+mS5BGZRBdMQWoCZb/oKpuGGWbnYDDAKrq0iTzx+jjmqq6qL0+H5jVXm/Zwtd6dNOBT5qg1tOAXeiC6IeA1wKnAOdCd/8qcENV3TrG/s+squsm6GNpZOD1ZI9tUtslOYA22rrKOhtMU7mSJEnS8s8gumK6jG6E8l5J1gE2Bq4CtgP+Msa+GWP5sDsGXt8NrNleHwHsXVUXJ3k1MGeCdk4DDgQeCbwXeHvb59S2fg8mDrOD7mLRlPM1prDfWLYBftZeH8Hkjm1S21XV4XRTf1l9w81G+8WBpEl4xL6HzHQJ0r3mHbLnTJcgSfcL3iO6Yvpf4EFJXgndvZjAfwNHjDOyOOJ04KVtvycCT5pi3w8GfpdkNWC/SWx/NvA04J6qup1uKvA/0gVUmML9oc0CuqANMOEDmMaT5Bl0I5afa4vGOrab2zom2E6SJEkSBtEVUlUV8ALgJUl+DlwB3A68axK7fxrYoE3JfQcwH7hpCt3/K124/AHwf5Oo9Q7gV8BZbdFpdEHukhagN6uqCdsZ8BHgdUnOANafwn4j9mkPHbqC7ny9qKpGRkTHOrajgLe3P/uy6TjbSZIkSQLSZRap08LfalV1ewtV/wtsXlV3zkAtOwEvr6oD++67b6tvuFlt+KpDZ7oMSdJSWuDUXEm6V5Lzq2r70dZ5j6iGPYjuSbSr0d0v+rqZCKEAVXU63VRhSZIkSSsQg6gWU1U3A6P+1kKSJEmSpoP3iEqSJEmSemUQlSRJkiT1yiAqSZIkSeqVQVSSJEmS1CsfViQtB5600bqc5yP/JUmStJJwRFSSJEmS1CuDqCRJkiSpVwZRSZIkSVKvDKKSJEmSpF4ZRCVJkiRJvTKISpIkSZJ6ZRCVJEmSJPXKICpJkiRJ6pVBVJIkSZLUK4OoJEmSJKlXBlFJkiRJUq9SVTNdg7TSS3IzcPlM16Fptz5w3UwXoWnndV0xeV1XTF7XFZPX9f5jk6raYLQVq/ZdiaRRXV5V2890EZpeSc7zuq54vK4rJq/risnrumLyuq4YnJorSZIkSeqVQVSSJEmS1CuDqLR8OHymC9Ay4XVdMXldV0xe1xWT13XF5HVdAfiwIkmSJElSrxwRlSRJkiT1yiAq9SjJ7kkuT3JlkrmjrE+Sj7f185NsOxN1amomcV33a9dzfpIzkmw9E3Vqaia6rgPbPTnJ3Ule3Gd9WjKTua5J5iS5KMllSU7pu0ZN3SS+D6+b5DtJLm7X9TUzUacmL8kXklyb5NIx1vsz0/2cQVTqSZJVgE8BewBPBF6W5IlDm+0BbNb+HQD8v16L1JRN8rpeAzyjqrYC/g3vbVnuTfK6jmz3n8BJ/VaoJTGZ65pkPeDTwPOq6q+Bl/Rdp6Zmkl+vbwB+WlVbA3OA/07ywF4L1VQdAew+znp/ZrqfM4hK/dkBuLKqrq6qO4GjgOcPbfN84MvVOQtYL8mGfReqKZnwulbVGVV1Y3t7FvConmvU1E3m6xXgTcAxwLV9FqclNpnrui9wbFX9EqCqvLbLv8lc1wIenCTA2sANwF39lqmpqKpT6a7TWPyZ6X7OICr1ZyPgVwPvf92WTXUbLV+mes3+HvjeMq1I02HC65pkI+AFwGd6rEtLZzJfr5sDf5VkXpLzk7yyt+q0pCZzXT8JbAH8FrgE+Kequqef8rSM+DPT/dyqM12AtBLJKMuGH1s9mW20fJn0NUvyTLogutMyrUjTYTLX9VDgHVV1dzfIovuByVzXVYHtgF2BNYEzk5xVVVcs6+K0xCZzXZ8DXAT8DbAp8IMkp1XVn5dxbVp2/Jnpfs4gKvXn18DGA+8fRfeb2aluo+XLpK5Zkq2AzwN7VNX1PdWmJTeZ67o9cFQLoesDz01yV1Ud10uFWhKT/T58XVX9BfhLklOBrQGD6PJrMtf1NcAh1f3dwiuTXAM8ATinnxK1DPgz0/2cU3Ol/pwLbJbkMe0BCX8HHD+0zfHAK9uT4J4K3FRVv+u7UE3JhNc1yaOBY4FXOKpyvzHhda2qx1TVrKqaBXwTeL0hdLk3me/D3wZ2TrJqkgcBTwF+1nOdmprJXNdf0o1yk+ThwOOBq3utUtPNn5nu5xwRlXpSVXcleSPd0zVXAb5QVZclObCt/wzwXeC5wJXArXS/wdVybJLX9b3AQ4FPt9Gzu6pq+5mqWROb5HXV/cxkrmtV/SzJ94H5wD3A56tq1D8foeXDJL9e/w04IskldFM631FV181Y0ZpQkq/TPeF4/SS/Bt4HrAb+zLSiSDdDQZIkSZKkfjg1V5IkSZLUK4OoJEmSJKlXBlFJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSYtJUkm+MvB+1SR/THJCe//wJCckuTjJT5N8ty2fleS2JBcN/HvlGH18M8ljl+ExLEiyfnt9y1K29fkkT5zktkckefEU2p6VZN+B97OTPHcJatyg/f1TSbpfWHWmC5AkScudvwBbJlmzqm4Dng38ZmD9B4AfVNVhAEm2Glh3VVXNHq/xJH8NrFJVV09v2UsvSej+zvo9I8uq6h+WYZezgH2BI9v72cD2wHcn20CSVavqj0l+l+TpVfWTaa9SkqaZI6KSJGk03wP2bK9fBnx9YN2GwK9H3lTV/Cm2vR/wbYAkq7RRxEuTXJLkrW35vCQfS3Jqkp8leXKSY5P8PMkHRxpKclyS85NcluSAiTpO8vYk5yaZn+T9bdms1sengQuAjYf2mZdk+7FqHcWzkpyW5Iokew0c54cH+v7Htu0hwM5t9PgddCF/n/Z+nyRrJflC2+/CJM9v7b06yf8k+Q5wcmvruHZuJWm554ioJEkazVHAe9t03K2ALwA7t3WfAo5O8kbgh8AXq+q3bd2mSS4aaOdNVXXaUNtPZ1GwnQ1sVFVbAiRZb2C7O6tqlyT/RBdctwNuAK5K8rGquh7Yv6puSLImcG6SY9ry+0iyG7AZsAMQ4PgkuwC/BB4PvKaqXj/OORmv1kGzgGcAmwI/TvI44JXATVX15CSrAz9JcjIwFzioqkYC6x+A7avqje39fwA/qqr9W3/nJPlh62dHYKuquqG9Pw+4N6RL0vLMICpJku6jquYnmUU3GvrdoXUntfs7dwf2AC5MsmVbPeHUXLoR1T+211cDj03yCeBEFo3uARzf/nsJcFlV/Q4gydV0o5bXA29O8oK23cZ0QXPUIArs1v5d2N6v3bb/JfCLqjprgrrHq3XQN9rU3p+3Wp/Q+t1q4P7RdVvfd07Q527A85Ic1N6vATy6vf7BQAgFuBZ45ATtSdJywSAqSZLGcjzwEWAO8NDBFS0AHQkc2UZNdwHOn2S7t9EFKqrqxiRbA88B3gC8FNi/bXdH++89A69H3q+aZA7wLGDHqro1ybyRdscQ4ENV9dnFFnaB+y8TFT1BrYttOsr70I0OnzTU95wJug3woqq6fGi/p4xS8xp051aSlnveIypJksbyBeADVXXJ4MIkf5PkQe31g+mmoP5yCu3+DHhc23994AFVdQzwr8C2U2hnXeDGFkKfADx1gu1PAvZPsnbre6MkD5tsZ1Oo9SVJHpBkU+CxwOWt79clWa21tXmStYCbgQcP7Dv8/iTgTe0hSiTZZpwSNwcunezxSNJMckRUkiSNqqp+DRw2yqrtgE8muYvul9qfr6pz28ji8D2iX6iqjw/tfyLdKOsPgY2ALyYZ+eX4O6dQ4veBA5PMpwt7406traqTk2wBnNly3S3Ay4G7J9nfZGu9HDgFeDhwYFXdnuTzdPeOXtBC5R+BvYH5wF1JLgaOAL4EzG3n8EPAvwGHAvPbfguAvcbo95l051aSlnupGp49IkmStOy0Bwv9GHh6VU02BGoCSU4Fnl9VN850LZI0EYOoJEnqXZLnAD+rqqlM6dUYkmxAF+yPm+laJGkyDKKSJEmSpF75sCJJkiRJUq8MopIkSZKkXhlEJUmSJEm9MohKkiRJknplEJUkSZIk9er/B3ZDRSyVLOEWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 936x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "### As of 2021-07-21, SciKit-Learn's example throws lots of convergence warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "# To use this experimental feature, we need to explicitly ask for it:\n",
    "from sklearn.experimental import enable_iterative_imputer  # noqa\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "from sklearn.linear_model import BayesianRidge\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.ensemble import ExtraTreesRegressor\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "N_SPLITS = 5\n",
    "\n",
    "rng = np.random.RandomState(0)\n",
    "\n",
    "X_full, y_full = fetch_california_housing(return_X_y=True)\n",
    "# ~2k samples is enough for the purpose of the example.\n",
    "# Remove the following two lines for a slower run with different error bars.\n",
    "X_full = X_full[::10]\n",
    "y_full = y_full[::10]\n",
    "n_samples, n_features = X_full.shape\n",
    "\n",
    "# Estimate the score on the entire dataset, with no missing values\n",
    "br_estimator = BayesianRidge()\n",
    "score_full_data = pd.DataFrame(\n",
    "    cross_val_score(\n",
    "        br_estimator, X_full, y_full, scoring='neg_mean_squared_error',\n",
    "        cv=N_SPLITS\n",
    "    ),\n",
    "    columns=['Full Data']\n",
    ")\n",
    "\n",
    "# Add a single missing value to each row\n",
    "X_missing = X_full.copy()\n",
    "y_missing = y_full\n",
    "missing_samples = np.arange(n_samples)\n",
    "missing_features = rng.choice(n_features, n_samples, replace=True)\n",
    "X_missing[missing_samples, missing_features] = np.nan\n",
    "\n",
    "# Estimate the score after imputation (mean and median strategies)\n",
    "score_simple_imputer = pd.DataFrame()\n",
    "for strategy in ('mean', 'median'):\n",
    "    estimator = make_pipeline(\n",
    "        SimpleImputer(missing_values=np.nan, strategy=strategy),\n",
    "        br_estimator\n",
    "    )\n",
    "    score_simple_imputer[strategy] = cross_val_score(\n",
    "        estimator, X_missing, y_missing, scoring='neg_mean_squared_error',\n",
    "        cv=N_SPLITS\n",
    "    )\n",
    "\n",
    "    \n",
    "##### NEW ADDITION HERE #########\n",
    "# This is the piece of code that adds imputations with cmfrec\n",
    "from cmfrec import CMF_imputer\n",
    "estimator = make_pipeline(\n",
    "    CMF_imputer(k=10, user_bias=False, lambda_=[0,0,10,10,10,10],\n",
    "                verbose=False),\n",
    "    br_estimator\n",
    ")\n",
    "score_simple_imputer[\"MF (cmfrec)\"] = cross_val_score(\n",
    "    estimator, X_missing, y_missing, scoring='neg_mean_squared_error',\n",
    "    cv=N_SPLITS\n",
    ")\n",
    "##### END OF NEW ADDITION #########\n",
    "\n",
    "\n",
    "# Estimate the score after iterative imputation of the missing values\n",
    "# with different estimators\n",
    "estimators = [\n",
    "    BayesianRidge(),\n",
    "    DecisionTreeRegressor(max_features='sqrt', random_state=0),\n",
    "    ExtraTreesRegressor(n_estimators=10, random_state=0),\n",
    "    KNeighborsRegressor(n_neighbors=15)\n",
    "]\n",
    "score_iterative_imputer = pd.DataFrame()\n",
    "for impute_estimator in estimators:\n",
    "    estimator = make_pipeline(\n",
    "        IterativeImputer(random_state=0, estimator=impute_estimator),\n",
    "        br_estimator\n",
    "    )\n",
    "    score_iterative_imputer[impute_estimator.__class__.__name__] = \\\n",
    "        cross_val_score(\n",
    "            estimator, X_missing, y_missing, scoring='neg_mean_squared_error',\n",
    "            cv=N_SPLITS\n",
    "        )\n",
    "\n",
    "scores = pd.concat(\n",
    "    [score_full_data, score_simple_imputer, score_iterative_imputer],\n",
    "    keys=['Original', 'SimpleImputer', 'IterativeImputer'], axis=1\n",
    ")\n",
    "\n",
    "# plot california housing results\n",
    "fig, ax = plt.subplots(figsize=(13, 6))\n",
    "means = -scores.mean()\n",
    "errors = scores.std()\n",
    "means.plot.barh(xerr=errors, ax=ax)\n",
    "ax.set_title('California Housing Regression with Different Imputation Methods')\n",
    "ax.set_xlabel('MSE (smaller is better)')\n",
    "ax.set_yticks(np.arange(means.shape[0]))\n",
    "ax.set_yticklabels([\" w/ \".join(label) for label in means.index.tolist()])\n",
    "plt.tight_layout(pad=1)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (OpenBLAS)",
   "language": "python",
   "name": "py3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}