Numerical Solution For Finding Langrange Multipliers-Corrected Gradient.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Solution For Finding Langrange Multipliers\n",
    "\n",
    "In this code we're looking at how to find the values of $\\lambda_i$.\n",
    "\n",
    "We start by defining a set of random c vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(suppress=True, precision=5)\n",
    "\n",
    "\n",
    "cs = []\n",
    "\n",
    "for i in range(10):\n",
    "    cs.append(np.random.randn(7,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note 'ls' is the vector of $\\lambda$s:\n",
    "\n",
    "Here are functions to calculate M, L, dL/dl, one to check the gradient (which we call at the end with a random set of ls).\n",
    "\n",
    "Note: The functions are written to minimise $\\log |M|$ (useful for using the scipy minimize function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value:\n",
      "-6.65549620755\n",
      "Approx\n",
      "[-0.71781 -0.35232 -6.99802  0.58674 -0.39791 -0.53393 -2.68435 -1.35602\n",
      " -0.68122 -0.54406]\n",
      "Analytical\n",
      "[-0.71781 -0.35232 -6.99802  0.58674 -0.39791 -0.53393 -2.68435 -1.35602\n",
      " -0.68122 -0.54406]\n",
      "Difference\n",
      "[-0. -0. -0. -0. -0. -0. -0. -0. -0. -0.]\n",
      "Ratio\n",
      "[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(suppress=True, precision=5)\n",
    "def calcM(ls,cs):\n",
    "    \"\"\"\n",
    "    Find the covariance matrix, M, as the lambda weighted sum of c c^T\n",
    "    \"\"\"\n",
    "    d = len(cs[0])\n",
    "    M = np.zeros([d,d])\n",
    "    ccTs = []\n",
    "    for l,c in zip(ls,cs):        \n",
    "        ccT = np.dot(c,c.T)\n",
    "        #print c,ccT,l,M\n",
    "        M = M + l*ccT       \n",
    "        ccTs.append(ccT)\n",
    "    return M\n",
    "\n",
    "def L(ls,cs):\n",
    "    \"\"\"\n",
    "    Find L = -log |M| + sum(lambda_i * (1-c^T M^-1 c))\n",
    "    \"\"\"\n",
    "    M = calcM(ls,cs)\n",
    "    Minv = np.linalg.pinv(M)\n",
    "    t = 0\n",
    "    for l,c in zip(ls,cs):        \n",
    "        t += l*(1-np.dot(np.dot(c.T,Minv),c))[0,0]\n",
    "\n",
    "    return (-np.log(np.linalg.det(M)) + t)\n",
    "    #return t\n",
    "\n",
    "#ignore this function (it was for testing)\n",
    "def shortened_L(ls,cs):\n",
    "    M = calcM(ls,cs)\n",
    "    return -(-np.log(np.linalg.det(M)))  \n",
    "\n",
    "def dL_dl(ls,cs):\n",
    "    \"\"\"\n",
    "    Find the gradient dL/dl_j\n",
    "    \"\"\"\n",
    "    M = calcM(ls,cs)\n",
    "    Minv = np.linalg.pinv(M)            \n",
    "    grads = np.zeros(len(ls))    \n",
    "    for j in range(len(cs)):        \n",
    "        grads[j] = -np.trace(np.dot(Minv,np.dot(cs[j],cs[j].T)))     \n",
    "    return np.array(grads)+1\n",
    "    \n",
    "def checkgrad(ls,cs):\n",
    "    approx_dL_dl = []\n",
    "    d = 0.0001\n",
    "    for i in range(len(ls)):\n",
    "        delta = np.zeros_like(ls)\n",
    "        delta[i]+=d\n",
    "        approx_dL_dl.append(((L(ls+delta,cs)-L(ls-delta,cs))/(2*d)))\n",
    "    approx_dL_dl = np.array(approx_dL_dl)\n",
    "\n",
    "    print \"Value:\"\n",
    "    print L(ls,cs)\n",
    "    print \"Approx\"\n",
    "    print approx_dL_dl\n",
    "    print \"Analytical\"\n",
    "    print dL_dl(ls,cs)\n",
    "\n",
    "    print \"Difference\"\n",
    "    print approx_dL_dl-dL_dl(ls,cs)\n",
    "    print \"Ratio\"\n",
    "    print approx_dL_dl/dL_dl(ls,cs)\n",
    "\n",
    "ls=np.random.rand(len(cs))\n",
    "\n",
    "checkgrad(ls,cs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The checkgrad function suggests the gradients are all computed correctly!\n",
    "\n",
    "In this part I try two ways of optimising,\n",
    "\n",
    "   1. using the scipy minimize function\n",
    "   \n",
    "   2. using my own optimiser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.99802014997\n",
      "Iteration limit exceeded    (Exit mode 9)\n",
      "            Current function value: -9.63990359477\n",
      "            Iterations: 101\n",
      "            Function evaluations: 1402\n",
      "            Gradient evaluations: 101\n",
      "Lambdas:\n",
      "[ 0.9614   0.13897  0.97277  0.      -0.       0.97788  1.00864  0.99754\n",
      "  1.00427  0.96032]\n",
      "7.02177262942\n",
      "Solutions to c^T M^-1 c (should all be 0<=x<=1)\n",
      "[[ 0.99912]]\n",
      "[[ 1.00029]]\n",
      "[[ 0.99096]]\n",
      "[[ 0.29621]]\n",
      "[[ 0.76422]]\n",
      "[[ 1.00011]]\n",
      "[[ 0.99121]]\n",
      "[[ 1.00072]]\n",
      "[[ 0.99533]]\n",
      "[[ 1.00058]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lionfish/anaconda2/lib/python2.7/site-packages/ipykernel/__main__.py:29: OptimizeWarning: Unknown solver options: xtol\n",
      "/home/lionfish/anaconda2/lib/python2.7/site-packages/ipykernel/__main__.py:26: RuntimeWarning: invalid value encountered in log\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def maxcMinvc(ls,cs): #this allows us to enforce a bound on the max(c^T M^-1 c) <= 1 (over all c_i)\n",
    "    M = calcM(ls,cs)\n",
    "    Minv = np.linalg.pinv(M)\n",
    "    maxv = -np.Inf\n",
    "    for l,c in zip(ls,cs):\n",
    "        v = np.dot(np.dot(c.transpose(), Minv),c)\n",
    "        if v>maxv:\n",
    "            maxv = v\n",
    "    if np.isscalar(maxv):\n",
    "        return maxv\n",
    "    else:\n",
    "        return maxv[0,0]\n",
    "\n",
    "print maxcMinvc(ls,cs)\n",
    "\n",
    "def findLambdas(cs):\n",
    "    ls = np.random.rand(len(cs))+0.5\n",
    "    rank = np.min([len(cs),len(cs[0])])\n",
    "    \n",
    "    #other contstraints\n",
    "    #cons = ({'type': 'eq', 'fun' : lambda ls: np.sum(ls)-rank},{'type':'ineq','fun':lambda ls:np.min(ls)})\n",
    "    #cons = ({'type': 'eq', 'fun' : lambda ls: np.sum(ls)-rank},{'type':'ineq','fun':lambda ls:np.min(ls)},{'type':'ineq','fun':lambda ls:1-maxcMinvc(ls,cs)})\n",
    "    cons = ({'type':'ineq','fun':lambda ls:np.min(ls)},{'type':'ineq','fun':lambda ls:1-maxcMinvc(ls,cs)})\n",
    "    #cons = ()\n",
    "    #res = minimize(L, ls, args=(cs), method='BFGS', options={'xtol': 1e-8, 'disp': True})\n",
    "    #res = minimize(L, ls, args=(cs), method='nelder-mead', options={'xtol': 1e-8, 'disp': True})\n",
    "    res = minimize(L, ls, args=(cs), method='SLSQP', options={'xtol': 1e-8, 'disp': True}, constraints=cons)\n",
    "    #res = minimize(L, ls, args=(cs), method='TNC', options={'xtol': 1e-8, 'disp': True}, constraints=cons)\n",
    "    #res = minimize(L, ls, args=(cs), method='TNC', options={'xtol': 1e-8, 'disp': True}, constraints=cons, jac=dL_dl)\n",
    "    ls = res.x \n",
    "    return ls\n",
    "\n",
    "ls = findLambdas(cs)\n",
    "print \"Lambdas:\"\n",
    "print ls\n",
    "print np.sum(ls)\n",
    "\n",
    "M = calcM(ls,cs)\n",
    "Minv = np.linalg.pinv(M)\n",
    "print \"Solutions to c^T M^-1 c (should all be 0<=x<=1)\"\n",
    "for l,c in zip(ls,cs):\n",
    "    print np.dot(np.dot(c.transpose(), Minv),c)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's my own gradient descent method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambdas:\n",
      "[ 0.95975  0.14205  0.96294  0.       0.       0.97757  0.99976  0.99822\n",
      "  0.99957  0.96014]\n",
      "Solutions to c^T M^-1 c (should all be 0<=x<=1)\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 0.29717]]\n",
      "[[ 0.76416]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n"
     ]
    }
   ],
   "source": [
    "def findLambdas(cs):\n",
    "    \"\"\"\n",
    "    Gradient descent to find the lambda_is\n",
    "    \n",
    "    Parameters:\n",
    "        cs = list of column vectors (these are the gradients of df*/df_i)\n",
    "    \n",
    "    Returns:\n",
    "        ls, storedls = vector of lambdas, and history of lambdas (for plotting)\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    ls = np.ones(len(cs))*0.7\n",
    "    \n",
    "    #learning rate\n",
    "    lr = 0.01\n",
    "    Nits = 10000\n",
    "    storedls = np.zeros([Nits,len(ls)])\n",
    "    for it in range(Nits): \n",
    "        delta_ls = -dL_dl(ls,cs)*lr\n",
    "        ls =  ls + delta_ls\n",
    "        ls[ls<0] = 0\n",
    "        storedls[it,:] = ls\n",
    "    return ls, storedls\n",
    "\n",
    "ls, storedls = findLambdas(cs)       \n",
    "print \"Lambdas:\"\n",
    "print ls\n",
    "M = calcM(ls,cs)\n",
    "Minv = np.linalg.pinv(M)\n",
    "print \"Solutions to c^T M^-1 c (should all be 0<=x<=1)\"\n",
    "for l,c in zip(ls,cs):\n",
    "    print np.dot(np.dot(c.transpose(), Minv),c)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pleasingly it matches the optimiser above (which uses approximations to the gradient)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAD7CAYAAABZqT4/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm0HOV5p5+fJIQEQgK0oQUJGxBE2BZhkSEe4ws2YTGD\nmDnM2DjHHkMyJh5sMyGxcWJ8jBNmxskcx7tjkyE4ZILxDJ6DAbMvwgcwGBBildg3CekKCQkEkkDL\nO390t9RqVd/uvl3d/X7V73NOHfXydVXpd7uefuurr6pkZgRBEAR+GNHrFQiCIAh2JsQcBEHgjBBz\nEASBM0LMQRAEzggxB0EQOCPEHARB4IxR3VyYpBibFwRB05iZ2vn8vpKtbb75S2Z2QDvLywt1cxyz\nJPucfTfXeU5lVa7zy5OJrB7y/RsvfpBTLj6qS2uzA2+ZTWQ1V1z8Ep+5eHavVwXwkc8/XLyaz188\niYms6fWqbGfqm91dl90mtC9mSfbtJtv+Oe0vLy+S78oYZEqvV6Eua5jU61XIxFtm3nLylM8aJvZ6\nFbYzON7PuhSd5MUMvjakWrxJp4K3zDayR69XYSc85RNy7j8KIWbwtSHVUk/OBw1M7/Ka7IynzOYO\nTHL3I9bLfI4a2PmHKuTcXxRGzOBLNLVkSefgHosZ/GQ2d2Ay4G8Po1f5HD2w6x5EyLl/KJSYwY9o\nsvAmnQreMvOWk6d8Qs6tM6PJyROFEzP42pBq8SadCt4y85aTp3xCzsWnkGL2jjfpVPAkH/CXk6d8\nQs7FprBi9rQRZeFNOl7xlpOn71XIubgUVszgayPKwpt0wGdm3nLylFHIuZgUWszgayPKwpt0wGdm\n3nLylFHIuXgUXszgayNKBY+ZhZzrE3KuzwFNTllI2l3S/ZIelvSYpG9ktPmUpEfK092S3t/uOveF\nmMHXRlSLN+FU8JiZt6w8ZRRyzh8zewc43sx+HzgcOEXS/JpmzwPHmdk84BLgH9tdbt+IGXxtRLV4\nE04Fj5l5y8pTRmuY6EbQBZLzhvLD3SldkdNq3r/PzN4oP72PHIZF95WYwddGVIs34VTwmJm3rLxl\nFHLOD0kjJD0MrARuNbMHhmj+J8CN7S6z78QM/jaiarwJp4LHzLxl5S2jkHM+mNm2clfGTOCDkuZm\ntZN0PHA2cGG7y+xLMYO/jagab8Kp4DEzb1l5yyjkXLJp1vQy8IuqqRFm9iZwJ3By7XuSPgBcCpxu\n1sq1+bPpWzGDv42oGm/CqeAxM29Zecso5JzNH1C6OH5lykLSJEkTyo/HAicCS2vazAJ+CXzazJ7L\nY936Wsze8SacCt7EA/6y8pZRyHnYTAPulLQYuB+42cxukHSupM+V23wd2Bf4cXlY3e/aXWjyt5bK\nAw+3EhqKRreo6hUec/OWlbeMvNyqqpnbVOV1a6llTbadSdxayhXeqptavFWDFTzm5i0rbxlF5ZwG\nIeYy3jagWrwJp4LH3Lxl5S2jkLN/oiujBm+7nrV421Wv4DE3b1l5y8h7t0ZeXRnW5O+iVkVXhlu8\nVTe1eKsGK3jMzVtW3jKKytkvIeYMvG1AtXgTTgWPuXnLyltGnuQcgt5B18W8iqndXuSw8LYB1eJN\nOBU85uYtK28ZeZEzRPVcISrmIfC2AaWCx9xCzkMTcvZFT8ScStXsHW+yqcabeMBfXt4yCjn7oWcV\ncypy9rbx1OJNNtV4zM5bXt4y8iTn3Jje5OSInnZlhJzzwZtsqvGYnbe8vGVUSDknRs/7mEPO+eBN\nNtV4zM5bXt4yCjn3lp6LGULOeeFNNtV4zM5bXt4yCjn3DhdiTglvG08t3mRTjcfsvOXlLaOQc29w\nI+ZUqmbwt/HU4k023vGWl7fvV/Jyntrk5Ag3YoaQc554k00Fr7l5y8tbTsnLeZhIminpDklPSHpM\n0peGaHu0pM2S/n27y3UlZkhLzt7xJpsK3qRTwVte3nLqUzlvAS4ws8OAY4HzJB1a20jSCOBbwM15\nLNSdmCEdOXvbcLLwJpsKXrPzlpe3nPpNzma20swWlx+/BSwBZmQ0/SJwNeRzCUGXYk4JbxtOFt5k\nU8Frdt7y8pZTv8m5gqQDgMMp3WKq+vXpwBlm9g9ALpcNdSvmVKpm8LfhpITX7ELOQ9NvcpY0jlJF\nfH65cq7mu8CF1c3bXd6odmfQSVYxlSkM9no1mmKQKe4uhF7NGia5u3B8Ba/ZecvMW05rmOjmYvtD\nMi375YUrYOHKxh+XNIqSlP/FzH6V0eQo4CpJAiYBp0jabGbXDneVu34HkzPs5y1/LhU5g7+7VNTi\nSTS1eM3OW2becuqUnGdqTT53MPlsk21/ln0HE0lXAKvN7IImlnc5cJ2Z/b/W1nRn3HZlVBPdGvnh\nbRe9Gq/ZecvMW05F7taQ9CHgj4ATJD0saZGkkyWdK+lzGR/JpdJ13ZWRKt52OWvxtotejdfsvGXm\nLadkujVaxMzuAUa20P6cPJabRMUMaVXNKeCtCqzGW0VYwVtm3nIqcuXcbZIRM6QlZ28bTRbeRFON\n1/y8ZeYtp5BzPiQlZgg554030VTjNT9vmXnLyZ2c41oZ3SHknC/eRFNNCvl5wFtO7uScGEmKOTW8\nbTRZhJxbw2Ne3nIKOQ+fhmKWdJmkQUmP1nn/I5LWlYeRLJJ0Uf6ruSspVc3gb6PJwqNsKnjMz2Ne\n3nIKOQ+PZirmy4GTGrT5jZkdUZ4uyWG9miLknD8eZVPBY34e8/KWU8i5dRqK2czuBtY2aJbLhTuG\nQ2pyDtrDm3Qg5NwMIefWyKuP+VhJiyX9WtLcnObZNCnJ2dsGk4VH0VTjMUOPmXnLqWdynt7k5Ig8\nxPwQMMvMDgd+CFyTwzwLjbcNJguPovGOx8y8fdeicm6Otk/Jrr4EnpndKOnHkvY1s9ez2i+5+Ort\njycNzGXyQD4FdkpXogN/p9Rm4e005Gq85ucxM29Z1Tt9+96Fm/ntws09WCN/NHV1ufIFoq8zs/dn\nvDfVzAbLj+cD/8fMDqgzn2FdXa4VUpIz+LtSWBbeRFON1/w8ZuYtq0bX1sjt6nLfa7Lt+dlXl+sF\nzQyXuxK4F5gj6WVJZ9dcWelMSY9LepjSBaM/0cH1bUhK/c3gb1czC4+76BW85ucxM29ZRbdGfZK4\nHnOrpFY1g79qJguPVWAFr/l5zMxbVvUq59wq5iaVo7MSqphTJLWqORU8VoEVvFWDFTxm5i0rz5Vz\noxPsym0GytdqflzSnXkst5BihvTk7G1jqYdH0VTwmqHHzLxl5VjOQ55gJ2kC8CPgNDN7H/Af8lho\nYcUMIedO4VE0Fbxm6DEzb1l5lHMTJ9h9CvilmS0vt8+l76rQYk4RbxtLinjNMOTcGI9ybsAcYF9J\nd0p6QNKn85hp4cWcWtUM/jaWLDxKphqvGXrMzVtWicl5FHAEcApwMvB1SQflMdPCk9rJJ+DvpIAs\nPJ5MUY3XDD3m5i2rkpxzuodgndOtFz4MCxe3PfdllO6gvQnYJOk3wDzg2XZmWsjhcvVITc6eNpSh\n8CaZWrzm6DE3T1nN09P5DJe7q8m2H8keLtfgBLtDgR9QqpZ3B+4HPmFmTw5/rfugKyNlvO1i1sPj\n7nk1XnP0mJvXrHpFoxPszGwpcDPwKHAfcGm7UoY+q5ghvaoZfFUxQ+GxAqzGa44ec/OQlZeKuRf0\nXcUcBwM7h8cKsBqvOXrMzWtW/ULfiRlCzp3Eo2RSwGNuqXznikhfjMooCt6OnNfD46iDCp4z9Jib\n57yaxtlF8JuhLytmSLNqDvLBcyUYlXMAfSxmSFPOqWwkHgVTjeccPWbnOa8i0tdihpBzJ/EomGo8\n5+gxO895FY2+F3OqpLKReBRMNank6IXIqzuEmEmzaoZ0NpKQ8/DwmpvXvOqxZUpzkydCzGVCzp3F\nq2QqeM3Ra25e8yoKIeaga3iVTAWvsvGam9e8ikCIuYqomjuPV8lU8Jql19y85pU6IeYaQs6B1yxD\nzv1DiDmDkHNn8SqYarxm6TU7r3mlSpySXTBSOYXW4+nHtXjN0mt2XvMaHN/sHVFyujB/DkTFXIdU\nq+aU8Fr9pYDX7IpYOUs6WdJSSU9LujDj/fGSrpW0WNJjkj7b7jJDzEOQqpxT2ji8CqaC5yy9Zuc5\ns1aRNAL4IXAScBhwVvmuJdWcBzxhZocDxwPfltRWb0SIuaCktHF4FUwFz1l6z64AzAeeMbOXzGwz\ncBWwoKaNAXuVH+8FrDGzLe0sNMTcgFSrZvAtlFq8C8Zzlh6z85xXi8wAXql6vqz8WjU/BOZKehV4\nBDi/3YXGwb8mSPEu2xW8HpDJwutBrQqes/SYnee8AO5duJnfLtycx6xOAh42sxMkHQjcKukDZvbW\ncGcYYg6CFvAsm5BzNmvIHpVxyEBpqvCdbz6d1Ww5MKvq+czya9WcDfwPADN7TtILwKHAg8Nb4+jK\naJro0ugOHnfLa/Gcp8f8POfVBA8AB0maLWk08Eng2po2LwEfA5A0FZgDPN/OQkPMLRBy7g4e5VKL\n5zw95uc5r6Ews63AF4BbgCeAq8xsiaRzJX2u3OwS4A8kPQrcCnzFzF5vZ7kys3Y+39rCJDvDft61\n5XWCVPuaK/R6t7IVvO2WZ+E5T4/5tZLXPD2Nmamd5UmyR2xO15aXF1Ext0jKVTOkVbl4rPxq8Zyn\nx/w85+WJEPMwSF3OKeFRLrV4lo3H/Dzn5YUYldGHeDhS3goeRxukhMf8uvkdbP6HIHNURk+IinmY\npF41R9WSL5Fn60Rm9Qkxt0HIuXt43CWvxXOeXvPznFkvCTEHyeBVLtV4Fo3X/Dxn1itCzG0SVXN3\n8SqXajxn6jU/z5n1gjj4lwMpX0sD4mBgJ/Ccqdf8OpWZ1x+joYiKOQDSq1hS2Ng8Z+o1P8+ZdZMQ\nc06k3qUB6W0UXuVSjedMvebnObNuEWLOkSLIOcgfz6IJOfskxBzsRGobhFex1OI511Qy7CdCzDlT\nhKrZs0SyCLG0T2ToixiVEWTieVRBFl5HGlTjPdMUMhwOqRUaEBVzRyhC1ZwiKVR93iWRQobdRtLJ\nkpZKelrShXXafF/SM5IWSzq83WWGmDtEEeTsXSJZpCAW77mmkGG3kDSC0s1WTwIOA86SdGhNm1OA\nA83sYOBc4CftLrehmCVdJmmwfHX+em1y/bUI/OBdIqniPdeQ83bmA8+Y2Utmthm4ClhQ02YBcAWA\nmd0PTCjfYmrYNFMxX07p1yKTTvxaFIUiVM3gXyK1pCIV77mmkmOHmQG8UvV8Wfm1odosz2jTEg0P\n/pnZ3ZJmD9Fkp18LSRMkTTWzdM9RDpInlQNZ3g8IFplnFr7Kswtf7fVqZJLHqIx6vxYhZtK/jkaF\nFAUScm6fVDIcinqV/74Dk5g/8IHtz2/+5qKsZsuBWVXPZ5Zfq22zf4M2LdH14XJLLr56++NJA3OZ\nPDC326vQdULOvSMVsXjOtlsZPrJwHY8sfKPjy2mRB4CDyr0GK4BPAmfVtLkWOA/4haRjgHXt9hjk\nIeaWfi1+7+Izc1hk0Cs8C6QeqcjZM93IcN7A3swb2Hv78//9zVeGaN0dzGyrpC8At1A6JneZmS2R\ndG7pbbvUzG6QdKqkZ4G3gbPbXW6zYlZ5yiL3X4siUpSqOegc3n/0+vUHzsxuAg6pee2nNc+/kOcy\nmxkudyVwLzBH0suSzpZ0rqTPlVfoBuCF8q/FT4H/kucKBv7wPpogi1RGGHjPNpUcU6eZURmfaqJN\nrr8WRaVIVbP36i6LVCo+79mmkmPKxLUygmHjXSBZpCIV79mmkiP43wvJIk7J7jJFOekkZVLZHfcu\nlFRyTJEQc9AW3uVRj1Sk4j3fVHJMjRBzDyha1exdHqnjPd+Qc/6EmIO+JSWhhJz7izj41yOKNEID\n/B+sqkdKB7G84zXLFPdQo2IOcsN7VVePVKq9VPMNWifE3ENS/CVvRKryCDnnQyo5eifEHASJEXIu\nPl0X86YVY7q9SNdE1eyHlITiPeOUsvRI18V8+2Gn8PiX5/HOa7t3e9FBF/EujnqkJBTvGaeUpTe6\nLuYTHruJrRtGcduhp/Lk19/Pu+t26/YquKOIVXPKpCSUkHNjVjOxqckTXRfz2Bkbmfejhxh46BY2\nLR/LbQd/nKcumcvmN2PkXtHwLo2h8CCUZvGec0pZtoKkfSTdIukpSTdLmlCn3Z9JelzSo5L+VdLo\nRvPu2cG/PQ94myP+6Xccd8/trH9qPLceeFpJ0G/0ZwVd1KrZuzSC7lBQOX8VuM3MDgHuAP6ytoGk\n6cAXgSPM7AOUzh35ZKMZ93xUxrg56znqX+7juHtu562n9+LWgz7O0r85rG8FHfghJZmk8AOYUp5N\nsgD45/LjfwbOqNNuJLCnpFHAHkDDO8D2XMwVxs1Zz5FX3M9x99zO28+M49YDT2PpXx8WfdAFIAVp\n1CMlmaSQc0p5NsGUyt2azGwl7PoHMLNXgW8DL1O65d46M7ut0YzddexWBP3WM+N46pLDuO2g03jv\nF5/mvec/zei9N/d69TpK0U7TribVU7bB76nGWaSQs5c81y18hDcWPjJkG0m3wk79jAIMuCijuWV8\nfm9KlfVs4A3gakmfMrMrh1yu2S7z6hiS7Az7eUufeeuZcTz13w5j5XXTOeBPnufAP3uKMftt6tAa\n9p6iihlwL4xGeJBJs6SQdaM8/1B3Y2b17jXaFJJsjg0t3wpPa15Ly5O0BBgws0FJ+wF3mtnv1bQ5\nEzjJzP5z+fmngQ82uuuTm66Meow7+C2O/Nn9DDx0C1s2jOT2uaew+E+P4u3nxvV61TpCUQ8CQhq7\n2kUhsu4K1wKfLT/+T8CvMtq8DBwjaYwkAR8FljSasXsxV9jzgLeZ94NFfGzpDYye9A53ffBEHjjr\nWNYt3rvxhwM3pCyM1PpHvWedWp4Z/C1woqSnKAn3WwCSpkm6HsDMfgdcDTwMPEKpK+TSRjN235VR\nj81vjuLFSw/kue8cwvj3v8FBFzzF5BNXorZ2fHxQ5O4MSGM3eyiiSyNf6uXpvSujkyRTMdey2/gt\nHPwXT3Hi89cz4z++zON/cTh3vO8UXvjpgWzZMLLXq9cWRe7OAP+VXCNSqvRSyDqlPLtFshVzLWaw\n+s4pPPe9Obx+7yRm//HzvOe8Z9lj/w0dWV6nKXrVDGlUc0MRlXO+1OaZV8U8ceuyptquGTkzKua8\nkWDyCas45ld3c9xvb2PrppHcOe8kfnfmhxi8eT9sW6/XMAh6R1TOaVGYijmLzW+OYtmVs3nx0gPZ\nvHY0s//4eWad8zxjp6cx3C6qZv+kVDVDGnlXMo2KuaDsNn4L7/nT5zh+0S3Mv/oeNr6yB3ccdir3\nLfg3LP/lTLZuKvR/PwlSqOSGIrUqL4W8U8u0E/SNmfY+ci2H//RBTnrlWqYtWM4LPz6Ym6YvYNE5\n83nt9qnYVhc/lEGCpCaSkLN/3J2S3WlGjdvC7HNeYPY5L7Bx+ViW/2IWj39lHu+sGMt+C5YxbcFy\nJh2/ipG7R6d0t0jhNOIg6CaF7mNuhbee3osV18xgxa9m8OYTE5jyhyuZdvpyJn9ssKengLfTz/zu\n+t0YvHcaG1ftwZxPL81xrTpD6nKO/uZ8OUvX5NLHzEtNXmNn9m5u+pj7rmKux7g56zn4K0s5+CtL\n2TS4Oyuvm8Hyq/fn0S8eyZgZG5h8wiomf3SQiR9+jdH7vtvr1c1kw+AerLp/KivumsnK38xg7ZJ9\nmHzkKmad9kKvV60v8HJxnmaJPRW/RMXcgG1bxBsP78Nrt0/ltdunsvb+iew+dRP7HP06e89fwz5H\nvc5ec9/smKyzKuZtW8T6F8ez9sl9Wb1oCqsfmsLqRVPYunEUk48eZL/jljPtuOVMPnqQUWO2dmS9\nOkURRJGSnMFv5v1cMYeYW8S2ivVL92Lt7yay7oF9WffQvqxfuhcjdt/GXoe+ybhD17PHAW8zdsYG\nxszYyJjpG9l9yiZ2G7+FEaMb91tv2yK2rB/Fu2t2Z9OrYxn96rtsWLEnby8bx5vPTmDdU/vw1ovj\n2WPa2+x96FomHrGKyUeuYtIRqxg3e33yp6R7lUSrhJzbJ8TcrYUVQMxZmME7K8ewful41i8Zz8aX\n9mTjq2PZtLw0vfPa7mxZvxsaYYzaawsjx21BI2z71VvNxNYNI9myfhTb3h3BqL22MHrfdxkzfSN7\nT1/LHtPfZo/pbzPhoHVMOGQt4w98I7lKuBU8SqJVUhMz+Mu9n8Ucfcw5IMGYaZsYM20Tk4/P/nKb\nwbZ3RrBl/W5seWvUjktql78GI/fYUpL22K27VL39cKJJNUXo+0ytvznwRYi5S0gwcsw2Ro55h90n\nv9Pr1Qm6QGpyLsIPYiYr0tNc35xgEqRFCidBNENqJ0oUJfduIOlMSY9L2irpiCHaTZD0fyUtkfSE\npA82mneIOQiCnQg5N81jwL8D7mrQ7nvADeXbTs2jSHcwCfqPoggitaoZipN9JzGzp8zsGbYfKdoV\nSeOBD5vZ5eXPbDGzNxvNO8QcuKYoggg59y3vAVZLulzSIkmXShrb6EPp9YoHQaKkdjCw8Dy0EBYt\nHLKJpFthp1sKidKYqq+Z2XVNLGUUcARwnpk9KOm7wFeBbzT6UBC4prCjBRKgENkvr/P6fgNw6sCO\n55d9c5cmZnZim0tfBrxiZg+Wn18NXNjoQ9GVEQRdJLo0CktmP7OZDQKvSJpTfumjwJONZhZidk6/\nnVxSjyLJIeRcDCSdIekV4Bjgekk3ll+fJun6qqZfAv5V0mJKozL+e6N5R1dGEARNUYhujRwxs2uA\nazJeXwGcVvX8EeDoVuYdFXOQDEWq2lKsmqFYfwPPRMUcBD0iRml0iebuxeqKqJgdE/3Lu1K0ii3F\nyrlofwOPNCVmSSdLWirpaUm7DPWQ9BFJ68oDqBdJuij/VQ2CwAsh587SUMySRgA/BE4CDgPOknRo\nRtPfmNkR5emSnNez74hquT5Fk0KKVTMU7+/giWYq5vnAM2b2kpltBq4CFmS0c3GB6SBIkZBzUE0z\nYp4BvFL1fFn5tVqOlbRY0q8lzc1l7fqUqJYbE0IIikxeozIeAmaZ2QZJp1Aa2zcnq+GSi6/e/njS\nwFwmD4TDqwkp9y+pjtLIa3zzkwtf48mFHfj/r8x/lp2mGTEvB2ZVPZ9JzdnnZvZW1eMbJf1Y0r5m\n9nrtzD5y8Yd3eh7D1XcQUm6NIp7w0M9ynjswmbkDk7c//+U3n2p3tZKlma6MB4CDJM2WNBr4JHBt\ndQNJU6sez6d0k9ddpJzFFAZ3mvqVfv6/BzsT/c1Bw4rZzLZK+gJwCyWRX2ZmSySdW3rbLgXOlPR5\nYDOwEfjEcFeoVlCrdrriXjEJKQ+fIlbNKRN/j3xoqo/ZzG4CDql57adVj38E/CjfVSuRJa2iyDqE\nHNQj1S6NIB+SPCU7dVmHkPOlqFVaqnIu6t+jmyQp5iy8yzpkHPQTruSc4LUyCiPmLIaSYaelHSIO\n8iDVqhmcybkDSPo74N8C7wDPAWfXu9Fq+QzqB4FlZnZ6o3kXWsxDEeIsFkWWQMjZLbcAXzWzbZK+\nBfxlecrifEp3LhnfzIzj6nJBEATDwMxuM7Nt5af3UTrHYxckzQROBf5Xs/MOMQeFocjjaFMd2wzF\n/rtUcQ5wY533vgN8mdLdtZuib7sygiA1oksjZ1YuhMGFQzaRdCvsdEBKlAT7NTO7rtzma8BmM7sy\n4/MfBwbNbLGkAZq82FuIOQiCrtAzOa+o98YATBmoev7NXVqY2YlDzVrSZyl1U5xQp8mHgNMlnQqM\nBfaSdIWZfWao+UZXRlAoir7bnHKXBhTr7yPpZEpdFKeb2TtZbczsr8xslpm9l9LlLO5oJGUIMQdB\ncqQu5wLxA2AccGv5zk0/BpA0TdL17cw4ujKCIOgqLvubh4GZHVzn9RXAaRmv3wXc1cy8o2IOCkeR\ndpfrkXrV3A9/o3boesVc+aWMP0wQtEfKozSgi5Xz8sZNvNGzinkqqwqxOxMEQZA3Pe/KqAg6JB3k\nSb/skUWXRjHpuZirCUEHQf8Rct4VV2KuEFV0EDRP6lUzhJxrcSnmakLSQdCYkHOxSGocc4zoCIKg\nZVb2egVax33FnEVU0UEz9NsPeFTNxSFJMVcTkg6CYhFyLoCYqwlBB/1OEapmCDl3vY+5cqZSJ79A\n1XLu9z9w0H+kfkZg0MODf90QNISkgyBIj56PyuiWoCEkHfQPUTVXsbbXK9A6bvqYJ7K6q1+k6I8O\ngqAdJP21pEckPSzpJkn7ZbSZKekOSU9IekzSl5qZtxsxV6gIuluSjlEdQVEpyoFAx/ydmc0zs98H\nfg18I6PNFuACMzsMOBY4T9KhjWbsTszV9KqKDkkHRSHk3DnM7K2qp3sC2zLarDSzxVXtlwAzGs27\n533MzdDNfugK0R8dBEEjJF0CfAZYBxzfoO0BwOHA/Y3mm9SF8nshaAhJB2kTBwKX1Xn9t+WpPpJu\nBaZWvwQY8DUzu87MLgIuknQh8EXg4jrzGQdcDZxfU2lnL9fMGrXJDUn2cztjp9faFV0vd9VC0r6J\nLqkdpCjm/6pLMTO1Mw9JBq802Xr/YS9P0v7ADWb2/oz3RgHXAzea2feamV/P+5jb7dPtdj90NdEn\nHaRC9DXnj6SDqp6eQan/OIt/Ap5sVsrgQMwVUhY0hKSDoA/5lqRHJS0GPgacDyBpmqTry48/BPwR\ncEJ5WN0iSSc3mrG7g3/tXtqzWs69qhKiTzrwSPQ154uZnVnn9RXAaeXH9wAjW52322tl5HHt5V4d\nLKymtoIOUXeH2HMJUsb9tTLyqD49CLpCVNNBL+nPqnl5r1egZXreldGKNIvQzVFNSLozRLUcpE7P\nxVyhm4JudXndILo88iGk3Jj+rJrTwt0JJv0u6Aoh6iDoX3pWMXsVdLPL7DbR7dGYqJabJ6pm3/Rg\nVMYa1jBx+3Nvgm51mb0gqumdCSEHRaMnFfNE1gCEoHOin0UdUh4+/VM1v9jrFWiZnh78S03QzS63\n12TJqoiyDikHRaX7B//eLMl4cPwOGXdL0EPNr1lSqaJrKZKsQ8hB0em6mEetgi1T6gu6Ws6Qr6Cb\nmV+zpFZXR13PAAAEFklEQVRFZ5GSrEPGnaF/ujPSoiddGaPK21iWoLOqZ/Ar6OEs2zP1BNhtYYeI\ng36m+2J+tfzv9B3VM4SgvdOsKFvNMQTce6Jq9kdTYi5fpu67lC4TepmZ/W1Gm+8DpwBvA5+t3Odq\nF17d+WllBYoo6FaWXxRCtIE/OnutDEl/DvxPYJKZvZ7xfkN/1tLwesySRgA/BE4CDgPOqr3Lq6RT\ngAPN7GDgXOAndWe4GFjFDkG/WppGrdrRxQElQVckDSVBVyRdTaNrILd6nea8r6k81PJfXPhSbstJ\nnchiZyKPNJA0EzgRyPyDNePPLJq5UP584Bkze8nMNgNXAQtq2iwArgAws/uBCZKmksUgJRmvoiTp\nakHjT9B5Sbqy/Op1eCk2vu1EFjvT7Tz6bc8uR74DfHmI95vx5y40I+YZ7HzTrGXsevvt2jbLM9oA\n8Pc/A24GHmZXQVcmOi9or1V0EARpIOl04BUze2yIZs34cxe6fvDvAjP+XuICKN17dpCd70FbzfSd\nR3BASdCNxkCXZt24z3giq1uqFKayKtfRCWPZ0PI6BEGniIOAuzLEXbIvAv6KUjdG9Xv5YGZDTsAx\nwE1Vz78KXFjT5ifAJ6qeLwWmZszLYoopppianRr5qQl/vdjC8la2MN/3ASuB54EXgM3lZU1p1Z9Z\nUzMV8wPAQZJmAyuATwJn1bS5FjgP+IWkY4B1ZjZYO6N2b0UeBEHQCmZ2QIfm+ziwX+W5pBeAI8xs\nbU3TZvy5Cw3FbGZbJX0BuIUdwz2WSDq39LZdamY3SDpV0rOUhsud3eT/LwiCoAgY5a4MSdOAfzSz\n0+r5s9HMVC6vgyAIAic0MyqjZSSdLGmppKclXVinzfclPSNpsaTDO7EeHmiUhaSPSFonaVF5uqgX\n69ktJF0maVDSo0O06ZfvxpBZ9Nt3I9hB7mLO/YSUhGlhcPlvzOyI8nRJV1ey+1xOKY9M+uW7UWbI\nLMr003cjKNOJijnfE1LSptnB5X1zUNTM7gZqD5BU0y/fjWaygD76bgQ76ISYcz0hJXGaHVx+bHm3\n/deS5nZn1dzSL9+NZonvRh/S0zuYBAA8BMwysw3l3fhrgDk9XqfAB/Hd6FM6UTEvB2ZVPZ/Jrpd3\nWg7s36BNEWiYhZm9ZWYbyo9vBHaTtG/3VtEd/fLdaEh8N/qXToh5+4BqSaMpDai+tqbNtcBnAIY6\nIaUANMyiuv9U0nxKQxh3uXRgwRD1+0775btRoW4WffrdCOhAV0ackLKDZrIAzpT0eUqndG4EPtG7\nNe48kq4EBoCJkl4GvgGMps++G9A4C/rsuxHsIE4wCYIgcEZHTjAJgiAIhk+IOQiCwBkh5iAIAmeE\nmIMgCJwRYg6CIHBGiDkIgsAZIeYgCAJnhJiDIAic8f8Bo/2pXY+Wq4gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f26059a4590>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEACAYAAACNlawWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXu4HGWVr99fCEEQiAQkBEIAQwIGA+EWElHTgEjiqOCI\nclFwcIAcBwYHZh5F1GHz6HjEOXMURA8T5CI6GG4OcpVrGgYN4ZaQmDsEAoQkCiEgcktgnT+qeqd3\n775U7/66u6p7vc9TT6qq1/6+r9mVXxar1vqWzAzHcRyn9Qxq9wIcx3G6FRdgx3GcNuEC7DiO0yZc\ngB3HcdqEC7DjOE6bcAF2HMdpEzUFWNJISfdJWihpgaSzKthdLGm5pHmSJoRfquM4TmcxOIHNRuAc\nM5snaWvgMUl3mdmSgoGkacBoMxsj6RDgUmBSc5bsOI7TGdT0gM1sjZnNi89fAxYDu5SYHQ1cHdvM\nAYZKGh54rY7jOB1FXTFgSbsDE4A5JR/tAjxXdL2K/iLtOI7jFJFYgOPwww3A12JP2HEcx2mAJDFg\nJA0mEt9fmtlvy5isAnYtuh4Z3ysdxzeecBwnMWamRn5+B8le6n97pZnt3si4oUgkwMAVwCIzu6jC\n5zcDZwDXSpoErDezteUM2735T09PDz09PW1dQwh6thc92zc+zvLljY/RKBcDZVNrMsYYs455vqD9\nf1ekhrQXgJeAy0runQa7lZlrKvBjoqjA5WZ2Ycnn2wK/AkYBmwH/YWZXNbq+mgIs6VDgi8ACSXMB\nA84j+hJmZjPM7HZJn5T0JPBX4JRGF+a0hjFj0iHCjtMuJA0CLgGOAF4AHpH02+JMLyIHc6GZfUbS\nDsBSSb8ys42NzF1TgM3s90SKX8vuzEYW4rSPMWOiP12InS5lIrDczFYCSJpJlNlVLMAGbBOfbwO8\n1Kj4QhdWwuVyuXYvIQi5LcOPWRDiVnNIe6ZtCp3yfEFnfZcalGZxPU//LK5LgHGSXgCeAL4WYuKk\nMeCOoVMeqtxWzRm3Hd6wC3A66ZTv8h76580OgKOAuWZ2uKTRwN2S9m00I6zrPGAnGe3yhh0nNIcQ\nveQtHGVYRfRyrUC5LK5TgN8AmNlTwNPA3o2uzQXYqciYMS7ETlfwCLCnpN0kDQGOJ8rsKmYl8HGA\nuMp3LLCi0YldgJ2auAg7nYyZvQOcCdwFLARmmtliSdMlnR6bfQ/4sKT5wN3A181sXaNzq5V5uZKs\n3XnAHcPYxnMkB4JnSlRmjD/bQZHUcCGGJFtWcm8sjRd4hMI9YKcu3Bt2nHC4B5xV2uQBF+PecF/c\nAw5LKA/YSpwGLe/vAdeqhIttcsCPgM2BP5vZYY2sDdwDdhrAvWGnEyiqhDsK2Ac4QdLeJTZDgZ8C\nnzKzDwGfDzG3C7DTEJ4p4XQAvZVwZrYBKFTCFXMicKOZrQIwsxdDTOwC7ATBRdjJMEkq4cYCwyTN\nkvSIpJNCTNx1lXBO8/A9JZwOZjBwAHA48F5gtqTZZvZko4M6TlB8hzUnTeR3hfyfqpokqYR7HnjR\nzN4E3pT0ALAf0JAAexZEVklBFkQSukmIPQsiLMGyIE4ouffrvlkQkjYDlhJtR7kaeBg4wcwWF9ns\nDfwEmApsQbS9xHFmtqiR9XkMOKsc1O4FJMNjw07aSVIJF+8NfCcwH3gImNGo+IJ7wNnlxCLH4NH2\nLaMeOt0bdg84LK3ygNuJe8CdgHvDjpNJagqwpMslrY03oSj3+RRJ6yU9Hh/fDr9MpyYHkQkh9rxh\np+VMKDlSRBIP+EqiCpFqPGBmB8TH9wKsyxkoGRBhcBF20oWkqZKWSFom6RtV7A6WtEHS34aYt6YA\nm9mDwMs1zFIRT3Fi3Bt2nMQkKUUusvsB0cu4IISKAU+WNE/SbZLGBRrTaZQMiDC4EDttJ0kpMsA/\nAjcA1bOK6yBEIcZjwCgze13SNOAmorK9svT09PSe53K5juk7lVoKIpyBTAkv4Ohu8vk8+Xy+HVOX\nK0WeWGwgaWfgGDM7TFKfzxohURqapN2AW8xs3wS2TwMHltst3tPQAnLiAKI+GRDhAlkUYk9DC0uo\nNLRZp0O+qHnQBff0K8T4HHCUmZ0eX38JmGhmZxXZXAf8HzN7WNKVwK1mdmMja4PkHrCoEOeVNNzM\n1sbnE4lEveFWHU4TcG/Y6UJyx0Gu6PqCe/qZJClFPgiYKUnADsA0SRvMrLR3XF3UFGBJ1xCtf3tJ\nzwLnA0MAM7MZwLGSvgpsAN4AjmtkQU4LOIjMiDC4EDtNp7cpJ1Ep8vFAn/INM/tA4Tz2gG9pVHwh\ngQCb2Yk1Pv8p0UbFTpZwb9hxgKgUWVKhFLnQEWOxpOlscjT7/Eioub0UOasMJAZciQyIcIE0C7HH\ngMMSrBT53pJ7R3gpspMmMpI3DJ6u5nQWvh+wswmPDTsdiKWs/LgY94Cdvrg37HQhtUqRJZ0o6Yn4\neFDS+BDzugA75cmQCLsQO42QsBR5BfAxM9sP+B5wWYi5XYCzSiv+t8q9Yac7qFmKbGYPmdkr8eVD\n9G/aOSBcgLNMq7bXy5AIuxA7AyBJV+RiTgXuCDGxv4TrBCYA85o8h+cNOxnlivkjeTT/etGdgRfq\nSjoMOAX4SKPrAhfgzqHgCbdCiDMiwuBC7MDBua04OLdV7/WlF/QT4CSlyEjaF5gBTDWzWlv0JsJD\nEJ1GK8ISHht2OoveUmRJQ4hKkfuUGUsaBdwInGRmT4Wa2AW4U/HYcC8eG3aqkaQrMvAdYBjwM0lz\nJT0cYm4vRc4qP6yjkrLZYQnIRFgCmhuS8FLksIQqRX7C+m5Pvp+WeSmy00LcG+7FvWEnTbgAdwse\nG+6Di3D3sIhxfY5yJGnKKeliScvj9mtB/ja5AHcb7g334t6wA8kq4eJ2a6PNbAwwHbg0xNwuwN2I\ne8N9cBHuepI05TwauBrAzOYAQyUNb3TimgIs6XJJayXNr2IT3DV3ahBC3Nwb7sW94a4mSSVcqc2q\nMjZ1k8QDvpLINS9Ls1xzJwEhvEz3hvvgIuy0kiQtiR6MeyVVoo9rLmlocaNOpwWEqE5rVTlzBtLV\nvIqus7g9/15W5ldWM0lSCbcK2LWGTd2EKEWu5Jq7ADeJDRvgqWf3Yu9RSzfdDLFXQyvKmX1PCafF\n7J7bjd1zm3zIBy54sNSkZlNOosq4M4BrJU0C1odwMlu+F0RPT0/veS6XI5fLtXoJmedPf4Ij/uVe\nVl03sv+H7g0Hxb3h1pHP58nn8y2fN0lTTjO7XdInJT0J/JVoQ56GSVQJF//LcIuZ7Vvms0uBWWZ2\nbXy9BJhS7l8Hr4QLw6pVMHG/VeUFuJgQAudVdL3UEmGvhAtLqEq479h5fe59V9/PXCWc4qMcNwMn\nA4R0zZ0AZClTIgMv6TxTwglNzRCEpGuAHLC9pGeB84EhNNk1dwKRldgwZCos4SGJ7LCAIO3bmkKS\nLIgTE9icGWY5TmI2r9M+S7FhSL0Qe2y4e5C0HXAtsBvwDPCFovZEBZuRRNlgw4F3gcvM7OJaY3sl\nXIaxCXW23M5K3jBkIiQBHpLoEs4F7jGzvYD7gG+WsdkInGNm+wCTgTPKNPbshwtwRnnjzS1Z9/J2\nQJ0iDB4bDozHhjueo4FfxOe/AI4pNTCzNWY2Lz5/DVhMgko5F+AMsvPOcNJxV/OhQxfyX9efiJl7\nw47TRHYsJBaY2Rpgx2rGknYn+psxp9bAviF7RrF1Ys6jE5l+zgx22nENP/33f2D0HisAUL1x2qyk\nq0G6Y8PL/NkOSag0tENnfYcX84t67y294MZ+40q6myh+23sLMODbwFVmNqzI9iUz277CfFsDeeC7\nZvbbmutzAc4mti56fjZsGMyP/t/Z/PDib3D6l/+Tb/7T/2abbV4DOliI0yrCLsBBCSXAY+2JPveW\nab+6xpW0GMiZ2VpJOxHVPXywjN1g4FbgDjO7KMnYHoLIKAuGRW1WNt98I18/69954oF9WfXCSD44\naQm/+PXJvPuuPDbsOGG4Gfi7+PzLQCXP9gpgUVLxBfeAM8t89uo9H79uWe/5Q48cwj+ddxFm4t++\nfR5HTLkXyb3hluAecFBS5AEPA64j2oxnJVEa2npJI4jSzT4l6VDgAWABUejCgPPM7HdVx3YBzibF\nAgx9Rfjdd8W1/30cPRdewC4jVvHd877NoYf8oXNFGNIhxC7AQUmLADcTD0F0CAuGje0NSwwaZJzw\nuZks/MM4vvSFX/Kl6f/Fx4+5h9v+PI139q3jufNMCcdpKu4BZ5RSD7iYYm8Y4O23N+e6m77A//3Z\nP/Pmm+9h+imXcuLYX7Pjdn9OPqF7w7VxDzgooTxgnnq3783Rg1LjAbsAZ5SZfBaAcSyqaFMqxGaQ\nfzDHlf/1FW6589N8ZNKDfGrcrRx10J3sMeKZZBPH4vaXt7Zm1lOHceeyo7hr2Se46eRj2Genymvp\nQ6fGhl2Ag5IWAU5SilxkO4jo6XvezD5Tc2wX4GxSEGCoT4QLvPbae7n5d5/hd/dM4667P8EWQ95i\nwuh5jN9jAbvssIrh261lyy3ewEy8vXEIa9btxOqXRrD0+b2Y98cJPP/KSCbvNpujxt7JJ8bcxb4j\n5qN6/qp0ojfsAhyUFAnwhcBLZvbDuGX9dmZ2bgXbs4EDgW1dgDuYYgEuMBAhhuil3dMr9+CJWyaw\n8Jl9WL1uBGtfHs4bb22JZGw+eAM7bbeGEduvZs+dn2T/MXPZa9elDJ77TuNfpJO8YRfgoKRIgHv3\nOI/zgPNm1m+fh3hDniuBfyPaF8IFuFP5V75VVnAHKsIFPFOiAVyAg5IiAV5XUgnX57ro/vVE4jsU\n+OckAtzylkROOBYxrp/gLmIcUF6IC1kS1YS4ULyRWIh9v2En7TyUhzn5qiY1SpFL6fcvraS/Adaa\n2TxJOSo3sOj7c+4BZ5N/5Vt9rt0brpNmCLF7wEEJ5gHPLPm9HF/fuElKkSV9H/gS0baUWwLbAL8x\ns5OrjZ0oD1jSVElLJC2Lg9Cln0+RtF7S4/FR7l8NJyClu/wXPN/Se+XuQ9+84Ur4DmuOAyQoRTaz\n88xslJl9gKir8n21xBcSCHCcVnEJcBSwD3BChY2GHzCzA+Lje7XGdRqnnAhXEuKKY9QQYfD9hp2u\n50LgSElLgSOAHwBIGiHp1kYGrhmCiBttnm9m0+Lrc4n6wV1YZDMF+Bcz+3SNsTwEEYjPMrPP9XgW\n9LmuFH7wsEQFQqzTQxBBSUsIopkkCUHsAjxXdP085Xd6nyxpnqTbJFV2uZwglHq1qfaGQ4Qlmo17\nw04bCLUXxGPAKDObQBSuuCnQuE4VSkV2AePTGRsGjw077WNOyVEnkraTdJekpZLulDS0gt1QSddL\nWixpoaRDao2dJA1tFTCq6HpkfK+XuAdS4fwOST+TNMzM1pUO1tPT03uey+XI5XIJluCUsmzFeMZ+\nIAo7lKajLWB8n5BEpdS0cmlsvWMMG1szJGET6gxJhEpZ887MHUk+nyefz7d7GeUoNOUsVMJ9M75X\nykXA7Wb2+Xhz9q1qDZwkBrwZUAg+rwYeBk4ws8VFNsMLPZMkTQSuM7Pdy4zlMeBAaMWm/44FIYb+\nIuux4QapZ50eAw5KsBjw2SW/lx/VnYZWsxJO0rbAXDMbXc/6aoYgzOwd4EzgLmAhMNPMFkuaLun0\n2OxYSX+UNBf4MXBcPYtwBsAjm56fZSs2hR3aGRtuS8pas/HYsJOsKecewIuSroxTcWdI2rLWwF6I\nkVF0bdHFwe4Np8Ibdg84KK30gBttyinpQOAhYLKZPSrpx8ArZnZ+1fW5AGcTnROfFML8AxThcjaV\n7vWOkUCEoQuF2AU4KMEEePwseCW/6eazFzSjEm44MDsuxEDSR4Bv1EzNdQHOJvoYm/7XuPhdq3vD\n7RNhF+CgBBPgj5b8Xv6n7hjwhcA6M7uw2naUku4HTjOzZZLOB7Yys36Vw31+xgU4m2hzYHJ8USrE\nLfSGwYW4FxfgoKRIgGs25Yzt9gN+DmwOrABOqbRxe+/YLsDZRFoOg8dUFmFwbxhaK8IuwEFJiwA3\nExfgjCL9HJgSiTC4N1yLVgjxNf5sh8QFOPRkLsBB6RVhcG84Cc0WYRfgoLgAh57MBbgptMsbrnYf\nulCIXYCDEkyAB5f8XjbWHQNO1JQz7gf398C7wAKiGPDbVcd2Ae4M+nnDEAmxe8PlaYYIuwAHJUUC\nXLMpp6SdgQeBvc3sbUnXAreZ2dVVx3YB7izcG66TkELsAhyUFAlwklLknYHZRPWZfwH+G7jIzO6p\nOrYLcOdR0xuGhoXYRbgMLsBBSZEAJ23KeRZRU87XgbvM7KSaY7sAdy7uDddJo0LsAhyUYAI8aBZY\nftNN618JF6AU+X3AjcDngVeAG4DrzeyaqutzAe5s3Buuk0ZE2AU4KMEEmNJnbWwzSpGPBY4ys9Pi\n65OAQ8zszGpjh9qQ3UkpZqdiNga4HzYuj27OJhKrgmAVNqmussNa6cbvxVTbYc2bgjodQM2mnMCz\nwCRJ75Ekou17F5ex64N7wF1EK7zhcjaV7vWO0SnesHvAQUmRB5y0FPl8oo7IG4C5wKlmtqHq2C7A\n3Ufi2DD0CnGxCENfQfXYcIwLcFDSIsDNxAW4S4lEGNJWzpxpEXYBDooLcOjJXIBTx0DKmd0broAL\ncFDCCfBlJXdPS40AJ3oJJ2mqpCWSlsWVIOVsLpa0XFFr+tS+xkhp07+6CfU9zE4F7qf3Jd1s+r6k\nK+4kG7+kW7ZifMU2SJU6M1drg/RI/vV+nyVtg1QXIVoL1ZizU54v6Kzv0giSCi3X3pF0QBW7mjpZ\nSk0BljSIqNX8UcA+wAmSSqtApgGjLXrdPh24NMnk7aBTHqqQ3yPKlIiFeOPyTUIM5TMlioS4QK1M\niYJNuXt35N9bdl1ZzJTolOcLOuu7NMgC4LNEnkpZkuhkOZJ4wBOB5Wa2Mn6jNxM4usTmaOBqADOb\nAwyNW3Q4GSIN3nA5suoNO52BmS01s+VExRmVSKKT/UgiwLsAzxVdPx/fq2azqoyNkwGa5Q2XE+JS\nPG/YyTBJdLI/Zlb1AD4HzCi6/hJwcYnNLcCHi67vAQ4oM5b54YcffiQ9aulTAv16psy4a8rY3Q3M\nLzoWxH9+ushmFmV0LalOljsGU5tVwKii65HxvVKbXWvYpObNo+M43YGZ7Z7Q7sgGp0qik/1IEoJ4\nBNhT0m6ShhBVetxcYnMzcDKApEnAejNbm2TVjuM4GaKSE5lEJ/tRU4DN7B3gTOAuYCEw08wWS5ou\n6fTY5nbgaUlPAv8J/EOir+I4jpNyJB0j6TlgEnCrpDvi+yMk3QqVdbLm2F4Y4TiO0x66bje0pEnV\naWYgCd9pRNLlktZKmt/utTSCpJGS7pO0UNKCeGPuTCJpC0lzJM2Nv8v57V5TJ9N1AkyCpOo0M9CE\n75RyJdH3yDobgXPMbB+iYu4zsvo7MbO3gMPMbH+iBLtpkia2eVkdS9cJcMKk6jQzoITvNGJmDwIv\nt3sdjWJma8xsXnz+GtE+sJnNgzezQm34FsBgotQtpwl0nQB3AANL+HZagqTdiTzHOdUt04ukQZLm\nAmuAu83skXavqVNJkgecOar0d/qWmd3SnlU5nY6krYl6gX0t9oQziZm9C+wvaVvgJknjzKzy1nXO\ngOlIAQ6QVJ1mBpTw7TQXSYOJxPeXZlauZU3mMLNXJc0CpkKVvUOdAdPtIYgsxoEHlPCdYkQ2fw+l\nXAEsMrOL2r2QRpC0g6Sh8fmWwJHAkvauqnPpOgGulFSdFQaa8J1GJF0D/AEYK+lZSae0e00DQdKh\nwBeBw+P0rcclTW33ugbICGCWpHlEcew740Irpwl4IYbjOE6b6DoP2HEcJy24ADuO47SJJC2JEpVZ\nZqUnnOM4TlpIkoZWKLOcF+c5PibpLjPrfTNa3BNO0iFEPeEmNWfJjuM4nUGS7SiTlFl6TzjHcZw6\nqSsGXKXM0nvCOY7j1EniSrgQZZaSPOfNcZzENNrGbAfJXup/e2XSVkXNJpEAJyizTNQTDqDdecc9\nPT309PS0dQ0h6DlS9BwJzGv3SirwaG2TAj0vQc/2zVtKy1hmHfN8Qfv/rkiNF0i+BFxWcu802K3M\nXFOBHxNFBS43swtLPt8W+BXRNgCbAf9hZlc1ur6kIYhaZZbeE65dpDXfpNF2747TIhLusX0GsNDM\nJgCHAf8RO6YNUXOAojLLBfEWdQacR/SviJnZDDO7XdIn455wfwUyWVKaWSaQTk+4IMJ1eMOO0wZ6\n99gGkFTYY7t4DwwDtonPtwFeMrONjU5cU4DN7PdELnctuzMbXUwryOVy7V5CEHIfKLlR8ITTKsRV\nRDi3ZctW0nQ65fmCzvouNSi3x3ZpF5BLgJslvQBsDRwXYuKuq4TrlIcqN7rCB2kOSVQIS+S2aulK\nmkqnPF/QWd8lAEcBc81sZ2B/4KdxYkJDdJ0AdwUTSLcQO04LeQ/wRNFRhiR7bJ8C/AbAzJ4CngYa\n7vvXkRuyOzEeG3YcDomPApf0N+ndYxtYTbTH9gklNiuBjwO/j4vMxgIrGl2bC3Cnk+HYsOO0AjN7\nR1Jhj+1CGtpiSdOJEw2A7wFXSZof/9jXzWxdo3O3dD9gSdbuPOCO4YcDyJFMowgX6AQhXubPdkgk\nNVyIIcmWldwbS+MFHqHwGHA34bFhx0kVLsDdSJpF2IXY6SJcgLsV94adLmHMmL5HOSRNlbRE0jJJ\n36hgk4t7/v0x7hbdMP4SrtvxTAmnyykqRT4CeAF4RNJvS/Y8Hwr8FPiEma2StEOIud0Ddtwbdrqd\n3lJkM9sAFEqRizkRuNHMVgGY2YshJnYBdjaRZhF2IXaaR7lS5NL9zMcCwyTNkvSIpJNCTOwhCKcv\nnjfsdBj51yH/RsPDDAYOAA4H3gvMljTbzJ5sdFDH6Y/Hhp0OIfcxyBVdX/DrfiZJSpGfB140szeB\nNyU9AOwHNCTAHoLIKq34X3KPDTvdQW8psqQhRKXIN5fY/Bb4iKTNJG1FVN28uNGJ3QPOMq3yBt0b\ndjqYJKXIZrZE0p3AfOAdYIaZLWp0bi9Fzir3lVRStkqE0ijEBdotxF6KHJRQpchWsq2Ofp2hUmRJ\nl0taW7QJRennUyStl/R4fHw7/DKdmrQqUyCtIQnwsISTOZKEIK4EfgJcXcXmATP7TJglOQ3RikyB\ntGdKQPu9YSc9lDoN/V/CtY2aHrCZPQi8XMMsFe68E+PesHvDTl0kKUWO7Q6WtEHS34aYN1QWxGRJ\n8yTdJmlcoDGdRvFMCRdipyYJuyIX7H4A3Blq7hAC/BgwKm7XfAlwU4AxnVC4N+wi7NQiSSkywD8C\nNwB/CjVxw2loZvZa0fkdkn4maVil3eJ7enp6z3O5nDf+axUeG47w2HBqyefz5PP5dkxdsyuypJ2B\nY8zsMEmlHZMHTFIBFhXivJKGm9na+HwiUWpbxVYdxQLstJhW5g1DeoXYRTiVlDpkF1xwQZBx809B\nvuHubfwYKI4NB3nvVVOAJV1DVMm3vaRngfOBIWzqlXSspK8CG4A3gONCLMxpIq0SIS/gcFJA7riS\nUuR7+pkkKUU+CJgpScAOwDRJG8ystGKuLrwQI6uUFmIMlG4v4Aj5/b0QIyjBCjHuLbl3RN9CDEmb\nAUuJ9gNeDTwMnGBmZUuNJV0J3GJmv2lkbeB7QWQWC/XSq1UvqNL6ks4zJboeM3sHKJQiLwRmFkqR\nJZ1e7kdCze0ecEaxdZscA4XyLt0bbgz3gIPSKg+4nbgH3AG4NxwI94adFuMecEYp9oCLcW84EAP5\n/u4BByWUB/zuS33vDdq+vwcsaSpRpkNhN7QLSz4/kU1ZEH8BvmpmCxpZG7gH3HG4NxwI94a7hoSV\ncCuAj5nZfsD3gMtCzO0C3IHYhEBC3MoqujQLsdPp1KyEM7OHzOyV+PIh+veMGxAuwB2Me8OBcG+4\n00nSlLOYU4E7QkzsHTE6nIIINxwb9io6r6JzkHQYcArwkRDjuQB3CTYh0As6b4MU4UKcGa6YP5JH\n868X3em3U0KSSjgk7QvMAKaaWa0tehPhWRAZZT57MX7dsgH9rGdKBKT4v4FnQQQlVBbEEza2z739\ntKzuSjhJo4B7gZPM7KFG1lSMx4AzzIJhY2sblSFobNi3unQyTsJKuO8Aw4CfSZor6eEQc7sHnFHm\ns1efa/eGU8A1/myHpFUecDtxD7hDcG/YcbKHC3AHsWDY2AEJcbC8YfC8YcepAxfgDqQRIQ6Ce8NO\niljEuD5HOZI05ZR0saTlcf/LIE+fC3AH496w49QmSSmypGnAaDMbA0wHLg0xtwtwRqn0L3kpXeUN\nuxA7AyNJU86jgasBzGwOMFTS8EYnrinAki6XtFbS/Co2wV1zpzZJRRi6xBsGF2FnICQpRS61WVXG\npm6SVMJdCfyEWP1LKXbNJR1C5JpPanRhTjIKIjyORTVtCyJcb8paJqvoIN0pa07LWJT/M4vyL7Z7\nGWWp6QGb2YNAtbK7prjmTn24N1wG94YdgNwUxvV8rvcoQ5JS5FXArjVs6iZEDLgprrlTmddeg1lX\nHNHvfrW3vKV4bNhxenkE2FPSbpKGAMcDpd2ObwZOBpA0CVhvZmsbnbjlm/H09PT0nudyOXK5XKuX\nkHleeQWu/84JHPaVe8t+vohxiUISEAnxQEISEDAs0YoqurRu7uP0ks/nyefzLZ/XzN6RVChFLnTE\nWCxpevSxzTCz2yV9UtKTwF+JdkRrmESlyJJ2I2rDvG+Zzy4FZpnZtfH1EmBKuX8dvBQ5DKtWwQcn\nvsqMVV+uaZtUiGFg5czBSpkh++XMXooclFClyN+x8/rc+66+n7lSZMVHOZrimju1SRJy8NhwGTwk\n4aSEmiEISdcAOWB7Sc8C5wNDaLJr7iSnVsihVZkS4Bu/O+ljAePbvYSK1BRgMzsxgc2ZYZbjDJQk\nIluvEA8kJBE0Zc1jw04KkLQdcC2wG/AM8IWi/nAFm5FE2WDDgXeBy8zs4lpjeyVch5Ek5NCKTAlv\nCup0EOcC95jZXsB9wDfL2GwEzjGzfYDJwBllOiv3wwU4gwwZAq+/siWzHji87OdJY8P1CPFA8Niw\n0yEcDfzDVVmPAAAMVUlEQVQiPv8FcEypgZmtMbN58flrwGISpOO6AGeQ978fDvzVQzz2xclcfepX\neGzdgWXtQr6kc2/Y6WJ2LCQWmNkaYMdqxpJ2J3qK5tQa2JtyZpSdj1nF+w9fy6Jv7cu9+0xj7b+P\nYNoXb0VlclVCv6Rre2wYurcpqFM3f84v4sV89Wdb0t1E8dveW4AB3y5jXjHfUNLWwA3A12JPuPq8\n3pIom+zF/F7BXPfQ9sw/80AGDXmX8Rc9zscOvr/iz9US2WbnDUMHt0HyPOCghMoDHmtP9Lm3TPvV\nNa6kxUDOzNZK2omo7uGDZewGA7cCd5jZRUnG9hBEhil4rsMmvcSUh+9it9OeYs7RH+WqvzuVv6ze\npurPVBszU7Fh3/jdaT43A38Xn38Z+G0FuyuARUnFF9wDzix70Xd30ILnuuHVwSz7/jieuWw0e/yv\nJ/nMv/yGLbd7s+wY7g0PgGrrdQ84KCnygIcB1xFtxrOSKA1tvaQRROlmn5J0KPAAsIAoRGHAeWb2\nu6pjuwBnk1IBhr6C+frKrVj63X1YfdNIPvxPDzDxrNlsse1bVX+mEkmFuO0iDK0R4krrdQEOSloE\nuJm4AGcUrYj+O479wIJ+nxUL5mvLt2bJBR9i7R0jOOgrD3PwP85m6KhXqv5MOdwbLkPpel2Ag+IC\nHHoyF+BgFAQYaoswwF+feS8rfjKGF64axa4ffYZ9T57LmL9ZyuAt3qn4M+UoZ7PmuZ3YfqcX2Xzz\njb332i7C0HohdgEOSigB5ql3+94cPajeEETNSrgi20FET97zZvaZmmO7AGcTXQsc3Pe/ZRIh3vDq\nYF64YVfW/fL9/Gn+Tuzx8ScZfdSTjJryNNvt8TIaZFWF2Ax2eOFFFj++Dw/fO5nZdx3Kuj9tzy9+\nfwK77/V0P/u2C3ErRdgFOCgpEuALgZfM7Idxx+TtzOzcCrZnAwcC27oAdzC6Nj4ZgAgX2G3VSlbc\nOYan7tyT52eP4q1X3sP791nLNru8yq7Dn2PIlm9jBhvfHsz6NduxfvV2rF66C2aw+/5PM2XKLD58\n1IPsvf8iNtvs3bJzFOiKrS6/7s92SFIkwL1b7MZpaHkz61dmHO8HcSXwb0RlyS7AnUqvABdoQIjH\nE9m+sW5L/rxwR/6yehv+unZrNr6xOcO1ls0238j7dlrP+0a8zE57rma7ndf1Fnx4bLgIF+CgpEiA\n15nZsErXRfevJxLfocA/uwB3MDonPjmk6GYAES5HyEwJ6GBv2AU4KK0U4BqVcFeVCPBLZrZ9yc//\nDTDNzM6UlCMS4E/XXJ8LcDbpFWDoK8LQNiHuem/YBTgowQT4O7NgUX7TzRsvCF4JJ+n7wJeIdkXb\nEtgG+I2ZnVx1bBfgbKKP0b8KzL3h+gkpxC7AQQkmwDNLfi/H1zdu/BJunZldWOslXGw/hYQhiESl\nyJKmSloiaVm8gH4TSlov6fH4KLeBhROaR+krIMV7Lz2i6IhZtmI8y1b07QxQqex4AeMrdhHwNkhO\nF3IhcKSkpcARwA8AJI2QdGsjA9f0gOO8tmXxxC8QtXA+3syWFNkkUnz3gMOhzeOTyfGfdXjD0N8j\ndm+Yxr1h94CDkhYPuJkk8YAnAsvNbKWZbQBmEm1QXEoqvlDXsHF59Ofs+LoObxhwb7gc7g07LSaJ\nAO8CPFd0/Tzld3qfLGmepNskJf9b5wycYhEuFuICc+gvxEXUG5aoROgd1rpi43fHIdyG7I8Bo8zs\ndUnTgJuAsn+Lenp6es9zuRy5XC7QErqN+4Epm0R48JhIhCezSYQLQjKHTSGJgggXhSWWrRjfLyRR\nbhP3ggiXC0skbQqaNCSRio3fW1VF5wCQz+fJ5/PhB67Zl6I6SUuRJQ0Ffg58iKgx51fMrOrsSWLA\nk4AeM5saX59L1JL+wio/8zRwoJmtK7nvMeCASD+Pz6ZEfwwes+nDQLFhKC+qaY8NQxtS1jwGHJRg\nMeCzS34vPxpQFkTNUmRJVwH3m9mV8ebsW5nZq1XHTiDAmwGFt3+rgYeBE8xscZHN8ELPJEkTgevM\nbPcyY7kAB6afCMMmIZ5cZFgsxCnJG05qAykQYagtxC7AQUmRANcsRZa0LTDXzEbXtb4kgihpKnAR\nUcz4cjP7gaTpRJ7wDElnAF8FNgBvAGeXc71dgJtHIm8YKguxe8PJqCbCLsBBSZEA1yxFlrQfMANY\nBOxH9KR8zczeqDq2F2J0Du4NV6fp3rALcFCCCfCxs+D5/KabD/WvhAtQinwg8BAw2cwelfRj4BUz\nO7/q+lyAO492esPQxeXMLsBBCSbAHy35vfxP3R5wklLk4cBsM/tAfP0R4Bu19oPwppwdiNmp8Vnc\nHXnj8v55w9BQ3jBUTj9rZcraQMhcU1Cn3dRsyhm/A3tOUuGhPAJqexPuAXc4DXvD4Jv7JOFR3AMO\nTIo84JpNOWO7/YjS0DYHVgCnlEtX6zO2C3B3UFGIUy7CSW0gBSJ8uD/bIUmLADcTF+Auwr3h6jQs\nxC7AQXEBDj2ZC3AqcG+4Mg2JsAtwUIIJ8OCS38tGF2Cnzbg3XJ0BCbELcFDSIsB1lCKfDfw9URny\nAqIY8NtVx3YB7m66wRuGFm116QIclBQJcM1SZEk7Aw8Ce5vZ25KuBW4zs6urju0C7Lg3XJ3EQuwC\nHJQUCXCSUuSdiZI8JwB/Af4buMjM7qk6tguwU8C94cokEmEX4KCkSICTdkU+i6gr8uvAXWZ2Us2x\nXYCdYtwbrk5VIXYBDkowAR40Cyy/6aY1pRT5fcCNwOeBV4AbgOvN7Jqq63MBdspRtzcM3hTUBTgo\nwQSY0t/x2GaUIh8LHGVmp8XXJwGHmNmZ1cb2UmSnLGXLmaF/9w1vCup0PjVLkYFngUmS3iNJRKXI\ni8vY9cE9YKcmrdj4PaveMBR5xO4BByVFHnDSUuTzgeOJtuWdC5xqUR/NymO7ADtJyPJWly0TYRfg\noKRFgJuJC7BTF94GqTIa5s92SLpBgBPFgCVNlbRE0rI4EbmczcWSlsedkT1C1qEk2uqyjtgwJN/q\nMkRsuNlbXTpp5P6Soz4kHSvpj5LekXRAFbuaOllKTQGWNAi4BDgK2Ac4QVJpEvI0YLSZjQGmA5cm\nmbwdNKXrahto5/cwOzUW4qIHutJLugJz6C/EBR7K1/2SrhJJBLYeEa5XiDvl+YLO+i4NsgD4LFXU\nO4lOliOJBzwRWG5mK+OA8kzg6BKbo4GrAeJecEPjHeJTR6c8VGn4HsG84Tn53ltZ94bT8HsJRSd9\nl0Yws6VmtpwoN7gSSXSyH0kEeBfguaLr5+N71WxWlbFxOpAg3vCq/iGJUN5wqJS1gXjDTleRRCf7\n4XnAThDSGhuu9DOln3tsuLuRdLek+UXHgvjPqj3dGsbMqh7AJOB3RdfnEjWbK7a5FDiu6HoJMLzM\nWOaHH374kfSopU8J9OuZMuOuGeBYs4ADBqqT5Y7B1OYRYE9JuwGriRKNTyixuRk4A7hW0iRgvUVN\n6vqQltQPx3G6AzPbPfCQlTQsiU72o2YIwszeAc4E7gIWAjPNbLGk6ZJOj21uB56W9CTwn8A/JPkm\njuM4aUfSMZKeI/Jyb5V0R3x/hKRbobJO1hzbCyMcx3HaQ9e9hEuaVJ1mBpLwnUYkXS5praT57V5L\nI0gaKek+SQvjlzdntXtNA0XSFpLmSJobf5fz272mTqbrBJgESdVpZqAJ3ynlSqLvkXU2AueY2T5E\nBdlnZPV3YmZvAYeZ2f5E3R2mSZrY5mV1LF0nwAmTqtPMgBK+04iZPQi83O51NIqZrTGzefH5a0Tb\nEGY2D97MXo9PtwAGE2UOOE2g6wS4AxhQwrfTGiTtTuQ5zqlumV4kDZI0F1gD3G1mj7R7TZ1KkjS0\nzFGlvci3zOyW9qzK6XQkbU3UiuZrsSecSczsXWB/SdsCN0kaZ2bJt5JzEtORAmxmR7Z7DU1kFTCq\n6HpkfM9pI5IGE4nvL82sXMeEzGFmr0qaBUyFOvbydBLT7SGILMaBexO+JQ0hSvi+uc1ragSRzd9D\nKVcAi8zsonYvpBEk7SBpaHy+JXAkUWWr0wS6ToArJVVnhYEmfKcRSdcAfwDGSnpW0intXtNAkHQo\n8EXg8Dh963FJU9u9rgEyApglaR5RHPvOuNDKaQJeiOE4jtMmus4DdhzHSQsuwI7jOG3CBdhxHKdN\nuAA7juO0CRdgx3GcNuEC7DiO0yZcgB3HcdqEC7DjOE6b+P+85iK+m8LjgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2605869410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "cs = []\n",
    "\n",
    "for i in range(3):\n",
    "    cs.append(np.random.randn(2,1))\n",
    "ls, storedls = findLambdas(cs)       \n",
    "\n",
    "other = ls[2]\n",
    "xs = np.arange(0,2,0.02)\n",
    "ys = np.arange(0,2,0.02)\n",
    "X, Y = np.meshgrid(xs, ys)\n",
    "Z = []\n",
    "for xlist,ylist in zip(X,Y):\n",
    "    zlist = []\n",
    "    for x,y in zip(xlist,ylist):\n",
    "        zlist.append(L(np.array([x,y,other]),cs))\n",
    "    Z.append(zlist)\n",
    "Z = np.array(Z)\n",
    "\n",
    "CS = plt.contourf(X, Y, Z, np.arange(-4,4,0.2))\n",
    "cbar = plt.colorbar(CS)\n",
    "plt.plot(storedls[:,0],storedls[:,1])\n",
    "\n",
    "\n",
    "storedls = np.array(storedls)\n",
    "plt.figure()\n",
    "\n",
    "for dim in range(2):\n",
    "    plt.subplot(2,1,dim+1)\n",
    "    Z = []\n",
    "    for xlist,ylist in zip(X,Y):\n",
    "        zlist = []\n",
    "        for x,y in zip(xlist,ylist):\n",
    "            #try:\n",
    "            grad = dL_dl(np.array([x,y,other]),cs)\n",
    "            zlist.append(grad[dim])\n",
    "            #except:\n",
    "            #    zlist.append(np.nan)\n",
    "\n",
    "        Z.append(zlist)\n",
    "    Z = np.array(Z)\n",
    "\n",
    "    CS2 = plt.contourf(X, Y, Z, np.arange(-1,1,0.1))\n",
    "    cbar = plt.colorbar(CS2)\n",
    "    plt.axis('equal')\n",
    "    plt.plot(storedls[:,0],storedls[:,1])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top figure is the values of L (over two values of $\\lambda$).\n",
    "The bottom two figures are the gradients of $L$, over two values of $lambda$. The first is the gradient wrt to $lambda$ on the x-axis, and the second wrt to the one on the y-axis. So you can see how it follows the gradient down (in this case) to the maximum of $L$, but with the constraint (artificially imposed) that the lambdas must be greater than (or equal to) zero.\n",
    "\n",
    "So in simple cases it seems to work!\n",
    "\n",
    "Below are examples of $c$ from real covariance functions, to investigate what structures they might have which cause a problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Problem c values\n",
    "\n",
    "The problem values of $\\mathbf{c}$ previously seemed to occur when I used 'realistic' values. But it may be the new code solves the previous problems. In particular I must have improved things as the new code doesn't require the sum of the lambdas to be enforced (to equal the rank).\n",
    "\n",
    "The vectors are the columns of:\n",
    "\n",
    "$$K_{*f} K^{-1}$$\n",
    "\n",
    "So we can produce some;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(suppress=True, precision=5)\n",
    "\n",
    "\n",
    "def k(x1,x2):\n",
    "    return np.exp(-(x1-x2)**2)\n",
    "\n",
    "X = np.arange(0,60,3.9)\n",
    "X[-1] = X[-2]+1e-10\n",
    "X[0:5] = X[10:15]\n",
    "Xstar = np.arange(50,55,0.2)\n",
    "\n",
    "Kstarf = np.zeros([len(Xstar),len(X)])\n",
    "for i,xstar in enumerate(Xstar):\n",
    "    for j,x in enumerate(X):\n",
    "        Kstarf[i,j] = k(x,xstar)\n",
    "        \n",
    "K = np.zeros([len(X),len(X)])\n",
    "for i,x1 in enumerate(X):\n",
    "    for j,x2 in enumerate(X):\n",
    "        K[i,j] = k(x1,x2)\n",
    "        \n",
    "K = K + np.eye(len(K))*0.0001\n",
    "C = np.dot(Kstarf,np.linalg.inv(K))\n",
    "\n",
    "cs = []\n",
    "for i in range(C.shape[1]):\n",
    "    cs.append(C[:,i][:,None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambdas:\n",
      "[ 0.       0.       0.5      0.5      0.50015  0.       0.       0.       0.\n",
      "  0.       0.       0.       0.5      0.5      0.49985  1.     ]\n",
      "Solutions to c^T M^-1 c (should all be 0<=x<=1)\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.00001]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 0.]]\n",
      "[[ 1.]]\n",
      "[[ 1.]]\n",
      "[[ 1.00001]]\n",
      "[[ 1.]]\n"
     ]
    }
   ],
   "source": [
    "ls, storedls = findLambdas(cs)\n",
    "print \"Lambdas:\"\n",
    "print ls\n",
    "M = calcM(ls,cs)\n",
    "Minv = np.linalg.pinv(M)\n",
    "print \"Solutions to c^T M^-1 c (should all be 0<=x<=1)\"\n",
    "for l,c in zip(ls,cs):\n",
    "    print np.dot(np.dot(c.transpose(), Minv),c)"
   ]
  },
  {
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}