jupyter_notebooks/CR_softmax (3).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import gc, argparse, sys, os, errno\n",
    "%pylab inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import h5py\n",
    "import os\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import scipy\n",
    "import sklearn\n",
    "from scipy.stats import pearsonr\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/chenxupeng/projects/DIP\n"
     ]
    }
   ],
   "source": [
    "cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [x] implement  https://github.com/Silver-Shen/Causally-Regularized-Learning\n",
    "- [x] test on training data\n",
    "- [ ] multi class\n",
    "- [x] numba accelarate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes:\n",
    "$$\n",
    "\\begin{array} { c } { P ( Y = k | X = x ) = \\frac { \\exp \\left( \\beta _ { k , 0 } + x ^ { T } \\beta _ { k } \\right) } { 1 + \\sum _ { l = 1 } ^ { K - 1 } \\exp \\left( \\beta _ { l , 0 } + x ^ { T } \\beta _ { l } \\right) } } \\\\ { \\text { for } k = 1 , \\ldots , K - 1 } \\\\ { P ( Y = K | X = x ) = \\frac { 1 } { 1 + \\sum _ { l = 1 } ^ { K - 1 } \\exp \\left( \\beta _ { l , 0 } + x ^ { T } \\beta _ { l } \\right) } } \\end{array}\n",
    "$$\n",
    "\n",
    "$\\beta$ should be a matrix with **feature_num** * **class_num shape**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "outputs": [],
   "source": [
    "def oneHotEncoding(y, numOfClasses):\n",
    "    \"\"\"\n",
    "    Convert a vector into one-hot encoding matrix where that particular column value is 1 and rest 0 for that row.\n",
    "    :param y: Label vector\n",
    "    :param numOfClasses: Number of unique labels\n",
    "    :return: one-hot encoding matrix\n",
    "    \"\"\"\n",
    "    y = np.asarray(y, dtype='int32')\n",
    "    if len(y) > 1:\n",
    "        y = y.reshape(-1)\n",
    "    if not numOfClasses:\n",
    "        numOfClasses = np.max(y) + 1\n",
    "    yMatrix = np.zeros((len(y), numOfClasses))\n",
    "    yMatrix[np.arange(len(y)), y] = 1\n",
    "    return yMatrix\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    #z = np.dot(X,W)\n",
    "    exp_x = np.exp(x)\n",
    "    sum_exp_x = np.sum(exp_x,axis=1)\n",
    "    sum_exp_x= sum_exp_x.reshape((x.shape[0],1))\n",
    "    sigmoid = exp_x / sum_exp_x\n",
    "    return sigmoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_num = 20\n",
    "sample_num = 500\n",
    "class_num = 10\n",
    "eps=10e-3\n",
    "X = 2*np.round(np.random.rand(sample_num, feature_num))-1; # 1000 samples and 20 features\n",
    "beta_true = np.random.random([feature_num, class_num]);\n",
    "Y = sigmoid(np.dot(X,beta_true))\n",
    "lambda0 = 1; #Logistic loss\n",
    "lambda1 = 0.1; #Balancing loss\n",
    "lambda2 = 1; #L_2 norm of sample weight\n",
    "lambda3 = 0; #L_2 norm of beta\n",
    "lambda4 = 0.001; #L_1 norm of bata\n",
    "lambda5 = 1; #Normalization of sample weight\n",
    "MAXITER = 1000;\n",
    "ABSTOL = 1e-3;\n",
    "W_init = np.random.rand(sample_num, 1);\n",
    "beta_init = 0.5*np.ones([feature_num, class_num]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7faf533b2d68>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIcAABEkCAYAAABixr4wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3WmAXFd1r/1VVV1V3V09zxpsWbasti1bnm08YMZAgARCgkOAQIDADRcIvAETuIGQhEACgRBCnIFAHLj4JTFjGE0Ag2w825ItD5LLlizZktXz3F3dNfb9gFV1qvustRohw0ny/O6HW2evfXqXxF9b1s4++8SWl5cFCBP/RX8BRBfhgIpwQEU4oCIcUBEOqAgHVIQDKsIBFeGAquEX/QWO2rlzZ4OIbBSRw+eff37pF/19/js51t/byIRDRDaJyD4RefrOnTsP/6K/zH8zG0XkxyKyRUT2r/WmKIVj3ZP//49/od/iv7d18osOx+Dg4CtF5H0ikhSRT2Sz2b9fw21DIiL9//y30jA7U2184qr3y4aPfaB6ffkdQ+YPyZeK7kCt6Sa3T66YN+vlSmVV23987TPyay99Q/W6KZl2x/m1rjPdPg8WJ8z6/ZOPmfXe3m7553/+sMiTv8drddzDMTg4uEFEPiQi54tIXkRuHRwc/FE2m93j3FoWEWmYnZGG6an6Lxm4Hh4eM3/IUqngfsdcOuP2WSgumfVSpRzaPjQ8Wv2cSTa648yXZ90+k4Upsz40PmrWA8K/tOKp+NfKc0Xkh9lsdjKbzS6IyJdF5GVPwTh4ij0Vf62sl/rpa0hELlrrzU9c9f5VbY998G+qn2/8Gb7Yz8Mdt33jF/0VjpunIhxxEQluL4uJyOq/oBUbPvaBur9GHvvg38im9/1B9fqkGw+Z96/lr5X2p+ivlTtu+4ZcfMmLq9dr+WvlVb3nu33uKdh/bdw9vs+srxvok6997dPuOCs9FX+tHJbavzxERAZE5MhTMA6eYk/FzPEDEfnTwcHBXhFZEJHfEJH/9RSMg6fYcQ9HNpt9YnBw8L0i8iMRSYnIZ7LZ7J1rvf9X7p2VsZHp6vVXReT8O2rX7+691Lz/X+YfcMfY1rzB7ZOOJcz67lz4Ot2mtr7q59lCzh3nmuHb3T49TW1mvaPR/muybQ3/dA/zlKxzZLPZL4jIF56Kn42fH/4Pb1ARDqgIB1SEAyrCARXhgCpK+zlERGR0YUaG5ifr2o4Erq9NZM37P9dwqjvGC8bvPrYvF3BW56bQ9p5UbU3i8Vn7/4IsItLTbK9hiIiMLEyb9T9ad4VZb2r3xwjDzAEV4YCKcEBFOKAiHFARDqgIB1SEA6rILYJJLCaxWGxFU+16ZNFeELqn09+3qT1WELS+pcus/2pyo9u+U+y9nSIip7X4G48OJJJm/dqFvWa9t6VbfskdZTVmDqgIB1SEAyrCARXhgIpwQEU4oCIcUEVuEaw93Sz5xpa6ts7A9eTinHn/1oJ/5FW6wV5UEhFJxOw/NzdXJle1vWBF+wmtve44lTW80mTReTi8WLF/zanCsT3xxswBFeGAinBARTigIhxQEQ6oCAdUkVvneGH7oMwX1tW1vazzrOrn7zQ8bN7/0aR9oOtajS7OmPWhhfCzQW8ar228ObVtXWifoGc39Lt9JhrttZ2HpuzTwDP5RXeMMMwcUBEOqAgHVIQDKsIBFeGAinBARTigitwi2LdnHpLRqfHq9etF5EtT91Wvf9Rzinn/v5Rb3TFuL9sLaSIiy2Jvwgl7U5OISLFc23gzVvBftPOFin2MlYjIbNE+JruxIWXW17K5KQwzB1SEAyrCARXhgIpwQEU4oCIcUBEOqCK3CPbuzHYpttQvHn245YLq5+2P//BnHqO/pdPtM7Zg7wTb0hG+yyvYfmB2xB3nNwfOcvt4C3KfHrnDrK/lqbowzBxQEQ6oCAdUhAMqwgEV4YCKcEBFOKCK3CLYB2Z3yvB07cV5O+QP5O3Tt1av377BftndF2cecMdYebZ6mOf1bzfrd87sD20fX6ot4J24hmOf7i6Oun12TT5q1q/oOd2sd3R1uGOEYeaAinBARTigIhxQEQ6oCAdUhAOqyK1zNMQTkown6tqC1/84fJt5f2dTi1kXEeloaHb73D5tv58tEQ//cxVcQ+loyLjjfGlz+JNzQRfm7F/T7VOPmPX+tL/eEoaZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBPtI03YpZ+qfePuXzIXVz3+8HL7J5qh9M0PuGMlYwu3TmbYXnsrL4YtXmYbG6ufs7BPuOOfcbx/ZJCJyRftWs37D1B6zHl/D5qbQ+47pLvyPQDigIhxQEQ6oCAdUhAMqwgEV4YAqcotg78rtlpH52hNvN4jI6+bvrF5/puUi8/7XxPyjlkZz9pFOIiKb2+wX82lHOg0HXgz4nv7L3XGunrrb7XPXwkGz3p6yd5y1JpvcMcIwc0BFOKAiHFARDqgIB1SEAyrCAVXk1jkWC3lZKCzVtQWv/6hsvxPtzI4T3THOS9lrGCIiBbGfRBtZnA5tb06mq58fFvvdbCIihcA74TTTSxNmfWNrj/szjgUzB1SEAyrCARXhgIpwQEU4oCIcUBEOqCK3CPaFjnNluVC/GeeHXU+rfn5F/iHz/oOz/nHRf9qx3u3zivldZr23qT20Pfik3HVDd4b2CWps8J94847jftz5NRf906dCMXNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotgP04lZTFdWxh6iYh8P3A9PDkVclfN2d2b3TFel9vt9lksFcx6rpx322PiH7d0Rbf9fjYRkStinWb9o1N3mfW2tH/WexhmDqgIB1SEAyrCARXhgIpwQEU4oIrcOscBWZTZFU+KBZ8c68vY72W/d9x+77uISHoNG2wa4vYR2CPz4estwfbTu05wx/ncGfNun8t32Zt5ZvP2k3WZ/JJZ1zBzQEU4oCIcUBEOqAgHVIQDKsIBFeGAKnKLYOOVRZkqL9S1jQSuJxfnzPvXssCViPt/Jja19Jr1zanwo5ZeOHBu9fOOyb3uONvvaHT7vLf1PLP+xVT403dHdXbbm4U0zBxQEQ6oCAdUhAMqwgEV4YCKcEBFOKCK3CLYvVMHZHhirK7tjolHqp8v6j7VvD+/7J8lvnN8n9vnzMZ1Zv1waTa0faqyWP18Vvsmd5w7xuyz3EVEftBmv5NuuBB+DvtRpaK9q03DzAEV4YCKcEBFOKAiHFARDqgIB1SEA6rILYKd0jogHYv1u7lOa99Q/Xz35H7z/o5G/6Bv71FHEZFduUNmPRUP/62bLtV2rV3ZuMUdZ7TdXuASEVmWZbPeELN/PQ3HOAcwc0BFOKAiHFARDqgIB1SEAyrCAVXk1jn2zQ7J8HT9Zp+904ernze0dJv3l5cr7hjNybTb5zXNp5n1j4zfFtr+2Hztu78k3e+O8+ezI26fw/PjZt17T1x5ic0+OM4IB1SEAyrCARXhgIpwQEU4oCIcUEVuEezU9vXSla8/CumMztoZ4vMl+xzvmcKCWRcR6W+yz08XEfmbSfu9aXHl/W3B9n2FFnectbx/rVyxF/ae0WUv2LW3+7/eMMwcUBEOqAgHVIQDKsIBFeGAinBARTigitwiWCIWX/UEV/D68bmxlbfUeUX/Be4Y+0v+U2YHZ+2X7sVi4YtgxUq5+vm3p292xzmz0z8aavfEAbP+0rK9yJWstLljhGHmgIpwQEU4oCIcUBEOqAgHVIQDqsitczTFU9Icr3+CK3h9Ybd9Ws7nh+5wx7hn49lun19yTgha1xj+3rTgKUQHF+y1EhGRe8btk4pERFpSTWb9XXP2xqT+5h75nDvKaswcUBEOqAgHVIQDKsIBFeGAinBARTigitwi2F80xyXeUtvcMycinwpcn/PYIyF3/XQuH/Xfq5YvFc36xOJcaPsDU49XP3tHVImIXNl1jtvnP+cfNuvpeNKst6f9J+/CMHNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotgb56Zk4mp2hNp14jIlYHrp/fYRxzdOmkvGImIXNa51e3zyJJ9JnlPMvwpsgsCO9UenHk8tE9QojX8ybm6PjH7z/C+6SNmfaGp5I4RhpkDKsIBFeGAinBARTigIhxQEQ6oIrfOcXB+VIZn60/veWR2qPo50W7nuWsNG1u+O3yv2+fZ/WeZ9ZvHHwptv2fy0erntRxd/cJF/8/nD5x1jot7B816V1f403keZg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHst3rOkYXl2bq2N/ZfXP386RH7WKdi2d/Y0t/iLwrdOLbHrKcbwp8ya0jUns5rTdrHNYmIvH7JX5C7qv18s74/VjDrLQ3t7hhhmDmgIhxQEQ6oCAdUhAMqwgEV4YCKcEAVuUWwO/JDMrE0Wb1+jYjcuHSoev2uvkvt+ytT7hg7Zx51+zQ1pMx6vhx+LFRwEe7KjL1DS0Tk6nH/rPZ/XHjQrA8tTJr1gVivvNIdZTVmDqgIB1SEAyrCARXhgIpwQEU4oCIcUEVuEexwbkJG5usfhzw4X3uR3sXLJ5r3f8I5J1xE5Fd6trt9rp98wKzHK+F/ruKBRxfvKE+445zdcZLb585x+7x371iomPhHS4Vh5oCKcEBFOKAiHFARDqgIB1SEA6rIrXNUlitSXq7UtQWvX7e027z/pT3+O9OKUnH7tKXsI5uubT49tP3LbU+rfv61sVvcceIxfw1Ce7ruKO8pv2VZdscIw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2C5UslyZeKK9pq1ydkesz7D5fn3DFuGc+6fT7Z/XSz/s7i6p9xrYi8c3l/9fqyXvt9dCIiN4/udfu0puxFsDPb7Q1QXR2cfY7jjHBARTigIhxQEQ6oCAdUhAMqwgFV5BbB0g0Nq3Y+Ba+9XU0jxVmzLiLSlLSPdBIR+W7C/jmpUvhvXSpWa29Yw5+9tbwksCvdatbXN7SZ9faE/1LEMMwcUBEOqAgHVIQDKsIBFeGAinBAFbl1jpjEJL7iJJrgdUMssfKWOnunDpl1EZELe051++xceMysb27qC23vStbWFB7Lj7vjvLL7XLfPxLL9Drevjuw06wPxXnmLO8pqzBxQEQ6oCAdUhAMqwgEV4YCKcEBFOKCK3CJYf1OHJDLlurb1ma7q51w5b96fjPu/pD0z/kJZcMwwM6Wc235G03p3nC9O3ef2eVnnWWb9gu4tZr2rkyfecJwRDqgIB1SEAyrCARXhgIpwQEU4oIrcIlh7Q7OUG+oXujobarurys655W/u2+aO8aGpO9w+6bh91NL+uWG3fbZp0R3n5Ey/2+eaYfv7djRmzPpczj4bXcPMARXhgIpwQEU4oCIcUBEOqAgHVJFb53h0YVhG5uqfFMvOPVH9vK3tBPP+j8zc7Y5RrvjveMtOP+H2CVMIvG/tiXn/Xfb3/f2z3T7dbzlg1sdyM2a9YTHtjhGGmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwTLFfOyUFyqawteH1gcNe/PNDS6Y0wu+u+Be/HAeWb9+xMPhrY3JGrHUl3R5b/j7c8+4C+2Xeg80fZEftKs97Z2u2OEYeaAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBepvbRTL1T2gNZGrHFm1psp8Qu2fmoDtGbMXZ6mEeXBwy66e3bXTbvzfiH+l00xreN3dux2az/lcNp5v1hoT9DjgNMwdUhAMqwgEV4YCKcEBFOKAiHFARDqgitwg2nV+QiaX6nVrB6zdnzjTvv6WYdcf4RN8z3T5XF+yfs28+fJEs2P6b6y50x7lu6E63z+3jD5v125bt77ou1Sdfkze646zEzAEV4YCKcEBFOKAiHFARDqgIB1SRW+d4Yec2mSvVH2N0Zc+51c/vHbrRvL813eSO8d3YlNtnKGf3ScYToe2JWO3P2xbxv0tjg7/ZZ7B9g1mfLy2Z9Z42nnjDcUY4oCIcUBEOqAgHVIQDKsIBFeGAKnKLYN+feUhGJmtnn79FRL4+ef+a7w+ePa75zvA9bp9Uwn7H20IlfOFpemmh+vnvJu9yx/lc2yVunwvXj5j10/YcNuvFFv8JvzDMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcIpgnEf/Z89yS8ndovahnu1nfmTsU2n5y+0D1cybhv3jv1TO3un1aFu3z3G/ttb/rcle7+K88XI2ZAyrCARXhgIpwQEU4oCIcUBEOqCK3znFh2ykyk69/QuvpnYPVz3/dnjfvf97QmDtGOm5v5BER+e7kA2Y9Vwz/Ho/N1d5Bt5ZjtDubWtw+pUrZrF80Ym9eGoj1ydfcUVZj5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsFuGLtfhoZqC0l/LCLfHNpZvf72sJ3ntWwGKq7hqbj2xoxdT4fXg+0v7bCP4hYRuWb4drdPT7P9jrbJt51v1ouZdnnUHWU1Zg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhEslUhKesV54MHrTNJ+imypXHTHqCwvu33Obz/ZrO8YfTC0fXKx9j66v35vqzvOtX/q/08QPE89TM/f7TLrAwO98pVfdYdZhZkDKsIBFeGAinBARTigIhxQEQ6oIrfOsbm1X9oW659I29q+vvr5yOKkeX8qvoZ1g6T/Z2Lfon2k9PqWLrf97r+edceJxfyn4rx1jou7t5r1zs5Od4wwzBxQEQ6oCAdUhAMqwgEV4YCKcEBFOKCK3iJYsku6Uom6tlNTvdXPO55tb/Y5+Zv2+85ERM5oO8Hts2tyv1kvV8IPjD4yX1uke0nqXnecUtk+0klE5Im5CbM+sjBt1geSvWZdw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2CfXfsvrpjn94jIl8duqt6feO37SOQehvb3TGG8lNunzM67IWyZzaG139//eXVz3PL/vFSnx3yj31qTKbM+gu6zzLrbd3+70kYZg6oCAdUhAMqwgEV4YCKcEBFOKCK3DpHR2OLFJqX6tq6A8c7zxeWVt5SZzafc8fY3nmS2+eB6cfN+o0vW7128LCI/Nmv1J5y6/zn3e44Z3dvdvscnB816z+Y2mPW+9O98r/dUVZj5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsGWSgXJFfN1bcHrDufda9syG90xtKOpgy7tPc2sX/Kl1U+hfe519e1x57gmEZHfbjjJ7fPRBntz0vaWE816R1uHO0YYZg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHsus5BiRUHqtfLInJT4ImuV+WHzPt/OHK/O8aZXZvcPvdOHzDri6VCaPvD00eqn/sz/uLTtaWDbp/1jeHnrB+1a9b+rn1NPe4YYZg5oCIcUBEOqAgHVIQDKsIBFeGAinBAFblFsDctPiFjudpuqq+IyK/nHqtezxUWzftP7/LPNb+k0d8t9vDMEbO+ua3fbX98bswd54vpU90+z5u1z1BfXLFzbqX0UpM7RhhmDqgIB1SEAyrCARXhgIpwQEU4oIrcOsfU0ryML87WtQWvl5RNNkdNL827YxxM2scoiYi8rf9Ss/7pybtD28eWZqqfn9dnHzstIvKs0fCfE/T+vqeb9Z0yZ9bbezneGscZ4YCKcEBFOKAiHFARDqgIB1SEA6rILYLFYjGJxWKr2o5qiCfM+9vSze4Y65vsJ8hERPZUZs16OpF02/9+wD+H/dxpfyPOX03eYdbzpaJZH4j3ylvcUVZj5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsEyyUZpTdUvDAWve1NtK2+p84eJk90x3jZ3l9tnoMPePVWqlN325xywF9JERCYX7V1cIiLn92wx6wcWRsx6Jpl2xwjDzAEV4YCKcEBFOKAiHFARDqgIB1SRW+eYzedkemmhri14vbGx27z/DZM3u2M0r+Hf/YXl8HWMo7T3twXbPxk/xR3nNc0Lbp9Di+NmfVOmz6x3N/ubm8Iwc0BFOKAiHFARDqgIB1SEAyrCARXhgCpyi2AN8YQkVzzVFrzelOw077932X7fmYh/RLaISCpmP1l3SXv4Bpxg+1tyj7jjdKdb3T6/3GxvYLp66BazPpDyn7wLw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqCC6CxVcd7RS8/t7EA+b96Ybw45iClor2+ekiIrmK3acx4f/WPTE/4fbpaMy4fd69bcis/+0T9q61irOrTcPMARXhgIpwQEU4oCIcUBEOqAgHVJFb55gv5Fdtxglel5UTdY7SnkQL6mm2TwcSEblldK9ZXzz0o1Vt9z8xL1+882PV66YTnuWOs5aNR1vvPGzWW50jvTMp3mWP44xwQEU4oCIcUBEOqAgHVIQDKsIBVeQWweIxkfiKd7wFr1POkU0v6z3PHeP/H7LfmfaTMe0/N5kTn7Oq7fbbvl7X3rKGxaeV77MLs1SyNx5NfvRXzHoh3SL73VFWY+aAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBmpLpVWeTB697Gu1dXJ8/cps7xroW/yzwjzduN+uvnb41tD0ZeBLujb0XuuP868Qut8/F3VvNeve7vm3WBwZ65ctfebk7zkrMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcItgZrRtl3WJLXds57Zurn2fK9jnea3kcsjGRcvu8Y+k+s/7SvnPd9i/N7nHH+b2uC9w+/3fW/i4l5xHRcqXijhGGmQMqwgEV4YCKcEBFOKAiHFARDqgit86xe/qgDE+O1bXdOVl7V9rze84y778/9pg7xomNPW6f6dKCWb9pbt+qtnesaH9h62nuOLeVx9w+MwV7bWdDa7dZ723pcMcIw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2CXdq5VWZKfXVtz+nZVv28c8Fe5GpJNbpj3DW1egFrpcaE/a6401s3hrZvbqp998+N+MdLdTa2uH2+2mI/Ofe6woNmnc0+OO4IB1SEAyrCARXhgIpwQEU4oCIcUEVuEeymib0yPFrbHfV+EfnP0doTX96CzsntA+4YY0szbh9tkeuo2XL4i/uC7c/s3RbaJ+j7w7vdPv+74RGz/g8pe3dcIum/8DAMMwdUhAMqwgEV4YCKcEBFOKAiHFBFbp0j9uT/W9l21D/0PtO8/7EGf2PLR2ZucvvMl5fM+kPTh9328aZZd5wT2/rcPsO5KbP+p22Pm/XuWJd82B1lNWYOqAgHVIQDKsIBFeGAinBARTigIhxQRW4RrDGZkqZk/fHTwet/LB8w73946og7RiLu/5loTqTNekM84ba/ps1+T5yIyHULe90+aefpu5elTjLrjWz2wfFGOKAiHFARDqgIB1SEAyrCARXhgCpyi2AvbD9D5gv1T6Rd2XVO9fPYsr1D6zXdm826iMj/Gfux22fnxH6zXigXQ9sXi/nq53+avMsdZ7FYcPv0NNuLWP9ZHjHrnZWCPMMdZTVmDqgIB1SEAyrCARXhgIpwQEU4oIrcOsdXxu+R4ZHayT7/S+qPif6N/vPN+3Oy7I7hvftdRKSzyT52OlcM/3PVlKxtErqq8yJ3nE/M7HT7DM/bT7w1JVJmvah8Vw8zB1SEAyrCARXhgIpwQEU4oCIcUBEOqCK3CPa9TadIPNNdvZ4RkQdOOa16fcmhPeb9tyT9d7y9Zf3lbp+rj9gbglpSTaHtwSfe/mT4Rnec1nT4zwnyjuzePzNk1hebeMcbjjPCARXhgIpwQEU4oCIcUBEOqAgHVJFbBHvDxKyMj01Xrz8vIi8OXJ/Xbj/R1hn3F8G+l3vU7ZNxFtPisZjbHlP6BH22+QK3z8umbzbry8v27rflNeyOC8PMARXhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgo3lZ2Vksf7xv+HA9amN9gvzvnTkTneMT/c92+3ztuVbzPq7O58W2h58BPJvZ3a547w2d7fbp7up1ayP5+wXDcaPcQ5g5oCKcEBFOKAiHFARDqgIB1SEA6rIrXOc3rpRBhYzdW1nt59U/XzD+IPm/Qnl3WtBu5PhR1MHtaWbzfoPK2Or2p6/on1i0V5/ENHfFRf0kbaLzfq7Sreb9aaU/b46DTMHVIQDKsIBFeGAinBARTigIhxQEQ6oIrcIdt/0QRmerF9gumtyX/Xzs7rPMO//SNp/uuucQ7e6fc7q2mTWbxwNP34q2H5e9ynuODvH97l93jy2w6yf1rnRrPc0dbpjhGHmgIpwQEU4oCIcUBEOqAgHVIQDKsIBVeQWwe7+i5dIqpirXu8VkdFPvbJ6vfmtXzXv376G3VedjfbL/URE7p98zKxXlsPPEw+2758fdsfZ0rHe7ZOv2DvX9k4dMuvr0nl3jDDMHFARDqgIB1SEAyrCARXhgIpwQBW5dY5L//jbMjI8Xr3+3o5flxPe/KXq9a92nWXe39rp/5Iuy/t93rxsnxDU2JAKbe9vqW2smV5acMcZW5px+8zmc2Y9HrP/jMecuvpzj+ku/I9AOKAiHFARDqgIB1SEAyrCARXhgCpyi2Bntp4gGxbrj3MOvtftK+P3mPcXyyV3jP+bSLp9vPezLZUKbntjwxrGEf89cKc6G4Juucw+oqrY3imH3VFWY+aAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBDuUnZWxpoq7twFLtGKhU3P7KHemMWRcR6U7Z70wTETkwN2LWtd1VwUWtZ3ae5o5z9/xBt0+xUjbr/3bXCWY93dcm57qjrMbMARXhgIpwQEU4oCIcUBEOqAgHVJFb5xjJTcvwwmRd25HA9d90X27e/y/lx90xchX/pJstbevM+tmp/tD2lwSeyPu3Ef899amE/z+Bt87xgeJOs97X3COfdUdZjZkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEe3H3dpmr1B+F9Iq+C6qf3zL0o595jFeuu9jt8/2Zh8z6yNL0qra3iMj3ZrPV66Zk+NFQQecEnubT7M/Zx2Q/q/VUs97a1u6OEYaZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBPvq2L0yPDxavX6ziHx++I7q9fk9W8z7Jwpz7hj35f13r3k6UuFP1gXb/zl1pvtzXjF1l9unt8lexPrm5P1mvT/VI691R1kMjeZvAAAgAElEQVSNmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwTrbMxIsan+WKbuwPXO8X0/l+8xkOk069mpJ9z2vQP2Di0RkUTO//P5rw0nm/X3pqfMemeL/WvRMHNARTigIhxQEQ6oCAdUhAMqwgFV5NY5JhZnZXSh/mmy4HVcOVb6qPQa3qu2WPSPfbq09RSzfqhpJrT9gsBmpL9bsDfhiIjM5RfdPh/qWDDrHys7x3Uv+0d+h2HmgIpwQEU4oCIcUBEOqAgHVIQDKsIBVeQWwU5s65PmxURd2ykdtXPID89PrLylTrlSccd4/fpL3T6fCzxlF2b8d85Y1faIiPznr6Sr153/Oraqz0obW3vcPjvGHjTrl5RLZn1dQ5983R1lNWYOqAgHVIQDKsIBFeGAinBARTigIhxQRW4RbLaQk+lC/c6n4HXRWfDpaPR3Pe0ujLp91jlPvMVaG932RDwR2uenFZOYWW9OhX+Xo5qSabOuYeaAinBARTigIhxQEQ6oCAdUhAOqyK1zLJUKklvxRFrwelmWzfunl+ynw0RE7s0fcPuc2Npr1vv+/t5VbTe/sr795QMXuuPcuuB/lxZnHeOhZw2Y9WJ7p4SfQ2Rj5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsF+sPUESXS3VK9HRWTvts3V6zP2HDTvz5eK7hjbu05y+2Rn7WWjhkT4Rp5g+zfGd7vjrOWYqjPbTjTrl95qH2/d098on3RHWY2ZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBKuUYhIrxVa1HbVULNj3L/vHPn2gssHt82ulg2Zd252VDDzltuh8VxGRfNlftLuv8phZXygu2WNk7N1zGmYOqAgHVIQDKsIBFeGAinBARTigitw6x3P2HZTh4dqx0DeLyGl7a++vj8XsU272nbbNHeP0h291+8Sdcc5r3+y2/2jkAXecq9Y/w+3z6cm7zXpjQ8qsp526hpkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsE29TaK5mV73hrrx1rlJ2yn0Tbstd+H5qIv5C2FtoCV7A9HvP/7H3w7g+6fT61+flmfS6fM+sthUV3jDDMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcIlhhuSyFSv173ILXLxo417x/98Ihd4xzMye4fb41svps86CeTFtoe39L7d1wazmH/Rtn/bHbx9uV9lRh5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsE6GzKynKw/Lqk72Vr9fO/C4+b9qbj/S7phYo/b5wX9Z5v1m6ceDm0Pnr3e3dQa2ifoA8v+CwBf3X2eWf/63ENmvbelyx0jDDMHVIQDKsIBFeGAinBARTigIhxQRW6dY7FckFw5X9cWvD48N27en0r470zb2r7e7TNVtp8Sa06m3fbpvL/Z51fbTnf7eAd2x8TeDOTVNcwcUBEOqAgHVIQDKsIBFeGAinBARTigitwi2EeT3SLJ+q91TbKv+vk5qbGVt9RpTzW7YzTF/bPADy1NmHXt3WzB9pNb+t1xvjC+y+3TmW4x68MLU/YPyB3b/8zMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcIthLxx+QodHR6vXtIvK00Z3V6w2t3eb9c2s453tjS7vbJxlLmPUDi6Oh7enATrSrYxvccX6pZJ/lLiIyE7N3lLWmm8x6JtXojhGGmQMqwgEV4YCKcEBFOKAiHFARDqgit85xSe9WmV7uq2t7dv+Z1c9HCtPm/blS3qyLiOzK2acDiYj0pzrM+nhu1m2/qtHemCQi0tMUfkx2UFvS3sB0bpO9ntLW5q/rhGHmgIpwQEU4oCIcUBEOqAgHVIQDKsIBVeQWwRolKU2x+qObgtezxZx5/1req1aqlN0+44vhi1xHLcuy2/5LSX+zT1/Cf9/cRxd2m/Vvjd9n1vsbeuUN7iirMXNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotgPxrfI8MjtafJ/lhEvjNyb/U6HrPP8T69Y6M7xoNT/k6w9nTGrGtnn3c21o5omo95p5aLfGLiNrePNtZRd2/aYtYrPR0y746yGjMHVIQDKsIBFeGAinBARTigIhxQRW6doy3VKIvp+ie82gPXmaR9Ss0js0PuGJva+tw+1yROMet/ngzfDLStpbbOcv3CPnecDmc9RUSkomwsOupphw6Y9f5ij3zeHWU1Zg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhFsfXO3pDP1G3pOzPRWPx/O2e9eO6PDf4LsvsmDbp+3dNjvgXt47Eho+01je6ufP9b3DHecD8/f4/ZpbbCPr96WsTc4dbbaR1hpmDmgIhxQEQ6oCAdUhAMqwgEV4YCKcEAVuUWwUqW86lim4PWLO7aZ939p4l6zLiLSELff3yYisn9m2KxrT6EF29899mN3nG+2X+z2eXX+AbM+2WDvjouV7SfmNGsOx+DgYJuI3Coiv5LNZg8ODg4+V0Q+LiJNInJdNpt935P9zhGRz4hIm4jcJCJvymazpWP6dviFWtNfK4ODgxeLyM0isvXJ6yYRuUZEXiIip4vIhYODgy94svu1IvLWbDa7VURiIvLG4/2l8fOx1v/meKOIvEVEjv4fFC4SkUey2eyBJ2eFa0XkysHBwU0i0pTNZm9/st9nReTK4/h98XO0pr9WstnsG0REBgcHjzatF5HgNu8hEdlotOO/oGP9D9K4SN1++ZiIVIz2NfvoZz+wqu2a71y95vvf9NMM9hS44cbrjvvP/NZx/4lrc6zhOCwi6wLXA/KTv3K09jV712vfLxOjk9Xra75ztbz+hW+tXl/UaB/8upZ/rRTK/n8fLy/bz4qkG5Kr2m648Tp5zjNeXr3Ol4ruOMfjXyvrm7rMendfl/zlNX/ijrPSsa5z3CEig4ODg1sGBwcTIvJKEbk+m80+JiJLg4ODlz3Z79Uicv0xjoFfsGOaObLZ7NLg4OBrReQrItIoIt8RkS8/WX6ViHz6yX/67hKRT/40P/vzpyaloae20eawiHzvjNr1ptvuNO/vamp1xzi3Y7Pb5+bRvWZ9sCV8BjuppfY0XXbaf0/9O+WQ2ydftmeg3RP2E2/rkv6R32F+qnBks9mTAp9vEJGzQ/rslp/8awb/xbF8DhXhgIpwQEU4oCIcUBEOqAgHVJHb7POiPVMyOjJevf6miJx7b+264ixr33tBtzvGmXf6C0/aO9yO2qccLxVsLziLVyIib07ax0uJiLy9aC+mxWP2n/GYU1d/7jHdhf8RCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhFsU1OvtDXXf60tzQPVzyc09Zj3X7xrxh3jthNPdPuc9+h+s15c8VTeUcFFunUt9t5OEZF3T9/u9mlssI+g8vbENsRZBMNxRjigIhxQEQ6oCAdUhAMqwgEV4YAqcotgv1PpkVKlftHnrZX+6udXTtxo3t+kHMcU9NsT/lng3iOIv9y36mE/ERF5Xu9Z1c/fHNrpjvPs/rPcPrdPPWLWT++wT7nobvMX48Iwc0BFOKAiHFARDqgIB1SEAyrCAVXk1jm+05CTmcAxRe8RkS8FrlOJ1Qe1BSXW8HTXvZP2MUkiIjGJmfUHcuFPoQXb13IE1Y/H7eOlRER+1HmBWf9o0l6TaW9od8cIw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2C7Zh6SIbHx6rX7xGR68fvq16vy3Sa988Wcu4YC5Ult8/Ii+zjmC768Zz7M7y3GYiIPLj0uNvnjRX72KdDE+NmvT/ZK281e4Rj5oCKcEBFOKAiHFARDqgIB1SEAyrCAVXkFsFObh2Q9sX6J94G22sv23tW2n666xNHfuyOcWnvaW6f591eMOuPz4257S/pP88dZ7Jl3u1TqNjHOiWcY528uoaZAyrCARXhgIpwQEU4oCIcUBEOqCK3znFFcr3kUi11bc9P1Y6j/qeJu8z7G+IJd4xHcyNun/5G+/SfZuUEoWD710d2ueN0Nra4fT6TPMOsv3zZPkHIey+ehpkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEG5eizEn9RpuRwHVPY5t5/5md690xnr3sH4P09/msWV8qhW8GCravZUFuoeg/ffeykr3wl3TGicfsI6zU+47pLvyPQDigIhxQEQ6oCAdUhAMqwgEV4YAqcotge4pjMlGYrGu7tzBc/fyGJvtptcU1rPd8aNbfoVVerpj1UqXstp/XfZI7znRpwe1zYMbeueYto7Xk/YW2MMwcUBEOqAgHVIQDKsIBFeGAinBAFbl1jkO5cRmZrz8158B87d/5OzrszT7TpUV3jPHcrNunozFj1uPKu+SC7aOFGXecdWn7uG4RkYvXbTLrSed9dK29vOMNxxnhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgjUlUtLcUH+kUvB6x/ge8/61PN11Y9dFbp/3JOwFrEPJidD2DS3d1c+5Ut4dZzTmL5Q9MGO/B+6M9hPMelcpfGOSh5kDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEmysuykwhV9cWvP63tkvM+z+V8heVPuQ+IybSJEmzPlMIf1It2F5UnooLWsuxT4Wy/Y63+6ceM+sD6ZxZ1zBzQEU4oCIcUBEOqAgHVIQDKsIBFeGAKnKLYCdl+qS1tf5rndq6rvr5VbO3m/ef1r7RHWP3xAG3T0+z/djlYjH87PNg++c6L3PH+e3JH7t93j9whVn/8NitZp2zz3HcEQ6oCAdUhAMqwgEV4YCKcEAVuXWO4aVpGV2sP976icB1R9o+juk5KX+d47HGUbfP1NK8WU8nwjcDpRK139Jb0v5mn3O6Nrt9Pjpur+1c3m0f+d3R2eGOEYaZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBJsvLMpsvv4JreB1usF+Eu3RZf/prpm838c7+7wt2Rza3ttUO2c8W/bPWD+0OO726WxsMeuvKneZ9WTZ3rikYeaAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBipWSFMrFurbg9bIsm/f/zcCcO8b14/4vu1ypmPXScvgur2C7911FREYXpt0+7c6C3Jcz9lFX7YmYvN0dZTVmDqgIB1SEAyrCARXhgIpwQEU4oIrcOsf53Vtkqtxd13ZZ3+nVz0fyU+b92x60j3oWEWlNNbl9NjR1m/X7Jg+Gth+ZDzydl7TXJ0REzu72n3h7bdLu86GZnWa9r7HHHSMMMwdUhAMqwgEV4YCKcEBFOKAiHFARDqgitwh239RBGZ4Yq2u7a2Jf9fOv959n3v9oyl4kExG5d9I/3roxkTLr2hNxwfbt6QF3nMNl+3gpEZGrRm8068/tO8usd7Ry7BOOM8IBFeGAinBARTigIhxQEQ6oCAdUkVsEu7R7q0yX++rant27rfr5hzNZ8/7Pp890x3hd84TbpzdlH5W0Lt0Z2r4lU3sf3b8N3emOk0k1un2aGuwFuR+M3m/WB+J98k53lNWYOaAiHFARDqgIB1SEAyrCARXhgCpy6xwPLjwho/P1Rz7fM197iu281pPM+6+N++9VG5qadPuckdlg1vcsPBHafmip9t03tfWF9gnKJNJuH2/j0SNzR+wxUv4YYZg5oCIcUBEOqAgHVIQDKsIBFeGAinBAFblFsFK5JMVyqa4teH3D2APm/Q3xhDtGWzr8/WxBF8bt96bdUdgX2j5XWKx+LjX4C3KXZk5y+8wsF8x6W/vJZr2TJ95wvBEOqAgHVIQDKsIBFeGAinBARTigitwi2Gx+UaaXFuragtdP7zvDvP+BucfdMebyi26f7+TtM9RzxbzbfvDFJ7jj9H3lLrfP8rL9rrhTO9ab9Xwh5o4RhpkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEO61zo/QW63dqbe8+qfr59omHzftLFX/31Vo8Mjdk1puT4Y8YBtt7v/yIO85dA/ZZ7iIi95dazfrrx+2z0ecbS2Zdw8wBFeGAinBARTigIhxQEQ6oCAdUkVvniEtMEisyG7w+pc1+b1p5ueKO0Zu0j64WEXl0ccSsLxSXQtvjsdrGmhcNnOuO87bygtvnUN7ewLS+xX46rzfT7o4RhpkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEy1eKslipP+YoeJ2OJ837s7PhZ5IHjcSn3T7rmsPf4XaU9sRbTGqLYN8d3e2O86K+c9w+l7XYZ6jvKNobk7rSHPuE44xwQEU4oCIcUBEOqAgHVIQDKsIBVeQWwcrLlVW7uYLXs6WceX9rqskdoztlP0EmInJwfvSYxkk31BbpTmmxd62JiLyo7O9Ke8/83Wa9M9Vi1uMlXgCI44xwQEU4oCIcUBEOqAgHVIQDqsitcxxemJDhubG6toNztTWHdRl7E87TWu33nYmI7M4ddvuc2rrOrJ+WDt+A88y2rdXP+4uT7jjvnvOPt76o/RSz/v3R+836YrN9PLaGmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwT7ZMt5Um6brWv7t7ZLqp//JGa/e+27Y/aCkIhIIu7/mZiKz5v1+ydXf493isiXhmqLWsmE/9v7/t7L3D5Xz9lPznHsE37uCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHszbN3y/B0befXbXKVXDl9c/V6a8d68/6WVKM7xvSSf9649ztz5boL3fYvDt3pDvOGsw65ff7xNvvP8GhuxqzHF1PuGKH3HdNd+B+BcEBFOKAiHFARDqgIB1SEA6rIrXPEY/XvSftJW+16KDdl3t+RzrhjzBfC388WtLXNXk8JW8N454r2m7oudscZ+P4dbp/gkdlh/rr/mWY93eWfHhSGmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwR7ee+5Mi/1T7z97sDTqp8/PXSbeX+hUnLHyCT9454PzI+Y9VQi/F1zwfbnzuxyx+lq8o/a3pzpN+vPTNjvrFuOi/i/K6sxc0BFOKAiHFARDqgIB1SEAyrCARXhgCpyi2AHK/MyXa5fBNsXuD6xtde8v6XBf+LtBelNbp/Pz9rHR3Wmw9+rFjybPVfKu+OM52bdPhe22ee5XzH+kFnvb+iVL7ijrMbMARXhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgu3PjcjYwkRd20MLR6qf/6h5u3n/D+Jz7hj/PLnT7bPykcyVRhbCd189Plt7eeEv95/tjnN7ZZ/b53sj95n17d0nmfXujH02uoaZAyrCARXhgIpwQEU4oCIcUBEOqCK3zpFOpKQxUX8cc/B6T0PRvH/HpL3xRUSkVCm7fRriCbO+LMtu+4kJ/wiqIxl785KISFfafiquIWZ/10Ts2OYAZg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHs0Ny4DM+O1rUdmK0dwXTtkv2E2N2b/afZXjhsH5MkInJ4xYajlTa29rjtuwqjoX1+Wo/ODJv1tnSzWe/P2AuHGmYOqAgHVIQDKsIBFeGAinBARTigIhxQRW4RrCWVltZUU11b8NrbxbVl74PuGO2N/g6tVNz+rTk8N+627/3qH7jjtP3S+9w+vZl2s54r2sdLlZfDd615mDmgIhxQEQ6oCAdUhAMqwgEV4YAqcusc8/klmc3n6tpWXltevf4St88Xhu90+3gn+2SS4cdoB9vf+zvfdceJOeOIiORL9madH7Sfaf+AVnudRMPMARXhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgvVl2iXWUqprW99SO555qWwvCH1l/B53jMt7T3f7LFXsce6ZfDS0vVCuffdLl/w/ew/1neX2uX36EbN+ydjdZn1dok++5o6yGjMHVIQDKsIBFeGAinBARTigIhxQEQ6oIrcI9vyWU2Vhsb+u7aVtZ1Q/P15ZMO9/tGAf1yQisjnR5vb5/NgdZl07Gz24s+t1c7e746xr7nT7nOicj/5Q8bBZj8f93Wah9x3TXfgfgXBARTigIhxQEQ6oCAdUhAOqyK1zfGFslwwP194H/2oR+dTwbdVr7d1qR1XWcIrNq/tPdvt0OKf/vK/tgtD2D/dcXv38zpEd7jhHFibdPt4Tb+9ad4VZb+7x13XCMHNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotg5Upl1RHWwetkwv7K5UrJrIuIvGd4h9vnteueZtaHpRLeHq+1e8dSi4hsaR5w+zTHU2Z9R3HIrHeVluQF7iirMXNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotgl/Rulenlvrq2Z/fXzvZ+R8FeWHpz+WF3jLXsFvv+vH3U0vDC1Kq235B3ysdHbq5eb+88yR3n0dyI22exVDDr3tnw6xr6zLqGmQMqwgEV4YCKcEBFOKAiHFARDqgit85x78xjMjI1Vtd2x9S+6ueXO/f3N3W4Y2Qa0m6f/XPDZv35fdvd9nWxJnecvbP2qTwiIomY/Wf499dfbtYzvTzxhuOMcEBFOKAiHFARDqgIB1SEAyrCAVXkFsEaE0lpaqh/wit4/dvtZ5v3fyNnb9IRETmUG3f7FFc8dbfSzrmDbrv3M0REzu3Y7PZ5ZkO/Wf9b5yju/kSv/Jo7ymrMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcIthsPifTS/XvcQtef60ha97/2OyoO8aHBp7h9vmLSXth6Yq2U932LeLvBPuLoRvdPo9m7KfiUs5RWCnlfXQeZg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHMM1u0jzj6lx5/gettzgKXiEhz0n5k8kBx9bFPK9u3JP1FsLXIlfJmvT1lv6ywOdl4TOMyc0BFOKAiHFARDqgIB1SEAyrCAVXk1jm8d7wVyvY73B5K+kdXzxcW3T7eBpqdE/vd9seax0L7BL133Ro2Hg3dZP+MzovNemMTxz7hOCMcUBEOqAgHVIQDKsIBFeGAinBAFblFsHQyKY3J+mOfgtcLhSXz/rblmDvGmV2b3D4PTj5u1r/WGX7e+Fc7Lqt+ft3iPe44Hwm8E07T4Dyx9oXiQbPeXeqSp7mjrMbMARXhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgr2q8xxZKM3Wtb2pp7bT6e9GbjHv/0bpiDvGfMleSBMR2dZ1oll/U/6BVW1fF5E3FWrtc2vYcXZax0a3T3m5YtZb4vYTbc1x/4WHYZg5oCIcUBEOqAgHVIQDKsIBFeGAKnLrHAeWczKzPF/X9siKa8sdY/bx1yIip3ed4PbZM3XIrC8vhz9ZN7IwXesj/tN3t77vXLdP/7uvN+vZ0hNmfaChzx0jDDMHVIQDKsIBFeGAinBARTigIhxQEQ6oIrcItnvhcRmdG69ru3Pu0ernC7q3mPc/Mj/kjpGOJ90+3+uwnxG7cvG+0PaOxtpR0xVloSyo7R3/4fZJxu3/mS7rPc2sd3Z3umOEYeaAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBZvM5mV5aqGsLXt+Ze8S8f/rf3+KOse0N/+72eUX5QbO+VC667YmY/2cvvoY+xYp93nt2wX7Kry9XcMcIw8wBFeGAinBARTigIhxQEQ6oCAdUkVvnSDckpbGhfjNO8Pqyjm3m/Zt+51/dMT7RcoHb51OxYbPeHE+Ftl/WNVj9vGPMXisR8Y+uFhHpbuow62O5GbOeWuJd9jjOCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHspR3bZKFYf7z1q7vPq34eqthHU2eS/oLPd5I5t8+usUfN+oZMd2j7gaWx6ueT2wbccYYXp9w+C0X715xO2E/wpZy6hpkDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsE+2HugIzNT1SvXyMi35p/uHo9nV8Iuaumr6ndHeNVS/47z75asd+rtqC8Jy7Y3t9o7+ASETkx0+v28Z6cu2/ioFlfKvLEG44zwgEV4YCKcEBFOKAiHFARDqgIB1SRWwSbyi/IxNJcXVvwulQpm/c/Mm0fgSQi8hvxUbfP9q6TzPqD04+HtgcX6bwzy0VELm852e3z3Sn7scoHNp9p1is9HeLvfVuNmQMqwgEV4YCKcEBFOKAiHFARDqgit87xjI6tMlvsr2v75e7av+N35Q6Z9w/l/CfI4hJz+zw6bx/7lC+FH28dbD8vc6I7zt2L9q9HRGRDc/jTdUf9Q97e4JQptMpvuKOsxswBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2CtUiDLK/4Wq2Ba2+R64LOU9wxXhDvc/u8d/hGsz7YucFt/8/x+91xLuo61e2za9o+guqh6cNmfUB6WQTD8UU4oCIcUBEOqAgHVIQDKsIBFeGAKnKLYN+ffVhGp8er128SkW9M1574SjovzHtOwl/g+of5B9w+Z3VtMusHF8KfmjuSm6x+fm6P/SSaiMi3hne5fc50vktLwj7vvauz0x0jDDMHVIQDKsIBFeGAinBARTigIhxQRW6dY1vrRtmw2FLXdm7bSdXPu2YPmPd/ZPw2d4w/7LnE7fPJqbvN+nM7z3Dbd8xk3XFia3j67uSU/cTbd8fsTUUDKf8I7TDMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcItih/KSML03UtR1YGqt+nssvmve/t/9yd4yPT97p9llSjnU66huj96xqe/eK9s7GllV9Vko3JN0+t8zsM+s7Os+xf0C7/967MMwcUBEOqAgHVIQDKsIBFeGAinBARTigitwiWEuiUfINTXVtrYHrziZ7YemDIz92x/jXzqe7fV43Zf+cK/svcNu/cOR2d5wr11/k9rlx5mGzftm4vWttXUOffM0dZTVmDqgIB1SEAyrCARXhgIpwQEU4oIrcOsdkcU7GC7N1bcHrDzfZG1t+d8E+llpE5EvJWbfP9s6TzPoNs6ufZvv/VrTHYv7TbEOlebeP93Makymznk76G4rCMHNARTigIhxQEQ6oCAdUhAMqwgEV4YAqcotgIwszMjw/Wdd2JHD9ptwt5v2z17/fHQMQDR4AACAASURBVGPgxR9x+yyWCmZ9c1t/aHtwY9J0YsEd55axvW4f72io9S1dZr091eyOEYaZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBHtmz+kys7yuru0F/WdXP986/Yh5f+sL/swdI5O034kmItKcTJv1JxYm3PZf7j3LHec7o7vdPpf0DJr1t5Tss80bkm3uGGGYOaAiHFARDqgIB1SEAyrCARXhgIpwQBW5RbDHCxMynq/fCbY/Xzv7vLK8bN4/kOl0x5gvLrl9mhrsRwzf3nZuaPv7ei6tfR7e4Y7zooHwnxM0XJoz638dP2LWu+KL8gF3lNWYOaAiHFARDqgIB1SEAyrCARXhgCpy6xxxiUtiRWaD1zN5+ymy6SX/GKUT2uzNMSIi44v20VDj8bLbnkr4xy1dP+Jv9tm//VSzftbex816f3PeHSMMMwdUhAMqwgEV4YCKcEBFOKAiHFARDqgitwj2xOKkjCyM1bU9FrhuT2fM+3NFf8Hnha32E2QiIt9cto9jyknFbS+WS+44X17D++a23n+rWS+Ui2a9pbDojhGGmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwT79/YtEivUdmqVRWRH1xnV60ucY5L6Mx3uGPcWxtw+S87Z59+b37eq7XdXtHc1tbrjvKvysNun2/k53hN8mVSTWdcwc0BFOKAiHFARDqgIB1SEAyrCAVXk1jlO/qtnSapSW2PYIyKD//LC6vWGVxw27x/OTbljjOZm3D6lSvgTbUdNKU/WHZwdqX5uS/vvVptYsk/tEREpVOxNQ+VK+Majo9ay6SgMMwdUhAMqwgEV4YCKcEBFOKAiHFARDqgitwj2tDd9RkaGa5txvr/j32XDr/919broLE4lYn7eSxX7CTERkVTC/q3RFska4onq579sv9gd5zoZdfvcOv6QWf9y+6VmPdHCO95wnBEOqAgHVIQDKsIBFeGAinBARTigitwiWHuqWQrplrq2zsD1delN5v2/Wxox6yIixWV7IU1EJBlLmPWD8+GLV+nAeed/vnCvO05w0UyTSTaa9ZdM3mTW16X75Ovye+44KzFzQEU4oCIcUBEOqAgHVIQDKsIBVeTWOU5Md0trY/2//U9u7Kt+vnRsp3l/2dkMJCIy/90/dfu0v/DPzXpPc/gGmkyqtiYxND/pjpN0NhWJiFSWl90+TwVmDqgIB1SEAyrCARXhgIpwQEU4oCIcUEVuESy/XJb8cmlFW+3a2xwTj8XcMc7+rU+5fdqdI5ve1nqO2/6+hR3uOCe29rp9tjT2m3XvV9zR6x/5HYaZAyrCARXhgIpwQEU4oCIcUBEOqAgHVJFbBHt0YURG58fr2rLzR6qfn95zmnn//kX/GKXH5/x3vL2j/3Kz/sGxW1e1/ZK8Qz44Xmv/k3XPdMf5wNCNbp+54qJZb0ykzHrvYt4dIwwzB1SEAyrCARXhgIpwQEU4oCIcUBEOqCK3CNbYkJSmhvpFneD1j0Yf/JnHeP3A09w+Hx2yj1Ia7NwQ2n5Sa+3RzWvm/e/6wOYz3T6XDx0w64dz42a9nDm2OYCZAyrCARXhgIpwQEU4oCIcUBEOqCK3zjGdX5CJpbm6tuD1bw1caN7/WGnWHeNzI3e6fd617gqz/ump8OOnRhanq58v69jqjrPtwC63zynt68x6zHnKr63RfnpPw8wBFeGAinBARTigIhxQEQ6oCAdUhAOqyC2CdaQzUmxcqmvrbmytfv7h7MPm/R9rPNsdY6J9zu3zNyO3mPUtysLUQFNn9fMPxh9wx/n1gQvcPjdM7TXrKef89ESMzT44zggHVIQDKsIBFeGAinBARTigIhxQRW4RrDvZKrFU/dnnfan26ud7J+2nv141/yN3jLnPvt7t0/+/vmDWj+TCX+4XbC+US6F9gr42bL/QUMTf6XVSW59ZzyQbzbqGmQMqwgEV4YCKcEBFOKAiHFARDqgit86xb25Ihmfqj5/eM3Oo+tnb2HJhr338tYjICW/6otvn1b32JpyNlWRo+zs7a0/k/d3sPe44rckmt8/h+Qmzno6Hf5dqPXZs/zMzc0BFOKAiHFARDqgIB1SEAyrCARXhgCpyi2BtqWZZSmfq2joC197Gl90zB90xKsvLbp9/Hb7DrCfiq/9cPV/eIX8ReMdbvlR0x9n/Ty9y+5x+1X+a9YemDpv1demCO0YYZg6oCAdUhAMqwgEV4YCKcEBFOKAiHFBFbhHs9zJnSL6l/vzyP2w5p/r5T6ftxamZpQV3jDO6TnT7HJgbMetLpfCFpWD76Z0nuONc8ydDbp+GeMKst6bt3WSZNE+84TgjHFARDqgIB1SEAyrCARXhgCpy6xwfGr9DhkdHq9e3yR/IVaM3rvn+wc6Nbp+D86Nun5f1nmfWty2Hry385cCzqp8/OnO3O84fTO1w+7xi3UVm/e64vdmnJ93mjhGGmQMqwgEV4YCKcEBFOKAiHFARDqgIB1SRWwQ7o2OjDBTqn3g7p2tz9XNZKub9e6YOmXURkc1t/W6fb0/Z72e7PuS9ac8WkY/N1I6rnlz03yV36KJT3T6n3bPbrHtPAZaX7M1CGmYOqAgHVIQDKsIBFeGAinBARTigIhxQRW4RLDvzhAxP1Z99fv/UY9XP3U2t5v1b2te5Yzwyc8Tts63Tfirueanwp9le31F7Ou/21nF3nHPut3dxiYh0NrbYPyNjf9f29g53jDDMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcIlhjMiVNyVRdW/C6Iva55d454CIif9X/TLfPtaWDZn1XZWpV28tWtN81uc8d5yW957h97ll8wqx/Z+Resz4Q65O3uaOsxswBFeGAinBARTigIhxQEQ6oCAdUkVvneGPb2bKUrz/e+qrO2rFH/zh/n3l/b6bdHeNLZXvdQERkaGn1OkbQg4uPh7bfOLan+nljS487zncnH3T7bG/fZNa7euzNQF1dne4YYZg5oCIcUBEOqAgHVIQDKsIBFeGAinBAFblFsGvm7pfRmdqTYs8XkU8EjlLKFfPm/amE/0vaM+sfDbW1db1ZXygshbanE8nq5z9Kn+GOc0PGfyfdLfOPmvWupL0I1lQOfx+dh5kDKsIBFeGAinBARTigIhxQEQ6oCAdUkVsEm1laWHVmePA6JvY5371N/k6wszPhRzYFrY+Hv+DvqLe3h/+Mf2h/WvXzexbtM8tFRI7MTbh9mpJps54IOYc9KFW2fy0aZg6oCAdUhAMqwgEV4YCKcEBFOKCK3DpHc6pRWlL1/y4PXi8v2yf7jOSm3TE2tG51+9xTGDPruWTXqrbTROT7ydomoJZSozvOmV3202wiIo8t2N+lO2kf+d3RkDHrGmYOqAgHVIQDKsIBFeGAinBARTigIhxQRW4RLB6LSTwWW9V21Ekt/eb9D0yFH8cU9BdX+Btstn/L7nMkP7mq7fdFZMd87UjrucKiO06zs5FHRH+67qgDCyNmPZcru2OEYeaAinBARTigIhxQEQ6oCAdUhAMqwgFV5BbBEhJb9QRX8Prc9Drz/j0x/0inLf/hvwfu1Bb72Kf9C8Oh7Yul2hFLF3dscce5b95ftDu5fcCsn9u80ay3tftPAYZh5oCKcEBFOKAiHFARDqgIB1SEA6rIrXPkKyVZKhfr2oLXOxb2m/f/3sAl7hhXP/Fjt89sU86sFyslt/2Wyaw7Tn7FrzVMLmkf6f1HGXsdZDlWlvBva2PmgIpwQEU4oCIcUBEOqAgHVIQDKsIBVeQWwZaKBVlc8R634PXJnX3m/f825R8p/dr1l7p9Hi/PmfWJRHg9+I63U9rsxSkRkdvH/IWy2by9IPdpsY91apEmeZk7ymrMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcItj2rpNkY6n+Ca0Le06tfp6t2Ecgvbv9AneMq+cfcPtoO71+GsN5/xz29rR/Lvm65k6z/sWp+8x6X6qHRTAcX4QDKsIBFeGAinBARTigIhxQRW6dYzg/LaNL9UdLHw5cD82vPlY6aJfYT8SJiPQ0t7l9xnIzZv3IM05Z1XZIRHZf2FO93njTAXectXh656BZH0jaJ/d0tnQc07jMHFARDqgIB1SEAyrCARXhgIpwQEU4oIrcItilmc0y19pV1/ac1q3Vz0OZBfP+bQl/wefzk/e4fRriCbN+yq2rj8jesaLd+xkiIqe1b3D7FJbtd7TdNLbXrA8k7KcENcwcUBEOqAgHVIQDKsIBFeGAinBARTigitwi2M6lwzKeq+38+n0RuTV3sHr99ObN5v2fndjpjnFSxl8U+uXOTWb983P3h7Z3rOEJtqCuBr//D8fsJ/QqyxWzvuzUNcwcUBEOqAgHVIQDKsIBFeGAinBARTigitwi2MTSnIwu1j+KGLz+3Mxt5v2nd57gjrF3evUurlV9xO5TVhaWRhZqRz19qP8Z7jj/Z/hHbp8NLd1m/fqujWa90tsh9isEw60pHIODg38iIr/55OW3s9nsHw4ODj5XRD4uIk0icl02m33fk33PEZHPiEibiNwkIm/KZrM/+wFb+Llz/1p5MgTPE5FzReQcETl/cHDwFSJyjYi8REROF5ELBwcHX/DkLdeKyFuz2exWEYmJyBufii+Op95a/ptjSETemc1mC9lstigie0Vkq4g8ks1mDzw5K1wrIlcODg5uEpGmbDZ7+5P3flZErnwKvjd+Dty/VrLZ7INHPw8ODp4qP/nr5e/kJ6E5akhENorIeqV9zb5w3dWr2m648bqf5kf8Qt1263/8VP3vkD/4mcdc/Jl/Qrg1/wfp4ODgNhH5toi8S0RK8pPZ46iYiFTkJzPRckj7mr3y5W+VkeGx6vUNN14nz3nGy6vX8wX7t2It/0G6b2bI7eMJ+w/S2279D7nk0l+rXkfqP0ivfqc7zkpr+qfs4ODgZSJyg4i8J5vNfk5EDovIukCXARE5YrTjv6C1/AfpCSLyHyLyymw2++9PNt/xk9LglsHBwYSIvFJErs9ms4+JyNKTYRIRebWIXP8UfG/8HKzlr5WrRKRRRD4+OFg9fuifROS1IvKVJ2vfEZEvP1l7lYh8enBwsE1EdonIJ3+aL/SyzrNkoTRb1/bantqR1f8wdIt5f3b6CXeMllSj22epVLQ7LNtlEZFvlYfdPpva/I1Hl2XsDU7XlFJmPVNqkRe7o6y2lv8gfbuIvF0pnx3Sf7eIXHQM3wURw/I5VIQDKsIBFeGAinBARTigIhxQRW6zz/cX9svYXO2Jt1eLyNfnHqpepxuS5v15b/FKRAplf3tJqWIftTTYEX5cU7D94OKoO8580X5nnYjIxy8eN+svuidm1ruKS8e0CMbMARXhgIpwQEU4oCIcUBEOqAgHVIQDqsgtgs0XF2W2kKtrC16/oPss8/41bNCS3Yv+E29T+Xmzrm1SDrZ/oufp7jhXTd3q9tl4w0GzfmJrr1lPlNLuGGGYOaAiHFARDqgIB1SEAyrCARXhgCp66xyFvMyteFg6eP2DqT3m/Yc//qvuGF2/77/jLR6z/9wMZDrd9r9cDD8C+6fVlm426+NLs2Y9kWedA8cZ4YCKcEBFOKAiHFARDqgIB/4fe/cdJ9dZ3v3/2tmd2dnei+Qm2bLHltyNe6MnIaGFOAktT4JD+QGhk5gQCDUhQEIgDUIzBEh4AtgJLZjm3nFvY0u2LFna3svs7uzs/P7AmrJzvve1VgSv87z4vPnnnPu6j+6V/NUBX5xzH4lwQIpdE6whkbCGRP26sfJ5eyrcEDrtz69116iz8BtiZmZNDeGtlAYbO93x+2b2uOu0pZrcOee0HR2sv2o5/J24+nS7u0YU7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE2ygudMSLdVbLm1u6S4d75kbW39JlecNnOqu4T05ZWb2rK4TgvWFYvT2Ut315Sbd5ubuyDlP9WeZW1sJ1l+38liw3p/vtc+6q9TizgGJcEAiHJAIByTCAYlwQCIckGLX55hZWajZVafy3PuM5nVTWXeNrsZWd05HXXgb7Z8tRj/I81CuvLPPtqYBd51Uwv9H8Egu/K2441o2B+tdzdEPJnm4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNsITVWaKu+o20yvN0fbg5dUpXxl1jdDW8dbWZ2TWL4Qdojkz3uuM/Hb3fXSdV7/8j2N5xRLDemQhv69Rex7ZPOMQIByTCAYlwQCIckAgHJMIBiXBAil0TbCG/bPMrS1Vjlee/0bUjeP23xvx9zVfXCu6cPxk4P1j/7MRtkeN3TpebZ88fOM1d5+i68JZNZmafHL4+WL+nLryN1WCiz97orlKLOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBGtJNlprKl01Vnl+e+6J4PW/03e6u8bXhm9153x19p5g/SU90dtLVY5/aegmd523b7rIneO9Ajq7shisNyd5EgyHGOGARDggEQ5IhAMS4YBEOCDFrs9RsKIVimvVYxXnuUJ4q+cfzz3srrH+jboom9Jdwfo3J+6qGXv9uvG/G3iGu062Lvz7MTNbcn7PnY3hB4a87+Ip3DkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7JlhfY7vVNa1WjQ02lRtSj82OBK+vT/h5P7w1esumSouF5WB9pbDqjt+SWHDXGTD/QZyGRH2wfmLzYcF6R5q9z3GIEQ5IhAMS4YBEOCARDkiEAxLhgBS7JthjsyM2ND1aNfbw9L7S8Tl94b3N9y1Nums01afcOXOruWD9L/qit4WqHP/e6lDknEpfH/M/WOg1wRpaw3/HG8x/8i0Kdw5IhAMS4YBEOCARDkiEAxLhgBS7PkeyvsFS677jVnn+vtXwgzpfaQ2/qWa2sZ19XrvpvGD9L/dfUzP2a/Y2+8vh8vjWjgF3Ha+HYeZvx33VxH3B+mCyz97grlKLOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBNvS2m9t7dU/1rHtm0rHl+YfCl6fWPXzfmr3VnfODxZ2BuvbOje543P58ANDZmbv7b/QnfOlXPiBoJHcdLCeqDu4ewB3DkiEAxLhgEQ4IBEOSIQDEuGARDggxa4JNlfI2cxq9TfLKs9/s/344PXfnn3QXeN3kke4c/41H/51zmw+PHL89Kby+BVzd7rrfCB3nTvHa2Ktrdsrfr3lVX9/9ch1D+oq/EogHJAIByTCAYlwQCIckAgHpNj1OaaXF2xiaa5qrPL86yu131arNLeygQdsFv3ewkf7wt+Yv9HmI8cLViwdF4vFyDmVntd/ijtn1/J4sF7v9EF62rvdNaJw54BEOCARDkiEAxLhgEQ4IBEOSIQDUuyaYD3pNltrqn44pb+po3S8vkG2ntcQMjP748Fz3DkfmQk/qPPtlm2R4+8slBtfP06l3XW+MxJu6pmZdTQ2B+uN67bJqrHi1AXuHJAIByTCAYlwQCIckAgHJMIBiXBAil0T7Nj0gPU1NVaNbW/aXDr+4fy9wevXNvD01ZWzD7hzJnPhZtu5i3fUjN1kZueOlcc/NBB+mszM7J+c34+Z2czKYrD+N21nBuvJpnZ3jSjcOSARDkiEAxLhgEQ4IBEOSIQDEuGAFLsm2P78jE3kp6rG9lacb+uI3nO8dP3ipLvG9PKCOyfdkArW1auOTRXXfWz6dnedyzrDDSwzs0anr/f6sWuD9cFkn33TLnXXWY87ByTCAYlwQCIckAgHJMIBiXBAil2f45GZIRueGq0ae2Bqb+nY28p5I/6n63x3zu8s1D7MU+n2o2rfeJtdN37KY+Fvs5mZXTZ8tTsnWR/+x9TuvBHXmmpy14jCnQMS4YBEOCARDkiEAxLhgEQ4IBEOSLFrgg22dFp9a6Fq7LDWntLxrPP213Ih767x/Jlb3Dkpp/F0zt7HasauWjf+Z/1+s+2Hq0PunNZEY7A+mAg3wdraO4J1hTsHJMIBiXBAIhyQCAckwgGJcEAiHJBi1wQ7vmmTbWqpfnLp5JbDS8dT6fAH/jY3tLlrLBUL7pyrJx8M1hN1de74Jyb8Zttiftmd473lN5qYDdZ7VlaCdYU7ByTCAYlwQCIckAgHJMIBiXBAil2f4/7FfTY6P141dsf846Xj4fmp9ZdUqRP9h0oNiXp3TtKZ84Kuk9zx78+GeyVmZgNNne4c1VM54B9twPkVeNgHhxjhgEQ4IBEOSIQDEuGARDggEQ5IsWuCHZ7usZam6sxubeovHTfXh9/+Gl4MN8k2ynsI58tDN9eMvX7deDoZ3iLbzOyB3/QaWGZHXfF4sH5pS3grrN61Rft7d5Va3DkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7JlhLImX5ddscVW57tGctvK1TYQN7o2/kSbCz+44L1nctDEeO9za3l47T9X4T7ITvRP86ld7bc26wfvnyzmC9oe7g7gHcOSARDkiEAxLhgEQ4IBEOSIQDUuz6HI/khm1scaJq7IHFfaXjhfxS8Pq39IV7AmZmfztygzvn7undwfopnVsix49pGSwd3z4R7j+YmT29b4c7Z1civDPP/sXJYL2wRJ8DhxjhgEQ4IBEOSIQDEuGARDggEQ5IsWuCHdbYY83p6swele4rHd8491Dw+mtXR9w1Lug93p3zk5F7g/XfTWx2x28qhH9WM7OE+dtUfXXsZ8F6e2P4G29NDeG3BBXuHJAIByTCAYlwQCIckAgHJMIBiXBAil0TbO/SuI3mqvc+fzRXbmy9fvP5weuvWdrrrtHb0OrOaU6lg/XbEos1Y+evG29s8N94u2s+vKWTmdnxHYcH63dOPBqs1y36P0cU7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE+xdqeNttXG2auxjjSeXjl+5/5rg9QMt/gf19iyMuXOanaenvjl6R83YW9aNb27pdtd5Vssx7pwfzD8crGc6DwvWe9v9nyMKdw5IhAMS4YBEOCARDkiEAxLhgBS7Pse7F+6ykblyH+LHZva6uVtK5yd1HxW8ftecv130SmHVnXNYR0+wvqWlP3L8pK7yz3f35GPuOh88pc2d85+31j5YVGlkYTpYn2/yf79RuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbAT2g63wVxL1dgpHVtKx5Or88HrP9np733+x2M/defsmQ8/ELSYX44cv2N8V+nY247JzOzE2/1tqjpSLcH6Jd2nBuutvR3uGlG4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNsDsmH7Xh8dGqsZvGsxu+/o+L/rZPqfqkO6dQXAvW10S9cnz/HV901+k55RXunKlcuPF3XX14W6fepW57mbtKLe4ckAgHJMIBiXBAIhyQCAckwgEpdn2OpmTKmpPVu+pUnnv9h560/wbZkeled86+5clg/eqtXTVjo2a2+9Ty9+OOP+cN7jpLqyvunE2t4Z159i6MB+v5nP8duSjcOSARDkiEAxLhgEQ4IBEOSIQDEuGAFLsmWGsybYuppqqxtorzJufBlvnVJXeNs5J97py76uqD9XeP1n4n7s1m9u7R8vbaT8yFt6U2M6szv0F1emt4q6vvzt8ZrC/l/UZbFO4ckAgHJMIBiXBAIhyQCAckwgGJcECKXRNsZnnBppaq3/CqPB9x9i1vTaXdNX6y5L8V190Q3mrpmrmdNWNvXjf+kk1Pc9f5xtBt7pyrRu8N1r/bdUGwnuhod9eIvO6grsKvBMIBiXBAIhyQCAckwgGJcECKXZ8j3RB+460hEX4IZ1vrJneN3Yuj7pwnbCJYV9tbTyzNlY6f3VT7QNB6/1ksunO83k02Fd6pqDGVtIP5mj13DkiEAxLhgEQ4IBEOSIQDEuGARDggxa4JdnTLgHW2Vr/Vdlzr5tLx7RO1D9lUen7qCHeND0w86s5pcRpPhbXo7acqx98+fZO7zua2HnfORW3bgvUHLLohd0Br3Yr5jx3V4s4BiXBAIhyQCAckwgGJcEAiHJAIB6TYNcGGl6dtbKl63/F9FecdjeE30b4wf7+7xkBLpztnZS38Zp3S2FB+KuvtPWe783PmPwn2dyPXB+vbO8ONv56VJXupu0ot7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE2xyac7GcjNVY5XnX207J3h9m/nNq99fCm+jZKZfdzxgeTUfOb6wUt57/Ws5f+/z3bP+q5lrziuTD04/EawPpsO/F4U7ByTCAYlwQCIckAgHJMIBiXBAil2fo7OxxfLp6u+09aTbSscvn7s5eH1TQ/gbcGZmJ7Yf6c4ZWZkJ1nfODEWO19WVv9n2P4f520ofOzPszjm1e2uwfqfzBp96O8/DnQMS4YBEOCARDkiEAxLhgEQ4IBEOSLFrgp3cfLjNtLZVjT2tdUvp+PrCI8HrF/JLwbqZ2Q1jD7lzvD3WN+IF+xfdOdu7/G2qvnl4+B/T0eFt2g8adw5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCfbA0pCNLVZ3de5a3Fs67k93BK+/8ZXhp6bMzLZ/ba87ZyI3F6yrJ84qx++f2uOus5Fm29NXw/ujr/9g4nrppP90XBTuHJAIByTCAYlwQCIckAgHJMIBKXZ9jkvTx9ly02zV2FubdpSO3z0Z/m7awOf9N8g2tXT9r+fMrkQ/yJOqL/+Rbm7ZHDmn0msbM+6c902F3/JbKYR3M8o7dYU7ByTCAYlwQCIckAgHJMIBiXBAIhyQYtcE+4/8HpvIl7/pdq6ZfS5f3tYo7WzrtLXJb3D9VdHf9ukFU9cF66n6ZOT4QsW22FNL+911Gvv9JljlltlRKreailJwtsdWuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbA9C2M2PDdWNbZrrvx0V9J5Q2wg6e83/tLJ29w5g86TYIVi9H7iHY3NpePf6tjurvOmoavdOVs7BoL1P02dEKynuv0/kyjcOSARDkiEAxLhgEQ4IBEOSIQDUuz6HKuFgq0WCjVjB/xp77nB6z8ze7e7byinawAAIABJREFUhveGmJnZjtbDg/Vrxx6MHJ9eWigd/08i667zzIET3Tk3Tj4crF/THN5Gu70+aTuCM6Jx54BEOCARDkiEAxLhgEQ4IBEOSIQDUuyaYKd0bbUj851VY2f1HFs6/tTU7cHruxpb3TX+cOBkd85Hhq4J1pvEltIN9eWHkcYXZyPnVPpAeps75zznjbev7Q9vC7Wp2G+vclepxZ0DEuGARDggEQ5IhAMS4YBEOCARDkixa4LtWhi20fnxqrGH5veVjluS6eD1j874e5//ZANvxXU3tQXrL+s6NXL80r6zSsefH7vVXef3lx9z52xx3nj7w5bwm3VNfbzxhkOMcEAiHJAIByTCAYlwQCIckGLX5+hrbLdEuvqNtMF0+eGfB6efCF5fZ+Gtns3MdiTD34Y3M7trJdx/6CpG7zBUOV7cwLbSLQ3h79CbmR2ZCv+8/7b4ULDeu9Rjz3BXqcWdAxLhgEQ4IBEOSIQDEuGARDggEQ5IsWuC7ZwdtuHp0aqxysbXaT1HB6+/6tfD34AzM9v6n/7WUOmG6G+4HfChiDfifsveVjV+Xt/x7jq3TjziztlpQ8F6e8WW2lFSK03uGlG4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNMM8d47uC9U3/7v+WPtFzgTvno0v3B+und2yNHL+4v/z22Z0zu911KreJUtL14Ybcn7c/LVhvbOONNxxihAMS4YBEOCARDkiEAxLhgEQ4IMWuCVZfV2f1dYl1Y+Xz/+48J3j9dzbw0NNbR6535zynL/xhvqsnoj8AeNtUuUnXkPAbXP1NHe6cPbNjwfqfLof3aR9M9NsVdqm7znrcOSARDkiEAxLhgEQ4IBEOSIQDUuz6HIVi0QrFtXVj5fPfnL4xeH3rYnj7azOz5dUVd849C+HtpZYLeXe8K+1/b273zIg7R31P7oCxH304WF+xBtvprlKLOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBPP8uOvsYP07jeG3w8zMPr4//HCMmdlxzYPB+vDCVOT4RvZerzTQ2uXOmcrNB+sdz/izYH1wsN+u+Na/PqWfy4w7BwIIByTCAYlwQCIckAgHJMIBiXBAil0TrLGhwRrX7Tteef70ifCTYHV1fhNqS8eAO2ckPxusn9y9xR1/bMF/ymt6acGd8+7Bi4L1o1bDv+dkN9s+4RAjHJAIByTCAYlwQCIckAgHpNj1OS7tOs1yq9U9hrf1lnfz+e7K3uD1q8WCu0ZHffibaGZme5bGg/UTmjdHjvfUt5SO788vu+u0pvw39D4ze1ewPr4Y7skMNrCzDw4xwgGJcEAiHJAIByTCAYlwQCIckGLXBBu3VZu36i2VRirOs7P7gtfnNtB4Gmzx3zJbXA3/OjcuPxw9PlUeP73raHedB+fC20uZmT29/bhgva8jFay39PGwDw4xwgGJcEAiHJAIByTCAYlwQCIckGLXBPuvqXtteLz8TbNXm9lXx39WOl8thJ/0esvmC901PjV0gzunqSHcWFq/P/sBq2vln+/RRf+Nt4UNNO2aLPytuLNWwv8Yk/mD+8fMnQMS4YBEOCARDkiEAxLhgEQ4IMWuz+FZv+vPeu+/1P8tffKDRXdOfSL896ZQiO5zVPIeGDIz27SBB4+uWXg0WP/e6lKw3p/rtS+4q9TizgGJcEAiHJAIByTCAYlwQCIckAgHpNg1wYrFohWLxZqxAxKJ8FbOx/7dne4aZ/Ye684ZbGgN1v/x8Nrvrg2Z2cOnbimdH33HLnedmQ1sb93d1Basz6+Em2DpZX+bqyjcOSARDkiEAxLhgEQ4IBEOSIQDEuGAFLsmWGdji+XT1U2dnnS5CfTEXHhP8qWGfLBuZtbSEn6bzczs28PhZtprl8+vGes1s3t39ZfOV9eit4aq9NFNz3TnfGIuvPf51FptQ65S5Vt4TwV3DkiEAxLhgEQ4IBEOSIQDEuGARDggxa4JNpdfspmVxaqxyvNNrd3B69uT/lNP+1am3Dnn9mWC9ZdM3Vozdo2ZvWSuPH5C1xHuOu8aucadk6gLP/12++AZwXqxr8P8lzcj1j2Ia/ArgnBAIhyQCAckwgGJcEAiHJBi1+c4q+Nom16u7mVc2FXuOdy1sCd4/e65UXeNgeZOd85sIhesb20bcMcPS/lbOp0xeJg756rZh4L150xng/WBxl77krtKLe4ckAgHJMIBiXBAIhyQCAckwgGJcECKXRPsHOuwpXUPtzy9rtxMmm9eCV6/Ulh111gt+m+AzeQXg/XZlej6/sXJ0vGe+bHIOZXSzrfkzMx2tIUfGppeDW8d1dMUfkBK4c4BiXBAIhyQCAckwgGJcEAiHJAIB6TYNcE+NnGLDY+Wm0c32Fvt3aPXls5/s//U4PVbWvqDdTOzdF34I4JmZkvF8PZRI4vTkeO51XKT7pX9Z7rr/GhhpzvnxvHwk2BHtYV/z01r4cahwp0DEuGARDggEQ5IhAMS4YBEOCDFrs9x1dZjrL6tp3Q+bWb3H3t86fyEh8PbTq//PlyU9kZ/95/ppfCW0cn66D+6tWJ5D53P7b/BXecbPU9351zWFt4me+f0/mB9ocl/ACoKdw5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCfasXY/Y0HB566ZbzCyTvb90fnzX4cHrR3LRD+FU6m5sc+f0ptuD9We1HB05/qrBc8rrFP0/3tfPhL/fZmY2tjgTrJ/bf3yw3t3jbz8VhTsHJMIBiXBAIhyQCAckwgGJcEAiHJBi1wQ7ry9j08XqfcWfNXBS6fgnI/cFr69P+HmfWwnva25m1pCoD9a/FtFsu9TMvjZ+R+l8e7v/jbcp54kzM7OpN4S/4db/L3cH64PJPneNKNw5IBEOSIQDEuGARDggEQ5IhANS7PocN45lbWik+jttPx65t3T8p5suDl5/Zc7fKWcj204/r/fk8Dojd0SOL+aXS8d7l8bddfIb2I572xceCdb/uu/CYL2xJ/zgksKdAxLhgEQ4IBEOSIQDEuGARDggEQ5IsWuCndh9lG1erW7anNG7rXR81cre4PUfTmwL1s3M/rzV3xrqrtwTwXpPU/Rbc5Xjs8vh78SZmXWLX6eS9+v8Y+6BYL1vqcfOclepxZ0DEuGARDggEQ5IhAMS4YBEOCARDkixa4I9vjBqI3PVT1DtnBsqHV/clQlen7U6d41fSxzjzvni2G3B+op4gmt6aaF0vJG3705pP8qdc/Xo/cG6t9VV3ZL/Tbso3DkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7JlhzQ6O1JtNVY5XnX738hcHrO37j/e4aram0O8drYBUL0U+TFa083pVudde5fvwhd06d09hbyq8E68v5vLtGFO4ckAgHJMIBiXBAIhyQCAckwgEpdn2O2ZWcTS8vVI1Vnh/zkk8Er+9v7nDX6G30t0Eayk0F64m66N5Boq789+2iNv/tu9ckw9tom5k9Z/qWYH1TS/gbbn0t/p9JFO4ckAgHJMIBiXBAIhyQCAckwgGJcECKXROsLdVkucbmqrH2ivPh+XBzanNrt7uG1+AyM8u0bg7W31QYiBz/csd5peOXDV/nrvNEb/gNPjP/99TR0Bystzt1hTsHJMIBiXBAIhyQCAckwgGJcEAiHJBi1wSrr0tYQ13101GV583JxuD1o4sz7hrtjX5TaHglvJXSK+ZqP8p3g73dXjF9fel8I0+l3TSedeek6sP/mO7/wNOC9ZV0mz3qrlKLOwckwgGJcEAiHJAIByTCAYlwQIpdn2N1rWD5teqtoyvPD3bHnUr/mj7VnfOapbvCv0b3he74Zbnwr2FmtlZcc+f0pMPfgev48/8J1jcN9tuVV7zSXWc97hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE2w8N1fzwE7leX1dOM8DzZ3uGi+dvcmdc3zHYcH605pr35rLrR/PucvY4W297pz122Ct9xsD4aZeZ5//ZxKFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBLNi0YrF4rqh8nlrY/j7bM9uO85d4vK5G905D8/uD9YvXBirGbvKzC4c2Vk6/7XuHe46d+aecOc8sye8h/pkMfyNt8RB3gO4c0AiHJAIByTCAYlwQCIckAgHpNj1ORqTSUsnU1VjledHtvQFr//SkP8gz+9tPtudc93czmD999tOjBx/VfcZpeMrFvxde/YvTLpzfqP5mGD9S/vDv+dN1m9vclepxZ0DEuGARDggEQ5IhAMS4YBEOCARDkixa4Id0dxrza3VmT26tfw9tV1zw8HrT+nZ6q7x9f23uHMaG1LB+n/YfTVjv29m/zFXHp/IzbnrPL13uzvn/xTDb7x9pbk9WO9Mt7prROHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBHp7ZZ0OTo1Vj900+Xjo+sfuo4PVzq/5eS82p8FtzZmandx4drH+2rfaPbs7Mftxf3i7qgqHd7jo3Tj3szjm3sBqsH9Ea3jqqI9XirhGFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBPPk1sJbHN19xVvdX6P1GX/qzrlj+tFg/dzZ2j+6q8zs3H3l63Kr4Z/VzKx1Aw255dV8sL5rZihYX2wquGtE4c4BiXBAIhyQCAckwgGJcEAiHJBi1+fY3NpjyXUvcB3Z3l86HlmcDl5/3Av+xl3j3P7j3TlNdclg/fHl8cjxnnRb6bi5vtFdJzu9z53j+XTfM4L1VE/4jTiFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBGtI1FsyUV81Vnm+rX1T8PoHpva6a3Qn/W2QJovzwfrR6X53PLsYfgjHzOykrvAbfGZmp6UGgvU/GbkuWB9M9tk37FJ3nfW4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNsCfmJmx4tnrbp90V58ViMXh9X0uHu8b+3AY+ute1I1gfWVuMHC9a+edbXgu/qWZmdsfELnfOnRZ+++7uI08O1tf6Om3ZXaUWdw5IhAMS4YBEOCARDkiEAxLhgBS7PkeizixRV7durHzeUB/+kb0+iJnZi7pOcud8c/LuYP2yrrMix59V11M63pOccNfp6/LfRssXwzvzPH86/FBRX+OK/ZO7Si3uHJAIByTCAYlwQCIckAgHJMIBiXBAil0TbEfnkbZppa1q7LTu8vfWLkiG33j7/OTt7ho7CzPunDWnmfbBiZtqxp61bry70X+zrivtz+mpbw7Wfzh2b7BeaD24ewB3DkiEAxLhgEQ4IBEOSIQDEuGARDggxa4J9tjCiI3MVe8r/vDc/tLxz1bCb4htbQ9vkWRmdsPYQ+6ctlRTsL4kvt+WL6yWjufzS+46709Gbx9V6VVLDwTrH+m7KFhvZO9zHGqEAxLhgEQ4IBEOSIQDEuGAFLs+x4s6d9hCfrZq7BU9p5WOvzB2W/D6qZXwttRmZp3pFndOc0P4+2wLG+hhHN7c4875ctrf/eecxqOD9U/M3xOs9y32WPT7eWHcOSARDkiEAxLhgEQ4IBEOSIQDEuGAFLsm2NfG7rDh4bHS+R+Y2WeGy2+R1VldxFVlnc1+g+uxmRF3zkpqNVh/88D57vgXp+9013lheqs758a18WB9c7o7WO9u7HTXiMKdAxLhgEQ4IBEOSIQDEuGARDggEQ5IsWuCndy1xQ7PV3+n7Wk920rHD8zuDV7/6Mywu8Yz+k9059w3H14nLZpxlePTSwvuOg9b9LfiKu3KjQbrld+Vi7K6wrZPOMQIByTCAYlwQCIckAgHJMIBKXZ9jq5EsyXqq3fN6a0vP8Azv+K/aeYZzvvbW7c0pIP1f576Wc3YC9aNt6bCv4aZ2R8s+XMeTYV35tm9GO6D5MQuRB7uHJAIByTCAYlwQCIckAgHJMIBiXBAil0T7Kdj99vQcLmp8x4z++5w+c2xSzadGbx+em3ZXeOHI+FtkszMfn3glGD9lcmOyPF/ajmjdPzyiWvddX5v7VZ3Tkcq/I23qVx4q6v0Uvh6hTsHJMIBiXBAIhyQCAckwgGJcEAiHJBi1wQrPvmf9WMHfGuk9gmsSpMff767Ru87w99MMzO7eWZnsP5gqnZ7qa+b2Z/ny9+PO6rd/37bnrkxd854YTZYP7K9L1jvawtvC6Vw54BEOCARDkiEAxLhgEQ4IBEOSIQDUuyaYM2ptLWmmqrGKs87G8N7m7/2Y0PuGk0NKXfOWjG8ldLmxi53/LcTm9x13jV/nTtnR9eRwfqj8+GtrlIrTcG6wp0DEuGARDggEQ5IhAMS4YBEOCDFrs9xQtvhNpir7mWc0rGldHz9aPhBnW8sTLtrdKb978DNreSC9eGV6HUqx7/lt1PsPf0XuHO+ld8TrJ/aEf5OXFcb33jDIUY4IBEOSIQDEuGARDggEQ5IhANS7Jpge5fGbTQ3XjX2aG6kdNzs7Ceey/vbPl27+Qh3zkunwt9e2780GTk+v1rem31bs998eu/Q1e6c9Q8/rbdSWA3WB5PhN+IU7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE2y5kLflQr5m7ABvH/CedJu7xjOG9rlzntueCdafn4p+C+217aeWjn+8Gn4Tzcys6LxZZ2bWWJ8M1vubovdhP6C3JfrtPA93DkiEAxLhgEQ4IBEOSIQDEuGAFLs+x+zSYs03yyrPP997cfD6W1L5YN3M7Ctj4S2yzcx+trQ//GtMP1Ez9lv2NvvwcPm7bid0Hu6uU1dX584ZX5wJ1s9qPzpY72zkjTccYhu6c2QymQ+Y2e+YWdHMPp/NZv8uk8k828z+zsyazOzr2Wz2L56ce6qZfc7M2s3sWjN7XTabDT/Hhlhy7xyZTOZiM3ummZ1sZk8zsz/JZDKnmNkXzOyFZnaCmZ2ZyWR+48lLvmJmb8xms8eZWZ2ZvfoX8YPjF88NRzabvcbMnvHk3/5++/ndptPMHslms489Of4VM7skk8kcZWZN2Wz25icvv9zMLvmF/OT4hdvQf61ks9l8JpN5v5m9w8z+08w2m1nlno5DZnZ4YHzDrrjiszVjN9/0Xxu+fscG5rzqKfw8T9VNN175C/zVf7k2/G8r2Wz2LzOZzN+Y2bfN7Dizqi/m1JnZmv38ThQ1vmEvfvGrq75IffNN/2XnnPvC0vkv699WtrSGP6TzYMS/rdx045V27nkvKp1v5N9W7p8K771h5v8/t88bPC1Y7+zvtLf869vcddbbyP/mOP7J/5Fp2Wx20cy+ZWZPN7PKHVgHzWy/mT0hxvH/oI38q+zRZvbZTCbTmMlkUvbz/xH6GTPLZDKZbZlMpt7MXmZm389ms4+b2VImkzn/yWtfaWbf/0X84PjFc/9rJZvNfi+TyZxlZneaWcHMvpnNZv8jk8mMmdk3zSxtZt8zs288ecnL7edhajezO8zsU0/lBzq37zibLlbf0p85cGLp+NKRa4LXNyTq3TWe13+KO+e/h+8I1t+56aLI8XcMXlg6/vLsPe46LcnwG3xmZjs6wttbLxTD/1WacurKRv8H6fvM7H3rxn5sZjV/ytls9m4zO+ugfhrECh1SSIQDEuGARDggEQ5IhAMS4YAUuyfBHs2N2tjiRNVYdrH8/+Xt6A43hGby4e2azMyun3nEnXNuX/iNtx+t1P5/Ky9ZN97fGH4Tzczs6OYBd87jubFg/Yim3mA9Xyy4a0ThzgGJcEAiHJAIByTCAYlwQCIckGLX5xhbnLHhhamqscrzEQt/w6052eiusZEHgrLz4acbL+44LnJ8S7K8i05znf/HO7S24M55UfsJwfrlY7cH6wNpv/cThTsHJMIBiXBAIhyQCAckwgGJcEAiHJBi1wTrTrfZalP1d9p6m9pLx23J8PfOTmvyX17ebuEtss3MPj17d7C+vzDnjj8wu9ddZ6DJ35Lp9tyjwfr67cDXW3HqCncOSIQDEuGARDggEQ5IhAMS4YBEOCDFrgm2rXnAelqqn+Y6oeWw0vFt07uC12/kjbd7Uq3unPHcbLA+KvYjv3W8/DbdkW197jqPzY64c7ytoR47M7z3+Wp3l40GZ0TjzgGJcEAiHJAIByTCAYlwQCIckAgHpNg1wUbzszaer37lcaji/PTOcMPn7pnd7hrP7tjqzslO7QvW1fZT27uOKB0/Nuc3uNINKXfOcW2bg/XXPR5u6nUsttib3VVqceeARDggEQ5IhAMS4YBEOCARDkix63PsnBmyocnqR1MenCy/ObY3NR68vq/J31L6M8M3uXM+2/f0YP0149Hfmntgqvyzfq4n+jtwlV47eb075/bxncH6cwZODtZXn9rXW0u4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNsIZEvSXrq3+synPv4ZixXPSbaJVW1/xvnk07f20+PHCxO/4nE36z7YROf5uqmdXwW3y7lsLvs/Uu8403HGKEAxLhgEQ4IBEOSIQDEuGARDggxa4JdnTboHXkqhtdmY7ytk/T+fAH82aW/Q/qHdnub8f0L0sPBevnNx9VM/ZMM3uwLlc6X7Oiu87vp7a4c/5++a5gPZUI/2NcWWPvcxxihAMS4YBEOCARDkiEAxLhgBS7PseF6cNsoamtauw5TeWeQr4p3Dv477qsu8bU8rw754jm3mD930durxl7w7rxLe397jof3MADQUe1hvsyC6vLwXqd1blrROHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBvjv/iI3NTpTOX2pm/zn7QOn8ipYtwes/szDprrH+jbqD0dbY5I7vnvW/rHZmzzZ3zpLzsM7Iuu3A12vJL7lrROHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBxnOzNrJY3dSpPD9zrvYJrEqNG/hmWv8G9kfPzoS/8Tbxl8+oGXvIzB677OzSed8HovdHr3T75C53Tnc6/A237W1HBOtdrV3uGlG4c0AiHJAIByTCAYlwQCIckAgHpNj1OQ5r7bXGtuo3tI5qK785tuw8+DK66G9vvW9+wp3T5fQWOv7yRzVjtzz3TVXjo7/pP8gz8F2/z1Eohr/RlnZ29mlM1LtrROHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXB9s2P2/DcWNXY43P+m2MHvGjgdHfOlSN3uHPq6sJbJR3Tsckdf9Gt/h/vRt6+m14Kb9n9t93hJlexod5W3FVqceeARDggEQ5IhAMS4YBEOCARDkiEA1LsmmBndB9jU6s9VWPn9R5fOr5tamfw+o00uJ7Zd6I75+bpR4L1+rrov1eLhfI+5Kclw/unm5nNduTcOXsWxoL1dzmfcOtcNXuLu0ot7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE+y28Z02NFr95Nd1o+W9z+sT4TyvFcMfCDQzu2rkbndO2tk+6ojm6Nclu5Pl8S+N/8xdZ987z3LnHPv34Vc8xwvhJ8XWCv5WWFG4c0AiHJAIByTCAYlwQCIckAgHpNj1ORqTSUsnq/+9vPL8XT3nBq8/c3nVXeM1q1l3Tk+qLVgfW47uPUznyz2H+RX/QZ7NH73FnbOYXw7WR+angvVNDf3BusKdAxLhgEQ4IBEOSIQDEuGARDggEQ5IsWuCLa+u2vJqft1Y+fxvp24LXu81jMz8LZ3MzJ7WclSw3tsQ/bDPCS2HlY7Pat3qrjPqPKhjZrZnaTxYP7H5sGC9s6/TXSMKdw5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCdaQSFjDug/WVZ4v5JeC1xfWwh/LW//rKVdN3Besr1n0m3U3TJafMmtuaHTXGd/ABwubU+lg/Rs3fDJYX1mrs4f9h9JqcOeARDggEQ5IhAMS4YBEOCARDkix63O0pZos19hcNdZRca62lT6gpSHcEzAzu6B5iztnthjeM3p7XUvk+Nt6zykd/8Nk+MEkM7OE8/sxM2tydhk6/KzXBOsDA732pf/4lLvOetw5IBEOSIQDEuGARDggEQ5IhAMS4YAUuybYjtbD7bDF6rfJTmvbUjr2tqbeSFOpKB7UqZQrrATrdyaSNWMvMLOvLjxUOve2yDbzm3pmZp9qPiNY/0z9aLDe1dblrhGFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBHtOsdNWitWZfWGxu3T8WOfm4PU7p4fcNTobop/iqjSaC7+JVkhGv1m3WChvOzW24L/NthEr0TtMldwy8UiwPpjsO6h1uXNAIhyQCAckwgGJcEAiHJAIB6TY9TmW6syW1+0+nas4z68Vgte/YfMF7hrfX9jpzlkphL8Vp7bRHl+cLR1vZBvtEzoPd+e8ff72YP2Y9sFgvbe1O1hXuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbC/Hr/FhkfLb3DdZG+1Px29pnR+ft/xwes/PXyju0Zz0t922ttq6ZvNp0SOX9VZ3vbpnxvDb82ZmX179C53TqYz/A23/bnJYL24fHD/mLlzQCIckAgHJMIBiXBAIhyQCAckwgEpdk2wtWLBCuue9qo8zy7sD17f29TurjG8MOXO+WT/M4L1P8rVfgPuG2b2R6vlt88uTe9w1/nWWviJMzOz7PS+YP35/acG6+1tHe4aUbhzQCIckAgHJMIBiXBAIhyQCAckwgEpdk2wnqZ2K7RUP0HV39JZOr6k48Tg9VctPuquMZfMuXNmEuH90dWe5ZXjn5+/312nK+3s6WRm3Y1twfr3J+4N1geTffb/uavU4s4BiXBAIhyQCAckwgGJcEAiHJBi1+coFNesUFyrGTvg8rHwFkgf7TrXXeODef8ts2uKE8H6XD66V1I5vlzIu+vMr/g9l/V/Huu9ufecYL25y38AKgp3DkiEAxLhgEQ4IBEOSIQDEuGARDggxa4JVvfkf9aPHbCQXwpe/1eL97hrnNG2xZ1z/dTDwXqmLXorpmNayvuQ3zvzuLtOR6P/vTmvCXaJzQXrRau38I7x0bhzQCIckAgHJMIBiXBAIhyQCAckwgEpdk2wlmSjtSbTVWMRbJSPAAAgAElEQVSV5+9rOyN4/ZV14+4aD+WG3TnHtIY/qnfreHSTrHL8uQPR+6NX+sGw/1Raa6opWH/udPiDhv2NvXa5u0ot7hyQCAckwgGJcEAiHJAIByTCASl2fY4zmo6wuZbqN7TOb9laOr5s/Kb/9Rqv6T3LnfMfc7XbV1d64WB0v6Vy/L+Gf+au8+rN57tzPjcU/m6d9wBUaikdrCvcOSARDkiEAxLhgEQ4IBEOSIQDEuGAFLsm2LWzO210uvzAzpvN7AfTD5bO0w2p4PX5DXwz7QuTfnOqJx3eUjqqwfXn68b/ZjD8nTgzs/+b3+vOaUmGm1h/1hPe9indx7ZPOMQIByTCAYlwQCIckAgHJMIBiXBAil0TbGZ50aaW5qvGKs9P7t4SvH7n3JC7Rs8Gvqt2WGN3sP6hnkzk+Nd6Li4dv3EmvE+7mdlkLrxlk5nZCwZPD9Y/PnVrsD7Q2GcXB2dE484BiXBAIhyQCAckwgGJcEAiHJBi1+c4pn3AOpaSVWOZzvJW0hvpY3i2N21259w6+2iwfunqrpqxq83s0tmbS+ffbg3vQmRm9orEg+6c74+Gt+ze3nVEsN7THO7ZKNw5IBEOSIQDEuGARDggEQ5IhAMS4YAUuyZYbm3FFgvLVWOV54v55fWXVCkWi+4aN0yHv99mZjXfmVvv+b3RW1dXjl8y5X9vbmZ50Z1zVWf4jbZPJcPbPnU0dLhrROHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBhhembHh+ompsX8X5wydEv2l2wPOH5oN1M7Pp/II7p78x3Dh6aHnUHa+rCzfSzMwaEvXunBct3BmsF+bXgvXBZJ+90V2lFncOSIQDEuGARDggEQ5IhAMS4YBEOCDFrgmWqk9aY33165CV56c88ljw+vbGZneNuZWcO6cz2RKsPzgVvWf5PZO7S8dHdwy663ymKbylk5nZS6evD9bfPXBhsN7Uxd7nOMQIByTCAYlwQCIckAgHJMIBKX59jkS9peqrf6zKc+8NsT/qO9Nd4/Jxf9vpQjH8AM2W9gF3fP/CpLvORzv3u3P+tveiYP3dEzcH6wPJPvO/NleLOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBDul/SibXqp+2+ysjmNKx9dOPhS8/p+HbnDXaGxIunPGl2eD9XM7tkWO72gu76u+e3bEXWf9FldRLkjOBOvzzsNLbSvhbaEU7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE+ymyUdseHysauyn4w+UjpsaUsHrB9o63TVWiwV3zkI+3Dj6eGu+Zmx23fh3NrAP+9Ja7a+z3u8uhj96eHL3lmC9p5MPAOIQIxyQCAckwgGJcEAiHJAIB6TY9Tl+rWeHza5VP9zy4v7TSsffnbg3eH1LfaO7xvDStDvnY23hN+fO3FP7ltmPzezMPbtK553pVnedsaXwgzxmZpubw32KiZW5YD2RD/eG5HUHdRV+JRAOSIQDEuGARDggEQ5IhAMS4YAUuybYj6cespGKh33eZmbfGb+ndO59E+36S7rcNfo+t8+d8576u4P1ObH9VOX4MZ2b3HUem/Hfinto+Ylg/bcGw1tkdzT7D0BF4c4BiXBAIhyQCAckwgGJcEAiHJAIB6TYNcGe03WCza5urhp7Qe8ppePZ4krw+r7P3ROsm5n1NXe4c6aW5oP1jnT0N+Aqx3fPjrrrvHfTxe6cH64OB+s/m98drPe39LprROHOAYlwQCIckAgHJMIBiXBAIhyQYtfn+PbYPTY8XO4PvMXMvj58W+m8U/QXDpj4wLPdNXre+yN3ztHt4e/QPzob3XtYzJe3qx5o8R+yef/Q1e6cLufNudZkU7BeX3dw9wDuHJAIByTCAYlwQCIckAgHJMIBiXBAil0T7F/an2aFpervq32j84LS8Z3pcJ673vNDd42OxmZ3Tncy3HjaKbauXqsY97bINvMbXGZm7anwz/u481DRakudu0YU7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE+zDxcdsfG2idP7vZvaOtYdL5yPj4X3LM52HuWvsnvffRHt0Mbwd04ldR7rjD8/ud9fxvllnZvalhqOD9bf3tgfr3d3+VlhRuHNAIhyQCAckwgGJcEAiHJAIB6TY9TnGc7M2sljdy6g8f3HvqcHrr5nfFaybmZ3dtc2dM7e2HKzfMR69zr2Tj5eOX735PHedr43f4c559vRtwfrp3ccE6/UHeQ/gzgGJcEAiHJAIByTCAYlwQCIckAgHpNg1wZbyK1VbJ5lVb6XkbW99RLrHXePq0fvdOW/fdFGwvtAV/TZbpqv8sNHXJ+5y11lY8d+KO6q9P1jfkgxvL9XeEH4YSOHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXB/r7rLFvNV2/79NWeckPq4/kngtcXbM1d49cGTnHn/MPoTcH66lohcvzh6fJbbok6f7ul5mSjO2d5LR+sZ1fGgvWefPTP6uHOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBLpu/y0Zmy02dn5jZq2dvKZ2/vO+M4PX/Nna7u8behnF3zhGtvcH6Menop7Oe3X9S6finY/4TZ22p8Mf7zMw+nD45WH/V+NXB+qbkgrtGFO4ckAgHJMIBiXBAIhyQCAckwgEpdn2OVKLeUvXVP1bl+ReHb1l/SZUvdF0QrJuZ/cHENe6c3Gr4zbpdM0M1Yx80s6tG7i6df2jw6e46H5kM/37MzF6/fGOw/rJNZwfrbX0d7hpRuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbDp5QWbzM1VjVWetzoPx7zgrf4bZGvv9t+Ke1bPjmD9xplHIse7m9pKx389cbO7zvxKzp2TrA//Y/rpXPTPckB/c6+9yl2lFncOSIQDEuGARDggEQ5IhAMS4YBEOCDFrgm2vetIG8i3Vo2d2nt06Xg6H35765SP+vuNX9F9sTvnd0dvCNZP7toSOb6lufwm3F2Tj7nrnNV3nDsnO7cvWFdbUG20rnDngEQ4IBEOSIQDEuGARDggEQ5IsetzrKzla7Zzrjx/dsu24PU/Wxl21/h0w4w752UDZwbrN+X2RI7PFfyHdyo9MLPXndOSSgfrI/NTwXp9LvWUfqYDuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbDx5TkbXapuUlWeP5qcXX9JlfunoptTlYpdR7hz2lPhxtFIbtod/6PBc9x1Prc//FCRmdmRrX3B+ube7mC9pytcV7hzQCIckAgHJMIBiXBAIhyQCAckwgEpdk2w2aVFm8rNV41Vnt+x5r9F5jk+FW4qmZldMxveSunZXdvd8W9M3h05p1JPc7s757G5kWD9a23hZlv9mr9GFO4ckAgHJMIBiXBAIhyQCAckwgEpdn2O4pP/WT92wJbW6G/IH/DlZv97ZucN3+PO8b5l/4OJ+2rG/mzdeF1dnb9OU3gdM7Pscnhnnw8lwvXuxKJ92F2lFncOSIQDEuGARDggEQ5IhAMS4YBEOCDFrgnW1thki40tVWMdFee75sLbOj1tyt9GqTPd4s55YmEiWP9o17nu+BtHr3bXeWADP+99W8Pfmzt9TzZYH2xacteIwp0DEuGARDggEQ5IhAMS4YBEOCARDkixa4LNLi/a9FL1G2+V5+/YHP4+2/kb6Pf87oy/1dIL+k8L1j+59HDN2L+Z2SdXyuPNyUZ3ncX8sjtn+6P3BusdjeGmXqr+4P4xc+eARDggEQ5IhAMS4YBEOCARDkiEA1LsmmD9zZ1mratVY5tay3t3/3Ql/Orfp+fDdTOz1bWCO6ehLvz3ZjgX/dG9yvFTO7e669w49pA754WDZwTrvXXhZltrh/+KaBTuHJAIByTCAYlwQCIckAgHJMIBKXZ9jpaGtLUlm6rGKs/vmNgVvj6Zdtd42aaz3TlXjN8VrG/viP5O3LFtm0vHD8z5b7PVJ/y/n/fl9gfrn0hsC9YTywd3D+DOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBPlR3uBXr2qrG/qVuS+n4efVjwevnV3LuGm11/m+7sLYWrBeLRXd8YQNvsx3bsdmds7K2Gqw/f/LaYH1TY7/9l73WXWc97hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE+wNSw/aSG68dP4/ZnZJrvzBvrN7jg1ef8vEI+4a35qu/Xjfeid2Hhmszxei95eqHH9Wb3jPcjOz2TW/UdZZH366bTa/GKz3NLW7a0ThzgGJcEAiHJAIByTCAYlwQCIckGLX57j+w79pqYp/b3/QzB7/l0tK571//JXg9flC+MEYM7PjWw9z59w8Xrt9dSW1O9Aj0+W30yaa5tx1ZpfDPQozs9ZUuM8xtW478PUacv4221G4c0AiHJAIByTCAYlwQCIckAgHJMIBKXZNsEs+dJNNjE6Wzj/73RfbBW//aelcvWlWqlu4bmZ298xud87z+k8J1o9KNEeO/8nmC0rHV8z5W1e3bmCbqhe2ZoL1fxwKf7OuztmqW+HOAYlwQCIckAgHJMIBiXBAIhyQCAek2DXBHpzea0MTo1Vjd088Vjqu/N5blM3pcN3M3z/dzOxLL08G611/e33N2MvtnfYP+8vjO7rDb82ZmT02N+LO+dTsdcH6YEtXsN7b1BasK9w5IBEOSIQDEuGARDggEQ5IhANS7PocdVZndXV11WMV5xO58Ftk8/noHXcqNSX9N8B6P3FbsJ5OptzxDxSPctd5SX6POyfhPKxzS6YvWF/t7rLx4Ayx7kFcg18RhAMS4YBEOCARDkiEAxLhgEQ4IMWuCTbY2mX1bdXfVzu8rbd0vKM5vGVTNjfsrjG+NOvOqW+IbnIdkFjXqDsgmagvHV+ZXnHXad5AQ+6vus4N1jP3hN94GxxcsK+7q9TizgGJcEAiHJAIByTCAYlwQCIckAgHpNg1wcYXZ21kYbpqrPJ8PBduYJ3U6T99NZKbduec2bUtWL9r5rFg3czs95b8P94vreTcOe+fuz1Y/62+8BZV7d0d7hpRuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbAX955ic8WZqrFXDJxZOi44e5v/bGl/sG5mlkr4v+2WuvC2T4W1NXf89YXwRwTNzPpbOt05nh9OPhCsD6R67bUH8ety54BEOCARDkiEAxLhgEQ4IBEOSLHrc1wxfrcNj4yVzl9nZl8ZKW/B1JoKfxPtru2b3TV23Bfdo6j004lw7+C49uh1jmkbdH/tSr2pdnfOI3Ph3s38Snirq5aV5af0Mx3AnQMS4YBEOCARDkiEAxLhgEQ4IBEOSLFrgl3Qk7HpwkDV2LP7Tiwdv2At/D2zjw2vumt0pCbcOaOLM8H6vtykO96aDDfszMz2z0f/OpXW7wW/Xr2zN3q9c73CnQMS4YBEOCARDkiEAxLhgEQ4IBEOSLFrgl0/kbXhsfKTYO8zsx+N3Vc6v0F8eO+AuWV/G6WN7Dfe1dQarI8tRDfJJis+ULj7v97prtP23Pf4P0s6/LN4TbIW5+k5hTsHJMIBiXBAIhyQCAckwgGJcECKXZ+jJdlY81Zb5XlfOrxV80u6w1s9m5l9a/Ied86z2zPB+n2NI5HjJ3WXt9c+/Pkfcdd5weDp7pwbZnYG67PLi8F6i/NGnMKdAxLhgEQ4IBEOSIQDEuGARDggEQ5IsWuCza4s2cy6pk7l+fraepc0hb/NZma2kPebQvctRze5Dnhw+gl3/OK+7e46dy7sdedMLc0H65/suzhYT3X7W0tF4c4BiXBAIhyQCAckwgGJcEAiHJAIB6TYNcFO6jrSDsu3VY2d0XNM6XgiP7f+kiqfn/Wf8qozfxuk+6b2BOundm+NHD+pq/wk2P0L0Y2ySkurK+6cb3eeF6y/afmhYL13pcc+5a5SizsHJMIBiXBAIhyQCAckwgGJcECKXZ/j3snHbWh8tGrs9vHyG1/erjyLef97Zt5OOGZm923dEayf9UT0W2gPV3yP7UU9/tt3/72Bt+9ek38wWD+qqS9Y70q2BesKdw5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCZZI1Fl9ojqzlef/0HFO8Pp/LPpvkE3mw2+QmZl9Phf+lpz6rlrl+Pem73fXaUs2uXOemBsP1ufy4S29B5r9Lb+jcOeARDggEQ5IhAMS4YBEOCARDkiEA1LsmmDHtm+27qXqvc+3dx5ROu4uFILXr9Stums8q+UYd86XZ+4O1tXWUZXjCdEoq1RoWHPn1Cfqg/XG+mSwnkoc3D9m7hyQCAckwgGJcEAiHJAIByTCAYlwQIpdE6whUW/JdU2byvNXLdwVvL492eyuMdAYbhqZmR3TMhisn9jREzn+hwNnl46/PHqbu85yIe/OOaw1eq0DvCfFkrl0sK5w54BEOCARDkiEAxLhgEQ4IBEOSLHrczw0/YQNT45Vjd0zubt0/MaB8FbPjxfD34AzM/uroWvcOaf3hB8IunzklpqxV68bP6qt311nYdX/3pzXx3jLpguD9ZY+vvGGQ4xwQCIckAgHJMIBiXBAIhyQCAek2DXBzu451qYL1Xt5X9y3vXT8+cnbg9evf1AoyuFtve6cuyYfC9abGlKR45VvnxWt6K7zN03+/uhvKYR/zw+tzQbrncWDuwdw54BEOCARDkiEAxLhgEQ4IBEOSIQDUuyaYHfO7LaRqeonwW6eeqR0fHbXtuD1140/5K6R2MAHAC/qOyFYv3Nmd+R4Q8UWTZ9NbHXX+e25//1bcdml4WC9d9l/qy4Kdw5IhAMS4YBEOCARDkiEAxLhgBS7PsfJ7UfY1FJb1djTOo8uHd8ztyd4/Yv6T3PX2G7+7j+fmAz3H45pjd75p3L8t2bvcNfpSPk/i9pK+4DF1eVgfWl1xV0jCncOSIQDEuGARDggEQ5IhAMS4YBEOCDFrgl20/jDNjQyWjX205H7Ssfv2Hxx8PpG8x/ked/Q1e6ct2+6KFj/wfLjkeMrxfI35jbS4JpcmnfnNDjfeBtZmA7WE4vRb+d5uHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbCt7YPWtlT9Dbbjug4rHf/T2M3B61cKq8G6mdkRbX3unE+O3BCsb20fiBxfXiu/XTaQ7nTXOa8tvMe6mdnzCm3B+ntydwfrfS0d7hpRuHNAIhyQCAckwgGJcEAiHJAIB6TY9Tk2N3Zac2P12JGN5W+575weCl7f1dTqrrFWXHPnJOrCf2+mVxbc8ZYm/xvy3xkL9yjMzN532PZgPZELP+DkP/4kft2DvA6/AggHJMIBiXBAIhyQCAckwgGJcECKXRPszundNjK5bnvryfL21q/ZfF7w+k/vu95d4882P92dc33TSLB+01g2cnx8sfy9tT9rO91d552F8LfkzMxO2XNPsF6fCP8dr8slg3WFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBEsl6i1VX/1jVZ7/ZDHcNOpIt7hrfG3+fnfOFxqOC9Zf1xH9xtvWivEPTN/irrOj+0h3zq658DfcjmvfHKz3dHS7a0ThzgGJcEAiHJAIByTCAYlwQCIckAgHpNg1wZqSjdaSrH6NsPJ87/x48PrCmv+qY3djeBslM7PfXrgzWD+n69jI8aPT/aXjtWLRXeecxsPcORemjwjWvzf/cLCeXPVfy4zCnQMS4YBEOCARDkiEAxLhgEQ4IMWuz/Gu5HG2mpqtGvtoakfp+BXz1wWv976HZmb22ubwNkpmZu+ZvzZYv3os+oEhNa7saRp05+xbmQrX5yeC9ULrwd0DuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbA3Td1sQ+OjpfNb7e320vGrS+dNycaIq8raG5vdNS4b+ok755kDJwXro/nZyPFMR/nhnYmVOXed8dV5d87YSvRaB1zf42wv1cU33nCIEQ5IhAMS4YBEOCARDkiEAxLhgBS7JtinO862wnJ10+db3ReXjp9Ihp/0+ljuXneNzrT/kcCbpx4J1nP55cjx+6f2lI438lTayMK0O6e3uT1Y/8PV8EcRe1eX7W/dVWpx54BEOCARDkiEAxLhgEQ4IBEOSLHrc3xobZeNFcpvcP1fM3tb4aHS+Wq+ELz+5h3+Vs4n373kzmluCD9UpHYQStWXv6fmfXvNzOzX+v23734wGv7G28zyYrC+2OzvdhSFOwckwgGJcEAiHJAIByTCAYlwQCIckGLXBGtvaLZ8srpJ1Zksf7ftI2vh7aCf+/Dj7hrphpQ7py8VfsDmntzuyPGVQr50/Kze8FtzZmY/HLvPnbO6Fm78dTSGt69ONySDdYU7ByTCAYlwQCIckAgHJMIBiXBAIhyQYtcEuzC5yRZTLVVjz02Vv3H2mvkHg9dvb97srpFK+L/t/bnJYL2vOXorpcrxn4z4Da5LNp3pzvnBZHg/9URdXbhu4bq+DhAIByTCAYlwQCIckAgHJMIBKXZ9js+O3mZDw+XtrX/H3mEfHyp/183bvnq5adVdY2gp/M00M7N0ffiBoHc0nRg5/het5a2mP5V8KHJOpa/vv8Wdk6gL/x1+16aLgvWmjvCDS3Ldg7oKvxIIByTCAYlwQCIckAgHJMIBiXBAil0T7E39Z1uurnp76/duKm9v/bXcw8Hr753b667Rl/abQrtmhoP1L6Zq36w7x8y+WCiPdyf9bbQ3tfrbVJ3YekSw/tdD14bXKPbbc+3V7jrrceeARDggEQ5IhAMS4YBEOCARDkiEA1LsmmBfmLnbRibHS+fPM7O/n7ytdN7Z2BJxVdnS6oq7xt75cXfOM/t2BOsPLO6PHJ9YmSsd/0PD8e4670wuuHMeXNgXrPe1RL99d0B3c5u7RhTuHJAIByTCAYlwQCIckAgHJMIBKXZ9jlx+xRbXfSe+8vzXu8P9h2sK4W/Qm5nNreTcOUckwv2UH8yNRY7vrRh/vD/6rbhKDfkNfO9+Mfy9e2+77twGej9RuHNAIhyQCAckwgGJcEAiHJAIByTCASl2TbBCsWiF4tq6sfL514duDV7fkPCbShtx+Uh4O6aJV2yvGdu5brznKz9113nZprPdOZtTncH6NWMPBOvLNMFwqG34zpHJZD5uZr3ZbPYPM5nMqWb2OTNrN7Nrzex12Wx2NZPJHGlmXzGzfjPLmtnLs9ns/C/g58YvwYbuHJlM5llm9n8qhr5iZm/MZrPHmVmdWekt3X82s3/OZrPHm9ntZvaeQ/iz4pfMDUcmk+k2sw+b2V89eX6UmTVls9mbn5xyuZldkslkkmZ2kZl9o3L8EP+8+CXayH+tfMbM3m1mB/YB2GxmQxX1ITM73Mx6zWw2m82urht/Sq741r/WjN1045VP9Zf5hdupxt/wkdLxLW/45fwsvyjBcGQymT82s73ZbPbHmUzmD58cTphZsWJanZmtRYzbk+NPyYt/+zU2XLFJ7U03Xmnnnvei8i9YDP+Sh+rfVjwjL8vUjO18w0ds2z9dVjrv+Ur43yLMNvZvKyNri8G6928rg4N99s1v1v6l83h3jt8zs02ZTOYuM+s2s1b7eQA2Va5tZvvNbNTMOjKZTH02my08OSf65Q78PyH4vzmy2exzstnsidls9lQze6+Z/Xc2m/0jM1vKZDLnPzntlWb2/Ww2mzez6+zngTIz+wMz+/4v6OfGL8HBNsFebmafzWQy7WZ2h5l96snx15vZlzKZzF+Y2R4ze+lT/YW70i2Wb6p+Q6un4vzjzacFr7908rpg3czs2A7/O3CXprYF6ydfeU/N2LfeYHbylSOl81R90l3nR7NZd4735Jq3N7pXVzYcjmw2e7n9/N9ALJvN3m1mZ0XMedzMnn5QPwlihw4pJMIBiXBAIhyQCAckwgGJcECK3ZNgM8sLNrVU/QhI5flfNjwYvH4j/9/K+PKsO+fz8v9a+7nN6eg9yyvHv5is/f9f1vujvN8ESzq/p2PbNgXrPW3+/upRuHNAIhyQCAckwgGJcEAiHJAIB6TY9Tl+vfckmy3OVI399sAZpeOHlkfXX1JlsXE5WDcz6075317btzgRrLenmiPHh5anSsffbOt317G8P2WtuP7R3HW/RLEQrK86dYU7ByTCAYlwQCIckAgHJMIBiXBAIhyQYtcEO9aaLLeuM7Tdyg2nyWR78PqZ1fBLx2Zm6UR4r3Azs/mVpWB9djl6nSfmyt+P+8JCeM9yM7PLBi5w53x65s5g/aGpJ4L1TY1s+4RDjHBAIhyQCAckwgGJcEAiHJAIB6TYNcG25ou2ulL95NNxFecfnwt/4K9xA1st7a1oVCn1B7lVUuUWS40N/s/yxfn73DlvbQ9vdXVlc3hfvu7uLneNKNw5IBEOSIQDEuGARDggEQ5IhANS7PoclyfGbDJRfmvsk2b2qcTwhq+fWV5w51zUd4I75965vcH6RC56d6DKrzq8pOdUd51/H73dnfPBpZuC9c+2nxusNxTDD0gp3DkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7Jtg9U7tteGKsauz2ifJW04e19oSvf/eZ7hpdl/nfJfQe9jmqPXpLp8rx78+Gt+I2M1sp+Ps+pevCb+hdthJ+YKgv32OfdlepxZ0DEuGARDggEQ5IhAMS4YBEOCARDkixa4I1JOprvtNWef74bHjv823vv95d49Ture6cVzQcGawv1kXvR/7q5u2l468u73LXmV7yn1zzvlu3uBre7z23yrZPOMQIByTCAYlwQCIckAgHJMIBKXZ9jrO7ttl0vrdq7OKe8htqy8XV4PWPLPpvxz2xFP5+m5nZJxLh790fka596Og5Zva9Qnn93mSbu87uRLhvY2bW2dgSrHcnw9+s60mzsw8OMcIBiXBAIhyQCAckwgGJcEAiHJBi1wTblmi1hfq1qrET6svbFn156q7g9Yv58IMvZhvbunosPxOs59eim3GPLo6UjieX5t11/qLvfHfOXRb+de5bGgrWl9f8t+qicOeARDggEQ5IhAMS4YBEOD9ee0MAACAASURBVCARDkiEA1LsmmBfHP2ZDQ+Xn456qb3TPjV0Q+n8iLbeqMtKBtKd7hp7F/xvvLWk0uF6Q3S9cny4MBU5p9KbX+v//Rz8q3uD9ZzT+JtL0wTDIUY4IBEOSIQDEuGARDggEQ5IhANS7JpgDYlEcNunuXwueL23LZT9/+3de5hdV3nf8ffMmTNzZs7cZ85oJNnIF+FjbAvbCNnxDZsQHi6loaSkgdZxTLGBUkiwE2hT/LTGSUpMSUkhhCcQCIGUUkKKkxhc8A3bsmVblm35PvgmCdma0dzvM+faP0DnMrN/650oarUp388/nL3W2lpn5J+WrZe91jazk7o2uGP2/tYZwf6hP7wvsv3lxcnq54vW8aLBgd+72x1TqpSD/ZlUuGDXlmp154jCygGJcEAiHJAIByTCAYlwQCIckGJX59ic6bfWVaclbenMVj9vS28K3r8ndcCdY3w5fKSTmdmJn3og2L9ciD4yur79jc1+PaXSH31Mdr350nKwf//CWLC/ueno1gBWDkiEAxLhgEQ4IBEOSIQDEuGARDggxa4INrUyt6ZIVX9981y44FMsl9w5btjwOnfM749HP8xzxIHzt65pe3lV+wn33+nO09fmn4/e5BxTNbuyGOzPrISLaHLeo7oLPxcIByTCAYlwQCIckAgHJMIBiXBAil0RLF8uWb7UeK54/XVq1W641Tpb2tw5bimNumN+Lbs92H/H82t3kZ1uZnc8v7l6PZjxd9+NL/pPpXmFsi1dg8H+bOfalxWuBysHJMIBiXBAIhyQCAckwgGJcECKXZ1juZBf8562+usPbAq/E+0dS+FTcMzM3jkbfk+cmdlE+1ywf1dl7U61r5vZHxSfrV4vFaN3xdVLJf1/BBUL74rzdvAlVzjZB8cY4YBEOCARDkiEAxLhgEQ4IBEOSLErgu0Y2GpT5cb3uF0yWDtqel9pPnj/m6cfc+f4Ye8Od8wnU+H3pt0zNRzZXn+8tXcstZnZyjoKZR5vl9+SOKLKw8oBiXBAIhyQCAckwgGJcEAiHJAIB6TYFcEOLI/b2NJEQ9sLS7Udam/seGXw/rduOMed47LR3e6Yt2w4O9i/eldeVPvGTK87z7gl3DEX9+WC/Xvnw+e9D6zjaKkorByQCAckwgGJcEAiHJAIByTCASl2dY7xxVkbXZhuaKu//m756eD9uUz4HXBmZpWI3Wqr7ZnbF+wviwd56tsf2OHXF4Zu90//2T37QrB/a2ZjsL8v7ddborByQCIckAgHJMIBiXBAIhyQCAckwgEpdkWws3tPtlcUehrazuuvPeDzxGz4wZazmvvcOcrZ090x+5bCxSl1jHZ9+zn3T0aOqffqvpPcMU9N/zjYn+lsCfa3N6XcOaKwckAiHJAIByTCAYlwQCIckAgHJMIBKXZFsMem99nI5FhD2+7J56qfu1vbg/f/zcwT7hxtyXDRyMxsdHE62H9p9ozI9td0n1z9vHv6eXeeO98R/nnMzDb+j/Cf4TtGHw/f3xR+B5zCygGJcEAiHJAIByTCAYlwQCIckAgHpNgVwYrl8pqzvOuve1oywfuvaz7NneOqyZ3umOamZLBfFbjq22/sOt+dZ+gb97hjetLhn7mrYyDYP9DV784RhZUDEuGARDggEQ5IhAMS4YBEOCDFrs7RlDBrSiRWtdWuC877zD4896A7R2dr9G61ep/JbA/239Qc/a65N/TWHgK6qeLvePPez7aeMd675NQRVR5WDkiEAxLhgEQ4IBEOSIQDEuGARDggxa4I1pPusHzbckNbX907yrpT4R1iB+bGgv1mZtv7T3XHfHwlvHPu2uS2yPZLKrXv+tXCi+48GzI97pgh5+zyZ2dfDvY3Lfs7/CLvO6q78HOBcEAiHJAIByTCAYlwQCIckAgHpNgVwT7R8iorts42tH2u9ezq52+mFoL37634had9i/5L9/KlYrD/+vwDa9puNbPrp2vtyYT/Z29b1xZ3zFhhNthftvALDStOv8LKAYlwQCIckAgHJMIBiXBAIhyQYlfn+G+VH9tE3U6xr5jZDXW1i8JyePdXtr3bncOrYZiZXd5/brD/wfxoZHuuc3P181g+XJ8wM3t4+gV3jHqf3BGX9b0q2N/T4z9QFIWVAxLhgEQ4IBEOSIQDEuGARDggEQ5IsSuCJRNNax6Sqb++IH1C8P7lNv+Io9tmh90xtyw8F+x/YWYksn33+LPVz1ds/AV/nuLT7pi5/FKw/9bF8DvehpJZ+4g7y1qsHJAIByTCAYlwQCIckAgHJMIBiXBAil0R7OvnJC01U/ta+8zs1vNr11tufSh4/1JhxZ0j6by/zcw/L7xSid5FVt/+zcN73Hnamv0jma7vDRfTPj3/aLC/L90Z7FdYOSARDkiEAxLhgEQ4IBEOSIQDUuzqHG9+cMoOj45Xr//ezLbdXbteKRaC9/e3d7lzzK4sumP628I759S77jd11t4bf2jef8fb3pNz7pgLD+4O9v9GX3h3XqbT/z2JwsoBiXBAIhyQCAckwgGJcEAiHJAIB6TYFcEmF+dsbGGmoa3+uiedCd7/4rWvcec48dPhopKZ2aLz0NBSMR/ZProwXf3clmp15znjuSfdMargdsShSvi7dlWiv6uHlQMS4YBEOCARDkiEAxLhgEQ4IBEOSLErgr13cIctJhrPDL924yXVz1+b2Ru8f9On1r57bbVCOXx+uplZ0Rmznh1v63nHm1fUMzO7ve+Vwf4rlqPPYT+iPx9+ek5h5YBEOCARDkiEAxLhgEQ4IBEOSIQDUuyKYN+cedxGJ8eq1//CzL40WTvq6eSODcH751Lhc8LNzNqa/KOWyhZd5DpidHk6sr23raP6+RXtWXeeRyf8FwBeNhE+q/2LmR3B/mZjOySOMcIBiXBAIhyQCAckwgGJcECKXZ0jmWhas8Or/vrAwtjqWxq8+JFz3Dm2fOaRo/tydTpb2iLb25tru9wem9zn/jqDmR53zNTSfLD/N1fCD0AN5gfsy+4sa7FyQCIckAgHJMIBiXBAIhyQCAckwgEpdkWw5WJ+zZFL9dfL4rilIzb/UfgdcGZmJef9bWZmf9J7UbD/T8sHItsHW2pnplcy4QeGzMxGF6MfGqp3StdQsP/KdHhHXDrDwz44xggHJMIBiXBAIhyQCAckwgGJcECKXRHshEy/pTsTDW0ndQ5WP5/aGt5F9r2RR905Ks5uNjOz3114ONg/tjgT2f7wxPPVzyd1hXfnmZldt+ESd4z3bW8Yuy/YP9SctYvsKnee1Vg5IBEOSIQDEuGARDggEQ5IhANS7OocL8wetpHpww1tP5p+ufr5mcrB4P2n9Wx259g3d9gdky8Vg/2vHzzLbV/P8db/Zfx+d8z1fRcE+9+WPTvY39XXHexXWDkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7IlgykVhTPKq/7mhJB+/vSIb7zfz3t5mZXd5/brD/vHz0b92V5YHq599Z/McfL2Vm9p3yoWD//Yd/FOwfskF7/1HMy8oBiXBAIhyQCAckwgGJcEAiHJAIB6TYFcGuHnitLVZmG9p+e0PtCKYvTYcLS8/NhwtGZma5dTwtdtPs08H+P52fXNN2v11rV47fVb0+pTt8XJPZ+p4We2T6xWD/pwYvDfa3DnDsE44xwgGJcEAiHJAIByTCAYlwQIpdneOh8qRNlaaq1+8ws3tLtfe6ndDeH7z/sYl97hwf7N7qjrl+eU+wvz3V6rafmh6MHFPv++s4iaiztT3Y/7eV0WB/byVv57mzrMXKAYlwQCIckAgHJMIBiXBAIhyQCAek2BXBRvLTNr7S+CDNS3XXB+bGVt/SQBWn6n108l53zEqxEOxPJaN/6wp1u+n2rYy786hfp97H+s4P9r+5MhfsryT835MorByQCAckwgGJcEAiHJAIByTCAYlwQIpdEeyyti22kOltaHtLpvbk1r/pDT9dtW3f4+4c6VSLO8YrThXE2ej17c/Wndmu/PrGcIHLzOyT47uC/f+xsBLs35gYtJvcWdZi5YBEOCARDkiEAxLhgEQ4IBEOSLGrc3x3dtgOT9cekrnCzL41XatdfGEp/GBLTzrjzpFN++88e2lhIth/UX8usv2ywTOrn3dNho+dNjPrN7/mcknf6cH+dCIZ7O/O9rhzRGHlgEQ4IBEOSIQDEuGARDggEQ5IhANS7IpgH2870wrtjcdbf7q99r619+fvC94/vbzgztHX2umOWSmFd7w9NnfAbe9saXPn+dtFv1A2nQ//TB/tfm2wP23+zxuFlQMS4YBEOCARDkiEAxLhgEQ4IBEOSLErgn1s+kEbmawd7XSvXWvvnbynen1e/yuD9w/P+7vM/mnGP/u87Iy5fWl/ZPtQuvbUVTblF5/um/CLYHmnIPfJ8gPB/g2tWbvQnWUtVg5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCHZhf86mSxsa2n4xe1b184sr4bPPb+jc7s7xv0r+meT3jD8T7H9b9uzI9lxr7Viqm0YfdudZj3RzeMvkUiEf7F92+hVWDkiEAxLhgEQ4IBEOSIQDEuGAFLs6x2mJTltoqjS0ndlUO6bptpnw8dV/ne0N9puZPb942B1TrHtXW5QH5l9027tb2915WpMpd0ypUg72f7jr3GB/ur/LnSMKKwckwgGJcEAiHJAIByTCAYlwQCIckGJXBHu2MmfT5ZmGtqfqrvvbwrvI7jr8lDvH9I1vdcdk/8MPgv2ji9Nue7lSiRxTz/t5zMwmnPPer1+6O9g/lBi0S+0qd57VWDkgEQ5IhAMS4YBEOCARDkiEAxLhgBS7IthDMy/a6FTjrrZ7p2pHI83nl4L397f5Tz1tuu52d0wm1Rrs70hFn2u+MdNX/fy55vCL+8zMfmX6XneM90TZH3WEzz5P9fIkGI4xwgGJcEAiHJAIByTCAYlwQIpdnSNfLNhKsfE45/pr75SbqeV5d47KOh7CSSQSwf602KlW3/79Nn+ev0i+zh1zTyp8Ms9VI87DPqlB+w4P++BYIhyQCAckwgGJcEAiHJAIByTCASl2RbBfGTjH5iuNO96u2HBe9fP1rwgfb33BU1PuHN3N/nFMT08fDPbvn4v+HvXtX5j23zf3tVTaHdPclAz2l5wjqsqVcL/CygGJcEAiHJAIByTCAYlwQCIckAgHpNgVwb4zsddGDtcKSe8zs68f3l29/ouRcEGnyXmCy8zs1zac4Y55tulQsP/mrh1u+5un73fnGVjHDr2WpvA/pi2ZbLC/v6cv2K+wckAiHJAIByTCAYlwQCIckAgHpNjVObpa2205nWlo66m7TibCeT4zc4I7x59P7nHHePWSq0vPrWn71qr29eysu+e0HnfMuU+OBPuXi+EdcfMZHvbBMUY4IBEOSIQDEuGARDggEQ5IhANS7IpgzYkmSyUad3jVX+fLxeD9t48+7s6xtWejO+arzeFi2nVNi5Htp6QHq59fmp9w57lo2N+hd0J7f7D/8cn9wf7lVcdorRcrByTCAYlwQCIckAgHJMIBiXBAIhyQYlcEG1mYskNz4w1tB+uuk84RSB/cdJE7x5+8dI875rLkYXdMlJ3jz1Q/b+32i23PTIWPlzIzO9Q0Geyf/Ej07rsjCu3d9oI7y1qsHJAIByTCAYlwQCIckAgHJMIBiXBAil0RrCfdYfn25Ya2/vba0Ug7uk4J3n/T7NPuHB/YfLE7prcS/q3ZW56JbP+l7FnVz7eMPurOsx7tqdZgf98f7w72Dw0N2nfe9g+fl5UDEuGARDggEQ5IhAMS4YBEOCDFrs5RKJcsX2rc1VZ/vXNyOHh/tq3bneOpor/LbMXZWferyc2R7b9ovdXPD7X732V8cdYds6EtfDTUWV3hI7J7+3qD/QorByTCAYlwQCIckAgHJMIBiXBAIhyQYlcESyYSa843r7+eyS/8o+fYP+vvZtvcGT5q6bb2tjVtrzOz26y2O61Y9s8c35Dxzz6/MbE12H9Tc/js886k/x65KKwckAgHJMIBiXBAIhyQCAckwgGJcECKXRFsKN1rqfbGF+fVn/3d29oRvP+s9ugntOp9aGVtAWu1dy48HOw/2Bz9NNnBfK19obDizrP6qbcoVxTCO9rMqQtuSGXtCneWtVg5IBEOSIQDEuGARDggEQ5IhANS7Ooc72neYoXmxh1av9V8au3z/H3B+/fNjrpz3HjyNnfMwtRysP+Jlej3qj1R97619pa0O88ruza5Y56ZCR+BveK8w20p79dborByQCIckAgHJMIBiXBAIhyQCAckwgEpdkWw35t92Eanx6rXd9g1ds30rup1e3P4qOcdveHdYWZmrxuJLmDV846UfkP/GZHtb9+4vfr5Vwv+TrN/PbXTHTPQHv512pPh7zrQFd69p7ByQCIckAgHJMIBiXBAIhyQCAckwgEpdkWwcqVspUq5oa3+OpsOF4TSCf9HWimFn5wyM2ttTgX7b514ck3bx1e1f/bMk9x5epYy7phFZ+fc1PJ8sL+0eHRrACsHJMIBiXBAIhyQCAckwgGJcECKXZ2jWC6tORa6/np46qXg/edtPMGdY0fPqe6YnePPBPvfvuHcyPZ/MvDq6uf3HwofO21m1t/a6Y758fx4sP/07vDP3N/Z584RhZUDEuGARDggEQ5IhAMS4YBEOCARDkixK4K1p1otk2o8Lqn+urMlfDT1t8bCx1KbmV0xuMMdM9l7YrD/mkLCbX/Hkr+z7vDCtDumJRl+8Mg7FmqojWOfcIwRDkiEAxLhgEQ4IBEOSIQDEuGAFLsiWKFcsny58b1n9dejTtHoBz1+gevK+WfdMVMr4V1k72lb+5TX183sPcXak2pTS+Ffw8zsI5succd8ZWJPsD+ZCP8ZTztFNIWVAxLhgEQ4IBEOSIQDEuGARDggEQ5IsSuCVSoVq1Qqa9qO2JDpCd7/oUp4u6SZ2Vntm90x+5sng/2rC3VRXp890x0zWvG3TJbK5WB/OVEJ9hcq4fsVVg5IhAMS4YBEOCARDkiEAxLhgBS7OsepnRutZ6nxnWX1xxo9MRXeRTZq/g6ywQH/3WtPTh4I9v/y0Gsi21+VHqp+vunQQ+48w6ed5Y65aTxcxygWwzWXfNE/zjsKKwckwgGJcEAiHJAIByTCAYlwQCIckGJXBHthfsRGZsca2oZnaw/weA++JBLRxzHVGyvOuWN60uF3r72lFF1Iq2//QUs6cky9780OumP606PB/oSFf+Zse7c7RxRWDkiEAxLhgEQ4IBEOSIQDEuGARDggxa4Idl7PqTad729oe13f6dXP/3vk0fAvEH5oyszMDq/MuGMml8KFsg/n71nTttOusQ9P1NrzJf8JrKv/8jJ3zO+85a5gf9nZ0VbM+IXBKKwckAgHJMIBiXBAIhyQCAckwgEpdnWOJxdessPz4w1tj8zXdrl5D/N4D76s16BzgtANmegdb5/trx1Xfc3Uve48X7nih+6Yrd0bg/0LxeVgf9b5WRRWDkiEAxLhgEQ4IBEOSIQDEuGARDggxa4I9oepnJVbGos+f96yrfr5d1vCxz7tXxgL9puZdTa3uWNSTclg/78dX/sAzi67pqH9goGcO88XCv775g7Mh3+mZue7pgvh3XsKKwckwgGJcEAiHJAIByTCAYlwQCIckGJXBLuttWzz6VL1+moz+7u665n5xeD9i4UVd45st3/2+chK+Ax1tcusvv1yG4ocU+8/5Q+5Y1qS4X9M6eaWYH+b06+wckAiHJAIByTCAYlwQCIckAgHpNjVOR5eOWQTy7X3yF9tZruWD1avx5dmg/ef1fsKd44nZsLvbzMzW3Hei9aWanXbf3/pMXeeb7f6DwTd2Bl+3/1dk88E+5eK4fsVVg5IhAMS4YBEOCARDkiEAxLhgEQ4IMWuCPbS0qSNrtq1Vr+L7a0Drw7ev78YfkjHzKzkHAdt5j9gs7Ur+iim07o2VT9vaPYfKvpUwj8C+6WCfxz3/w2sHJAIByTCAYlwQCIckAgHJMIBiXBAil0RzHNBpTPYf2Ey3G9m9pnWKXdMqin8W3NwccJt31c+7M6znt1o+XIx2F8sl4L9JadfYeWARDggEQ5IhAMS4YBEOCARDkiEA1LsimDndG+x6ZXuhrbze7dWP39s9IfB+8/s87dDLhbXcTRUOh3sV0+K1bero6Hq/VXLq9wxV5WeD/ZPLc0H+/PFcBFNYeWARDggEQ5IhAMS4YBEOCARDkixq3PcNzZsh0YbH5K5ffTx6uf+tvAusmyz/7DP40v73DFv7D0j2P93Y49Gtk8t12oOvekOd57COv58zuQXgv1bugaD/dnOfneOKKwckAgHJMIBiXBAIhyQCAckwgGJcECKXRHs1O4h615p3AV2eu8J1c8/mn45eP/D5RfdObb1n+SO+e54+Nzyy7Ovddu/MnK/O889Q0l3zCbrC/YnEolgf3sy+px2DysHJMIBiXBAIhyQCAckwgGJcEAiHJBiVwQ7tDRtowuNRyodrLv+0MaLgvd/eWy3O8e1TSe7Yz7SNBbs/8vRB9e0Xb2q/beHLnbn+a9ju9wxy84L/JoS4T/jQ63+Dr/IX/eo7sLPBcIBiXBAIhyQCAckwgGJcECKXZ3jk5mzrdg529D2xc7zqp+vnQk/hDOfX3Ln+IK95I6ZWwn/OhWrRLYX6o6i/tyY/7DPv89e6I65vRQ+Jvvu0SeD/eUKx1vjGCMckAgHJMIBiXBAIhyQCAckwgEpdkWwfzf/iI3M1h60+aGZvXe2VkwqlcNHRnu7v8zMTmjudscs9p4Y7P9Gx9pdaMtmtvfEc6rXvzwT3p1nZvY/l8NHV5uZPTdzKNjf6rwnriWZcueIwsoBiXBAIhyQCAckwgGJcEAiHJAIB6TYFcGWC3lbKjTu0Kq/Tli4yOXt/jIzG2zyj0F6RjzpdcTmHYtr2p5f1T7+97Nrxqy2UFh2x3jnvc+shM9GX8+75qKwckAiHJAIByTCAYlwQCIckAgHpNjVOSpWsUqlscZQf93fHv47/4EHv+jO0bHt3e6Yvrbwu+Jyt8ysafveR81yt4xUr0tHWV9YbXxx7Vz1vNpOuRyu2chf96juws8FwgGJcEAiHJAIByTCAYlwQCIckGJXBDuzd4sNFRoLUOcOnFr9vHci/A63zLZ3uXNclH2VO2bPVHgnWrEcfZRS/YM3V244353nryf3umPO79ka7H9NU0+wvz0bLhwqrByQCAckwgGJcEAiHJAIByTCAYlwQIpdEeyNqU221NLR0Pb2lldUPz+V/HHw/tVPkUX5XCrtjvFOJE8lo3/r6tsfK4y786xnx9vdE08H++91ngTb0Jy1N7mzrMXKAYlwQCIckAgHJMIBiXBAIhyQCAek2BXBPj++20ZGa2efv9Wutf88urN63dvWEXVb1c0dpwb7zcy2v/ywO+bSwTOC/Y/M7Itsb2mq/ZaO5qfdeQacI53M/PPcu1Pt4Tnae905orByQCIckAgHJMIBiXBAIhyQCAek2NU5su3dZpliQ9tQpvb39F/vCNcfLl96zp3jnw1td8c8vTIS7O9pzfyD2pX+lvDxUmZmz86G3/E2uhCup8ymi8F+hZUDEuGARDggEQ5IhAMS4YBEOCARDkixK4K1NqUsnUw1tNVff2FqT/D+yaU5d46OAX/Hm+ftHTm3/W9mn3J/ndPaB90x8xl/V1xItoOHfXCMEQ5IhAMS4YBEOCARDkiEAxLhgBS7Itj+uTEbmTnc0Pb8TO2prBM7B4L33+w8KWZm9vqpR9wxF/VHF7mO+NLYg2va/tWq9rbmFneemw8/6o5R56wf8d6NFwT7O9q63TmisHJAIhyQCAckwgGJcEAiHJAIB6TY1Tn2vucsSy3MVK+HzWzyN2s71Po+G37Y59LEhDtHZ0ubO2bn+DPumCjFUq0mMVWYd8d/cNNF7pi/mgjXZb46+kCwfyiRtXe7s6zFygGJcEAiHJAIByTCAYlwQCIckAgHpNgVwXZ87WkbGakdb33Xu8yGPr+3et3anIq6rWqpsOLOsZ5dcaf1bg72Pzv9cmR7oVw7YunGode78xywgjvmzX1nhr9LPlz46+tmxxuOMcIBiXBAIhyQCAckwgGJcEAiHJBiVwTLl0qWLxVXtdWum5z3nX29/zJ3jsvH73THTKz4hTLPDyrj7pjHpg+4Y0qVcrD/E52vDfa3JP33yEVh5YBEOCARDkiEAxLhgEQ4IBEOSLGrcyQTCUsmmla11a69v/NfV/CPlG5K+H8melrC72pTp+10tdbeK3/H6BPuPOuRbAp/3yu/fXGwP9/UYv6b79Zi5YBEOCARDkiEAxLhgEQ4IBEOSIQDUuyKYK/pO8Wmin0NbRcM1I6avn/iR8H7Dy1MuXO8b9OF7pgvHdoV7L80G32M9vbuU6qfH5vzH+TxHl4yMytXKsH+P7jq7mB/e7bL3vRnZ7vzrMbKAYlwQCIckAgHJMIBiXBAIhyQCAek2BXBHpvZbyNTYw1tu6dqzzH9Qv9pwfuHkuEnuMzMhovT7pit3RuD/Wc097jtO1f889N72zrcMSvFfLD/z6bC58FvaBmwN7mzrMXKAYlwQCIckAgHJMIBiXBAIhyQYlfnWCqs2GJ+uaGt/vq95Wzw/gebi8F+M7Nvjz3kjjlv4JXB/uHSrNueWMeDPLMri+4Y71327W2dwf7m5NH9Y2blgEQ4IBEOSIQDEuGARDggEQ5IhANS7IpgJ3ZmLb3YmNmTu4eqn68cvyt4/3oKT+8a2uGOGS6E35v2yOw+t/0z/eHjmMzMPjwW/nnMzLrrjpKKcnom/D663vboB5M8rByQCAckwgGJcEAiHJAIByTCAYlwQIpdEaw/1WmJlsanuQZbuqufL9x4UvD+vcuH3DnuX9zvjjmjbVOw/4XKSGR7ue5s9k8sPOLO07KOp7RSzpiF8kqwv7Uc3jGnsHJAIhyQCAckwgGJcEAiHJAIByTCASl2/eATSAAAAsNJREFURbD2ppStNLU0tGXqrr8//XTw/rn8kjvHYHu3O2a8tBDsVy8RrG9/V/c2d57/PrXXHXNi20Cwv1AObwH1tlMqrByQCAckwgGJcEAiHJAIByTCASl2dY5d4z+yQ4cPN7TdefiJ6ufBTHj3VsV5H5qZ2cTynDvmhHR/sP+UzAa3/Zx80p3nllb/eOuHx58P9r9r43nB/q4Wv64ThZUDEuGARDggEQ5IhAMS4YBEOCDFqc6RNDMbzK6tL2wcGqx+7m/vCv4irctt7kTNSb/+0NfVG+wv1W1earhvsHZfysLf1cxsoBiup5iZLaTDz2N0ZcN1jI7+6vHX/g9eJ7GeotH/C3v27LnYzO453t/j/3OXbN++fed6B8dp5dhtZpeY2SEzO7pHl6AkzWyj/eT3eN1is3IgfvgPUkiEAxLhgEQ4IBEOSIQDEuGARDggxalC2iCXy/1LM7vOzFJm9sfDw8OfP85fScrlcnea2aCZFX7a9P7h4eEHjuNXOiZiWSHN5XKbzWynmW03sxUzu8/M3j08PPzUcf1iEXK5XMLMDprZluHhYf+ltj9D4vqvlV8yszuGh4cnh4eHF8zs22b2zuP8nZTcT//3B7lcbm8ul/vQcf02x1Bcw7HJfvJ/wB1xyMxOOE7fxdNrZreb2TvM7A1m9oFcLvfG4/uVjo24/jdHk5nV//suYWbRD1AcZ8PDw7vMbNeR61wu92Uze6uZ3XrcvtQxEteV46D95P9iPmLIzF4+Tt8lKJfLXZzL5d5Q15Sw2n+Y/kyL68pxm5ldn8vlsma2YGb/3Mzed3y/ktRjZjfkcrkL7Sd/s/oNM/vA8f1Kx0YsV47h4eGXzOzjZnanmT1qZt8YHh5+8Ph+q2jDw8M3m9l3zewRM9tjZl/56b9qfubF8q+yiIdYrhyIB8IBiXBAIhyQCAckwgGJcED6Py7YuJh/m/lpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x5760 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=plt.subplots(1,figsize=(10,80))\n",
    "ax.imshow(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y.sum(axis=1) # each sample's 10 classes probs should sum to 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [ ] TODO: softmax more complex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "def sigmoid(x):\n",
    "    return 1/(1+np.exp(-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## J_cost\n",
    "- [ ] TODO: softmax have different lambda0 term\n",
    "\n",
    "Calculate the loss function without the non-differentiable part"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "W = W_init\n",
    "beta = beta_init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to implement softmax, change lambda0 term with softmax loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "def softmaxEquation(scores):\n",
    "    \"\"\"\n",
    "    It calculates a softmax probability\n",
    "    :param scores: A matrix(wt * input sample)\n",
    "    :return: softmax probability\n",
    "    \"\"\"\n",
    "    scores -= np.max(scores)\n",
    "    prob = (np.exp(scores).T / np.sum(np.exp(scores), axis=1)).T\n",
    "    return prob\n",
    "\n",
    "def computeLoss(x, yMatrix,wt,regStrength):\n",
    "    \"\"\"\n",
    "    It calculates a cross-entropy loss and gradient to update the weights.\n",
    "    :param x: An input sample\n",
    "    :param yMatrix: Label as one-hot encoding\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    numOfSamples = x.shape[0]\n",
    "    scores = np.dot(x, wt)\n",
    "    prob = softmaxEquation(scores)\n",
    "\n",
    "    loss = -np.log(np.max(prob)) * yMatrix\n",
    "    totalLoss = np.sum(loss) / numOfSamples\n",
    "    grad = ((-1 / numOfSamples) * np.dot(x.T, (yMatrix - prob))) \n",
    "    return totalLoss, grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- logistic regression\n",
    "$$\n",
    "\\begin{aligned} l ( \\beta ) & = \\sum _ { i = 1 } ^ { N } \\left\\{ y _ { i } \\log p \\left( x _ { i } ; \\beta \\right) + \\left( 1 - y _ { i } \\right) \\log \\left( 1 - p \\left( x _ { i } ; \\beta \\right) \\right) \\right\\} \\\\ & = \\sum _ { i = 1 } ^ { N } \\left\\{ y _ { i } x _ { i } ^ { T } \\beta - \\log \\left( 1 + \\exp \\left( x _ { i } ^ { T } \\beta \\right) \\right) \\right\\} \\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- softmax regression\n",
    "$$\n",
    "\\sum _ { k = 1 } ^ { K } y _ { k } ^ { ( i ) } \\log \\left( \\hat { p } _ { k } ^ { ( i ) } \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy_loss(W,X,Y,beta):\n",
    "    '''\n",
    "    part of J_loss first term\n",
    "    '''\n",
    "    n = X.shape[0]\n",
    "    Y_cap = sigmoid(X@beta)\n",
    "    #weight_term = (W*W).reshape(-1,1)\n",
    "    #print (weight_term.shape,(Y.reshape(-1,1)*np.log(Y_cap)).shape)\n",
    "    loss = np.sum((W*W)*(Y*np.log(Y_cap)))\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-382.6933796528916"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cross_entropy_loss(W,X,Y,beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lambda0*sum((W*W)*(np.log(1+np.exp(X@beta))-Y*(X@beta)))\n",
    "def J_cost(W,beta,X,Y,lambda0, lambda1, lambda2, lambda3, lambda5):\n",
    "    return lambda0*cross_entropy_loss(W,X,Y,beta) \\\n",
    "         +lambda1*sum(balance_cost(W,X)) \\\n",
    "         +lambda2*((W*W).T@(W*W)) \\\n",
    "         +lambda3*sum(beta**2) \\\n",
    "         +lambda5*(sum(W*W)-1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[24332.18264953]])"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J_cost(W,beta,X,Y,lambda0, lambda1, lambda2, lambda3, lambda5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## balance cost\n",
    "- [ ] TODO: it seems do not need to change, means than we can adapt this function into any gradient descent problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def balance_cost(W=None,X=None,*args,**kwargs):\n",
    "    m = X.shape[1]\n",
    "    f_x=np.zeros([m,1])\n",
    "    for i in np.arange(0,m):\n",
    "        X_sub=np.copy(X)\n",
    "        X_sub[:,i]=0\n",
    "        I=(X[:,i] > 0).astype('double')+eps\n",
    "        loss=( np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ) / (np.dot((np.multiply(W,W)).T,I.reshape(-1,1)))\\\n",
    "            -(np.dot(X_sub.T,(np.multiply((np.multiply(W,W)),(1 - I.reshape(-1,1)))))) / (np.dot((np.multiply(W,W)).T,(1 - I.reshape(-1,1))))\n",
    "        #print (loss.shape)\n",
    "        f_x[i]=np.dot(loss.T,loss)\n",
    "    return f_x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.38730847],\n",
       "       [0.41165626],\n",
       "       [0.39316231],\n",
       "       [0.2727617 ],\n",
       "       [0.31718681],\n",
       "       [0.44008927],\n",
       "       [0.37368856],\n",
       "       [0.21282965],\n",
       "       [0.19291022],\n",
       "       [0.16060995],\n",
       "       [0.44438877],\n",
       "       [0.15128009],\n",
       "       [0.29727955],\n",
       "       [0.33768393],\n",
       "       [0.20408237],\n",
       "       [0.35794777],\n",
       "       [0.30523414],\n",
       "       [0.38001066],\n",
       "       [0.2351875 ],\n",
       "       [0.28573664]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "balance_cost(W,X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## balance_grad.m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def balance_grad(W=None,X=None,*args,**kwargs):\n",
    "    n,m=X.shape\n",
    "\n",
    "    g_w=np.zeros([n,m])\n",
    "    for i in range(0,m):\n",
    "        X_sub = np.copy(X)\n",
    "        X_sub[:,i] = 0 # the ith column is treatment\n",
    "        I = (X[:,i]>0).reshape(-1,1).astype('double')+eps\n",
    "        J1 = (X_sub.T@((W*W)*I.reshape(-1,1)))/((W*W).T@(I.reshape(-1,1))) \\\n",
    "            -(X_sub.T@((W*W)*(1-I).reshape(-1,1)))/((W*W).T@(1-I).reshape(-1,1))\n",
    "        dJ1W = 2*(X_sub.T*((W*I)@np.ones([1,m])).T*((W*W).T@I) \\\n",
    "                  -(X_sub.T@(((W*W)*I)@(W*I).T)))/((W*W).T@I)**2 \\\n",
    "                  -2*(X_sub.T*((W*(1-I))@np.ones([1,m])).T*((W*W).T@(1-I)) \\\n",
    "                  -((X_sub.T@( (W*W) * (1-I) )) @  (W*(1-I) ).T ))/((W*W).T@(1-I))**2\n",
    "        g_w[:,i] = (2 * dJ1W.T @ J1).ravel()\n",
    "\n",
    "    return g_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 20)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "balance_grad(W,X).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## prox_l1.m\n",
    "\n",
    "- [ ] TODO: it is also universal?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def prox_l1(v=None,lambda_=None,*args,**kwargs):\n",
    "    x=np.fmax(0,v - lambda_) - np.fmax(0,- v - lambda_)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## main function\n",
    "\n",
    "- [ ] TODO: careful change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_CE(lambda0,W,X,Y,beta):\n",
    "    '''\n",
    "    part of J_loss's first term grad\n",
    "    '''\n",
    "    n = X.shape[0]\n",
    "    Y_cap = sigmoid(X@beta)\n",
    "    #original sigmoid(X@beta)-Y)*(W*W)).T@X\n",
    "    grad = (lambda0*((Y_cap-Y)*(W*W)).T@X).T\n",
    "    return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 10)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_CE(lambda0,W,X,Y,beta).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span class=\"burk\">TODO</span>: \n",
    "- [x] grad_beta \n",
    "- [x] grad_W\n",
    "- [x] 要根据grad_CE，对比cross_entropy_loss和原loss改\n",
    "- [ ] 改J_loss迭代时的判别条件，grad_beta.T@(z-beta) 现在不是scalar 是class_num*class_num\n",
    "prox_l1在控制什么？是L_1 norm of bata？\n",
    "- [ ] J_cost(z, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5) z改了之后如何修改\n",
    "- [ ] 2*lambda0*Y*np.log(sigmoid(X@beta))*W 导致grad_W维度变了，所以上面的z因为prox变了，维度不匹配，注意从原loss推导保证维度不变，，"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-363.45200780475614"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = X.shape[0]\n",
    "def sigmoid(x):\n",
    "    exp_x = np.exp(x)\n",
    "    sum_exp_x = np.sum(exp_x,axis=1)\n",
    "    sum_exp_x= sum_exp_x.reshape((x.shape[0],1))\n",
    "    sigmoid = exp_x / sum_exp_x\n",
    "    return sigmoid\n",
    "\n",
    "lambda0*np.sum((W*W)*Y*np.log(sigmoid(X@beta))) #NEW CE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "outputs": [],
   "source": [
    "#Original loss: lambda0*np.sum((W*W)*(np.log(1+np.exp(X@beta))-Y*(X@beta))) \n",
    "#Original   grad_beta =     lambda0*(((sigmoid(X@beta)-Y)*(W*W)).T@X).T  \n",
    "grad_beta = grad_CE(lambda0,W,X,Y,beta)\n",
    "\n",
    "#Original loss: lambda0*np.sum((W*W)*(np.log(1+np.exp(X@beta))-Y*(X@beta))) \n",
    "#Original  grad_W = 2*lambda0*(np.log(1+np.exp(X@beta))-Y*(X@beta))*W\n",
    "grad_W =   2*lambda0*Y*np.log(sigmoid(X@beta))*W\n",
    "grad_beta.shape,grad_W.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import KFold, train_test_split, GridSearchCV, cross_val_score\n",
    "from sklearn.metrics import accuracy_score, roc_auc_score, f1_score, recall_score, precision_score, \\\n",
    "    roc_curve, precision_recall_curve, average_precision_score, matthews_corrcoef, confusion_matrix\n",
    "\n",
    "def report_metrics(y_test, y_pred):\n",
    "    scorers = {'accuracy': accuracy_score(y_test, y_pred),\n",
    "           'recall': recall_score(y_test, y_pred,average=None),\n",
    "           'precision': precision_score(y_test, y_pred,average=None),\n",
    "           'f1': f1_score(y_test, y_pred,average=None),\n",
    "           'mcc': matthews_corrcoef(y_test, y_pred)\n",
    "    }\n",
    "    for metric in scorers.keys():\n",
    "        print('{} = {}'.format(metric, scorers[metric]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_num = 20\n",
    "sample_num = 500\n",
    "class_num = 10\n",
    "eps=10e-3\n",
    "X = 2*np.round(np.random.rand(sample_num, feature_num))-1; # 1000 samples and 20 features\n",
    "beta_true = np.random.random([feature_num, class_num]);\n",
    "Y = np.random.randint(0,9,500).reshape(-1,1) #sigmoid(np.dot(X,beta_true))\n",
    "lambda0 = 1; #Logistic loss\n",
    "lambda1 = 0.1; #Balancing loss\n",
    "lambda2 = 1; #L_2 norm of sample weight\n",
    "lambda3 = 0; #L_2 norm of beta\n",
    "lambda4 = 0.001; #L_1 norm of bata\n",
    "lambda5 = 1; #Normalization of sample weight\n",
    "MAXITER = 1000;\n",
    "ABSTOL = 1e-3;\n",
    "W_init = np.random.rand(sample_num, 1);\n",
    "beta_init = 0.5*np.ones([feature_num, class_num]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of training samples: 375, test samples: 125\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X,Y.astype('int'))\n",
    "print('number of training samples: {}, test samples: {}'.format(X_train.shape[0], X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy = 0.096\n",
      "recall = [0.15789474 0.25       0.         0.0952381  0.11111111 0.07692308\n",
      " 0.11764706 0.16666667 0.        ]\n",
      "precision = [0.2        0.06666667 0.         0.15384615 0.05       0.09090909\n",
      " 0.4        0.04545455 0.        ]\n",
      "f1 = [0.17647059 0.10526316 0.         0.11764706 0.06896552 0.08333333\n",
      " 0.18181818 0.07142857 0.        ]\n",
      "mcc = 0.0013368562728601418\n"
     ]
    }
   ],
   "source": [
    "model = LogisticRegression()\n",
    "model.fit(X_train,y_train)\n",
    "y_pred = model.predict(X_test)\n",
    "report_metrics(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy_loss(W,X,Y,beta):\n",
    "    '''\n",
    "    part of J_loss first term\n",
    "    '''\n",
    "    n = X.shape[0]\n",
    "    Y_cap = sigmoid(X@beta)\n",
    "    #weight_term = (W*W).reshape(-1,1)\n",
    "    #print (weight_term.shape,(Y.reshape(-1,1)*np.log(Y_cap)).shape)\n",
    "    loss = np.sum((W*W)*(Y*np.log(Y_cap)))\n",
    "    return loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mainFunc(X, Y, \\\n",
    "    lambda0, lambda1, lambda2, lambda3, lambda4, lambda5,\\\n",
    "    MAXITER, ABSTOL, W_init, beta_init,paras_save_path=None):\n",
    "\n",
    "    n,m = X.shape\n",
    "    W = W_init\n",
    "    W_prev = np.copy(W)\n",
    "    beta = beta_init\n",
    "    beta_prev = np.copy(beta)\n",
    "\n",
    "    parameter_iter = 0.5\n",
    "    J_loss = np.ones([MAXITER, 1])*(-1)\n",
    "\n",
    "    lambda_W = 1\n",
    "    lambda_beta = 1\n",
    "\n",
    "    W_All = np.zeros([n, MAXITER])\n",
    "    beta_All = np.zeros([m, MAXITER])\n",
    "\n",
    "    \n",
    "    # Optimization with gradient descent\n",
    "    for iter in tqdm(range(1,MAXITER+1)):\n",
    "        \n",
    "        # Update beta\n",
    "        y = np.copy(beta)\n",
    "        beta = beta + (iter/(iter+3))*(beta-beta_prev) # fast proximal gradient\n",
    "        f_base = J_cost(W, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\n",
    "        grad_beta = grad_CE(lambda0,W,X,Y,beta) \\\n",
    "                   +2*lambda3*beta\n",
    "        while 1:\n",
    "            z = prox_l1(beta - lambda_beta*grad_beta, lambda_beta*lambda4)\n",
    "            print ('beta',J_cost(W, z, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5),f_base , np.sum(grad_beta.T@(z-beta)),\\\n",
    "             (1/(2*lambda_beta))*sum((z-beta)**2))\n",
    "            #grad_beta.T@(z-beta)\n",
    "            if J_cost(W, z, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5) \\\n",
    "               <= f_base + np.sum(grad_beta.T@(z-beta)) \\\n",
    "               + (1/(2*lambda_beta))*sum((z-beta)**2):\n",
    "                break\n",
    "            lambda_beta = parameter_iter*lambda_beta\n",
    "            if lambda_beta<eps*eps*eps:\n",
    "                break\n",
    "\n",
    "        beta_prev = y\n",
    "        beta = z\n",
    "        # Update W\n",
    "        y = np.copy(W)\n",
    "        W = W+(iter/(iter+3))*(W-W_prev)\n",
    "        f_base = J_cost(W, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\n",
    "        print ('test1')\n",
    "        grad_W = 2*lambda0*Y*np.log(sigmoid(X@beta))*W \\\n",
    "                +lambda1*balance_grad(W, X)@np.ones([m,1]) \\\n",
    "                +4*lambda2*W*W*W \\\n",
    "                +4*lambda5*(sum(W*W)-1)*W\n",
    "        print ('test2')\n",
    "        while 1:\n",
    "            z = prox_l1(W-lambda_W*grad_W, 0)\n",
    "            print (z.shape)\n",
    "            print ('W',J_cost(z, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5),\\\n",
    "                   f_base , np.sum(grad_W.T@(z-W)),\\\n",
    "                    (1/(2*lambda_W))*sum((z-W)**2))\n",
    "            #grad_W.T@(z-W)\n",
    "            if J_cost(z, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\\\n",
    "                    <= f_base + np.sum(grad_W.T@(z-W))+(1/(2*lambda_W+10e-100))*sum((z-W)**2):\n",
    "                break \n",
    "            lambda_W = parameter_iter*lambda_W\n",
    "            if lambda_W<eps*eps*eps:\n",
    "                break\n",
    "\n",
    "        W_prev = y\n",
    "        W = z\n",
    "        \n",
    "        W_All[:,iter-1] = W.ravel()\n",
    "        beta_All[:,iter-1] = beta.ravel()\n",
    "\n",
    "        J_loss[iter-1] = J_cost(W, beta, X, Y,\\\n",
    "                              lambda0, lambda1, lambda2, lambda3, lambda5)\\\n",
    "                     + lambda4*sum(abs(beta))\n",
    "        print (J_loss[iter-1] , J_loss[iter-2])\n",
    "        if (paras_save_path is not None) & (iter%10==1):\n",
    "            #'output/models/crlr/somedir'\n",
    "            if not os.path.exists(paras_save_path):\n",
    "                os.makedirs(paras_save_path)\n",
    "            np.savetxt(paras_save_path+'/beta.txt',beta)\n",
    "            np.savetxt(paras_save_path+'/W.txt',W)\n",
    "            np.savetxt(paras_save_path+'/J_loss.txt',J_loss)\n",
    "        if (iter > 1) & ( abs(J_loss[iter-1] - J_loss[iter-2])[0]  < ABSTOL) or (iter == MAXITER):\n",
    "            break\n",
    "    W = W*W\n",
    "\n",
    "    return W, beta, J_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[13614.00065905, 13854.54173009, 13978.63253249, 13815.90896964,\n",
       "        13945.65751411, 13970.95611385, 13933.12020583, 13831.93626605,\n",
       "        13791.65313333, 13960.45349887],\n",
       "       [13854.52866687, 13645.01375407, 13820.53518615, 13877.85153554,\n",
       "        13854.17121754, 13958.31136596, 13897.7580181 , 13848.57365875,\n",
       "        13945.21157007, 13994.85571938],\n",
       "       [13978.62090025, 13820.51904867, 13510.89005773, 13975.09598579,\n",
       "        13920.5997876 , 13910.08264812, 13879.08504086, 13876.5030067 ,\n",
       "        13942.50709543, 13882.99129954],\n",
       "       [13815.91655597, 13877.86953868, 13975.10594116, 13656.51260558,\n",
       "        13876.82046914, 13788.79074326, 13931.21278304, 13928.37602901,\n",
       "        13907.21176414, 13939.0129163 ],\n",
       "       [13945.66233384, 13854.18279987, 13920.60218341, 13876.84387437,\n",
       "        13729.67621262, 13829.95981004, 13842.16692607, 13930.76536129,\n",
       "        13906.38613111, 13860.6303284 ],\n",
       "       [13970.96163116, 13958.32840462, 13910.07282788, 13788.79751598,\n",
       "        13829.96146685, 13627.27797046, 13944.01303835, 13939.11402549,\n",
       "        13845.15951072, 13883.17136055],\n",
       "       [13933.11571933, 13897.74690109, 13879.05569959, 13931.23904027,\n",
       "        13842.17645763, 13944.03254014, 13638.97458285, 13974.9314136 ,\n",
       "        13895.35653279, 13760.23087988],\n",
       "       [13831.95648533, 13848.57959749, 13876.48887361, 13928.37815785,\n",
       "        13930.76908089, 13939.11326391, 13974.93985457, 13536.56140373,\n",
       "        13879.68396312, 13950.42133856],\n",
       "       [13791.64571557, 13945.23066126, 13942.5026604 , 13907.20154861,\n",
       "        13906.38541239, 13845.16076032, 13895.33941381, 13879.67779227,\n",
       "        13674.6732816 , 13909.00887442],\n",
       "       [13960.44793166, 13994.84416317, 13882.9706366 , 13939.0273654 ,\n",
       "        13860.63898025, 13883.17138296, 13760.24673553, 13950.41764214,\n",
       "        13909.01361672, 13556.08038312]])"
      ]
     },
     "execution_count": 244,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_base = J_cost(W, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\n",
    "lambda_beta = 1\n",
    "z = prox_l1(beta - lambda_beta*grad_beta, lambda_beta*lambda4)\n",
    "f_base + grad_beta.T@(z-beta) + (1/(2*lambda_beta))*sum((z-beta)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prox_l1(v=None,lambda_=None):\n",
    "    x=np.fmax(0,v - lambda_) - np.fmax(0,- v - lambda_)\n",
    "    return x\n",
    "#prox_l1输出的最终会成为小循环中的新beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5.29118955, -1.60008736,  6.04240349, -3.0168948 ,  1.33163729,\n",
       "         1.08109643,  1.93694497,  2.24196092,  0.99816596,  1.27196264],\n",
       "       [ 3.36171466, -4.71676828, -0.42294623,  5.91425351, -2.22389619,\n",
       "         4.05181215, -1.02673809,  0.54782394,  2.93819288, -3.42344835],\n",
       "       [ 0.17863048,  3.47923778,  0.02881023,  1.3612402 ,  3.15001418,\n",
       "         3.13260724, -1.77160138,  2.20817822, -2.1069944 , -4.66412254],\n",
       "       [-4.6348411 , -1.13397937, -2.80130799,  2.18424242,  3.2029203 ,\n",
       "         3.50295381,  2.72290273, -1.74277643, -1.9783191 ,  5.67820472],\n",
       "       [ 6.82442529,  1.79846633, -3.88895982, -0.90876531, -1.80311503,\n",
       "        -2.87516321,  5.04877597, -2.88157351,  6.05271758, -2.36480831],\n",
       "       [-3.93808891,  6.5425408 , -0.33364548, -0.82070218, -0.89953817,\n",
       "        -3.19310437,  3.28831344,  8.89948704, -2.5932743 , -1.94798788],\n",
       "       [ 4.15585688,  3.40986765,  2.81348749,  0.78767159,  0.15056303,\n",
       "        -5.72769724,  0.72932893,  5.19218421, -4.87995266, -1.63530988],\n",
       "       [ 2.20613565,  1.13556997, -3.8050011 ,  1.667805  , -0.69076727,\n",
       "         0.39134425,  1.53786319,  3.16804616, -3.03737208,  2.42237625],\n",
       "       [-2.59014093,  3.74955543,  5.82403516,  2.71223536,  4.26207471,\n",
       "         1.57614999,  5.40202549, -5.71553797, -4.73983898, -5.48255826],\n",
       "       [-2.58134886,  2.28185101,  5.95074459,  6.93322488, -2.95973192,\n",
       "         5.09218773, -3.05080201, -3.31904584, -0.34687567, -2.9982039 ],\n",
       "       [-0.62172117, -1.83031176, -4.16774877,  2.69964229,  3.40198373,\n",
       "         2.27828496, -0.04142548,  5.39166231,  0.57488623, -2.68525235],\n",
       "       [-3.70647616, -3.42408052,  5.69377408, -6.11828451,  1.19329415,\n",
       "         3.13094112,  0.95652742,  2.38953671,  1.61755479,  3.26321292],\n",
       "       [-1.9360063 , -5.92634303, -3.85848136,  2.59753161,  0.22367575,\n",
       "         3.78655701,  7.48062926, -5.31783512,  0.81178205,  7.13649013],\n",
       "       [-3.50535968, -1.55195039, -0.2993995 ,  1.24763873,  3.73605169,\n",
       "        -1.21347004,  7.17510284, -4.65833695, -3.90515593,  7.97687922],\n",
       "       [-2.28383026,  2.26044372,  4.18788146, -0.86263955, -2.0898642 ,\n",
       "         5.37713779, -1.07407903,  2.07103233, -0.83984002, -1.74424224],\n",
       "       [-1.95912158,  4.52904024,  5.22213882, -4.04776707,  3.26443304,\n",
       "        -3.16800488,  2.7994337 , -2.58695633,  1.84319045, -0.89638638],\n",
       "       [ 2.3424484 , -3.46669663, -1.83614184,  2.64623697,  0.14899924,\n",
       "         2.21277719, -3.85441461,  6.79549635,  2.30747932, -2.2981844 ],\n",
       "       [-3.87765085, -2.17227133,  7.99257739, -3.75858345, -1.09123439,\n",
       "         0.74722276,  1.67679219,  1.55489604, -2.49534883,  6.42360047],\n",
       "       [ 3.68164066, -1.43294644,  6.73444178, -2.18338234, -6.15096746,\n",
       "        -7.03106016,  4.44733177,  4.54913531,  0.96888965,  1.41491723],\n",
       "       [ 2.73619489,  3.39060905,  3.56089649,  3.71277342, -0.88251245,\n",
       "        -2.06918143, -3.02759231, -0.29825481, -4.88363905,  2.7607062 ]])"
      ]
     },
     "execution_count": 218,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.fmax(beta - lambda_beta*grad_beta, lambda_beta*lambda4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f932e73d95a44307ab490440a3e0d3a7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>HBox</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta [[2389.84435352]] [[2389.84435352]] -4965479.659246893 248269.98645027477\n",
      "beta [[2389.84435352]] [[2389.84435352]] -2482739.8296234463 124134.99322513738\n",
      "beta [[2389.84435352]] [[2389.84435352]] -1241369.9148117232 62067.49661256869\n",
      "beta [[2389.84435352]] [[2389.84435352]] -620685.1041466783 31033.762980366013\n",
      "beta [[2389.84435352]] [[2389.84435352]] -310342.60890479916 15516.887173329018\n",
      "beta [[2389.84435352]] [[2389.84435352]] -155171.30445239955 7758.443586664509\n",
      "beta [[2389.84435352]] [[2389.84435352]] -77585.81896008384 3879.238466720659\n",
      "beta [[2389.84435352]] [[2389.84435352]] -38793.04633959888 1939.6329193160245\n",
      "beta [[2389.84435352]] [[2389.84435352]] -19396.640694598318 969.8282121379008\n",
      "beta [[2389.84435352]] [[2389.84435352]] -9698.32034729916 484.91410606895073\n",
      "beta [[2389.84435352]] [[2389.84435352]] -4849.1601736495795 242.45705303447522\n",
      "beta [[2389.84435352]] [[2389.84435352]] -2424.58008682479 121.22852651723764\n",
      "beta [[2389.84435352]] [[2389.84435352]] -1212.290043412391 60.61426325861842\n",
      "beta [[2389.84435352]] [[2389.84435352]] -606.1450217061976 30.307131629309424\n",
      "beta [[2389.84435352]] [[2389.84435352]] -303.0725108530984 15.15356581465467\n",
      "beta [[2389.84435352]] [[2389.84435352]] -151.53625542654882 7.576782907327292\n",
      "beta [[2389.84435352]] [[2389.84435352]] -75.76812771327876 3.7883914536640826\n",
      "beta [[2389.84435352]] [[2389.84435352]] -37.88406385663571 1.8941957268316743\n",
      "beta [[2389.84435352]] [[2389.84435352]] -18.942031928320446 0.9470978634160961\n",
      "beta [[2389.84435352]] [[2389.84435352]] -9.471015964158589 0.47354893170788465\n",
      "test1\n",
      "test2\n",
      "(375, 10)\n",
      "\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (20,10) (10,1) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-271-ffa9ff3f2893>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m W, beta, J_loss = mainFunc(X_train, y_train,\\\n\u001b[1;32m      5\u001b[0m         \u001b[0mlambda0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m         1000, ABSTOL, W_init, beta_init,)\n\u001b[0m",
      "\u001b[0;32m<ipython-input-270-40a4d7b8297f>\u001b[0m in \u001b[0;36mmainFunc\u001b[0;34m(X, Y, lambda0, lambda1, lambda2, lambda3, lambda4, lambda5, MAXITER, ABSTOL, W_init, beta_init, paras_save_path)\u001b[0m\n\u001b[1;32m     56\u001b[0m             \u001b[0mz\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprox_l1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mlambda_W\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mgrad_W\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     57\u001b[0m             \u001b[0mprint\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 58\u001b[0;31m             print ('W',J_cost(z, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5),\\\n\u001b[0m\u001b[1;32m     59\u001b[0m                    \u001b[0mf_base\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrad_W\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m                     (1/(2*lambda_W))*sum((z-W)**2))\n",
      "\u001b[0;32m<ipython-input-256-fdb25557ee4c>\u001b[0m in \u001b[0;36mJ_cost\u001b[0;34m(W, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\u001b[0m\n\u001b[1;32m      5\u001b[0m          \u001b[0;34m+\u001b[0m\u001b[0mlambda2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m          \u001b[0;34m+\u001b[0m\u001b[0mlambda3\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbeta\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m          \u001b[0;34m+\u001b[0m\u001b[0mlambda5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-255-8ed69c0408ef>\u001b[0m in \u001b[0;36mbalance_cost\u001b[0;34m(W, X, *args, **kwargs)\u001b[0m\n\u001b[1;32m      6\u001b[0m         \u001b[0mX_sub\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m         \u001b[0mI\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'double'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0meps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m         \u001b[0mloss\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mX_sub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m             \u001b[0;34m-\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_sub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m         \u001b[0;31m#print (loss.shape)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (20,10) (10,1) "
     ]
    }
   ],
   "source": [
    "W_init = np.random.rand(X_train.shape[0], 1);\n",
    "beta_init = 0.5*np.ones([feature_num, class_num]);\n",
    "\n",
    "W, beta, J_loss = mainFunc(X_train, y_train,\\\n",
    "        lambda0, lambda1, lambda2, lambda3, lambda4, lambda5,\\\n",
    "        1000, ABSTOL, W_init, beta_init,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((500, 10), (500, 1))"
      ]
     },
     "execution_count": 288,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.shape, W.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "balance (20, 10) (10, 1)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (20,10) (10,1) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-289-10d646755360>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mlambda_W\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mz\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprox_l1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mlambda_W\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mgrad_W\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mJ_cost\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbeta\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlambda5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-256-fdb25557ee4c>\u001b[0m in \u001b[0;36mJ_cost\u001b[0;34m(W, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\u001b[0m\n\u001b[1;32m      5\u001b[0m          \u001b[0;34m+\u001b[0m\u001b[0mlambda2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m          \u001b[0;34m+\u001b[0m\u001b[0mlambda3\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbeta\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m          \u001b[0;34m+\u001b[0m\u001b[0mlambda5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-286-bf20c72e618f>\u001b[0m in \u001b[0;36mbalance_cost\u001b[0;34m(W, X, *args, **kwargs)\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m         \u001b[0mloss\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mX_sub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m             \u001b[0;34m-\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_sub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mI\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m         \u001b[0;31m#print (loss.shape)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (20,10) (10,1) "
     ]
    }
   ],
   "source": [
    "grad_W = 2*lambda0*Y*np.log(sigmoid(X@beta))*W \\\n",
    "                +lambda1*balance_grad(W, X)@np.ones([X.shape[1],1]) \\\n",
    "                +4*lambda2*W*W*W \\\n",
    "                +4*lambda5*(sum(W*W)-1)*W\n",
    "lambda_W=1\n",
    "z = prox_l1(W-lambda_W*grad_W, 0)\n",
    "J_cost(z, beta, X, Y, lambda0, lambda1, lambda2, lambda3, lambda5)\n",
    "                   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {},
   "outputs": [],
   "source": [
    "def balance_cost(W=None,X=None,*args,**kwargs):\n",
    "    m = X.shape[1]\n",
    "    f_x=np.zeros([m,1])\n",
    "    for i in np.arange(0,m):\n",
    "        X_sub=np.copy(X)\n",
    "        X_sub[:,i]=0\n",
    "        I=(X[:,i] > 0).astype('double')+eps\n",
    "        print ('balance',( np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ).shape,\\\n",
    "               (np.dot((np.multiply(W,W)).T,I.reshape(-1,1))).shape)\n",
    "         \n",
    "        \n",
    "        loss=( np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ) / (np.dot((np.multiply(W,W)).T,I.reshape(-1,1)))\\\n",
    "            -(np.dot(X_sub.T,(np.multiply((np.multiply(W,W)),(1 - I.reshape(-1,1)))))) / (np.dot((np.multiply(W,W)).T,(1 - I.reshape(-1,1))))\n",
    "        #print (loss.shape)\n",
    "        f_x[i]=np.dot(loss.T,loss)\n",
    "    return f_x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_sub=np.copy(X)\n",
    "X_sub[:,0]=0\n",
    "I=(X[:,0] > 0).astype('double')+eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ],\n",
       "       [-0.23464869],\n",
       "       [-0.21352383],\n",
       "       [-0.09310426],\n",
       "       [-0.1121662 ],\n",
       "       [-0.02299837],\n",
       "       [ 0.09378211],\n",
       "       [ 0.02552632],\n",
       "       [ 0.06420116],\n",
       "       [-0.03598019],\n",
       "       [ 0.01770026],\n",
       "       [-0.02633627],\n",
       "       [-0.22665082],\n",
       "       [-0.0214175 ],\n",
       "       [-0.05105503],\n",
       "       [-0.00787399],\n",
       "       [ 0.0141895 ],\n",
       "       [ 0.08536466],\n",
       "       [-0.09877685],\n",
       "       [ 0.15218183]])"
      ]
     },
     "execution_count": 277,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "( np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ) / (np.dot((np.multiply(W,W)).T,I.reshape(-1,1)))\\\n",
    "        -(np.dot(X_sub.T,(np.multiply((np.multiply(W,W)),(1 - I.reshape(-1,1)))))) / (np.dot((np.multiply(W,W)).T,(1 - I.reshape(-1,1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J_cost(W,beta,X,Y,lambda0, lambda1, lambda2, lambda3, lambda5):\n",
    "    return lambda0*cross_entropy_loss(W,X,Y,beta) \\\n",
    "         +lambda1*sum(balance_cost(W,X)) \\\n",
    "         +lambda2*((W*W).T@(W*W)) \\\n",
    "         +lambda3*sum(beta**2) \\\n",
    "         +lambda5*(sum(W*W)-1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((500, 10), (20, 10))"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.shape,beta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = X.shape[1]\n",
    "f_x=np.zeros([m,1])\n",
    "for i in np.arange(0,m):\n",
    "    X_sub=np.copy(X)\n",
    "    X_sub[:,i]=0\n",
    "    I=(X[:,i] > 0).astype('double')+eps\n",
    "    loss=( np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ) / (np.dot((np.multiply(W,W)).T,I.reshape(-1,1)))\\\n",
    "        -(np.dot(X_sub.T,(np.multiply((np.multiply(W,W)),(1 - I.reshape(-1,1)))))) / (np.dot((np.multiply(W,W)).T,(1 - I.reshape(-1,1))))\n",
    "    #print (loss.shape)\n",
    "    f_x[i]=np.dot(loss.T,loss)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.dot( X_sub.T, np.multiply( np.multiply(W,W),I.reshape(-1,1) ) ) ) / (np.dot((np.multiply(W,W)).T,I.reshape(-1,1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 10)"
      ]
     },
     "execution_count": 254,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_base , np.sum(grad_W.T@(z-W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy = 1.0\n",
      "recall = 1.0\n",
      "precision = 1.0\n",
      "f1 = 1.0\n",
      "mcc = 1.0\n"
     ]
    }
   ],
   "source": [
    "y_pred = (sigmoid(np.dot(X_test,beta))>=0.5).astype('int')\n",
    "report_metrics(y_test, y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 417,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "230.724px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}