isotonic regression.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.collections import LineCollection\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.isotonic import IsotonicRegression\n",
    "from sklearn.utils import check_random_state\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "n = 100\n",
    "x = np.arange(n)\n",
    "rs = check_random_state(0)\n",
    "y = rs.randint(-50, 50, size=(n,)) + 50. * np.log(1 + np.arange(n))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100,)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100,)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #############################################################################\n",
    "# Fit IsotonicRegression and LinearRegression models\n",
    "\n",
    "ir = IsotonicRegression()\n",
    "\n",
    "y_ = ir.fit_transform(x, y)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 1)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[:, np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/sklearn/linear_model/base.py:509: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.\n",
      "  linalg.lstsq(X, y)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr = LinearRegression()\n",
    "lr.fit(x[:, np.newaxis], y)  # x needs to be 2d for LinearRegression\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.31347944])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lr.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #############################################################################\n",
    "# Plot result\n",
    "\n",
    "segments = [[[i, y[i]], [i, y_[i]]] for i in range(n)]\n",
    "lc = LineCollection(segments, zorder=0)\n",
    "lc.set_array(np.ones(len(y)))\n",
    "lc.set_linewidths(0.5 * np.ones(n))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[[0, -6.0], [0, -6.0]],\n",
       " [[1, 31.657359027997266], [1, 31.657359027997266]],\n",
       " [[2, 68.93061443340548], [2, 68.93061443340548]],\n",
       " [[3, 86.31471805599453], [3, 77.4581957130341]],\n",
       " [[4, 97.47189562170502], [4, 77.4581957130341]],\n",
       " [[5, 48.587973461402754], [5, 77.4581957130341]],\n",
       " [[6, 130.29550745276566], [6, 100.37627113452281]],\n",
       " [[7, 74.97207708399179], [7, 100.37627113452281]],\n",
       " [[8, 95.86122886681098], [8, 100.37627113452281]],\n",
       " [[9, 152.1292546497023], [9, 142.77028731284636]],\n",
       " [[10, 139.89476363991855], [10, 142.77028731284636]],\n",
       " [[11, 162.24533248940003], [11, 142.77028731284636]],\n",
       " [[12, 166.24746787307683], [12, 142.77028731284636]],\n",
       " [[13, 93.9528664807629], [13, 142.77028731284636]],\n",
       " [[14, 143.4025100551105], [14, 142.77028731284636]],\n",
       " [[15, 153.62943611198907], [15, 142.77028731284636]],\n",
       " [[16, 130.6606672028108], [16, 142.77028731284636]],\n",
       " [[17, 181.51858789480823], [17, 159.46232588434722]],\n",
       " [[18, 143.221948958322], [18, 159.46232588434722]],\n",
       " [[19, 187.78661367769953], [19, 159.46232588434722]],\n",
       " [[20, 183.22612188617114], [20, 159.46232588434722]],\n",
       " [[21, 141.5521226679158], [21, 159.46232588434722]],\n",
       " [[22, 131.7747107964575], [22, 159.46232588434722]],\n",
       " [[23, 185.9026915173973], [23, 159.46232588434722]],\n",
       " [[24, 182.94379124341003], [24, 159.46232588434722]],\n",
       " [[25, 121.9048269010741], [25, 159.46232588434722]],\n",
       " [[26, 134.79184330021644], [26, 159.46232588434722]],\n",
       " [[27, 196.61022550876018], [27, 180.64627212456048]],\n",
       " [[28, 187.3647914993237], [28, 180.64627212456048]],\n",
       " [[29, 199.05986908310777], [29, 180.64627212456048]],\n",
       " [[30, 168.6993602242573], [30, 180.64627212456048]],\n",
       " [[31, 187.28679513998634], [31, 180.64627212456048]],\n",
       " [[32, 206.825378073324], [32, 180.64627212456048]],\n",
       " [[33, 225.31802623080807], [33, 180.64627212456048]],\n",
       " [[34, 215.76740307447068], [34, 180.64627212456048]],\n",
       " [[35, 178.1759469228055], [35, 180.64627212456048]],\n",
       " [[36, 159.54589563221123], [36, 180.64627212456048]],\n",
       " [[37, 150.87930798631928], [37, 180.64627212456048]],\n",
       " [[38, 152.1780823064823], [38, 180.64627212456048]],\n",
       " [[39, 148.4439727056968], [39, 180.64627212456048]],\n",
       " [[40, 174.6786033352154], [40, 180.64627212456048]],\n",
       " [[41, 168.88348091416842], [41, 180.64627212456048]],\n",
       " [[42, 203.0600057846781], [42, 180.64627212456048]],\n",
       " [[43, 148.20948169591304], [43, 180.64627212456048]],\n",
       " [[44, 197.33312448851598], [44, 181.42419146615785]],\n",
       " [[45, 173.43206982445474], [45, 181.42419146615785]],\n",
       " [[46, 173.50738008550292], [46, 181.42419146615785]],\n",
       " [[47, 217.56005054539455], [47, 183.5004138990913]],\n",
       " [[48, 167.59101490553132], [48, 183.5004138990913]],\n",
       " [[49, 180.6011502714073], [49, 183.5004138990913]],\n",
       " [[50, 221.59128163621628], [50, 183.5004138990913]],\n",
       " [[51, 202.56218592907138], [51, 183.5004138990913]],\n",
       " [[52, 176.5145956776061], [52, 183.5004138990913]],\n",
       " [[53, 183.44920232821372], [53, 183.5004138990913]],\n",
       " [[54, 150.36665926162357], [54, 183.5004138990913]],\n",
       " [[55, 151.26758453675748], [55, 183.5004138990913]],\n",
       " [[56, 188.1525633917275], [56, 183.662501486821]],\n",
       " [[57, 206.02215052732095], [57, 183.662501486821]],\n",
       " [[58, 158.87687219528598], [58, 183.662501486821]],\n",
       " [[59, 192.71722811110502], [59, 183.662501486821]],\n",
       " [[60, 172.54369320866556], [60, 183.662501486821]],\n",
       " [[61, 235.35671925225458], [61, 194.14905688424437]],\n",
       " [[62, 161.15673631957662], [62, 194.14905688424437]],\n",
       " [[63, 199.94415416798358], [63, 194.14905688424437]],\n",
       " [[64, 216.71936349478185], [64, 194.14905688424437]],\n",
       " [[65, 190.48273710132125], [65, 194.14905688424437]],\n",
       " [[66, 161.23463096954828], [66, 194.14905688424437]],\n",
       " [[67, 225.97538525880535], [67, 204.21455447880095]],\n",
       " [[68, 202.70532522986298], [68, 204.21455447880095]],\n",
       " [[69, 219.42476210246798], [69, 204.21455447880095]],\n",
       " [[70, 198.13399385206577], [70, 204.21455447880095]],\n",
       " [[71, 174.83330595080275], [71, 204.21455447880095]],\n",
       " [[72, 210.52297205741954], [72, 210.52297205741954]],\n",
       " [[73, 247.2032546602085], [73, 212.9511496110062]],\n",
       " [[74, 256.8744056768155], [74, 212.9511496110062]],\n",
       " [[75, 166.53666701431655], [75, 212.9511496110062]],\n",
       " [[76, 181.1902710926842], [76, 212.9511496110062]],\n",
       " [[77, 266.8354413344796], [77, 220.28290585628667]],\n",
       " [[78, 221.47239262335108], [78, 220.28290585628667]],\n",
       " [[79, 181.10133173369405], [79, 220.28290585628667]],\n",
       " [[80, 211.72245773362195], [80, 220.28290585628667]],\n",
       " [[81, 254.33596236321267], [81, 222.26157242700558]],\n",
       " [[82, 245.9420303898299], [82, 222.26157242700558]],\n",
       " [[83, 239.54083994216566], [83, 222.26157242700558]],\n",
       " [[84, 178.13256282451584], [84, 222.26157242700558]],\n",
       " [[85, 240.71736481267536], [85, 222.26157242700558]],\n",
       " [[86, 220.29540593272918], [86, 222.26157242700558]],\n",
       " [[87, 176.86684072391034], [87, 222.26157242700558]],\n",
       " [[88, 250.43181848660697], [88, 223.73693616781551]],\n",
       " [[89, 226.99048351651325], [89, 223.73693616781551]],\n",
       " [[90, 253.5429753258425], [90, 223.73693616781551]],\n",
       " [[91, 191.08942885245202], [91, 223.73693616781551]],\n",
       " [[92, 196.6299746576628], [92, 223.73693616781551]],\n",
       " [[93, 276.1647391135002], [93, 232.61198322512337]],\n",
       " [[94, 235.69384458002705], [94, 232.61198322512337]],\n",
       " [[95, 201.21740957339182], [95, 232.61198322512337]],\n",
       " [[96, 257.73554892516916], [96, 232.61198322512337]],\n",
       " [[97, 192.24837393352863], [97, 232.61198322512337]],\n",
       " [[98, 264.7559925067295], [98, 246.50725090306705]],\n",
       " [[99, 228.25850929940458], [99, 246.50725090306705]]]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 1)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[:, np.newaxis].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXl4FFX2//86CQESAkHCIlsSRMTBQSOi4gKiOLKL2yhOPgrjEgV0RGdURuY3g45xFNxlxG8cUdF2xwWQcWFxZURWAUVAMQm7IUAgC5Dl/v6oSugk3enudHd6yXk9Tz+pvnWr7q0ueNetc889R4wxKIqiKNFLTKg7oCiKogQXFXpFUZQoR4VeURQlylGhVxRFiXJU6BVFUaIcFXpFUZQoR4VeiTpEJEVEikQkNtR9CTYiMlBENoW6H0p4o0Kv+I2I5IjIxX6e4zMRuSkQ/THG5BljEo0xFYE4XzhjjPnSGNM71P1QwhsVekVxg4g0C+fzKYq3qNArAUVEThSRz0WkUET2isibTvvOFZEV9r4VInKuXZ4FDARm2iaXmfXVt/d9JiL/FJGvReSQiHwiIu3tfWkiYqqEVUTaiciLIrJTRPaLyPtu+j7ePt8TIlIATLPLbxCRjfaxH4tIqtMxl4jIJruPz9rXflNDzicWT4jIryJyUETWi8hv7X0jROQH+1p3iMhf7PLBIrLdqT+/sX+bAyLyvYhc6rTvJRH5t4h8aJ9nuYj0bOi9ViIIY4x+9OPXB8gBLra3XwemYg0iWgLn2+XtgP3AdUAz4Fr7e7K9/zPgJqdzelP/Z+AkIN7+/rC9Lw0wQDP7+4fAm8BxQBxwgZvrGA+UA7fbbcYDY4CfgN/YZX8Dltn12wMHgSvsfXcAZVXX0YDzDQVWAW0Bset0tvftAgba28cB/eztwcB2ezvOPvd9QHPgIuAQ0Nve/xJQAJxlt+0A3gj1vx/9BP+jI3ol0JQBqUAXY8xhY8xXdvlIYIsx5hVjTLkx5nXgR2C0m/N4U/9FY8xmY0wp8BaQXvskItIZGA7caozZb4wpM8Z8Xk//dxpjnrHbLAVuBf5ljNlojCkHHgLS7VH4COB7Y8y79r6ngd1+nK8MaA2cDIhdZ5fT79pHRNrY17HaRd8HAIlYD7yjxpglwAKsh2QV7xljvrXbdrj6zZToQ4VeCTT3YI1Gv7VNBzfY5V2A3Fp1c4Gubs7jTX1nUS3BErnadAf2GWP2e9F3gG21vqcCT9mmkAPAPqzr62r3sbq+McYA22sd7/X5bGGeCfwb+FVEskWkjX3clVgPllzbPHSOi753AbYZYyqdyhrymylRhgq9ElCMMbuNMTcbY7oAtwDPisiJwE4skXMmBdhRdWitfZ7qe8s2oJ2ItPWyfu1+bANuMca0dfrEG2OWYZlTulVVFBFx/t6A82GMedoYcwbQB8ssdbddvsIYMwboCLyP9QZTm51AdxFx/n/dkN9MiTJU6JWAIiK/F5EqsduPJXSVwELgJBH5g4g0E5FrsMRsgV13D3CC06k81fcK2/TxX6wHznEiEicig3w4xXPAX0XkFPv6kkTk9/a+D4G+InKZPfE7CTi+oecTkTNF5GwRiQOKgcNApYg0F5EMEUkyxpRhzQtUujj3cqxR+j32dQ7GMnW94cP1KlGICr0SaM4ElotIETAPuMMYs9UYUwCMAv6MNSF4DzDKGLPXPu4p4CrbE+VpL+r7wnVYNu4fgV+Byd4eaIx5D3gEeENEDgIbsGz+2H35PTDd7mMfYCVwpCHnA9oAz2M9IHPtc85wuoYc+5hbgQwX5z6KJezDgb3As8D1xpgfvb1eJToRy6yoKIq/2CaT7UCGMWZpqPujKFXoiF5R/EBEhopIWxFpgeXWKMA3Ie6WotRAhV5R/OMcLH/+vVhmk8tsN0pFCRvUdKMoihLl6IheURQlygmLIEvt27c3aWlpoe6GoihKRLFq1aq9xpgOnuqFhdCnpaWxcuXKUHdDURQlohCR2qvHXaKmG0VRlChHhV5RFCXKUaFXFEWJclToFUVRohwVekVRlChHhV5RFMVbHA5IS4OYGOuvwxHqHnlFWLhXKoqihD0OB2RmQkmJ9T031/oOkFEnmGhYoSN6RVEUb5g69ZjIV1FSYpWHOSr0iqIo3pCX51t5GKFCryiK4g0pKb6VhxEq9IqiKN6QlQUJCTXLEhKs8jBHhV5RFMUbMjIgOxtSUzEIpKZa38N8IhZU6BVFaSwi1DWxBhkZkJPDK/94A3JyIkLkQYVeUZTGoMo1MTcXjDnmmhiJYu8PIXrYqdArihJ8Itg1MWCE8GGnQq8oSvCJYNfEgBHCh50KvaIowSeCXRMDRggfdir0iqIEnwh2TQwYIXzYqdArihJ8Itg1sUG4mnQN4cNOhV5RFP/xxpskQl0TfcbdpCuE7GGnQq8oin+o62RN6pt0DdHDToVeURT/CJY3SWP7nAeqvTD0MNJ49Iqi+EcwhK2xY78Hsr2UFOt4V+UhwuOIXkS6i8hSEflBRL4XkTvs8mkiskNE1tqfEU7H/FVEfhKRTSIyNJgXoChKiAmGN0kwfc5djdwD2V4Yehh5Y7opB/5sjOkDDAAmiUgfe98Txph0+7MQwN43FjgFGAY8KyKxQei7oijhQDCELVjmD3fzCa5G4A1tLww9jDwKvTFmlzFmtb19CNgIdK3nkDHAG8aYI8aYX4CfgLMC0VlFUcKQYAhbsHzO3Y3cY92MRRvaXph5GPk0GSsiacDpwHK76DYRWScis0XkOLusK7DN6bDt1P9gUBQlEnE2gUydCllZgRM2V28JcXFQVOTfZKm7EXpFRaOYWxzrHaQ9mcY4xtJ+envaT29PzP0xpD2ZhmN98CabvRZ6EUkE5gKTjTEHgVlATyAd2AU85kvDIpIpIitFZGV+fr4vhyqKEmrcmEDS1n8ZmPPXfktITgYRKCjwz4XT3Qi96i0kiOYWx3oHmfMzyS3MBTEUlBZQUFqAwZBbmEvm/Mygib1XQi8icVgi7zDGvAtgjNljjKkwxlQCz3PMPLMD6O50eDe7rAbGmGxjTH9jTP8OHTr4cw2KojQ2bkwg/Ra/Ebg2nM0fiYlw9Gid9nyeLK1vPiHI5papi6dSUlbidn9JWQlTFwcnwJk3XjcCvABsNMY87lTe2ana5cAGe3seMFZEWohID6AX8G3guqwoSshxYwJpVVjQqO35PFkawonSvEKnvhZ2gxW3gGMBrL7BdZ0A4s2I/jzgOuCiWq6U00VkvYisAy4E7gQwxnwPvAX8AHwETDLGVASl94qihAY3JpDipORGba9OeZiGYqishE4HRsOSB2DWGnhiG3z4HOztDeaYDKckBcfX3uOCKWPMV4C42LWwnmOygCYUlk5RmhhZWTUXGAEkJLB6yFgGNVZ7IpatPi3t2KRpYy6y8kBREaz6sQuf3wgffgh79nwAUgHdv4aL74He86H9j9XqmhCXQNaQ4MimroxVFMV3qoRz6lRMbh6SmgJZWeRsiQuO0NdoLxcRsSZl4Zigx8fXH2OmEcjLgwULYP58WLoUjhw5n6QkGDoUThuYy9Rt/WiXDPtL99Muvh2QTEHJPlLbppA1JIuMvsHppwq9oigNIyMDMjJ4ZdpbXD/taqts2ltBb6+4bQcSC/fW3FdSUlfkqwhijJnKSvj222Pivm6dVX7iiTBxIrQoWMoD/7mQuDiY9tmLSMF+fpi4i06JnarPMWfaW1w/+eqg9RE0qJmiRAYhSiodjvg84RvgGDOHDsHKjV3546Cf6RyXzznnwMNZ5Rx3dA+PPgo//ghbtsDjj8NveuQTF2cd996P73Feynk1RL6x0BG9ooQ7jR3gK8wpTkquO6IHy9e+tLTOvEEgFj3l5loj9vnz4bPP4OjR82jLfobxEaOZzzA+ol3eETg+G3rXvSdb929l3Z51PHaJT8uNAoYKvaKEO57imzcxVg8Zy6CPZtcV9KeesrZrzRs05DeqqLBMMvPnW2aZ9eut8s5phbQ8902OpjqISV7Gp7HlvJkA7eyu7Nvyf6Q8OZURvUawcMtCcsnjrukTKS0vBeCx/z1Gp8ROQbPFu0OFXlHCnTCMbx5KcvoOZNAVA9wLeu15Ay85eBA+/RSef/9M/vIs5OdbIXAGDoTHHoPKXvP4x/prqxc97XM6tqDVse3cwlxmrZxlfREoKD1matp5aCeZ8623scYUe7XRK0q4E8Kk0mFLgHzhf/kFnn4aLrkE2reHq66CNZu68LvfwWuvwd69lvfMXXfBzJ//VO/KVm8J5gpYd6jQK0q4E4bxzSOVigr4+muYMgXue3YoJ5wAd9wB27ZZfz//HJ65ex4OB1x7LbRte+zYQK5azSvMq55gv+7+sUGfYFehV5RwJ5DL9pug905hIbz1Flx/PXTqBOefb5li2rQ6zOMZq9jS5QI2bophxttpDNrmIDbGuDxPIFetpjRrVx0UTgh+nl0VeiWwNEEhaRQCYapoQkm89+xrxVNPwcUXWyaZa66xVqcOHw5vvmnZ35/r9yB3vjeIE3d+4VUEzn9e+E/EZZAA30iISyBrEcHLoOUCnYxVAoe6AYY3Uey9U14O//vfsYVLGzeOBKBPH8u+Pno0nHOONbnqWO8g/aWp5PXNpV1P6/h91Z4zJexLmEm76a9b5aX77BWsxyZVWzdvTdHRouryGnXsVa5VXjd5hXl19mUNySLjb9e5vpAgTbCr0CuBIxqFpCqfaF6eNfnZQHe9sCDCvHcc6x1MXTyVXPL4+5P31HFZrCxtw/4NZ3FzP6Fs0xBMSTuIPUpcj2Uw7D3anvo1e5JzmFG6jxeWtYNlllgLgsFYHjFO3jI1tp08ZZy3ASpMBa9c8YpLr5n6VrnW2JcytVETiKvQK4EjwoTEI9H2hpKS0qjiUhtn4b7rgRugtJR9LSppd0QgoRX7KotrjJ6dBbnaZbHgRNg0mYLNoyBvIFTGcTQhH3rNh5PmQ89PKGt5CIADAJb7eg2xNri2wXtLldeMX+6RboLCBWuCXW30SuCINjfA+t5QIhFP3jtBnF+pk13JFFPQshIjUNDSUFBZhOFY1iWwBbkiFnIGwccz4Jkf4Zkt8MnjUNwRznkMbjgP/nI8XD4eTpkLtsgHG789cJwm2JHgx8XXEb0SOBp5lOI3nswy0faG4ibiJBkZQX978ZRdqQalbeGnYbBpNPw0HA4fBzFHocdSOGsmnLQAjsvxu0/+EBAPHDtIW2OgI3olcIQwe4/PeOOBEmlvKP4k3Qjy24vHEfDek2DZXfDiUpieD3Nfh60Xw8nvw9VXwr3JcN0wOHtmyEU+mHHjg4UKvRJYnIUkK8sSinB0tfRG2CJpoZK/rpNBfnupMwKuaAa/DIaPHoOnN8HMTfDJY1DaDs5/BG4cYJlkLrsB+rwLLYoC0o8qqtwkk+OTSY5PRhCX25ia26lJqWSPzm70WDX+oqYbJSikrf8SZswO34lMb4StPlNHuOGvx1OQJ2r/ev5fufXtv1qmmE2jLdPMkbYQewR6LIEBT1kmmbY170vVhGxqUqpLl8Xaro3JCVXlBZa7pHFynRTYFw8pbVO9TvIxZ9pbXH/P1ce2gxw3Plio0CtBod/iN8Lb1dJbYXOVXCMc8XdE3oD5lRpeNNMnAlDAvmPbJftoW3QWRzcOo2TDRbAtH0wstNpDiz4fEJc2j6KTPybZlEBCKwoqikmWVtXeOCnFsWSdkEnGhGe9ugRnUbY6WHcOZs6WuIgVa39QoVeCgtvkEOEykRlpE8ee8HdE7uXbi2PWRKZuzSa3VQUCGKFmhMbKZhR8fxpsHgWbRnNg/4lWeae1MPAhWvT5lOdvvYXrThsHjKtx7i+uvN0OP1xpl1RAwsvQ5ryGDQ5cTXYGMwNWGKM2eiUoFCclu94RLhOZkTRx7A2BmE/wEGbBMWsimTtmkZtYAWKLPEBJO/ju/+DtN2D6XpizGFZMgOTNMHIC3NkdJpwOF/2dI8d/yf+31PUEb71vgYpf6IheCQpuk0OE04g5Uswy3tAI8wlTt2ZTkggYIP83sHm0ZW/ffo5lkkncBX3eht7z4YRF0Ny1O6U7D5ywfwuMYHREr/iPC7e+nL4Do2vEHGiCsTgpQDHaXXH0KOTuGQz/fRKe/gme/QEWPQJlCTAwC24+E+7qCmNuhpPnuRV5cO+DHvZvgRGMjugV/3Cz0CZt2A0w7ZnoGTEHkggJrZCfDwsXWoHCPv4YOLQImpVaXjLnzrC8ZJJ2+HTO+nzQI+ItMELREb3iH27c+votfqNx2o/EsMghDK3gWO8g7ck0xjGW9tPb0356e2LujyHtyTQmLJhIl3svYdxXa4hLW07HTpWMHw9zP97N0d/MgWtHwz3JkDEKzvx/kLQDscPGBMIH3fktsDHCAjQldESv+Icb+6lbe2sgaeyRcaAiWYYotEJVvJmSspJjnjLlzSHnYnI3j+K5zaPhQA8AyjuvhEEPQO/5mOPXcMQpGYcYy0yf6sH9sUE+6I0YFqAp4VHoRaQ7MAfohHV/s40xT4lIO+BNIA3IAa42xuwXEQGeAkYAJcB4Y8zq4HRfCTlu3PqKk5JJDHbbjRgWOaALwEIURbI63kxRB9gywppM/fkSONoampVYE6jn/8syybTZ5fY8RiA1KZWcaTlB7W9EYQ8CrsvNg5fuCbuFdd6YbsqBPxtj+gADgEki0geYAiw2xvQCFtvfAYYDvexPJjAr4L1WGg9PphE3bn2rh4wNft8acWQcUNe/QIZW8MJ0ZQysWwe58zPgP8vg0d3wwUuw/Wzo64BrR1mxZP4wBvo/X6/IVxHI/KkRj1P4icZIC9gQPAq9MWZX1YjcGHMI2Ah0BcYAL9vVXgYus7fHAHOMxTdAWxHpHPCeK8HHm/gpbvzRc/oODH7/6gs6FmDbfUBd/wLlw+/m/qSt/5LDh+Gjj+C226zLP+00YEkWVMbC4GmQ2Q/u6g6jJ0DvDyHusE9NBzJ/asQTAeGsfZqMFZE04HRgOdDJGFP16N+NZdoB6yGwzemw7XaZEq64E0Vv/wEH0a2vXtyNjEeMCHhu1IC7/gXiN6t1f3bTidkl1zBj3jDat7fyo86eDenp8Pzz8MziucTeci4M/id0WUOd9Kde5uOIxOiNfmH//7ju/rGuBw0REM7aa6EXkURgLjDZGHPQeZ8xxuD1P5Pq82WKyEoRWZmfn+/LoUogqW/UHu7/gN2NjBcuDPgIa/WQsWEXydLk5vEdp/IgUzmbb+jMbm5kNmvLT+X6661k2AUF8MEHcNNNcOkZZ1JpKmnTog0CJJdAcrE1uZp6ACasjiG1WTJirPKqfc7bkRq9saGkrf/Ss1kmAsJZe+V1IyJxWCLvMMa8axfvEZHOxphdtmnmV7t8B9Dd6fBudlkNjDHZQDZA//79/cvtpTSc+kbtIU495xWuVrdeF/jEyzl9BzLoigEhi2RZHUCsfA+33fAuRzcO5UibXDho/1frupz4Hn+jec/5bE9dx1utknlrA+xbUTe59QODH+COAXe49yKKibEe+rUwCGJyGuV6wwWvgvNFQNwkjyN624vmBWCjMeZxp13zOBaVaBzwgVP59WIxACh0MvEo4UZ9o/ZIisfuTLBGWI1honJhRpu5ZC43/P1rcp97Eqbv5dCLb3Bk1e+hy0q49Ab48/Fw8wBKL86isMc6iKE6JV/t9HwA9y25D8d6R/X1UFlZ83rc/E5uzVdRjFdzM57SAnoy/TQC3ozozwOuA9aLyFq77D7gYeAtEbkRyAWqnGQXYrlW/oTlXvnHgPZYCSz1jdojKR67MxEwwnKJbUYzJSWsJZ35uaNZcH1vVlT2B66EpFxIf8lKgp32GcQdaVAzHpNbu/n9Vg8Zy6AGtRiBVL3tuLNIuwln7fI89m8pELJV0B6F3hjzFXWnbaoY4qK+ASb52S+lsfAkipEY+CsCH1ClpbD4zsUsKHmMBYxiB90QKjm7cjlcdJ8VKKzjBvf/E32kXvdI5/SCTmadnC1xTUPoa4tzbXwZNDTiWo/60BAIkUAwl/m7mtAcNy58UwB6S6g8geqj1n3cOfNdlq46gUsvheRkGJ0/GwcZnM1yXmQ8uzme/3Euqaf+CzoFTuTBC/dId2adpoArca7CV1fYMHFoUKEPdxqSC9TXB0PtPK8vv+zSN9uvNpo6Dgfm5kxW5SYzzfydM3Ln0vX2K3hpQX/WrYMbboCPOl7PXtozl6sYz8t0xPJG++v3xwW0K03OPdJX3ImwiO8PvTDxyFGhD3d8XYzhb5Job4KU1bNQR6lJSQnMmweZtwrdSjfTn1U8wN9pyWEe4q8sTxzIL7/AzJkw9PGhtEioZU1NSODQpcMA6JzY2W3QsBrJrdU90j8CKc5h4tCgQc3CHXeji9xcayRd2/bsr03QmyBl9T4MnvHcRpjgKuepc7LpGomnnXOhOm3vK91HO0moznPa7ohQUdaLAz9eRPwvV3H0p4FUHG0OzUfR/KSPadVzAcWnLGST7GUTMDUB2s1of+xc9yZAaUz1uUiIoaDgdVrEtmDGJTNqCHSdHKnVF+ZwOe/yxbAbGDQ3cu5PyAjkZL6b+Y7GNoWp0Ic77rxiwPUMvq8PBi/bqxGkLJQRKwOEy0iONi63xc02UFBRAgX9YPNoCjaNht39AChtuxVOn2V5yaR+wdFmZRytOsapLzXOZYqhpb3d0kBlEQBHKo6QOd+61x5H4019MtVfAi3OYRCRU4yLhRGNTf/+/c3KlStD3Y3Q42oBC9QdXdQmNdWyHYIl5u4eDGCNTGpNJlWPbA/kkRrXjhHL97OwZyV5SdCuBBDYFw/tEiw/6n0lBVY5sC+Bmtt2nYKSfSQntKuz7W7E3JDtlKQURvQawcItC8k9kOeyDXdtO4trgziaAFsvtlLpbRkJRZ1BKqD7MkvYT1oAHTYGdAIV7KiRk3MAe0TvgyeUr/X9PX7OtLe4vleZZ+8n+9+9vx5S3vTP398g3BCRVcaY/p7q6Yg+XHAXWz072/pMnYrJzXWtG84jbFevnc7UMuPUHtnmlhcw64xj1QtaOW07j2hbeVHH3Qg4QNu5hbnMWjmr/jbqadtnCrvB5lGWuP9yEVS0hBaFcOJHlrif+BG0Cu5bTSRFjXQO7ezWhzxM/MyjHZ2MDRc82dZzcihOau/6WOdJohrukm5wejBUxyhX6lIpsP0sWPIAzFoDT2yDD2dBwUlw5iy4fgjc3QF+PxZOcwRd5CGyokZ6Fdo5EJEfw2DlabijQh8ueOFv63VgLR8eDJE0QmwUjrSCjZfBB/+Bx3bCf5bDl/dBi4Pwu7th0snwp14w7C44YQk0K2u0rkWaW6RX4QP89TP3NhZ8E38YqNCHC164dDnn1PQmjrk3D4ZIGiEGjQMptFpzN3GvLYLpBfDme8gPV0LaZySOzOC4SR2QP15A8sUvkpyyF5GakR9ruDIeFpJjEqtzpE7oP4HUpFT3rpBeuEt6k281qDRQJL0K7eyvK6M3bwQRkBgk2KiNPlzw1qXLh5AE3kRczBqSdcxG31SojIEdZ8HmUcjmSzF7+lIM9OoF6T2WMGHbDM4/soi478vhe+wJ7FeP/W7uJrydJsVr50h16wrphLs6XudbDQZ+2NBXDxnLoI9m1/9v2l9XRm/eCMIkDEEo0RF9uBCorEOuzltPKICMvhk8M/yYb7XzKNTXkafb7XoW8Lg6Z43ze1j846qvNftnjb7b7U+E76+gxdsvIo/ughf+B1/dS+/ifcy4djU//gibN8OojL1c+ML/EZfa1XUkQmi8Ze3hYG7ww4bu/Abq9rf0FPnRE968EYRJGIJQoiP6cCJEAcRObHciAAuuXcDIk0Z6rO/N6LQGXizgcT6n87aRGOt1uxb1xUavOj4nB+b3+wcL9p/LZwwGWpDAPq7kv4xiAcPMRxxXdACWpUJvp3N58ntujDj94eKN4q9IeuND7o+fuTdvBJGQVyHI6IheYcWOFQCc2fXM4DTgZtTmTV5Z9zHQTZ1RbkUFLFsGby3qS9++0KMH/Gn//eSSyu08w2dcwK90xMH/cS1vcBwHrAN9Hdk1xrL2cMlDGiaxWtzizRtBmIQhCCUq9Aordq4gJSmFjq061l/RH1NCA6MhupxQhupR7sGb/8w7d3zBuHFw/PFw3nnw0f9606EDPP44bO4ymB/5DY9yNxfwBc2oqNuIUzJxr67NX3ODN4SLuSESRNLTv63GuF9hjgp9IInQiI4rdq7gzC4eRvMh8lyoYee12UoPnuJP/I5PaF+ax++fHsT8+TB0KLz+Ojxz9wcsWQJ33gm9pt/s8kFRTa1k4l5fW7DD+IbLSLohIhkOcwu1acphl1GhDxz+Ro0MEQUlBWzdv9Wz0IfQlFAxNoOvXs3hXh7hFDbQk61M5il20JXJPMnnXMCvv8Krr8LYsdCqpZNve22hSk62Ps6iFYRk4n4TTiNpX0RSXRnDEhX6QBEuNlUfWbnTijHk0T7fyKaEwkJ46y34f++dRceOMHAgPM6ddGYXTzCZn+jJD5zCdO5lUGouzepzK3AWqr17rY+zaIWLmcSZSDU3ROj/g2hHvW78pSoQmbtAYmHuwrVipzURe0bnM+qv2AieCz//DPPnw/NzLuDGLCgvh8T4zlx+FYwaBUML3yNp8h8Dnws2XL0ywiDqoc+E40NT0RG9Xziba9wRarHwwIqdK+id3Juklkn1VwyCKaGiUvjyS7jnHpjy72GceKJlVz9Q1JI//xm++gqe+cs85syBq6+GpJuvDs4oN5zMJN7gqw28MW3m4TK3oNRAR/T+UF9uSfBZLDwlwqg3+YWvyTKcwvW2imuFY72j/uX1AYrRfeAAfPS3r5j/4l7+WzKE/f+EuNgKTkotYcr/14ZRo+CrOR9XryP4+dNaPvTBGOWGOjmE/VZ4XW4evHRP/W376l/f2P74gUzaoQQMFXp/cPM6agBJTfVJLLxKhOFr6F8v6xeXFXuX1KKBIrt5MyxYYJllvvyikorK82lPPqOZz2jmc0nzr1ib/nsG/claPPWVzy0EgFCZSXwVYl+X8zf28v9QPzQVl6jQ+4PbbEztSaxKBOIloQ4XXFJWwtTFUwMSNKusDL7++pi4b95slf/2t3B34nOMPvhfVuY2AAAgAElEQVQqZ7OcWCqtHaVEXBrCgOGrEPtqAw+FzTwS5xaiHBV6f3Dzmrp6yFifU7aFQ7hgf/qwbx989JEl7v/9r2Wiad4cBg+G22+3JlPT0oCY28BFSINISkMYUHwVYl8njsN1ollpVHQy1h/cBCLzuLTfxcKqcAgX7EsfjIFNm+DRR+GCC6BjR+vnWLwYLr8c5s61vBg//hhuu80WeXArMO5DHUQ5vk5e+jpxHGkTzUpQUKH3FnerXj1Eh3R5HhcLq7JajCBGQnc7vElqUVYGS5bAXXfBSSfBySfD3Xdbo/cpU+B//4Ndu2D2bLjiCmjd2sVJ3AjP6iFjw3NFZbDxVYh99a+PVH98JaB4VBYRmS0iv4rIBqeyaSKyQ0TW2p8RTvv+KiI/icgmERkarI43KoFc9erGJnvZYwvAQJsWbXwOCew2xK+X9etLalFQYK04veYaaN8ehgyBf/8bTjwRZs60nm3ffQcPPggDBljPwXpxtVI1Pp6B786E666rs6Iybf2Xvv/GoaT2w2rixPofXg0RYl+X8zfx5f8KYIyp9wMMAvoBG5zKpgF/cVG3D/Ad0ALoAfwMxHpq44wzzjBhTWqqMZbE1/ykplZXefkfb7rcroOIy3N9dCKGaZiPf/q4RvV6z+Vr217Ur6w05ocfjJk+3ZiBA42JibG62KmTMTfeaMx77xlz6JDHJrzj1VeNSUhw/dvan0NJ7X0+rTe/QVDw4npMQoJVT1ECALDSeNBXY4znEb0x5gtgn5fPjTHAG8aYI8aYX4CfgLN8eO6EJ4H0XHBje110Whuaxzbn/JTzfT+nnxw9atnWJ0+2siz16WMtYioqsl5Ali+HnTvhP/+Byy6DxMQANexpHQLQqnCv92acUJt+vLgeDQeghAJ/vG5uE5HrgZXAn40x+4GuwDdOdbbbZZFNID0X3HjqLO5/HOccn0ZCXD2RFgPI3r2Wd8zMt8/h9ifg4EFo0cIyzfzlL5aXTLduQe6EFw9Krxf5hEOiDm8f/BoOQGlkGjr7NwvoCaQDu4DHfD2BiGSKyEoRWZmfn9/AbjQSgfRccOGpUzDrMdYeyePiEy4OTH9dYAx8/z088gg8OPtCOnWC66+HzXntufpq+OADyx7/4Ydw662NIPLg24PS00g4HIJpeXs96tqoNDINEnpjzB5jTIUxphJ4nmPmmR1Ad6eq3ewyV+fINsb0N8b079ChQ0O60Xj4m8+1tscO1PDUWdq/PQbDkB5D6hzjjxniyBH49FP405+gZ09rwdKUKXC0rBl/+xusWAFP3jWf55+HSy+FVq18bsI/XD1ARdzXr28kHA7BtFxdT23UtVEJAQ0y3YhIZ2PMLvvr5UCVR8484DUReRzoAvQCvvW7l2GA41SYOhlyD0ByQhHsuIN991/nRVyZAtqVAFfDvgRIKcxlxJvXszDnDnLtekVHiwAYO3csDw15iIx1NNgMcbC4BS+/bK1I/fhjy87esiVcfDHce69lkln8/KfVsWR+WBDoX8oH3C2XdxcNtL6RcDgsDHJ1PSNGWPHuNRyAEkI8Cr2IvA4MBtqLyHbgH8BgEUnHWuKYA9wCYIz5XkTeAn4AyoFJxhgXudsiC7/j0DgN8nLbwqx+lVBeUKdeXmGeFXPm03gyvFkW73Bg7pvK6Xmt+ddTOczv8Ee+2XIpBujSBf7wBxg9Gi66yPNAM2S4Wy7va2CscAmmpcv/lTDEo9AbY651UfxCPfWzgKh6Nw14HJp6rBMlZSVMTS8h4zMXO20zxJEj8Pk/FjP/sYMsKF9KDj3gAPQ/uIrxv1nBba+O4PTT67eChDUNCYylwbQUxS0a68YLGjsOTZ6L0PC/0oEP241jwZXwySdQVDSEeM7hd3zKfTzESD6kS+Uuina2J7FfmE9ue0NDRsY6mlYUl6jQe0FKUgq5hfUkFwl0e3HJmPhS1pf2ZD6jWcAolnM2piCGrsstLRv9/0ZxEYuJ53CNY5tscDBFUdyisW68IGtIVuP4t5e1oMXWy+i1fhlpCb9yGuv4G1lUNo/n/ivXs3o1bNsGzz0HI1M31BF58CE4WKgXFymK0mjoiN4LMvpmUHy0mFsW3AIGkm1z/b4ELI+aqu0ES2RrZHAq2UdywrHsTylJKYzoNYKFWxaSV5hH2/KTObrxdxRvuBD55RKOHE1gWQL87nfw99EwciQcf3x63U75ESI5bf2XMGN2aBcXKYrSeHgTJyHYn7CPdWOM2bBng2Ea5vVZk+rGM6knfomrWDJr1hjzwAPGnHnmsVO0a1NsJkwwZuFCY0pKvOzUq69a8XZEjElONiY52VQiVlk98VQOJbX3GLtHUZTwh0DFumkSuAtB7MSOQ9a6ry7DfU9QXVpquVJPmGA5g5x+OvzjH1Zz/7zqO9Z2Hk7+wUSeXZjG8H0O4uO97HdVVMJXXrEaKSioEfnRnTnGrR1fl+YrSnTizdMg2J+QjuhdRRx0MUJ/cc2LhmmYnwp+8uq0O3cak51tzOm9t1efvlUrYy6/3JjZs43Zvdv7tj3iRXRNZ3REryjRATqi9xIvY6TsOGiP6Ft3cXkaY2D1anjgATjzTGvBUmYm5O1uyx//aKXZ27sX3n0X/vhH6NTJ+7Y94uPy/9VDxmrWIUVpQqjQeymSOw/t5LiWxxEfd8yuUlJi5Ui95Rbo3h3OOAOmTYNmzSzNXLcOHrvjQ2bOhKFDrVAEDWnbIz6mo8vpO1CzDilKE0K9bryMkbLj0A66tO7Cjh1WhMf5860Y7qWlVnz2Sy6xwg2MGGHlT61izVz/2/ZIQ5b/6+IiRWkyqNB7EMnKSssks8IxkkMbBtNtklUlLQ1uuskS90GDrFjugW7ba3T5v6Io9eGNIT/Yn5C7V9puilWuiUX/ed188IExN91kTOfO9lylVJiOJ282//qXMRs2WG6S3uAxrZ2zi6QHt0i/qXWdmtJOUSIbvJyM1RE9QEYG2y/I4G+3riLfnMGS2+DwYWjdGoYNgxEjKrlh4/HcPDSTKRc9GPC2G2XkHQ4ZmBRFCQlNVugrK2HlSmsydf58WLsW4AxOOMGaXB01yjLJNG8OOw/txjyeT9fWEZwVsT4PHxV6RYlqmpTQFxdbGZfmz7cmVPfssRYtnXuulWKvIvcjpswcVie8785DOwH3rpURQThkYFIUJSREvdDn5Vmj9gULYMkSK5Z7UpJlkhk1CoYPh2Q7DticaQddxnCv8qHv2iaCR/ThkIFJUZSQEHVCX1kJ3357zCSzbp1V3rOnFYJg9GgYOBDi4rw/Z1SM6MMlA5OiKI1OVAh9UVFNk8yvv0JsLJx3HsyYYY3ce/dueMalHYd2ECMxdGrVKbAdb0zUBVNRmiwRLfRffw3//CcsXQpHj1ommeHDrVH7sGHQrl1g2tl5aCfHJx5PbExsYE4YKnSRlKI0SSJa6MvL4Zdf4PbbrVH7eef5ZpLxlp2Hdka2x42iKE2aiBb6QYNg06bgt7Pj0A56Htcz+A0piqIEgYgOatZQm7uv6IheUZRIJqKFvjEoLStlX+m+yPa4URSlSaNC74FdRbuACPehVxSlSaNC7wFPCUcURVHCHRV6D1QtllIbvaIokYpHoReR2SLyq4hscCprJyKfisgW++9xdrmIyNMi8pOIrBORfsHsfKBwrHeQ9mQa4xhL++ntaT+9PTH3x9B+entumHcDAMMcw3Csd51sW1EUJZzxZkT/EjCsVtkUYLExphew2P4OMBzoZX8ygVmB6WbwcKx3kDk/k9zCXBBDQWkBBaUFGKztkjIrZMD2g9vJnJ+pYq8oSsThUeiNMV8A+2oVjwFetrdfBi5zKp9jx8T/BmgrIp0D1Vm/cTis1FAxMdZfh4Opi6dWi7knSspKmLrYx8TdiqIoIaahNvpOxphd9vZuoCoITFdgm1O97XZZ6KlKvJGbC8ZUJ97IK3QR0bEe8go1rK+iKJGF35Oxdjor4+txIpIpIitFZGV+fr6/3fCMm8QbKUW+xa9JSdKwvoqiRBYNFfo9VSYZ+++vdvkOoLtTvW52WR2MMdnGmP7GmP4dOnRoYDd8wE2CjayPK2gR611m74QyyGoxwrv2bDPRdfePrTYTKYqihIKGCv08YJy9PQ74wKn8etv7ZgBQ6GTiCS1uEmxkHExlZK+RAAhCcnwyyfHJiIHkYusjBlIPQPY8yHhkoee2nMxEwjEzkYq9oiihQCzLSz0VRF4HBgPtgT3AP4D3gbeAFCAXuNoYs09EBJiJ5aVTAvzRGLPSUyf69+9vVq70WM0/nJJjV5OQANnZXBb7Nj/u/ZEfb/uxepeRGEukayNiZTepj7Q019mcUlMhJ6dB3VcURamNiKwyxvT3VM9j9EpjzLVudg1xUdcAkzx3LwQ4Jd4wuXlI6rHEG2uevI9zu59bo3pxUjKJhXvrnseb1Huan1VRlDCiaa2MzciAnBxe+ccb1sg6I4OCkgLyCvM4/fjTa1RdPWSsNeJ3xtvUe+4eBpqfVVGUEBD9Qu/Cd96ZtbvXAtQR+py+AyE72zK3iFh/s7O9y9CUldXwh4SiKEqAiW6hd+M7n7b+y+oqa3avAeD0zqfXPd5+A6CysvoNwCsyMhr+kFAURQkwEZ1hyiNufOf7LX4DeAawhL57m+60T2gf2LY1P6uiKGFCdI/o3Ux+tiosqN5es2uN69G8oihKlBDdQu9m8rM4KRmwYtdsKthUxz6vKIoSTUS30LuZFF09ZCwA6/aso9JUqtArihLVRLeNPiMDx8Gvmbo1m7xWFbQ7IpAQQ0HFv0me/jqHyw8DMGnhJIrKisjoqzZ1RVGij6gWesd6B5n7X6YksQKAgpYGKotAoKD0mJ1+x6EdZM7PBFCxVxQl6ohq002DYs1rMDJFUaKMqB7R+xo7Pq8wF6Za8XAEjgUjA3WVVBQlYonqEX33pO6eKzmRUhTr0u+eqZpVSlGUyCWqhT7jt96PwhPiEsj6uML1Tg1GpihKBBOVQu+YNZHUu5vxr6/+hVRCckxizVjzTtsYITUplezR2WQcTHV9Qg1GpihKBBN1NnrHrIlk7phFSaL13QiUHinila4TyJjwbJ36c6a9xfWTr7a+ZOE6Zr0GI1MUJYKJuhH91K3ZlMTVLCuJs8o9osHIFEWJQqJuRJ/XyrWd3V15HTQYmaIoUUbUjehTimN9KlcURYl2ok7os07IJK7W4D2hzCpXFEVpikSd0GdMeJYhiaciBsRAalEs2W4mYhVFUZoCUWejB2iV1ove+UfZOGljqLuiKIoScqJuRA9W6IOUJPV9VxRFgWgW+jYq9IqiKBCFQn+4/DB7ivfoiF5RFMUm6oR+W+E2ABV6RVEUm6gT+qrQxCr0iqIoFlHndaNCryiNS1lZGdu3b+fw4cOh7krU0rJlS7p160ZcXJznyi7wS+hFJAc4BFQA5caY/iLSDngTSANygKuNMfv9accX8grzEIRubbo1VpOK0qTZvn07rVu3Ji0tDREJdXeiDmMMBQUFbN++nR49ejToHIEw3VxojEk3xvS3v08BFhtjegGL7e+NRl5hHscnHk+LZi3qr6gpAxUlIBw+fJjk5GQV+SAhIiQnJ/v1xhQMG/0Y4GV7+2XgsiC04Za8g1740DscVjji3FwEcyxloIq9ojQIFfng4u/v66/QG+ATEVklIlXBZDoZY3bZ27uBTq4OFJFMEVkpIivz8/P97MYxvFosNXWqpgxUFKXJ4K/Qn2+M6QcMByaJyCDnncYYg/UwqIMxJtsY098Y079Dhw5+dqP6nN4JvbvUgJoyUFGCj202JSYmYGbT2NhY0tPTOeWUUzjttNN47LHHqKysrPeYnJwcXnvtNb/bjgT8EnpjzA7776/Ae8BZwB4R6Qxg//3V3056S35JPofLD3sWenepATVloKIEFyezKSZwZtP4+HjWrl3L999/z6effsp///tf7r///nqPUaH3AhFpJSKtq7aBS4ANwDxgnF1tHPCBv530Fq9dK7OyrBSBzmjKQEUJPo1gNu3YsSPZ2dnMnDkTYww5OTkMHDiQfv360a9fP5YtWwbAlClT+PLLL0lPT+eJJ55wWy8qMMY06AOcAHxnf74HptrlyVjeNluARUA7T+c644wzTCCY+8NcwzTM6p2rPVd+9VVjUlONEbH+vvpqQPqgKE2NH374wfvKIsZYY/maHxG/+tCqVas6ZUlJSWb37t2muLjYlJaWGmOM2bx5s6nSm6VLl5qRI0dW13dXL1xw9TsDK40Xet1gP3pjzFbgNBflBcCQhp7XH3xaLKUpAxWl8UlJscw1rsqDRFlZGbfddhtr164lNjaWzZs3+1UvEomqlbF5hXkkxCXQLr5dqLuiKIorsrIsm7yz+SYIZtOtW7cSGxtLx44duf/+++nUqRPfffcdlZWVtGzZ0uUxTzzxhFf1IpGoinVT5XGjPr2KEqZkZEB2NqSmgoj1Nzs7oG/X+fn53Hrrrdx2222ICIWFhXTu3JmYmBheeeUVKiqsXKOtW7fm0KFD1ce5qxcNRNWIPrcwV2PcKEq4EwSzaWlpKenp6ZSVldGsWTOuu+467rrrLgAmTpzIlVdeyZw5cxg2bBitWrUC4NRTTyU2NpbTTjuN8ePHu60XDUSV0OcV5pHeKT3U3VAUpZGpb/Tdq1cv1q1bV/39kUceASAuLo4lS5bUqOuqXjQQFaYbx6yJpNzdjF+Lf+Xtb17AMWtiqLukKIoSNkS80DtmTSRzxyy2JVpP9MIWhswds1TsFUVRbCJe6KduzaakVojmkjirXFEURYkCoc9r5do2565cURSlqRHxQp9SHOtTuaIoSlMj4oU+64RMEspqliWUWeWKoihKFAh9xoRnye46gZblgIHUoliyu04gY8Kzoe6aoigucKx3kPZkGjH3x5D2ZBqO9f6HKU5MTGzQcQ899JBf7f79739n0aJFXtXNyckhPj6e9PT06s/Ro0eZN28eDz/8MADvv/8+P/zwg199coVYcXFCS//+/c3KlSv9OsfJM0+mb6e+vP37twPUK0VRvGHjxo385je/8aquY72DzPmZlJQdC4GQEJdA9uhsMvo2fBFVYmIiRUVFjXZcQ8jJyWHUqFFs2LDBbZ3x48czatQorrrqqjr7XP3OIrLKHEvj6paoWDBl7IQjI3uNDHVXFKVJM/mjyazdvdbt/m+2f8ORiiM1ykrKSrjxgxt5ftXzLo9JPz6dJ4c96VX7u3bt4pprruHgwYOUl5cza9YsBg4cyOuvv85DDz2EMYaRI0fyyCOPMGXKlOoVtaeccgoOh4PHH3+c2bNnA3DTTTcxefJkcnJyGD58OOeffz7Lli2ja9eufPDBB8THx9cQ5hUrVnDHHXdQXFxMixYtWLx4Ma1bt/bY55deeomVK1fyhz/8gXnz5vH555/z4IMPMnfuXHr27OnVdXsi4k03AAWlBZSWl2r4A0UJc2qLvKdyX3nttdcYOnQoa9eu5bvvviM9PZ2dO3dy7733smTJEtauXcuKFSt4//33efjhh6sTljgcDlatWsWLL77I8uXL+eabb3j++edZs2YNAFu2bGHSpEl8//33tG3blrlz59Zo9+jRo1xzzTU89dRTfPfddyxatIj4+Pg6/fv555+rzTaTJk2qse/cc8/l0ksvZcaMGaxduzZgIg9RMqKvCk/cPal7iHuiKE0bTyPvtCfTyC2sG6Y4NSmVz8Z/5nf7Z555JjfccANlZWVcdtllpKens2TJEgYPHkxVytKMjAy++OILLrvsshrHfvXVV1x++eXVMW6uuOIKvvzySy699FJ69OhBeroVXuWMM84gJyenxrGbNm2ic+fOnHnmmQC0adPGZf969uzJ2rXu33iCRVSM6H2KQ68oSsjIGpJFQlzN7G4JcQlkDQlMmOJBgwbxxRdf0LVrV8aPH8+cOXMCct4WLVpUb8fGxlJeXh6Q8zYWUSH02wq3ASr0ihLuZPTNIHt0NqlJqQhCalKq3xOxzuTm5tKpUyduvvlmbrrpJlavXs1ZZ53F559/zt69e6moqOD111/nggsuAKzAZmVlln/2wIEDef/99ykpKaG4uJj33nuPgQMHetVu79692bVrFytWrADg0KFDDXoY1A6dHCiixnTTIrYFHRI6hLoriqJ4IKNvRsCEvTafffYZM2bMIC4ujsTERObMmUPnzp15+OGHufDCC6snY8eMGQNAZmYmp556Kv369cPhcDB+/HjOOusswJqMPf300+uYaVzRvHlz3nzzTW6//XZKS0uJj49n0aJFPrt9jh07lptvvpmnn36ad955J2B2+qhwr7zmnWtYs2sNm2+PntRfihIp+OJeqTQcf9wro8J0k1eYpxOxiqIobogaoVf7vKIoimsiXujLKsrYdWgXKW1U6BVFUVwR8UK/49AODEZH9IqiKG6IeKFXH3pFUZT6UaFXFEWJcqJG6NXrRlGaLq781Z977rmArYz1lsGDB9O7d+/qeDbvvPMOYMWxASuC5WuvvdaofYIgLpgSkWHAU0As8B9jzMPBaCevMI/k+OQ6y6oVRWna3HrrrUE9vzEGYwwxMTXHyw6Hg/79a7q2L1u2DDgm9H/4wx+C2rfaBEXoRSQW+DfwO2A7sEJE5hljAh5Rf9vBbWq2UZQwYfJkCHTMrvR0eNK7KMU1mDZtGomJifzlL39h8ODBnH322SxdupQDBw7wwgsvMHDgQCoqKpgyZQqfffYZR44cYdKkSdxyyy0UFRUxZswY9u/fT1lZGQ8++CBjxowhJyeHoUOHcvbZZ7Nq1SoWLlxIamqqx75Uxb2fMmUKGzduJD09nXHjxnHnnXc24BfxnWCN6M8CfjLGbAUQkTeAMUDAhT6vMI+exwUunKeiKNFJeXk53377LQsXLuT+++9n0aJFvPDCCyQlJbFixQqOHDnCeeedxyWXXEL37t157733aNOmDXv37mXAgAFceumlgBWy+OWXX2bAgAEu28nIyKgOUbx48WKSk5Or9z388MM8+uijLFiwIPgX7ESwhL4rsM3p+3bgbOcKIpIJZAKkpDR8RJ5XmMeFaRc2+HhFUQJHQ0bejcUVV1wB1Awz/Mknn7Bu3bpqW3phYSFbtmyhW7du3HfffXzxxRfExMSwY8cO9uzZA0BqaqpbkQfXpptQE7KgZsaYbCAbrFg3vh7vmDWRv27N5mBiBa98OZOzc8o1T6yiKG6pCjXsHGbYGMMzzzzD0KFDa9R96aWXyM/PZ9WqVcTFxZGWlsbhw4cBquPVRxLB8rrZATi7wXSzywKCY9ZEMnfMYltiBQAHWhgyd8zCMWtioJpQFKUJMHToUGbNmlUdqnjz5s0UFxdTWFhIx44diYuLY+nSpeTm1k2W0hCCFYbYE8ES+hVALxHpISLNgbHAvECdfOrWbEriapaVxFnliqI0PUpKSujWrVv15/HHH/fquJtuuok+ffrQr18/fvvb33LLLbdQXl5ORkYGK1eupG/fvsyZM4eTTz45IP089dRTiY2N5bTTTuOJJ54IyDm9IWhhikVkBPAklnvlbGOM2xQyvoYpjpkmGHHRpoHKaaEPu6woTQkNU9w4+BOmOGg2emPMQmBhMM6dUhxLrm22qV2uKIqi1CQiV8ZmnZBJQlnNsoQyq1xRFEWpSUQKfcaEZ8nuOoHUoljEQGpRLNldJ6jXjaKEiHDIVBfN+Pv7RmzO2IwJz5KBCruihJqWLVtSUFBAcnIyIi4mzxS/MMZQUFBAy5YtG3yOiBV6RVHCg27durF9+3by8/ND3ZWopWXLlnTr1q3Bx6vQK4riF3FxcfTo0SPU3VDqISJt9IqiKIr3qNAriqJEOSr0iqIoUU7QVsb61AmRfKChwSTaA3sD2J1IoSled1O8Zmia190Urxl8v+5UY0wHT5XCQuj9QURWerMEONpoitfdFK8ZmuZ1N8VrhuBdt5puFEVRohwVekVRlCgnGoS+qcYmborX3RSvGZrmdTfFa4YgXXfE2+gVRVGU+omGEb2iKIpSDyr0iqIoUU5EC72IDBORTSLyk4hMCXV/goGIdBeRpSLyg4h8LyJ32OXtRORTEdli/z0u1H0NBiISKyJrRGSB/b2HiCy37/mbdqrKqEFE2orIOyLyo4hsFJFzmsK9FpE77X/fG0TkdRFpGY33WkRmi8ivIrLBqczl/RWLp+3rXyci/RrabsQKvYjEAv8GhgN9gGtFpE9oexUUyoE/G2P6AAOASfZ1TgEWG2N6AYvt79HIHcBGp++PAE8YY04E9gM3hqRXweMp4CNjzMnAaVjXHtX3WkS6An8C+htjfouVfnQs0XmvXwKG1Spzd3+HA73sTyYwq6GNRqzQA2cBPxljthpjjgJvAGNC3KeAY4zZZYxZbW8fwvqP3xXrWl+2q70MXBaaHgYPEekGjAT+Y38X4CLgHbtKVF23iCQBg4AXAIwxR40xB2gC9xorkm68iDQDEoBdROG9NsZ8AeyrVezu/o4B5hiLb4C2ItK5Ie1GstB3BbY5fd9ul0UtIpIGnA4sBzoZY3bZu3YDnULUrWDyJHAPUGl/TwYOGGPK7e/Rds97APnAi7a56j8i0ooov9fGmB3Ao0AelsAXAquI7nvtjLv7GzCNi2Shb1KISCIwF5hsjDnovM9YPrJR5ScrIqOAX40xq0Ldl0akGdAPmGWMOR0oppaZJkrv9XFYo9ceQBegFXXNG02CYN3fSBb6HUB3p+/d7LKoQ0TisETeYYx51y7eU/UaZ//9NVT9CxLnAZeKSA6WWe4iLPt1W/v1HqLvnm8Hthtjltvf38ES/mi/1xcDvxhj8o0xZcC7WPc/mu+1M+7ub8A0LpKFfgXQy56Zb441eTMvxH0KOLZd+gVgozHmcadd84Bx9vY44IPG7lswMcb81RjTzRiThnVvlxhjMoClwFV2tai6bmPMbmCbiPS2i4YAPxDl9xrLZDNARBLsf+9V1x2197oW7u7vPOB62/tmAFDoZBxM3UkAAACzSURBVOLxDWNMxH6AEcBm4Gdgaqj7E6RrPB/rVW4dsNb+jMCyVy8GtgCLgHah7msQf4PBwAJ7+wTgW+An4G2gRaj7F+BrTQdW2vf7feC4pnCvgfuBH4ENwCtAi2i818DrWPMQZVhvcDe6u7+AYHkW/gysx/JKalC7GgJBURQlyolk042iKIriBSr0iqIoUY4KvaIoSpSjQq8oihLlqNAriqJEOSr0iqIoUY4KvaIoSpTz/wN5NyWLgYxQ+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(x, y, 'r.', markersize=12)\n",
    "plt.plot(x, y_, 'g.-', markersize=12)\n",
    "plt.plot(x, lr.predict(x[:, np.newaxis]), 'b-')\n",
    "plt.gca().add_collection(lc)\n",
    "plt.legend(('Data', 'Isotonic Fit', 'Linear Fit'), loc='lower right')\n",
    "plt.title('Isotonic regression')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}