synthetic_experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This is the synthetic experiment used in the online adaptive pca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io\n",
    "import scipy as sp\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "font = {\n",
    "        'weight' : 'bold',\n",
    "        'size'   : 22}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### synthetic data generating functions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_cov_mat(eigval_use,dim_n):\n",
    "    \"\"\"\n",
    "    input: \n",
    "          eigval_use: the specified eigenvalues for the covariance matrix\n",
    "          dim_n: each data sample's dimension or # feature\n",
    "    output:\n",
    "          Cov_1: generated covariance matrix\n",
    "    \"\"\"\n",
    "    # This function is used to generate one covariance matrix\n",
    "    # with specified eigenvalues\n",
    "    H = rnd.randn(dim_n, dim_n)\n",
    "    Q, R = np.linalg.qr(H)\n",
    "    Cov_1 = np.dot(Q,np.dot(np.diag(eigval_use),Q.T))\n",
    "    #print np.linalg.eigvals(Cov_1)\n",
    "    return Cov_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_all_cov_mat(dim_n,rank_k,seg_num):\n",
    "    \"\"\"\n",
    "    input:\n",
    "          dim_n: each data sample's dimension or # feature\n",
    "          rank_k: rank of each covariance matrix\n",
    "          seg_num: number of different covariance matrices\n",
    "    output:\n",
    "          Cov_all: a list to store different covariance matrices\n",
    "    \"\"\"\n",
    "    # This function is used to randomly generate all different covariance matrices\n",
    " \n",
    "    idx_all = []\n",
    "    Cov_all = []\n",
    "    eigval_all = []\n",
    "    for i in range(seg_num):\n",
    "        idx_range = range(0,dim_n)\n",
    "        rnd.shuffle(idx_range)\n",
    "        idx_use = idx_range[:rank_k]\n",
    "        idx_all.append(idx_use)\n",
    "        eigval_gen = (1.*rnd.rand(1,dim_n)[0] +1.)/2. #scale to be [0.5,1]\n",
    "        eigval_use = np.zeros((1,dim_n))[0]\n",
    "        eigval_use[idx_use] = eigval_gen[idx_use]\n",
    "        Cov_gen = gen_cov_mat(eigval_use,dim_n)\n",
    "        Cov_gen = (Cov_gen + Cov_gen.T)/2.\n",
    "        #print np.linalg.matrix_rank(Cov_gen)\n",
    "        Cov_all.append(Cov_gen)\n",
    "        eigval_all.append(eigval_use)\n",
    "    #print idx_all\n",
    "    #print eigval_all\n",
    "    return Cov_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_data(Cov_all,dim_n,seg_num,data_num):\n",
    "    \"\"\"\n",
    "    input:\n",
    "          Cov_all: a list with different generated covariance matrices\n",
    "          dim_n: each data sample's dimension or # feature\n",
    "          seg_num: number of different covariance matrices\n",
    "          data_num: number of data samples generated by each covariance matrix\n",
    "    output:\n",
    "          data_gen: synthetic data as a matrix with shape (dim_n,data_num*seg_num)\n",
    "    \"\"\"\n",
    "    # This function is used to generate the data samples\n",
    "    # based on different covariance matrices.\n",
    "    # The order is: the first data_num samples are generated by\n",
    "    # the first covariance matrix. Then the second size of these\n",
    "    # is generated by the second covariance matrix, and so on.\n",
    "    mean = np.zeros(dim_n)\n",
    "    for i in range(seg_num):\n",
    "        #print i\n",
    "        Cov_use = Cov_all[i]\n",
    "        #print np.linalg.matrix_rank(Cov_use)\n",
    "        eigval = np.linalg.eigvals(Cov_use)\n",
    "        #print eigval.real\n",
    "        data_1 = rnd.multivariate_normal(mean, Cov_use, data_num)\n",
    "        for j in range(data_num):\n",
    "            if np.linalg.norm(data_1[j,:])>1.:\n",
    "                data_1[j,:] = data_1[j,:]/np.linalg.norm(data_1[j,:])\n",
    "        if i is 0:\n",
    "            #print \"initial data \"\n",
    "            data_gen = data_1.T\n",
    "        else:\n",
    "            #print \"current step \"\n",
    "            #print j ,\n",
    "            data_gen = np.hstack((data_gen,data_1.T))\n",
    "            #print data_1\n",
    "            #print data_gen\n",
    "    return data_gen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### synthetic data generation setups "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n = 20 #data dimension\n",
    "k = 2 #rank of the covariance\n",
    "seg_num = 3 #number of different covariance matrix\n",
    "data_num = 200 #number of data under each covariance\n",
    "T = seg_num*data_num #total time or horizon\n",
    "use_store_data = False # whether use already stored data\n",
    "store_new_data = False # whether to store the new generated data\n",
    "if use_store_data:\n",
    "    load_data = scipy.io.loadmat('arrdata.mat')\n",
    "    data_gen = load_data['data_gen']\n",
    "else:\n",
    "    Cov_all = gen_all_cov_mat(n,k,seg_num)\n",
    "    data_gen = gen_data(Cov_all,n,seg_num,data_num)\n",
    "if store_new_data:\n",
    "    # you need to change the directory path to your own one\n",
    "    scipy.io.savemat('/Users/Documents/online_adaptive_PCA/experiment/arrdata.mat',\n",
    "                     mdict={'data_gen': data_gen})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 600)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_gen.shape # this is the synthetic data shape (#feature, #samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the best fixed projection result: the batch PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the batch pca result\n",
    "Cov_data = np.dot(data_gen,data_gen.T)\n",
    "#print Cov_data.shape\n",
    "eigvec_cov,eigval_cov,eigvec_covh = np.linalg.svd(Cov_data)\n",
    "#eigval_sort_idx = np.argsort(eigval_cov)\n",
    "#Proj_mat = np.dot(eigvec_cov[:,eigval_sort_idx[-k:]],eigvec_cov[:,eigval_sort_idx[-k:]].T)\n",
    "Proj_mat = np.dot(eigvec_cov[:,:k],eigvec_cov[:,:k].T)\n",
    "Proj_mat = Proj_mat.real\n",
    "#get the step sum cost\n",
    "batch_step_cost = np.zeros(T)\n",
    "for i in range(T):\n",
    "    data_ob = data_gen[:,i]\n",
    "    batch_step_cost[i] = (np.linalg.norm(data_ob-np.dot(Proj_mat,data_ob)))**2\n",
    "batch_cost_step_sum = np.cumsum(batch_step_cost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xb13feda90>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3l30VhLCTEMK+BwgggopK3bViQcCFiigWK1q1Lk/V2sef1LbaWmlBQdzQioqiPsWFurBVFgmrBAhLICRsCUtCyJ7M/fsjQ6UKIZDlzPJ5XddcZ+acmcz3npl85sx9n8Wcc4iISOir5nUBIiJSNRT4IiJhQoEvIhImFPgiImFCgS8iEiYU+CIiYUKBLyISJhT4IiJhQoEvIhImanhdwIkiIiJcdHS012WIiASV1atXH3TONTvd/QIq8KOjo4mPj/e6DBGRoGJmyWW5n7p0RETChAJfRCRMKPBFRMKEAl9EJEwo8EVEwoQCX0QkTCjwRUTChAJfRMRDqUdymLlkB8t3HKr05wqoHa9ERMJFflExD7+/gY/X7QVg0rAODO7QtFKfU4EvIlLFfD7HfXPW8XnCfu4Y2p4xA6Po2LxBpT+vAl9EpApl5hby/+Zv4vOE/Tx+dTfuuCCmyp5bgS8iUkV2H8ph0j9Wk7D3KBOGtq/SsAcFvohIlfhq8wEmz1mLAdNu6sfVvVtVeQ0KfBGRSjY3PoVHPthAj9aNmHFrf1o3rutJHQp8EZFK9PWWAzw67zvO7xDBzHH9qVfLu9hV4IuIVALnHK9+s4s/fr6Fbq0aMuNWb8MeFPgiIhXO53M8/clmXv1mJ5d2bc4fR/amfm3v4/a0e9qaWSczW2hmh8wsy8y+MLMO/mXXm9l2M8szs0Vm1v6Ex91lZqlmlmtmH5tZ5e5RICISANKy8rhnzhpe/WYn44dE8/K4OCIa1Pa6LKBsh1Zo47/fk8BrwHBglpm1BN4BjgIPAf2BNwDMrC/wErDZ/7irgecrungRkUCy9UAW10z9N59t3M8jV3Tlt9d0p1o187qs/yjLb4xlzrmLjt8ws5uBHsBYoDbwjHNurpkNAG71r/3f5r/7b5xzq8zsGmCsmU10zuVVbBNERLy362A2N728kmoG8ycPpUfrRl6X9COnXcN3zhUcv25mcUATYAlwvPtmj3+a6p/GnGJZDSCynPWKiASc7WnHuOWVlfic4+07zwvIsIczOFqmmXUBPgZ2AZNPdhf/1J3JMjObaGbxZhafnp5e1nJERALCvDWpXPb8YjJyCnlj/MAqOSbO2SrTsLGZdQe+BvKBS5xz+8xsp39xW/+0jX+60385vmyvf1kR3/8K+A/n3ExgJkBcXNzJvixERALOsfwiJr21mqXbDhLX7lz+NLI3Mc0CN+yhDIFvZpHAIkq6ch4HBpnZIEoGbP8APGJmLYARwL+dczvMbDZwLzDFzL4AzgfmqP9eRELB0bxCbn9tFWtTMrjrohgmX9KJBgGw2eXplKXCDkAz//Vnjs90zpmZjQWeBZ4DVgLj/ctWm9kvgceAC4DPgPsrsG4REU+kHslh0ltr2LL/KH8f25cre1X9MXHOljkXOL0ocXFxLj4+3usyRER+JDOnkH9u2MuzCxIpKPIx/eZ+XNy1uddlAWBmq51zcae7X+D/BhER8VDK4Rwmz1nLupQMAGKa1efFm/vTpWVDjys7cwp8EZFTSDmcw82zVpKRU8CvL+tMv3bnMjimKWaBszPVmVDgi4icxO5DOYyY/g0FRT5mTxhI36hzvS6p3BT4IiI/kJlTyPjXv6XI5/jwl0MCetv6M1HmHa9ERMLB/sw8bpyxnN2Hc3jplv4hE/agNXwRkf9IOZzDTbNWcCS7kNduG8jgDqF1kF8FvogI8F1qJuNfX0VBUTFv3TGI2MjGXpdU4RT4IhLW8gqLuf/ddXyesJ/WjeryzsRBdGwefJtcloUCX0TCVkZOAY9+8B2fJ+xn0rAOTBjaPmBOVlIZFPgiEpb+lbCfxz/aSFpWPg9f0YW7h3X0uqRKp8AXkbDz0do9PDh3PVFN6vH2mFjO7xDhdUlVQoEvImEjM7eQZz7dzDurUhgc05SZ4/rTsE5Nr8uqMgp8EQl5B4/lM+mt1axOPoIDbhoUxW+v6U6dmtW9Lq1KKfBFJKRl5RUy4fVVJB7I4u5hHbm4azP6t2vidVmeUOCLSMjKKShiwuvxJOw9yku39Gd49xZel+QpHVpBRELSvsxcbnttFfHJh3l+dGzYhz1oDV9EQtDSben84s3VFPocfx3Tl2v7tPa6pICgwBeRkPLFpgP88h9riGlWn5fHxRHZpJ7XJQUMBb6IhIT0rHw+WJPKswsS6dmmEW+MH0DjerW8LiugKPBFJKjlFRbz2Icb+XBtKj4HF3SKYPrN/cJq+/qyUuCLSFByzrEg4QCPzttARk4hE4a2Z/SASDq3CM0Dn1UEBb6IBB3nHH/4fAszFifRpnFdXhjTl4s6N/O6rICnwBeRoOKcY8onm5n1752MHRjF41d3o35tRVlZ6FUSkaDh8zmemr+J15ft4rbzo3ny2u6YmddlBQ0FvogEvCPZBby2bBdLtqazLiWDO4a257Gruynsz5ACX0QC2tdbDvDAe+vJyCmkc4sGPHFNd24fEq2wPwsKfBEJSAeP5fP7Tzczb80eIpvU5S839uHiLs0V9OWgwBeRgLNsx0HuenM1WXlF/DS2Nb8f0UsDsxVAr6CIBIzMnEKe+1cib65IJqZZfV4fP5D+7c71uqyQocAXEc/lFhQz5dNNfL5xPwePFXBD3zY8eW0PGtXT3rIVSYEvIp7KzC05Qcnq3Ue4tGtzJg3rELYnKKlsCnwR8cyhY/lMfHM1G1IzmHZTP67q1crrkkKaAl9Eqlyxz/HS4h1M/Wob+UU+po7tq7CvAgp8EalSRcU+Hv5gA/PW7KFnm3P4zZXdOL9jhNdlhQUFvohUmdyCYibPWcuXmw/w4E86M/nSTl6XFFYU+CJSJYqKfdz9j9Us2prOUz/twbjB0V6XFHYU+CJS6VIO5/DU/E0sTExnyoie3DyondclhSUFvohUGuccn23cz+MfbeRwdgH3D++ssPdQtdPdwcymmtkBM3NmNv+E+Yv8845fMk5Y1s3MlplZvpklmtllldUAEQlMh47lc+fseO7+xxrq1arOZ/ddwH3D1WfvpbKu4b8D3HuS+ZuBp/zXC06YPweIAh4AJgFzzSzKOZd5toWKSPBIPZLDuFe/Zc+RXO4f3plfDIuhdo3qXpcV9k4b+M65e80smpMHfhrwiXMu6/gMM+sL9AGmO+emmVku8Aow0j8VkRBVVOxjxpIkpi3cTvVqxlt3DGJAtPaaDRSn7dI5jQuBo2Z21Mwe889r75/u8U9T/dOYk/0BM5toZvFmFp+enl7OckTEKz6f47531vHsgkSGdoxg/uShCvsAU57A/wC4BRgFpABPm9kFJ7nf8YNXu5P9EefcTOdcnHMurlkznYRYJFj9aUEin3y3j/+5siszx8XRrml9r0uSHzjrrXScc387ft3MWgFTge7At/7Zbf3TNv7pzrN9LhEJXD6fY+bSJF5avINbzoti4oUn/TEvAeC0gW9mVwM9/TcjzewOYBUlAf8hkAv8CvABq5xza81sAzDGzBIoGbTNouQXgYiEkD0ZuTz58Ua+3JzG8G7N+d21PXRGqgBWljX8h4CL/Nd7Ay9TEuLpwKNAIyAJGOecW+O/303ALOAvQDJwo3MuAxEJGZv3HeWml1eQlVfE767tzs/P13lmA11ZttIZdopFL5XymARg8FnWJCIBzDnHzCVJ/O3r7TSsU4O5vzifjs0beF2WlIH2tBWRMsvOL+KPn29h9vJkYiMb89fRsURHaHA2WCjwRaRM3li2i+cWJJKVX8Rt50fz5LXd1YUTZBT4InJas5Ym8fQnm7mgUwS3D2nPsC7NFPZBSIEvIqV6c0UyT3+ymat7tWLq2L5Ur6agD1YKfBE5pRmLd/DMZ1sY3q05fxndR2Ef5BT4IvIjWXmFPLsgkdnLk7mmdyueHx1LzerlPRKLeE2BLyL/JSOngHGvfsuG1Eyu6tWSF8aoGydUKPBF5D82pGYw4Y14jmQX8NfRsfw0trUGZ0OIAl9EAPhq8wHueXstTerX4sO7h9CrbSOvS5IKpsAXEd5ckcyTH2+kR+tGvHJbHM0b1vG6JKkECnyRMOac408LEnlx0Q4u7dqcqWP7Ur+2YiFU6Z0VCVNH8wr56xfbePWbndw0KIqnrutBDW2JE9IU+CJhZmFiGjMXJ/HdnkyO5Rfxs35tmXJ9Tw3OhgEFvkiYcM7x+cb93PfOOlo1rsN1sa25rk9rBrVvorAPEwp8kTCwZf9Rnp6/mX9vP0j3Vufw9p2DaFyvltdlSRVT4IuEsMPZBcxcksTMJTtoWKcmT17bnVvOa6e9ZsOUAl8kBK1IOsRHa/fwecJ+MnIKGdm/LY9f3U1r9WFOgS8SQnYfymHKp5tYkHCAhnVqMKRDBPdc0pGebbQTlSjwRUJCTkER0xZu5+UlO6lR3Xjo8i5MGNqeOjWre12aBBAFvkgQ25eZy+zlyby7KoXD2QXc0LcNj1zZlRbnaE9Z+TEFvkgQOpZfxHMLEnl3VQr5RcUM79aCuy6KoX+7Jl6XJgFMgS8SRHw+x+Jt6fz2442kHsllRGwbfjW8M1FN63ldmgQBBb5IENh1MJv34lP4aO0e9mbm0a5pPebeNZi4aK3RS9kp8EUC2LH8ImYuSeLFRdvxObigUwSPXNmVy3u01ICsnDEFvkgAyiss5v/W7eVPCxI5eCyfa/u05vGru2kwVspFgS8SQHYezGb6wu0s2ppOelY+vds2YtbP44iNbOx1aRICFPgiHvL5HDsPZbPtQBZfb0lj3po91KxejfM7NOW2IdEM6RBBNZ1PViqIAl/EA0eyC3hx8Q6WbE1ny/4sAOrUrMboAZHcN7yTzjgllUKBL1JFin2OZTsOsj4lgzeWJ3Mku4Durc/hsau60a/dufRq04haNXRQM6k8CnyRSnZ82/kXvtzGupQMAPq0bcTr4wfQo7WOcSNVR4EvUkmy84t47ZudzFuzh6SD2UQ0qM2zI3tzRc+WNKxT0+vyJAwp8EUqmHOOBQkHePqTTaQeySU2sjEvjInlyp6t1GUjnlLgi1Sg1CM5PDR3A8uTDtGpeQPeu2swA9trb1gJDAp8kQqQcjiHeWv2MGtpEg54+vqejBkQSQ2dWUoCiAJfpByy84t4+pNNzPk2BYDBMU3508jeRDbRwcwk8CjwRc6Cc46FiWk88VECezNz+fngdtx8Xjs6t2jodWkip6TAFzlDCxPTeHlJEst2HCKmWX3eu2swA3TUSgkCZepgNLOpZnbAzJyZzT9hfjczW2Zm+WaWaGaXnbBsiJlt8C9bY2b9KqMBIlXlcHYB97+7jvGvrWJ9SgZPXtudz++7UGEvQeNM1vDfAe79wbw5QBTwADAJmGtmUUA+8AGQC9wPPAa8b2adnHPF5a5apArtz8zjvfgUXvtmJ1l5Rdx3aSd+eXFHbWIpQadMge+cu9fMojkh8M2sL9AHmO6cm2ZmucArwEjgMNACeNg5N93MWgJPAMOAryqyASKVxTnH3NWp/L/5m8jKK6J/u3P5/YhedGmpfnoJTuXpw2/vn+7xT1P90xigUSnLFPgS8FKP5PCbDzeyZGs6A9s3Ycr1PenYvAFmOnKlBK+KHLQ9/p/gzmSZmU0EJgJERUVVYDkiZ2518hF+M+87Eg9kUa9Wdf73uh7cel47HaJYQkJ5An+nf9rWP21zwvzDpSz7L865mcBMgLi4uJN9WYhUusPZBTy7YAvvrEqhdaO6PHJFV67p3Urb00tIKVPgm9nVQE//zUgzuwNYDGwAxphZAiWDtlmUDNbmAWnAJDPLAiYAu4BFFVm8SHkV+xzvrNrNswsSycorYsKQ9tw3vJMObiYhqaxr+A8BF/mv9wZeBsYDNwGzgL8AycCNzrkMADMbBUwDXgASgDu1hY4EkrW7j/DbjxP4bk8m58U04amf9tSOUxLSyrqVzrBSFg8+xWOWAL3OoiaRSpWUfowZi5N4Nz6FFufUZurYvlzbu5UGZCXkaU9bCRtbD2Qxa2kS78WnUqOaMfHCGO69tBMNauvfQMKDPukS8talZPCHzzazatcRAG47P5oJQ9trQFbCjgJfQlZmbiHPLUjkrZXJRDSozR0XtOfOC2KIaFDb69JEPKHAl5BzLL+It1cmM3NJEoezC/j54GgevKyztryRsKfAl5BRWOzjtW92Mm3hDjJzCzm/Q1N+c1U3erbRicJFQIEvISAzt5DpC7fz/upUDmUXcEnX5tx7aSdiIxt7XZpIQFHgS9Aq9jnei0/huQWJHM4p4MqeLRkVF8mwzs20iaXISSjwJejsOpjN+tQMZi5JImHvUQZEn8sb1w5U143IaSjwJWjkFRYzbeF2Xly0gyKfo3WjOvxtbF+u0U5TImWiwJeA5pxjdfIRfv/pZjbvyyK3sJgb+rXh9iHt6di8AXVqVve6RJGgocCXgLX1QBa/+78Elu04RPOGtRkzMJLh3VowpGOE16WJBCUFvgSUtKN5vLw0iYWJ6WxPO0ajujV54prujIpryznajl6kXBT4EhBSDufw1spk3li2i6Jix9BOEVzRoyW3D21Pk/q1vC5PJCQo8MVTuw/l8NT8BL7akgbAdX1a88BPOtOuaX2PKxMJPQp88YTP5/jHymSe+WwL1c2YfHFHRg+Mok3jul6XJhKyFPhS5VbtOsyjH2xgR3o2F3SK4I8/601rBb1IpVPgS5VISj/Gkq3pbNiTybw1e2h7bl3+OjqWn8a21jb0IlVEgS+VJvlQNutTM3l7ZTIrkkrOa1+rejVuH9KeBy/rTH2deESkSuk/Tiqcz+eYvXwXT83fhM/BufVqcs/FHRnZvy1RTepRrZrW6EW8oMCXCvXlpgP8feF21qVkMKxLMx6+vCvREfWoV0sfNRGv6b9QKsR3qZl8sWk/0xft4Jy6NfnTyN6M6t9W/fMiAUSBL+WSlpXHF5sO8NuPEyj2OQZEn8srtw3QXrEiAUiBL2clcX8Wv3p3HZv3HQUgNrIxM27tT/OGtbVWLxKgFPhyxhZuSeOB99ZRs3o1HruqG30iG9MnshG1a+jIlSKBTIEvZZaWlcfXm9N47KONtGtaj9duG6BDIIgEEQW+nNbWA1lM/WobCxL2U1js6BfVmDduH0hD9dOLBBUFvpRqT0Yut8xayaHsAm6Ma8v1sW2Ii25CdW1LLxJ0FPhySpm5hYx/7VtyC4v59N4L6NKyodcliUg5VPO6AAlMuw5mc/vrq9h5MJsZt/ZX2IuEAK3hy39JOZzDWyuSefWbndSqXo3nR8dyfgedUlAkFCjwBYC8wmLeXZXCU/M3UexzjOrfloeu6ELzhnW8Lk1EKogCP8w553hzRTJ//GwL2QXFnBfThKev70XH5g28Lk1EKpgCP4xl5hZy9z9W8832Q1zQKYIb4yIZ3q0FdWtpByqRUKTAD0POOb7cnMaUTzaRciSXR67oyi8uitEhEURCnAI/zCSlH2PynLUk7D1Kh2b1mX37QIZ01KCsSDhQ4IeJnIIiZi3dyRvLduGAp6/vyegBkdSsri1zRcKFAj8MrNl9hCc+2kjC3qP0jWrMc6P60KGZBmVFwk25V+/MbJeZuRMu6/zzh5jZBjPLN7M1Ztav/OXKmXprRTJjZ65gX2YeL4+L48O7hyjsRcJURa3hLwFe9F8/YmZ1gA+AXOB+4DHgfTPr5JwrrqDnlFKkHslh2sIdzPl2Nxd2bsZfR8fSpH4tr8sSEQ9VVODvBD5xzmUBmNkIoAXwsHNuupm1BJ4AhgFfVdBzyklsT8viwffWk7C35MQkP+vXlj/8rJf66kWkwgJ/HPBzM0sH/gdo5J+/xz9N9U9jUOBXimKfY863u3nuX4lUN+MXF3XgpkFRtG5c1+vSRCRAVETgvwwkAnWAPwAzgEd+cJ/jG3i7Hz7YzCYCEwGioqIqoJzwcyS7gDtnxxOffER7yorIKZU78J1zU45fN7O+wAN8v0bf1j9t45/uPMnjZwIzAeLi4n70hSCl+3bnYR6dt4HUI7n85cY+jOjbRjtQichJlSvwzawX8HvgM//fGkfJQO1SIA2YZGZZwARgF7CoPM8n/+2D1ak8/MEGGtetyZu3D2RQTFOvSxKRAFbekbyDQHXgKUq6c5KBEc65vcAo4BjwAiXhP0pb6FSc177ZyYNz13NeTBMWP3yxwl5ETqtca/jOuX3AVadYtgToVZ6/Lz/mnOP5L7cx9attXN6jBVPH9qV2DR3sTEROT3vaBpGiYh8Pvb+BD9fuYVT/tjxzQy9qaHNLESkjBX6QOJpXyL1z1rIoMZ0Hf9KZey7pqMFZETkjCvwAV1js48M1e/jbwm3sy8hjyoie3DyonddliUgQUuAHsO1pWTw4dwPrUzLo2eYc/jwqloHtm3hdlogEKQV+ACos9jFzSRIvfLmNerWrM3VsX67t3UpdOCJSLgr8AJOwN5OH5m5g076jXN2rFb+7rgfNGtb2uiwRCQEK/ADy7qrdPPbhRhrXq8VLt/Tnip4tvS5JREKIAj9ALN2WzmMfbmRwh6b8bWxfGtfToYxFpGIp8APAM59uZsaSJDo1b8C0m/txTp2aXpckIiFIge+h5EPZPPPpFj5P2M/1sa353+t6KuxFpNIo8D2wPS2LZxcksjAxnepmjI6LZMqIntprVkQqlQK/ChUW+/hq8wHue2cdtWpUY0RsGx64rDMtzqnjdWkiEgYU+FWgoMjHd3syeeKjjWzad5T2EfV5d+J5NFfQi0gVUuBXoqJiH7+eu56P1u0FoGn9Wvzlxj5c2bMVdWvpCJciUrUU+JWgqNjHnFUpvL1yN5v3HeWmQVHERjbmyp4taahBWRHxiAK/gmXlFTLh9Xi+3XWY6Kb1eGFMLD+NbXP6B4qIVDIFfgVKO5rHL95azfrUTJ4b1Yef9dP5ZUUkcCjwK0BeYTHTF25n5tIkfD74+9i+XNmrlddliYj8FwV+OaVn5TNm5nJ2pGdzde9WPHRZF6Ij6ntdlojIjyjwz1J+UTEfrN7Dswu2kF1QzAtjYrmuT2t14YhIwFLgnyHnHPM37GPKJ5vZfzSP2MjGPDuyN51aNPS6NBGRUinwyyi/qJilWw/yzqrdfLk5jbbn1mXGrf25rHsLrdWLSFBQ4JfBnoxcbnxpOXsycqlVvRq/Gt6JCUPba5t6EQkqCvxSpGXl8fwXW/liUxr5RcXMGhfHeR2a0qC2XjYRCT5KrpPw+Rz/3LCX//3nJo7lFTGkY1N+NbwzfSIbe12aiMhZU+D/wOHsAsa/vor1KRk0rV+L/5s8hK4tz/G6LBGRclPgnyDlcA6jZyznYHYBz43qwzW9W1Gnpg5yJiKhQYHvt/VAFiNfXIbPwXt3DSZW3TciEmIU+MCChP08/P4GateszjsTz6NDswZelyQiUuHCOvAPZxdwz9trWLbjEF1aNOTlcXFENa3ndVkiIpUibAP/pcU7mPb1dnILi3ngJ525e1gHnVNWREJa2AV+Zk4hb67YxXP/2sqg9k14+Iqu9G93rtdliYhUurAK/ANH8xg9Yzm7DuVwSdfmzLi1PzW1Vi8iYSIsAj89K5+/f72NuatTAXjttgEM69JMx8ARkbAS8oE/Nz6Fp/65ibyiYq7r04a7Loqhs45sKSJhKGQDv9jnmL5wO3/+YiuDY5oyZURPYrS5pYiEsZAM/KN5hfz6vfX8a9MBru3Tmj+P6kOtGuqrF5HwFnKBX1js4+631rA86RCPXdWNOy5or756ERGgUld7zWyImW0ws3wzW2Nm/Srz+fIKi/n13PX8e/tB/nBDL+68MEZhLyLiV2mBb2Z1gA+AhsD9QAvgfTOrlKOR7TyYzQ3Tl/Hxur08dHkXRsVFVsbTiIgErcrs0rmSkpB/2Dk33cxaAk8Aw4CvKvKJFiWmcc/ba6lR3Xjl53Fc2q1FRf55EZGQUJmB394/3eOfpvqnMVRw4Ec3rU+/dufyzA29aNO4bkX+aRGRkFGVg7bHO9Pdf800mwhMBIiKijqrPxwdUZ/Ztw8sV3EiIqGuMgdtd/qnbf3TNj+YD4BzbqZzLs45F9esWbNKLEdEJLxV5hr+Z0AaMMnMsoAJwC5gUSU+p4iInEKlreE75/KAUcAx4AVKwn+Uc664sp5TREROrVL78J1zS4BelfkcIiJSNjregIhImFDgi4iECQW+iEiYUOCLiIQJc86d/l5VxMzSgeSzfHgEcLACy/GS2hKY1JbAEyrtgPK1pZ1z7rQ7MgVU4JeHmcU75+K8rqMiqC2BSW0JPKHSDqiatqhLR0QkTCjwRUTCRCgF/kyvC6hAaktgUlsCT6i0A6qgLSHThy8iIqULpTV8EREpRdAHflWfN7e8zGyqmR0wM2dm80+Y383MlvnbkWhml52wLODaaGadzGyhmR0ysywz+8LMOviXXW9m280sz8wWmVn7Ex53l5mlmlmumX1sZk29a8X3zGylvx05ZhZvZhf65wddW6DkFKP+z5Ezs7/75wXVZwzAzHb523D8ss4//5T1lvaeecnMGpvZbDPLMLNjZrbEP7/q2uKcC9oLUAfYT8kx9u+m5OxaSUB1r2srpeaplBw91AHzT5i/DjgM/BLYCGQCjQK1jZScqnIxcI+/TQ5YCLQE8oA1wGQgC1jif0xf//2+AB4GioDZXr8n/tqeB8YD/+Ova2uwtsVf3++BbH+Nfw/Gz5i/5l3+z9kY/+Xy0uot7T3z+gLMA4qBPwN3AK9WdVs8fxHK+QKO8H+gH/Lffsp/+1KvaztN3dEnBv4J4THNf/t2/+0JgdpGoNYPbh+i5BDY9/vrG+WfP9t/uwPff9EN8C9bAhQCdQLgPTFKdnwZ6A/KLUHclt5ALvDQ8cAPxs+Yv5ZdwOtAwxPmnbLe0t4zj9sR46/jLaAW/i/Tqm5LsHfplHbe3GBSWjsCso3OuYLj180sDmhCSeidaVtqAJGVWmzZNALSgZVAASVrYEHXFjOrBswCpgGrTlgUdJ+xE4wDjppZmplNIDiMIEhdAAACOElEQVTb0t0/HUDJCkW2mf2RKm5LsAf+D530vLlBqLR2BFQbzawL8DEla2KTT3YX/zTQ23IMuAy4l5Kf2U+d5D7B0JbxlPyCnM33pxVtBNT8wf2CoS0ALwM3ArdS8kU8g+/rOy4Y2lLbP60PjAa+oaQr8IfnJKnUtlTlScwrQ5nOmxsESmvH4VKWecrMugNfA/nAJc65fWZWWltOXLbXv6yI79dcPOOcK6KkP/4LMxsJXEzJaTohuNoSCTQD1p8w7xagtf96UH3GnHNTjl83s77AA3z/Gp+s3oalLPPSLv90qXNunpk1Ay7h+xCvmrZ43UdXzn6xOsAB/wswiZKfPjsJgMGmUmq+GniEkm/p9ZR0HXTyXz9xQO0o0DhQ20hJsKRREnKP8v2gWitKvgBW8/1A01L/Y/rz44HONwPgPbkceIWS/uzfUdIXvz9I29IdGOm/POmv8TN/vcH2GesF/JOSwcx7Kelyy6Hky+uk9Zb2nnncFgM2+P9n7gRW+D8zPauyLZ6+CBX0Ql4IfEfJz721QJzXNZ2m3kX+f8ITL7cBPYDl/jd4K3BFILeRkq10ftgO5192A7DD35YlnDDIxPdbIuT5/5kjAqAtA/wBmAtkULK10fHB2KBqyyneo+Nb6QTbZ6wV8CklR5DMAeKBy09Xb2nvmcftOf765/lf/5uqui3a01ZEJEyE2qCtiIicggJfRCRMKPBFRMKEAl9EJEwo8EVEwoQCX0QkTCjwRUTChAJfRCRM/H97reHRWpQYSAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(batch_cost_step_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.44810583e+02 9.16410797e+01 7.79436149e+01 6.27130162e+01\n",
      " 4.94685765e+01 2.33555255e+01 4.13723141e-14 3.29695387e-14\n",
      " 2.46543718e-14 2.23175864e-14 1.86092680e-14 1.53692641e-14\n",
      " 1.45109593e-14 1.11958310e-14 9.07576581e-15 7.82910881e-15\n",
      " 6.36242413e-15 5.86091365e-15 3.49415475e-15 2.12134494e-15]\n"
     ]
    }
   ],
   "source": [
    "print eigval_cov # Check the eigenvalues' information to explain the failure of batch PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixture decomposition function: Algorithm 2 in the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixture_decompose_lift(w_input,n,k):\n",
    "    \"\"\"\n",
    "    input: \n",
    "          w_input: the eigenvalues/singular values vector for the decomposition\n",
    "          n: data sample dimension or # features\n",
    "          k: the rank of our desired projection matrix\n",
    "    output:\n",
    "          r_all: all decomposed corners with n-k non-zero elements equal to 1\n",
    "                 Note: Here instead of the value of 1/(n-k) used in the algorithm 2 in paper, \n",
    "                       we lift it to 1 in order to avoid the small value resulted from 1/(n-k)\n",
    "          p_all: each corner's associated coefficient in [0,1] in order to do the sampling\n",
    "    \"\"\"\n",
    "    # Note: the input w_input is already lifted be in the range [0,n-k]\n",
    "    \n",
    "    #print \"start decompose\"\n",
    "    d = n-k\n",
    "    p_all = []\n",
    "    r_all = []\n",
    "    #w_input = 1./n*np.ones(n)\n",
    "    w_use = w_input\n",
    "    count_num = 0\n",
    "    #print \"mixture decompose: \"\n",
    "    p = 1.\n",
    "    #while np.sum(w_use)/d>1e-5:\n",
    "    while p>1e-7:\n",
    "        #print count_num ,\n",
    "        #print np.sum(w_use)/d\n",
    "        #print w_use\n",
    "        count_num = count_num + 1\n",
    "        diff_corner_val0 = w_use-np.sum(w_use)/d \n",
    "        if np.all(diff_corner_val0<=-1e-6):# the case when no components equal to w/d\n",
    "            #randomly select a corner\n",
    "            #print \"go in random\"\n",
    "            # to avoid selecting the one with 0 probability\n",
    "            idx_corner = rnd.choice(len(w_use),size=d,replace = False, p=w_use/np.sum(w_use))\n",
    "            w_no_select = np.abs(w_use-0.)\n",
    "            w_no_select[idx_corner] = 0.\n",
    "            r = np.zeros(n)\n",
    "            r[idx_corner] = 1.\n",
    "            s = np.amin(w_use[idx_corner])\n",
    "            l = np.amax(w_no_select)\n",
    "            p = np.minimum(s,np.sum(w_use)/d-l) \n",
    "            p_all.append(p)\n",
    "            r_all.append(r)\n",
    "            w_next = w_use - p*r\n",
    "            w_use = w_next\n",
    "        else: # the case when some components equal to w/d\n",
    "            idx_gen = range(n)\n",
    "            # get the needed corner\n",
    "            diff_corner_val1 = np.abs(w_use-(np.sum(w_use)/d))\n",
    "            idx_init_corner = list(np.where(diff_corner_val1<=1e-6)[0]) \n",
    "            if len(idx_init_corner)==0:\n",
    "                print \"error in decomposition \"\n",
    "            idx_corner_part = [idx_gen[ii] for ii in idx_init_corner] # initial index of the corner\n",
    "            # check if corner idx is enough\n",
    "            if len(idx_corner_part)<d: # the case when the # components equal to w/d is smaller than d\n",
    "                                       # we need to get more\n",
    "                diff_num = d-len(idx_corner_part)# get the additional # components needed\n",
    "                diff_corner_val2 = np.abs(w_use-(np.sum(w_use)/d))\n",
    "                idx_init_no_corner = list(np.where(diff_corner_val2>1e-6)[0])# get the unsatisfied index smaller than w/d\n",
    "                p_no_corner_part = [w_use[jj] for jj in idx_init_no_corner] #extract the unsatisfied index of w\n",
    "                pp_test = np.array(np.abs(np.array(p_no_corner_part)-0))\n",
    "                pp_test_asd_idx = np.argsort(pp_test)\n",
    "                pp_test_dsd_idx = pp_test_asd_idx[::-1]\n",
    "                if pp_test[pp_test_dsd_idx[diff_num-1]]<1e-8:# check whether the highest unsatisfied component \n",
    "                                                             # of w is very close to 0, then no need to further process\n",
    "                    break\n",
    "                # randomly select the remaining # of index to add to corner index based on their value in w\n",
    "                idx_no_corner_choice = rnd.choice(len(p_no_corner_part),size=diff_num,\n",
    "                                            replace = False, p=p_no_corner_part/np.sum(p_no_corner_part))\n",
    "                idx_init_no_corner = np.array(idx_init_no_corner)\n",
    "                idx_corner_add_part = idx_init_no_corner[idx_no_corner_choice]\n",
    "                idx_corner = np.hstack((idx_corner_part,idx_corner_add_part)) # the final corner\n",
    "                w_no_use = np.array(p_no_corner_part)\n",
    "                w_no_use[idx_no_corner_choice] = 0.\n",
    "                l = np.amax(w_no_use)\n",
    "            else: # the case when there are enough components equal to w/d\n",
    "                idx_corner = idx_corner_part \n",
    "                l = 0.\n",
    "            r = np.zeros(n)\n",
    "            r[idx_corner] = 1.\n",
    "            s = np.amin(w_use[idx_corner])\n",
    "            p = np.minimum(s,np.sum(w_use)/d-l)\n",
    "            p_all.append(p)\n",
    "            r_all.append(r)\n",
    "            w_next = w_use - p*r\n",
    "            w_use = w_next\n",
    "        if count_num>3*n: \n",
    "            print \"too many steps for decomposition, something wrong\"\n",
    "            print w_input\n",
    "            break\n",
    "    return r_all,p_all\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Capping algorithm function: Algorithm 3 in paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def capping_alg_lift(w,n,k):\n",
    "    \"\"\"\n",
    "    input:\n",
    "          w: the eigenvalue or singular value vector that needs to be capped\n",
    "          n: the data dimension or # feature\n",
    "          k: the rank of our desired projection matrix\n",
    "    output:\n",
    "          the capped w with each component in the range [0,1]\n",
    "    \"\"\"\n",
    "    d = n-k\n",
    "    d_upper = 1.\n",
    "    w_lift = w*d # we lift the input to make each component in the range [0,d]\n",
    "    if np.all(w_lift<1./d_upper+1e-7): # the case when the input w already satisfies the capping constraint\n",
    "        #print \"no capping \"\n",
    "        w_lift[np.where(np.abs(w_lift)<1e-10)] = 1e-10\n",
    "        return w_lift\n",
    "    \n",
    "    else: # the case when the input w not satisfies the capping constraint\n",
    "        idx_asd = np.argsort(w_lift)\n",
    "        idx_dsd = idx_asd[::-1]\n",
    "        w_lift_dsd = w_lift[idx_dsd]\n",
    "    w_tilde = w_lift_dsd\n",
    "    i = 1\n",
    "    while np.max(w_tilde)>1./d_upper+1e-7:\n",
    "        w_tilde = np.zeros(n)\n",
    "        w_tilde[:i] = 1./d_upper\n",
    "        sum_remain = np.sum(w_lift_dsd[i:])\n",
    "        for j in range(i,n):\n",
    "            w_tilde[j] = (d-i)*w_lift_dsd[j]/sum_remain\n",
    "        i = i +1\n",
    "\n",
    "    w_return = np.zeros(n)\n",
    "    w_return[idx_dsd] = w_tilde\n",
    "    #let the minimum of w_return to be bigger than 1e-10\n",
    "    w_return[np.where(np.abs(w_return)<1e-10)] = 1e-10\n",
    "    return w_return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### test the capping alg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "##test the capping alg\n",
    "#w_input = np.array([1./17,1./38,1./36,1./72,\n",
    "#                   1./36,1./33,1./36,1./43,3./48,1./54])\n",
    "#w_input = w_input/np.sum(w_input)\n",
    "#print w_input\n",
    "#w_input = capping_alg_lift(w_input,10,2)\n",
    "#print w_input\n",
    "#print 1./(10-2)\n",
    "#print \"start mixture decomposition \"\n",
    "#r_all,p_all = mixture_decompose_lift(w_input,10,2)\n",
    "#print \"done\"\n",
    "\n",
    "#print np.sum(p_all)\n",
    "#for i_w in range(len(p_all)):\n",
    "#    if i_w is 0:\n",
    "#        w_out = r_all[i_w]*p_all[i_w]\n",
    "#    else:\n",
    "#        w_out = w_out + r_all[i_w]*p_all[i_w]\n",
    "#print w_input-w_out\n",
    "#print r_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Follow the Leader algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def follow_the_leader(k,data_gen,T):\n",
    "    \"\"\"\n",
    "    input: \n",
    "          k: the rank of our desired projection matrix\n",
    "          data_gen: the synthetic data with shape (#feature, #samples)\n",
    "          T: time horizon\n",
    "    output:\n",
    "          cumulative loss or regret\n",
    "    \"\"\"\n",
    "    step_cost_ftl = np.zeros(T)\n",
    "    print \"follow the leader starts:\"\n",
    "    for tt in range(T):\n",
    "        if tt<k:\n",
    "            continue\n",
    "        else:\n",
    "            Cov_data_ftl = np.dot(data_gen[:,:tt],data_gen[:,:tt].T)\n",
    "            eigvec_cov_ftl,eigval_cov_ftl,eigvec_covh_ftl = np.linalg.svd(Cov_data_ftl)\n",
    "            Proj_mat_ftl = np.dot(eigvec_cov_ftl[:,:k],eigvec_cov_ftl[:,:k].T)\n",
    "            data_ob = data_gen[:,tt]\n",
    "            step_cost_ftl[tt] = (np.linalg.norm(data_ob-np.dot(Proj_mat_ftl,data_ob)))**2    \n",
    "    print \"follow the leader done\"\n",
    "    return np.cumsum(step_cost_ftl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### online PCA or static online PCA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def static_pca_online(n,k,eta,alpha,data_gen,T):\n",
    "    # although the function name is static_pca_online,\n",
    "    # it is also used in adaptive pca, since the later one augments the previous one.\n",
    "    \"\"\"\n",
    "    input:\n",
    "          n: the data dimension or # feature\n",
    "          k: the rank of our desired projection matrix\n",
    "          eta: the stepsize used in step (14a)\n",
    "          alpha: the fixed-share value in step (14b). For static PCA, it is 0. For adaptive PCA, it is non-zero\n",
    "          data_gen: the synthetic data with shape (#feature, #samples)\n",
    "          T: time horizon\n",
    "    \"\"\"\n",
    "    W_init = (1.*(n-k)/n)*np.eye(n)\n",
    "    W_cur = W_init\n",
    "    step_cost_static_pca = np.zeros(T) # one time step cost resulted from sampling of the corners\n",
    "    exp_cost_static_pca = np.zeros(T) # expected step cost\n",
    "    print \"online pca starts: \"\n",
    "    for tt in range(T):\n",
    "        if tt%100 == 0:\n",
    "            print tt ,\n",
    "        # get the eigen-decomposition result\n",
    "        if tt is 0:\n",
    "            eigvec_W_cur,eigval_W_cur,eigvec_W_cur_h = np.linalg.svd(W_cur)\n",
    "        else:\n",
    "            eigvec_W_cur = eigvec_W_hat\n",
    "            eigval_W_cur = eigval_W_hat_cap_lift\n",
    "        # call the mixture decomposition function\n",
    "        r_candidate,p_dist = mixture_decompose_lift(eigval_W_cur,n,k)\n",
    "        # randomly select the corner\n",
    "        idx_pick = rnd.choice(len(p_dist),p = p_dist/np.sum(p_dist))\n",
    "        r_corner = r_candidate[idx_pick]\n",
    "        if np.abs(np.sum(r_corner)-(n-k))>1e-6:\n",
    "            print \"wrong corner \"\n",
    "        Mat_corner = np.dot(eigvec_W_cur,np.dot(np.diag(r_corner),eigvec_W_cur.T))\n",
    "        # get the projection matrix\n",
    "        Proj_use_mat = np.eye(n) - Mat_corner\n",
    "        # get the data sample\n",
    "        x_cur = data_gen[:,tt]\n",
    "        # get the project loss\n",
    "        step_loss = (np.linalg.norm(x_cur-np.dot(Proj_use_mat,x_cur)))**2.\n",
    "        step_cost_static_pca[tt] = step_loss\n",
    "        exp_cost_static_pca[tt] = np.trace(np.dot(W_cur,np.outer(x_cur,x_cur)))\n",
    "        # some error tests\n",
    "        if np.any(eigval_W_cur<0.):\n",
    "            print \"negative eigval of W_cur!!!\"\n",
    "            print \"the most negative value is \"\n",
    "            print np.amin(eigval_W_cur) ,\n",
    "        if np.any(np.iscomplex(W_cur)):\n",
    "            print \"W_cur is complex matrix!!!\"        \n",
    "        if not np.all(np.isfinite(W_cur)):\n",
    "            print \"W_cur is not finite!\"\n",
    "            scipy.io.savemat('/Users/Documents/online_adaptive_PCA/experiment/W_cur.mat',\n",
    "                     mdict={'W_cur': W_cur}) \n",
    "        # this seems wrong, but the minus part is due to avoid log(0)\n",
    "        logm_W_cur_cap = sp.linalg.logm(W_cur) - np.log(n-k)*np.eye(n)\n",
    "        if np.any(np.iscomplex(logm_W_cur_cap)):\n",
    "            print \"logm_W_cur_cap is complex matrix !\"\n",
    "            scipy.io.savemat('/Users/liudingyi/Documents/online_adaptive_PCA/experiment/W_cur.mat',\n",
    "                     mdict={'W_cur': W_cur})\n",
    "        W_hat = sp.linalg.expm(logm_W_cur_cap-eta*np.outer(x_cur,x_cur))\n",
    "        # some error tests\n",
    "        if np.any(np.iscomplex(W_hat)):\n",
    "            print \"W_hat is complex matrix!!!\"\n",
    "            break\n",
    "        if not np.all(np.isfinite(W_hat)):\n",
    "            print \"W_hat is not finite!\"\n",
    "            scipy.io.savemat('/Users/Documents/online_adaptive_PCA/experiment/W_hat.mat',\n",
    "                     mdict={'W_hat': W_hat})    \n",
    "            break\n",
    "        W_hat_svd = (n-k)*W_hat/np.trace(W_hat)\n",
    "        W_hat_svd = (W_hat_svd + W_hat_svd.T)/2.\n",
    "        # add identity matrix to avoid computation precision near 0\n",
    "        W_hat_svd_off_1 = W_hat_svd + np.eye(n)\n",
    "        eigvec_W_hat, eigval_W_hat_svd_off1,eigvec_W_hat_h = sp.linalg.svd(W_hat_svd_off_1)\n",
    "        # subtract the added identity eigenvalue\n",
    "        eigval_W_hat_svd = np.abs(eigval_W_hat_svd_off1 - 1.)\n",
    "        eigval_W_hat = eigval_W_hat_svd/(n-k)\n",
    "        # the fixed-share step\n",
    "        eigval_W_hat = (1-alpha)*eigval_W_hat + alpha/n \n",
    "        eigval_W_hat = eigval_W_hat/np.sum(eigval_W_hat) #make sure sum is 1.\n",
    "        # call the capping algorithm\n",
    "        eigval_W_hat_cap_lift = capping_alg_lift(eigval_W_hat,n,k)\n",
    "        # error test\n",
    "        if not np.all(eigval_W_hat_cap_lift<=1.+1e-6):\n",
    "            print \"wrong capping\"\n",
    "        # construct the new density matrix\n",
    "        W_next = np.dot(eigvec_W_hat,np.dot(np.diag(eigval_W_hat_cap_lift),eigvec_W_hat.T))\n",
    "        W_cur = W_next\n",
    "    print \"online pca done \"\n",
    "    return np.cumsum(step_cost_static_pca),np.cumsum(exp_cost_static_pca)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### different algorithms's setup and test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n is \n",
      "20\n",
      "k is \n",
      "2\n",
      "D is \n",
      "1.8964892818408743\n",
      "L is \n",
      "213.48073311729274\n",
      "static pca stepsize is \n",
      "0.18872643671886297\n",
      "1e-05\n"
     ]
    }
   ],
   "source": [
    "# determine stepsize, since k<=n/2, we can use D = (n-k)log(n/(n-k))\n",
    "print \"n is \"\n",
    "print n\n",
    "print \"k is \"\n",
    "print k\n",
    "D_upper = (n-k)*np.log(n*1./(n-k))\n",
    "D_upper_static = k*np.log(n*1./k)\n",
    "print \"D is \"\n",
    "print D_upper\n",
    "L_upper = batch_cost_step_sum[-1]\n",
    "print \"L is \"\n",
    "print L_upper\n",
    "eta_static = np.log(1+np.sqrt(2.*D_upper_static/L_upper))\n",
    "print \"static pca stepsize is \"\n",
    "print eta_static\n",
    "eta_static = 1.\n",
    "eta_adaptive = np.log(1+np.sqrt(2.*D_upper/L_upper))\n",
    "alpha = 1./(1+T*(n-k))\n",
    "alpha = 0.00001 # the fixed-share value can be set heuristically\n",
    "print alpha\n",
    "eta_static = 1. # the stepsize for static online pca\n",
    "eta_adaptive = 1. # the stepsize for adaptive online pca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### run the two online pca algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "online pca starts: \n",
      "0 100 200 300 400 500 online pca done \n",
      "done\n"
     ]
    }
   ],
   "source": [
    "trial = 5 # the number of runs for each algorithm\n",
    "num_comp = 4 # the number of algorithms in the experiment\n",
    "cumsum_cost_pca = np.zeros((num_comp,trial,T)) # matrix to store all the loss results\n",
    "for tr_i in range(trial):\n",
    "    #print tr_i ,\n",
    "    cumsum_cost_pca[0,tr_i,:], exp_cost_static_pca = static_pca_online(n,k,eta_static,0,data_gen,T)\n",
    "    cumsum_cost_pca[1,tr_i,:], exp_cost_adaptive_pca = static_pca_online(n,k,eta_adaptive,alpha,data_gen,T)\n",
    "print \"done\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### get the follow the leader and batch pca results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "follow the leader starts:\n",
      "follow the leader done\n"
     ]
    }
   ],
   "source": [
    "cumsum_cost_pca_mean = np.zeros((num_comp,T))\n",
    "for num_comp_i in range(num_comp-2):\n",
    "    cumsum_cost_pca_mean[num_comp_i,:] = cumsum_cost_pca[num_comp_i,:,:].mean(axis=0)\n",
    "cumsum_cost_pca_mean[2,:] = batch_cost_step_sum\n",
    "cumsum_cost_pca_mean[-1,:] = follow_the_leader(k,data_gen,T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot the different algorithms' results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_total_loss(plot_list,T_list,Loss_list):\n",
    "    fig,ax = plt.subplots(1,1,sharex=True,figsize=(5,5),squeeze=False)\n",
    "    colors = ['b:','g--','r','k']\n",
    "\n",
    "    labels = ['Online PCA', 'Online Adaptive PCA',\n",
    "              'Best fixed Projection','Follow the Leader']\n",
    "    titles = ['Total Loss Comparison']\n",
    "    ax = ax.flatten()\n",
    "    for conf_i in plot_list:\n",
    "        ax[0].plot(T_list,Loss_list[conf_i,:],colors[conf_i],\\\n",
    "                       label = labels[conf_i],linewidth = 4,markersize=10)\n",
    "\n",
    "    ax[0].legend(loc = 2)   \n",
    "    plt.xlabel('Time Step',fontsize = 14, fontweight='bold')\n",
    "    plt.ylabel('Total Loss',fontsize = 14, fontweight='bold')\n",
    "    #plt.savefig('subspace_synthetic_online_comp.pdf',bbox_inches='tight',transparent=True,type ='pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### two online pca comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAFECAYAAACnCrwxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VNX5+PHPk0ASCHvYhLDJvgoaxA2sSvFLFWstKq4Fxa0uP7XydWuV1r0urYJad1wqBZSvCxartiIqYgFZBAQBgxB2gmJYErbn98e5k2SSmWQyZObOJM/79bqvmXvm3HufbE/Ovefec0RVMcYYUzUpfgdgjDHJyJKnMcZEwZKnMcZEwZKnMcZEwZKnMcZEwZKnMcZEwZKnMcZEwZKnMcZEwZKnMcZEoY7fAUSrefPm2rFjR7/DMMbUMAsWLNiuqi0qq5e0ybNjx47Mnz/f7zCMMTWMiHwfST07bTfGmChY8jTGmChY8jTGmCgk7TXPUPbv309eXh6FhYV+h2JqoIyMDLKzs6lbt67foZgEUKOSZ15eHg0bNqRjx46IiN/hmBpEVcnPzycvL49OnTr5HY5JADXqtL2wsJCsrCxLnKbaiQhZWVl2VmOK1ajkCVjiNDFjv1umtBqXPI0xJh4seVaz/Px8xowZQ4sWLUhPT6d///5MmzYtom3Hjx+PiPDGG28ArqXTp0+faott7dq1iEjx0qxZM0aNGkV+fn5xnSeffJKePXuSnp5OmzZtuPrqq4P2kZOTg4hwwgknVFtcpmbZtw/y82H37vB1XnkFnnsOnn4aJkyA/ftD1/vqK7jlFrjxRrj+enjppfD7fOghGD0aLr0ULrkEli8/rC+jcqqalMsxxxyjZS1fvrxcWbwNGTJEAb3ssst04sSJmp2drYD++9//rnTbu+++WwGdNm2aqqpOnjxZZ86cWW2x5ebmKqADBgzQ119/Xc8555ziWFVV77rrLgW0a9euOmHCBP3zn/+sxx9/fPH2a9asUUBTUlJURPT777+vttiSRSL8jsVTfr7qpEmqjz+u+oc/qP7lL+Hr3nijat26quCWv/0tfN0GDUrqgWpBQeh6f/97cL3zzgu/z1NOCa4bwZ9cSMB8jSAH+Z4Eo10iSZ6lv5EQ/pv1zDPB9a64Inzdinz88ccK6AknnFBc9v777yugJ598sqqq/uY3v1FAb7nlFm3btq1mZ2fr7NmzVbV88gS0d+/eQZ9de+212rVrV23evLlOnTpVVVWLior0d7/7nbZp00YbN26sI0eO1K1bt5aLL5A8zzjjjOLvV+AYu3bt0oyMDE1LSwtKigcPHix+f//99yugN998swL6yCOPRPeNSmI1KXnu3Kk6Y4bqo4+qHjoUus7SpcF/G337ht/fzTcH13344fB1GzcOrvvDD6HrTZkSXO/cc8Pv89RTg+t++GH4uhWJNHnaaXs1+vrrrwE4/vjji8sCp7dLly4Nqvv5559z9dVXk5eXx/jx4yM+xkcffcS1117Lzp07ue222wB44IEHePTRRxkxYgQ33ngjM2fO5Jprrgm7j/3797Nt2zbeeustANq3b8+yZcsoLCykZ8+etG/fvrhuSkrJr8jUqVNp2rQp9957L02bNmXKlCkRx20Sx549cO650KoVnHkm/O53sGpV6LqNGgWv79gRfr8NGgSv79oVvm6dMjdJHjgQWb1Dh8LvM6VMNtMYz6peo+7zTCbjx49n2LBh3Hvvvaxduzbi7W6++WauvPJKnn76aVZ5v/EzZswA4Jlnnimu98EHH4TdxwcffEDLli0BaNu2Lffffz/79u0Dwvcor1q1ikWLFnH22WezZcsWhgwZwttvv01ubq7d95hk6teHlSuh9F1Xc+dCt27l6zZsGLxe6vJ4OYHkmZrqtiubzEq76CKXxOvUcUtaWuh6ffvCww+7faamQpcu4fd5yy1uvyJu6d07fN3qYMmzGvXt2xeAuXPnFpd98cUXAOU6fpo1awZAnTp1OHjwYMTHKL3dIe/fsKpSp04dZsyYQWpqKkDxZ6EMGjSI++67j2bNmtGrVy/S09PZvXs3GRkZLF++nLy8PLKzs4v3k5KSUtzKfOutt4pbrOBao7feemvE8ZvEcMEF4J0oAfDFF66jpayGDeHii6FxY2jWDLKyXOsvVGK84Qa3pKW55FWRxx+PLM6uXV1SjMTpp0dWr7rU6OQZabP9yivdcrh+9rOfMWTIEGbPns0VV1zBgAEDeOihhwC46667Dv8AYYwYMYIFCxbw8ssvM3ToUJYvX05ubi7Dhg0LWb958+acdtppQWWZmZnceuut/PGPf2To0KFcf/31FBUVMX36dD777DOmTJlCZmYmr732WvE2F198sSXPJDVqFNxxB3ToACefDEOGhK6XmgqvvhrZPtPTqy++ZFCjk6cfpk+fzi233ML//d//8corr9CjRw8eeeQRTj311Jgd8/bbb2f37t1MnjyZt956i06dOnHVVVdVeT/jx48nKyuLp556iptuuommTZty1llnsWLFCpYuXcqZZ57J2WefXVx/yJAhzJw5k9WrV9OlovMpE3eFhe62nquvDt0K7NQJFi92p8V27390RGN9VTVGcnJytOxgyN988w09e/b0KSJTGyTL79iFF8LkyfDrX7skWvbapQlPRBaoak5l9ay33Zga5t13XeIEePNNOO648L3pJnqWPI2pQQ4ehP/93+CyNWsqftrHRMeSpzE1SEoKTJwIgwe79dRU+OQT6N/f37hqIuswMqYGEYHTToNTT4VZs2D+fBg0yO+oaqa4tjxFZK2IaKllkVd+oogsEZEiEflKRI6OZ1zG1DQicMopMG6c35HUXH60PGcDT3vvfxCRDOBNYC9wE3An8IaIdFXVyO8eN8aYOPLjmmcu8J6q/kNV/wUMB1oBT6nqU8ALQCfgZz7EdtgSeUi6gIkTJxYPS7dy5cqw9c4880xEpEqPj5a1Z88exo8fz6RJk4rLJk2ahIjwyCOPRL3f0myoPeMHP5LnpcBPIrJVRC7HJUqADd5rnvd6ZNkNReRKEZkvIvO3bdsWh1Cr7pxzzmHSpEmcddZZPPbYY+Tn53Peeefxn//8p8r7mjx5crUlmNKmTp1aPODH1KlTq33/pe3Zs4c//vGPQcnz5JNPZvLkyYwYMaJajzVgwABef/11TjnlFKZMmcL/et3Od999N9dddx0HDx7k0Ucf5aabbmLJkiXF23333XcsWLCAlJQU5s6dy7p166o1LlNDRTL0UnUtuFPykcDFuCR5APgdoMCFXp2rvPWxFe0roiHpxhPVcvQzR0c1lFWiD0mnqrphwwZNSUnRUaNGaZs2bbRPnz7FnxUWFurFF1+sjRs31hEjRujgwYMV0NzcXN26dav2799fMzMzNTMzU0866SRdunSpqqq+9NJLxeOCDhgwQLOysvRhbzyyDh06qPfzVEDvvvvu4voPP/ywTpw4UQGdPn26qqrOmDFDAX3ssccScqi9RB2S7oUXVD/6SPXHH/2OJPmRiEPSqep9qvqGqr4GTAFSKWlpZnuvbb3X3HjGVh2SYUi6adOmcejQIc4991zOOeccli5dynJvyO2//e1vvPbaa5x66qkMHjyYOXPmFG+XkpLCOeecw+OPP85tt93G4sWLufHGG4P2/f7773PVVVfRunVrxo0bx+LFi7n//vsB6NmzJ5MnT2bkyJFB25x//vnUqVOn+NLGG2+8QWpqKhdccIENtRehwkL3GObQodCkCfTs6UYsMrEVtw4jEekL3A/M9I57Ka6T6FNgK3CNiBQAlwNrgVnxis0Pfg1JN2XKFNLS0ujRowd79uxh4sSJTJ06lfHjxzNr1iwAHn74YTp37sw777zDZ599BkBRURHvv/8+X3zxReAsovifRcBll13GVVddRZ06dRg7diyffPIJF154IQAtW7Zk1KhRAJR+rLZ58+acfvrpzJgxg127dvHOO+9w2mmn0bp1axtqL0KLFgVPY1FU5IadM7EVz9727biW5p+A+sBy4E5V3Sgi5wJPAo8Dy4ArNAl72hN9SLr169czd+5cVJXepQY7nDJlSsjWbyBJAjzxxBPMmTOH6667jhEjRnD55ZdTUFAQMsbS20Uy4+RFF13Ee++9x6233sqOHTu4+OKLq/R1BdTWofb++9/g9YED/Ymjtolb8lTVTcAvwnw2G+hb7ce8O76DniT6kHRTp05FVbn99ts59thjAXjhhReYMWMGX3/9NaeccgpvvfUW48aN4/jjjw/6JxBIiLt27eLTTz8lLy+Pxo0bB+3/xRdfpF27djzxxBOICCeffDKNGjUiJSWF1atX8/e//52TTjqpXPy//OUvadCgAU8//TT169fnV7/6VZW+roDaOtRenz5w2WUuiS5bBt6P1sRaJBdGE3FJ1Angtm/frqNHj9asrCxNS0vTfv36FXfsqJZ0GM2bN09VVTMzM7VDhw6qGlmHUeCz3r17K97ETEVFRTpu3DjNzs7WevXqaa9evfTxxx8vF9uxxx6rIhLU6TJlyhQF9M477yzuMGrUqJGeccYZetJJJxV3GOXl5enAgQO1SZMmOnbsWO3Tp482btxYVUs6jMaOHVuuw0hV9dZbb9XMzEwF9NVXXw3qMAq45JJLFNALL7ywuCzSr6tsh1EoTzzxhPbo0UPr1q2rLVu21LFjx+o333yjgJ555plBdYcPH66Arlq1qtx+EuF3rCI//WSdRoeLCDuMbEg6c9gmTZrEmDFjePjhh7kl0mG/k5T9jtV8NiSdMcbEkA0MYg7b6NGjGT16tN9hGBNXNa7lmayXIUzis98tU1qNSp4ZGRnk5+fbL7mpdqpKfn4+GRkZfodSbOtWNwNm6SmETfzUqNP27Oxs8vLySNTn3k1yy8jIKL5PNBE88AD89a9u+LkOHdxTRklyd1WNUKOSZ926dZPmqRBjDsehQy5xgptie+1a92SRiZ8addpuTG3x1Vfly+zJoviy5GlMEurWzc2QefHFEHjQK8yDVyZGatRpuzG1RaNGMGqUWw4cgJ073WRvJn6s5WlMkqtTB7Ky/I6i9rHkaYwxUbDkaYwxUbDkaYwxUbDkaUyS2b0bqjB+tokRS57GJJmbboJ69aBrVzj9dIhiYlZTDSx5GpNkvvvOzVm0ejV88IFriZr4s+RpTJLJLTOvrD2R7A9LnsYkkUOHyrc0O3b0JZRaz54wMiaJpKTA5s3w00+uBbpuHTRo4HdUtZMlT2OSUKNGcNRRbjH+sNN2Y4yJgiVPY4yJgiVPY4yJgiVPY4yJgnUYGZNE/vAHaNkSjjzS3d/Zo4frgTfxZ8nTmCSxZw/ce2/JekqKmznTkqc/7NtuTJJYuzZ4vV07qFvXl1AMljyNSRr2WGZisdN2Y5JEly5wzz0uiebmwrHH+h1R7WbJ05gk0b07/P73fkdhAuy03RhjomDJ0xhjomDJ0xhjomDJ0xhjomDJ05gksGcP3H47rFzpdyQmwHrbjUkCl1wC06fDtGnu5vjBg+FPf/I7qtot7i1PEckQkZUioiIy0SvrKSJzRKTI+2xYvOMyJpHt2OFe16yBWbNsxsxE4Mdp+11AdpmyyUAP4GZgPzBNRBrHOzBjElVeXvD6lVf6E4cpEdfkKSL9gJuA8aXKBgBHAZNV9UngMaARMDKesRmTyM44A37xCzjiCKhfH375S78jMnG75ikiKcDzwJPAvFIfBZ7Q3eC9Bv7HHhmn0IxJeH/9q98RmLLi2fIcA3QEXgHaemWNgbLjwoj3qmV3ICJXish8EZm/bdu2WMVpjDGVimdvezugBbC4VNnFQBvvfeA6aCCxlhlDBlT1WeBZgJycnHLJ1Rhj4iWeyXMqsNR73xt33fN94PfAi8AoEVkGXAMUAG/GMTZjjKmSuCVPVV0OLAcQke1e8RpVXSAiF+Kuhz4GfA+cp6o/xis2Y4ypKl9uklfVWZRc20RVlwHH+xGLMYnu4EFITfU7ClOWPWFkTILr3BkaNoT+/d1y9dWQmel3VEZUk7PfJScnR+fPn+93GMbE1LZtbrbMgLQ0KChwryY2RGSBquZUVs8GBjEmgS1eHLzeu7clzkRhydOYBLZiRfD6UUf5E4cpz655GpPArrgChg93z7bn5UHHjn5HZAIseRqTwNLTXYdR585+R2LKstN2Y4yJgiVPY4yJgiVPY4yJgiVPY4yJgnUYGZOgNm2CTz+F7Gw3b9ERR0Ad+4tNGPajMCZBffEFnH9+yfovfgHvvedfPCaYnbYbk6DWrw9ezy4785fxlSVPYxJU2Unf2rXzJw4Tmp22G5Ogevd2E72tX+8Safv2fkdkSrPkaUyCGj3aLQFJOgBajWWn7cYkCZHK65j4seRpjDFRsORpjDFRiOiap4jUBxoAO1W1SETOAE4FFqvqK7EM0BhjElGkHUbPA+cDx4tIY+BdQAFEJEtV/xKj+IwxJiFFetp+DPCTqv4XOM8rW4ebAXNMLAIzpjZ7/XU44QQYORJuuAH++U+/IzJlRdrybAus9t73A75R1d4i8h3QMRaBGVObrVjhHs8MaNTIPZ5pEkekLc+DQIb3vhvwtfd+J6XmXzfGVI+NG4PX27TxJw4TXqTJcw3QVUSWAo2ABV55G2Bj2K2MMVHZtCl43ZJn4on0tP2vwCSgF/AD8KqI9AZaAB/GJjRjaq8nn4S1a10LdONG6N/f74hMWRElT1V9RUS+BroCn6rqZhFJAYZTci3UGFNNOna0mTITXcTPtqvqQmAhgIikA1nAZ6q6O0axGWNMworomqeI3CAi/xSRo0WkHfAtsAjYICLHxTRCY4xJQJF2GF0AnIZLmlcA7XC97I2Au2MTmjHGJK5Ik2cXYJ2q7gKOA/Jxj2f+BBwdo9iMMSZhRZo8G+Lu6QR3n+d8VZ2F6yxqHIO4jKm1vvsO5s6Fdetg3z6/ozHhRJo8twC9RORW3Cl74Cb5ZrhWqDGmmjz/PBx/PHToAOnp8OCDfkdkQok0eb6He8Lofm99hog0AbKB5bEIzJjaquwN8s2a+ROHqViktyqNA/bj7vN8V1Vni8hA4E3gnVgFZ0xtZI9mJodIb5LfDfy/MmXzcL3wxphq1K0bbN/ukuiWLdC2rd8RmVAivkneexzz90COVzQPuE9Vl8UiMGNqqwkTSt4fOGBzFyWqSEeSPwr4DKhPyShKRwIjROREVV0So/iMqdXq2Py2CSvSDqM/AplAAW4U+Xe995neZ8YYU6tEmjwH426I76mqv1TVX+JGWCrwPouIiHwpIgUiskdE5ovIEK/8bBFZLSKFIjJLRDpV9Qsxxph4ijR5NgA2qGrxTRSquhHIw91AH6k5wA3APUB/4HkRaQ38A5ecx+Gm/Hi5Cvs0xpi4izR5rge6i8j1ItJcRLJE5Aagh/dZpG7GnfL/GygCDuF67NOBB1R1AvB/wGAR6VyF/RqT9FThhRdgzx6/IzGRiDR5TvPq/hX3tNFWIDBj5pQqHK8xsA34EtgHjAUCp+gbvNc87/XIshuLyJXe6f78bdu2VeGwxiS+N96AsWNh2DAYMwYeeMAlVJOYIk2efwJm4XraSy+zcKfgkdoFDMOdumd4+y0r0Jtf7tdGVZ9V1RxVzWnRokUVDmtM4nvtNff6+ecwaRK8/LLdppTIIr1Jfi9wqogMBQZ6xfNU9aOqHExVD+Cm7fhQREYCpwAzvY+zvdfALcG5Vdm3Mcls7Vp4993gshtu8CUUE6Eq3UXmJcuPAEQkRUTu8Mrvr3BDV/903Jzvc3CDi5yAuwTwGnAvcKuItAJ+hRuhfk1VYjMmmW3aBKed5l43bYIWLeCaa/yOylRENMqLKt5UHHuBQ6paaRL2noV/CeiM6yxaCPyvqs4TkXOAh3Gtzy+BMZUlz5ycHJ0/f35UsRuT6A4csBvk/SIiC1Q1p7J61fHjieiqjPcsfJ8wn00HpldDLMbUCJY4E1+kHUbGGGNKseRpjDFRqPDkQES+rejjao7FGGOSRmVXVrrg7re0RGlMjBQWurmKRKBhVR52Nr6qLHnOIcTN6sbUJLm5sGIF/PijW445Bo49NnTd3/wG1q+H/fvdMm0atGtXvt7nn8OIEXDwIBw6BCeeCO+/H3qfo0bB22+7EeMvuADat7d7PJNBhclTVU+KVyDGVLeiIli1Cr791r2/IMy8B6+8AuPHl6z//vfhk+cXX7h9BuzeHbrewYPwww+V1wPX8gQ3cvyjj7pkbMkz8dkNEaZGmj3bPSNeVOTW27ULnzybNAle//HH8PtNTw9e378/dL3U1OD1gwfD77N16+D1o48OX9ckDkueJmnt3QtpaeUTFUCvXiWJE9yp9t69UK9e+bpVSZ5pacHr1ZE869eHBg0gJQX693etT5P4LHmapPPll25u82nTYMoUOP308nWaN4fevWFZqRm2Vq+Gvn3L1+3UybVSmzSBxo3dnOnhTJzohoyrW9ct3bqFrnfMMW4St9RUt9StG36fTz3lFpNcLHmapHLjjfD44yXrr70WOnkCDBniWpG9erkkV79++HpDhkR2/IoSa2l160JWVmR1TXKK+tl2v9mz7bXTjBmuFzsgM9NNz5uZWb5uUVH5a5Q1maqyfc928n7KY2PBRjYWbGT7nu3s2LuDXft2sffAXrfs38s1OdcwvOvwcvtYuX0l571xHqqKejfaRPr+1z1/zQNDHwgZW4+JPdh/KMw1jgp0y+rGzItmhvzs11N/zcJNCyvc/rcDf8stJ9xSpWPG89l2Y+LmF7+A7t1h5Uq3npHh3ofqZKnpibOgqIC5eXNZsGkBc9bPYd7GeWzetTmibYd3KZ84AQoPFLJkS3ST4W7eHf7Y3/3wXVTJs37dMKcLwMaCjeT+WPHIlT/s/aHCzw9H2OQpIvsi3Ieqag3/NTWJIiUFfvc7N+r6FVfA8OGhW5013ddbvubEF0+kYF9BVNsfOHQgZLkcxujLyXoWG62Knm2vU4XFmLgZOxb+9S8YObJ2Jk6APi37cHaPs6PePmzytIcJI1ZR4rsvblEYU4qqu11nzJjQnS42NYVrIT474llWbF/BvI3zisvr1alHxyYdaduoLW0btqVlZkuaZjSlUXoj6tWtR7069ahXtx79WvULud8uzbqw6KpFxccIJNNI3jfOaBw23m+u/SaqrzMtNS3sZ2+c+wZFB4vCfg7QJKNJhZ8fDuswMglnwgT3hE3z5i6JXnKJJcxwcn/I5c7/3MlJ7U8ip00OA1oPoG5qBfdFmUpF2mFUpeQpIu2ANkDQLcCqOqfKER4mS5410/LlcNRRbiT1gNtuczNJGhMP1drbLiItcfOpHxfiY410P8ZUJiXF9aDv2uXW69WDyy/3N6ZEsO/gPu6dfS/pqencOeROv8MxRD4Y8oPA8ZSfejiwGFMtevSAl14qWX/6aejSxb94EsVDnz3EPbPvYfGWxX6HYjyRJs/TcS3MwHx+y4E/ADuAC2MQl6nFRo50tyP95jdw6aV+R5MYPsp1s3wf1IMUHihk2+5tPkdkIk2eLYCVqvqMt75LVe8DtgEjYxKZqdUefBBefNE6igJW71gNwLbd2xj26jB+KIzdzd8mMpFeq9wD7Cv1vpOINAeygP+JRWCmdrPZI0vs3rebjQUbAfh03aekSAodm3T0NygTcctzA9Dee78SaA5swSXPnTGIyxjjmbZ8WtB6xyYdK7z/0cRHpMnzX8AWEekLTCC4s2hCjGIzptZTVe6ZfU9Q2YhuI8LUNvEU0cmRqt4M3Oytfi0i3+FuW1qiqmFmZjEmMpMmufl9zj8f/ud/Qg9YXFutzF/Jdz98V7yeKqlcd+x1PkZkAiJqeYrIHSJS3O+pqrNV9c/ADhEZFrPoTK3w8stuUONzzoFWrdxkaMb5bN1nQetDjxxKl2Z271YiiPS0/V7g6hDljwP/rL5wTG1z4ADMKfV8WkGBGwHeOCu3rwxaH9R2kE+RmLIqPG0XkTalVtNE5AhKborPxHUi2c0kJmrr1rk5ywNatrSb4kv7dse3Qevdm3f3KRJTVmXXPNd7rwoMAPJC1NlSrRGZWqVRI3jiCTc98KpV5Sdjq+0Wbw5+oqhbVphJk0zcVZY8A61KJXwL84XqC8fUNs2bw/XX+x1FYlq/cz3f7/y+eD0tNY0+Lfv4GJEprbLkeYX3+hywBveMe8AeYIWqVjyJiDEmKp+v/zxofWCbgWTUyfApGlNWhclTVV8AEJFTgFWBdWNM7JXtaT+p/Uk+RWJCifQ+z4tFpJ6IXAUExrmbB7yqqntjFp0xtdjybcuD1o/PjnDeYxMXkY7nmQV8CpTu6rsM+H8iMkRV82MRnDG12b8u/he5P+ayYvsKVmxfwbFtj/U7JFNKpMMv3Af08N4Hputr6JXdA/y2muMyNdxf/wrvvAMDB7pl8GB3g7wpUTe1Lt2yutEtqxtndT/L73BMGZHeJH8GsB8YoaqNVbUxcBZwELAHbU2VbNgAN90EH38Mf/4znHuuS6TGJJNIk2cr4FtVfS9QoKozcCMstYxFYKbmuuOO8mU5lc4YY0xiiTR55gOdReSoQIGI9Ae64EaTNyZi3bvDWWfBgAElZX3s9kWTZCJNnh8BGcB8EVkqIl/jetvTgA8j2YGIdBWRj0UkX0QKRORDEensfXa2iKwWkUIRmSUinaL5YkxyuOMON/jHpEkwbBjMnQt1bbZck2Qi7TC6ExiKO33vVap8s/dZJNrikvXdQDfgeuB5EbkA+AduXqRxwP3Ay8CQCPdrklT37vCvf/kdRWJ6dM6j7Nq3i7aN2pLdKJsT2p1Ao/RGfodlSon0Ps913mn69cBAr3geMEFVI322fY6qnhxYEZGLgN7ABUA68ICqThORgcAlItJZVddE+oWY5JOe7ncEieuZBc+waseq4vXFVy+mX6t+PkZkyqpsVKW7gDxVfdFLkr+P9kCqWjx2jojkAM2AN4HAKfoG7zUw+MiRuEdCjalVVJUNBRuCyrIbZfsUjQmnsmue44Gx1XlAEekOvA2sxbVky1XxXjXEtleKyHwRmb9tm029amqmnUU72bN/T/F6Rp0MmmY09TEiE0qkHUbVQkR6AZ8AB4BTVXUTkOt9HPjX2tZ7zS2zOar6rKrmqGpOixYtYh4hSc9hAAAVU0lEQVSvqX4HDvgdQeLL+yl45Me2DdsiNgdzwonkmme6iLSjgkGPVXVdZTvx9jELd7r+e2CQiAzCdRY9CNwqIq2AXwGf2fXOmmn4cNi1C7Kz3TJ6NBx1VKWb1SobfrJT9mQQSfLsjzvFDkcj3E9nINBcfKB4Y1XxetwfBh4BvgTGRLA/k2R27oSPPgouu/JKf2JJZGWvd7Zt1DZMTeOnSG9VOuxzBlWdFW4/qjodmH64xzCJbUNwTqBDB+jZ059YElm5lmdDa3kmokiS5wZstHhTDUIlT1NeuWue1vJMSJEkzzxV/WPMIzE13tChsG0bbNzoEmmGDYoeUrnT9oaWPBNRpKftxhw2ETdnUfPm0M/u9w7L7vFMDpXdqrQO2BSPQIwxjp22J4fK5jDqGKc4jDFA0YEitu/ZXryeIim0btDax4hMOHG9Sd7UbocOucWE9832b4LWWzdoTZ0Uu7qWiCx5mrh46SWYMMHvKBLfe9++F7TeqYmNzpioRLXcI+RJIScnR+fPn+93GCYCW7ZAr17uhviBA92o8e3auQ4kU96OvTtYlb+Kb/O/pWm9ppzZ7Uy/Q6pVRGSBqlY6t4GdD5iYe/992LEDHnywpGz1aujc2b+YElmzes0YlD2IQdmDSNbGTW1gp+0m5l5/PXi9TRtLnJGyAUESlyVPE1OqbuCP5s1Lyv75T//iMaa6WPI0MSXiphfOy4N//AOuu85GUTI1g13zNHGRng7nn+8WU97e/XtJTUklLTXN71BMhKzlaUwCmPDfCaTfm07zPzen79N9eXre036HZCphLU9jEsCmAvcUdP7efPL35lOwr8DniExlrOVpTALYuGtj0PoRDY7wKRITKUuexiSAQMsz4IiGljwTnZ22m5hZuBCefRbat3dLnz7W0x7Opl3BybNNwzY+RWIiZcnTxMzChfC3v5WsX3gh/P3vsTnW7n27KTxQiKIc0kOoKoqi6q1770OVZTfKJr1Oerl95u/JZ2PBxqBtItW+cXuy6meFjPPb/G/LlZdredppe8Kz5GliZl2ZOVXbtw9db9/BfeT+kMviLYtZsmUJm3dtZtuebWzfs52CogJ27dvF7v27qVenHmtvXBtyH9fNvI5JiyZFFeeiqxZxVOvyTeLJSydz/czro9rn8yOe5/KjLy9XvmzbMgY9P6jCbdNT02mS0SSq45r4seRpYub774PXA8lTVZm5eibzN85n5uqZfJn3JUrlz3Bn1Ak/b4ccxhyF4Y59OPs8HEc0PMIey0wCljxNzNxwAwwaBMuXwzffuJkyfyz8kbHvjOXNb96s8v4KDxRy8NBBUlNSy312WMkzzOAbfiWwE9qd4MtxTdVY8jQxM2CAW0o77ZVf85/c/0S9z937d9MovVG58gZpDWia0ZQUSUFEEAQRceve+3Bl4Z7qyaqXRe8WvcttE4lm9ZqFLK9ftz79W/cPu91RrY7i0WGPRnQM4y8bz9PE1bqd63h87uM899VzQTeCH9HgCHq26MnRrY+mc7POtKjfghaZLWic3pjMtEwy62bSIK0BDdIa2CmtialIx/O05Gl8sbNwJ0NfHcrQTkO5bMBldM3q6ndIxgA2GLJJcIf0EJ+M/oT6dev7HYoxUbHkaXzRtF5Tv0Mw5rDY45mm2m3c6Kbe+GjJcm56/+aQN4Ubk+wseZpq9+GHMHw4/Pzy2fz1y7/QfWJ3fv7qz/lk7Sd+h2ZMtbHkaard6tXem95Ti8s++u4jPl//uT8BGRMDljxNtVuzBhhyD3T6OKj8xHYn+hOQMTFgHUam2nXsvJ+Uzn+m9DAag9oOYnCHwb7FZEx1s5anqXZDL/uUQ3V2BZU9f9bzpIj9upmaw36bTbWb+N+JQeu/6PoL+rTs41M0xsSGJU9TrfYf3M8Haz4IKrt24LU+RWNM7FjyNNVq4eaF7N6/u3i9VWYrhncZ7mNExsSGJU9TreZtmBe0PrjDYBvIw9RIljxNtVqwaUHQes4RlY6vYExSsuRpqtXb84KTZ+f6x/gUiTGxFbfkKSJPiMgWEVERmVGqvKeIzBGRIhFZKSLD4hWTqV579+9lR+qyoLK2KUf7FI0xsRXvluc/QpRNBnoANwP7gWki0jiuUZlqMW/9Ykg5WFKw40gG9Ag9oroxyS5uyVNVbwD+UrpMRAYARwGTVfVJ4DGgETAyXnGZ6rN1cyqsPBMK3LS59X48hozwc7YZk9T8fjyzk/e6wXvN816PDFVZRK4ErgRoH24eW+ObM48eyBc3vMuaNbBw9Ua0216/QzImZvxOnmUF7mkJOTeIqj4LPAtuGo54BWUik5EBxx3nloto43c4xsSU373tud5rtvfatky5McYkpLi1PEXkDCDwgHM7ERkLfAIsAUaJyDLgGqAAqPqk3sYYE0fxbHmOAx703vcDngNOBC4EVuI6i9KA81T1xzjGZYwxVRa3lqeq/qyCj4+PVxwmNlTVHsM0tYrf1zxNDXHje7fR4sG2HP3Ucfx6yrl89N1HfodkTEwlWm+7SVILc9eyvWgj27dtZOG2L1n34QiGPu93VMbEjrU8TbXYULA+aL1pSjufIjEmPix5mmqxfV9w8uzcwpKnqdkseZrDduDQAQpkY1BZr+zsMLWNqRkseZrDtrFgI1pqrswW9VswdrQ91G5qNkue5rCt3xl8yt6+cXvq1fMpGGPixJKnOWzrfwpOnu0a2/VOU/NZ8jSH7dv8b4PW2zWy5GlqPkue5rDsO7iPP33yp6AyS56mNrDkaQ7L9f+8noN6MKhsYNuBPkVjTPxY8jSH5e2Vbwetd0oZQsbmk32Kxpj4seRpolZ0oIgtu7cEleU+9jLvvWcDhJiaz5KnidqmXZuCCwobwY8dsRlSTG1gA4OYqLXKbMWs38xi1eaNXHfHBooKhXr14Oyz/Y7MmNiz5GmiVq9uPU7ueDInd4Q1g+DBB+HBx6FFC78jMyb2LHmaanHTTdCwIVx/vd+RGBMfljxNtWjZEu64w+8ojIkf6zAyxpgoWPI0xpgo2Gm7icrizYt57qvnyG6UTXajbHq36M2AIwb4HZYxcWPJ00Rl4eZFPDnvyeL1bkWjWPrHydSt62NQxsSRnbabqDzz3+DZ3XKXtKWO/Ss2tYglT1NlCzctZO6mz4LKmqZmY9O2m9qk1rYVVmxfwbsr3yW9Tjo3DLohZJ2XF73M2h/XRrX/W0+6lYw65aei+HDNh8xZPyeqfV561KV0atqpXPnSrUt5c/mbUe1z6JFDObH9ieXKd+zdwYQvJ4TcZvwn48uVHXloeFTHNyZZ1arkuWvfLl5Z/Apvr3ybD9Z8AEDnpp3DJ8/FL/Px2o+jOtaNx90YMnl+sOYDHvnikaj2ObjD4JDJc9nWZSETWiTq160fNnlGus/2af34wzXdozq+McmqViXPi6ZfxDsr3/E7jBrn3ctepV8rv6MwJr5qzTXPeRvmhUyca35Yw8aCjSG2MJE45ohj6Nuyr99hGBN3tabl+d6q98J+NurMIxCgeXN4s9Slw0uPupTB7QcDsHARLF5c8ln//tD/qND7mzIVzhyeToo3G+/06ZCV5d4PPXIo9evWB2D9evj3f0q2a9cOTjs19D4/+jfcdnVH6hW69QcegBNOcO97t+zNXUPuAmD3Hpg2rWS7+vXhvHND73PhIvj7QycwY6dbHzMGRo9275tmNC3eZ+Br2ru3ZNvzzoOOzVtyUb+LEOspMrWRqiblcswxx2ikDh06pP2e7qeMp3jp+1Rf/dun05S6uxVUQbVNm/D7uOsuLa4HquPHh6/bqlVw3U2bQtd7++3geiNGhN/nGWcE13333dD18vKC68X7azIm2QHzNYIcVCtanku2LGHJliVBZW+c9waZhd1gv09BGWOSWq245rl069Kg9d4tetMtq5tP0RhjaoJa0fJc/9P6oPWfH/lzwF3jnDWrpDwtLfw+xoyBU0tdj+zQIXzd6dNhf6kWbbNmoeudcELw8QPXRUN54AEYN65kvU+f0PX8/pqMqS1qRfKcMTs4ee7Z7OYVT0+HkyOc6LFjR7dEItCRU5nmzSM/ft8IO7T9/pqMqS1qxWn7+oJ1QeuFW7N9isQYU1PUiuS57dC3Qeu92hzpUyTGmJqixifPov372ZvxXVDZ8IHWWWSMOTw1/ppnaqrw6eUf8/nKlcz6eiVrt2+hX/dGfodljEly4u4J9Z+InAg8DXQHlgFjVfWrcPVzcnJ0/vz58QrPGFNLiMgCVc2prF5CnLaLSAbwJtAQuAloBbwhIqm+BmaMMWEkRPIEhuMS5lOq+hTwAtAJ+JmfQRljTDiJkjwDg1Ru8F7zvFfrFjfGJKRESZ5lBYbpCbogKyJXish8EZm/bds2H8IyxhgnUZJnrvcauHu9bZlyAFT1WVXNUdWcFi1axC04Y4wpK1FuVZoJbAWuEZEC4HJgLTDLx5iMMSashGh5qmohcC6wC3gcl0jPVdWDvgZmjDFhJErLE1WdDdh8DsaYpJAQLU9jjEk2ljyNMSYKljyNMSYKCfNse1WJyDbg+ypu1hzYHoNwqluyxAnJE6vFWf2SJdaqxtlBVSu9FzJpk2c0RGR+JA/8+y1Z4oTkidXirH7JEmus4rTTdmOMiYIlT2OMiUJtS57P+h1AhJIlTkieWC3O6pcsscYkzlp1zdMYY6pLbWt5GmNMtbDkaYwxUagVyVNEThSRJSJSJCJficjRPsbyhIhsEREVkRmlynuKyBwvxpUiMqzUZ3GPX0S6isjHIpIvIgUi8qGIdPY+O1tEVotIoYjMEpFOpba7SkTyRGSviLwtIllxiPVLL8Y93nivQxIxTu+4Gd7PV0VkoleWUD9777hrvRgDy6LK4qno+x3DOJuIyCsi8qOI7BKR2XGLU1Vr9AJkAJtxY4P+Fjda/XdAqk/xPIEbOUqBGaXKFwE7gGuBpcBOoLFf8eOmQPkEuM6LWYGPgdZAIfAVcD1QAMz2thng1fsQ+F/gAPBKHL6nfwHGALd7x/w2EeP0jn0/sNs7/sRE/Nl7Ma31fv6jvOX0iuKp6Psd4zinAweBR4GxwIvxijPmvyx+L8CvvF/Ucd76n7z103yMqWPp5Fnqj/lJb/0yb/1yv+IH0sqs5+OGCrzJO/65Xvkr3npnSv4pDPQ+mw3sBzJiHKvgniI51ktMKxI0zn7AXmBcIHkm4s/eO9ZaYBLQsFRZ2Hgq+n7HMMYjvWO8BqTh/VOJV5y14bQ9GeZHqihGX+JX1X2B9yKSAzTDJZmqxloHaBfLWHGttG3Al8A+XAskoeIUkRTgeeBJYF6pjxLuZ1/KpcBPIrJVRC6vJB4/Yu3lvQ7E/dPcLSIPxSvO2pA8ywo5P1KCqSjGuMYvIt2Bt3EtkeurGE+8Yt0FDANuwJ2y/amKscQjzjG4M45XKJlmpjFQtwqxxPNn/xxwHnAJ7h/SM6WOH0k88Yg13XvNBM4HPsddhik7TnFM4kyYwZBjKKL5kXxWUYw7KvgspkSkF/AfoAg4VVU3iUhFsZb+bKP32QFK/rvHhKoewF2//FBERgKn4KZ2SaQ42wEtgMWlyi4G2lQQp28/e1W9L/BeRAYAN1Py/QkVT8MKPouVtd7rp6o6XURaAKdSkhBjG2esr534veBaIlu8b841uOZ6Lv51GJ0B3Ir7T7cYd4rZ1XtfutPgJ6CJX/Hj/ti34pLKbZR0HByBS6YLKLng/qm3zTGU74h5NcZxng68gLtGOB537XJzAsbZCxjpLXd7x5/pxZJoP/u+wLu4zpYbcJdE9uASfch4Kvp+xzBOAZZ4v6dXAHO9n2WfeMQZsy8skRZgCPA17vRjIZDjYyyzvD+c0stooDfwhfeD/Rb4Hz/jx/W2l41Tvc/OAdZ4sc6m1MV2Sno3C70/wOYxjnOgl3D2Aj/i7ggYmGhxhvneBnrbE+1nfwTwT9wwbnuA+cDplcVT0fc7hrEGvneF3vfuwnjFaY9nGmNMFGpjh5Exxhw2S57GGBMFS57GGBMFS57GGBMFS57GGBMFS54mrkTkZ6VG6vmZ3/EYEy1LnqZahBjCLNQyHncD+Jfe8pPPMY/0hoLb4Q1Nlyci/xaRS0vVsWRvQqoNj2ea+FiIe7IH3KNvgcfeFuFuRgbIU9WvgOPiHFs5IvJLYJq3uhn3qF8b3ON963HPoBsTlrU8TbVQ1V+p6nGqehxu9KCA4nJVfT5US05ExpcqGyYi33gDG08VkYYicpOIbBSRbSIyQUSK/+mLSJqI/MEbRLhI3ODNr4tINhW7wHv9BGijqkeramugC24oNryW8seltvnYizHwuYjItSKy2Gu57hSRd7wxAQLxjS71tf1SRD73BuFdJSK/iuJbbRKEJU+TaN7APbNcDzgX97zyfbjHL5vjBmceU6r+m7hRlLoAK71tLwA+F5GmFRwn8LvfHbhURLqIiKjqGlWd5X2WB3xTaptvcJcb1njrT+DG5OyHe3Z6LzACmCMioYY4mwJk4R4Z7AJME5G+FcRoEpglT5NorlHVHrjhxcANqPFz3OAp33tlpwGIm27jTK9suKr2w43LuB1oj3t+PZxncSOQt8a1NFcBm0XkRRFpD6Cqz5fZx2+9FvQ9ItIRN5AHwFWq2gvoACzDDTV3e4hjPuJ9bb1w40+mArdU9M0wicuSp0k073qva73XH1T1c1U9REnybOW9Diq13b9ERIEfcC1UqODaqqp+hBtU5EVKBsZtiWvVfioi9SuJcyAlQ5894x27EDdQRbhjT/GOnUfJP4c+lRzHJCjrMDIJRVUDPfAHvNfSPfKBUWykzCvAfyk/oO26So61EDeUHV5L8k+4wX/b46bH+DzctmWOvRiXOEvbWNGxKT+wsEkyljxNMvtvqfePqeoUcB05wGDcEHUhich1uBbnDFXdr6prReRTXPKEkqS9p9RmmaXez8MlawEmq+pDpfZ9DCWjnJd2LvC1iLQBjvfKllb8JZpEZaftJml5HTuBEeP/ISLfisjXuNknPwEqmqb3JNzMi7u9nvpluKkmwN1eFegoWoMbYBngFRGZKyIjVTUX+JtX/qCIfO/1uu/AjX9ZPH1wKTeLyDfevhsAh3CzPpokZMnTJLtf4UZmX4HrsMnGTTP7KG7g6XCew3UUrcFdQ+2Ou9/zVWCEuqk9UNV83Gjq64GmuOusrb19XOd9thh3vbQTsAl4GncXQFnn4UZlT/eOe76qLqnyV2wSgg2GbEwMicho4CVvtZOqrvUvGlOdrOVpjDFRsORpjDFRsNN2Y4yJgrU8jTEmCpY8jTEmCpY8jTEmCpY8jTEmCpY8jTEmCv8fyA9+zK4t+Z0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T_list = range(T)\n",
    "plot_list = range(num_comp)\n",
    "plot_list = range(2)\n",
    "Loss_list = cumsum_cost_pca_mean\n",
    "plot_total_loss(plot_list,T_list,Loss_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### all algorithms' comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAFECAYAAAB8sRYTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FFUXwOHfDSEJhJ7Qu4hSLCAIggJSBKVIVxQkoPSigDRRBD4hoBTpVRBQiIAUEQTpoiBVQQERaUJA6dJSSDnfH5NdMlvCEpLdTXLf55mH7J3Z2ZOwOZm9c++5SkTQNE3TUo6PpwPQNE1Lb3Ri1TRNS2E6sWqapqUwnVg1TdNSmE6smqZpKUwnVk3TtBSmE6umaVoK04lV0zQthenEqmmalsJ8PR1AaggODpYSJUp4OgxN09KZ/fv3XxaRvPc6Ll0m1hIlSrBv3z5Ph6FpWjqjlPrbleN0V4CmaVoK04lV0zQthenEqmmalsLSZR+rIzExMYSHhxMVFeXpULR0KiAggCJFipA5c2ZPh6J5WIZJrOHh4WTPnp0SJUqglPJ0OFo6IyJcuXKF8PBwSpYs6elwNA/LMF0BUVFRBAUF6aSqpQqlFEFBQfoTkQZkoMQK6KSqpSr9/tIsMlRi1TRNs7h9+zZ//+3SsNT7phOrG125coWOHTuSN29e/P39qVChAsuWLXPpucOHD0cpxddffw0YV0ePPfZYisV2+vRplFLWLU+ePLRp04YrV65Yj5k2bRply5bF39+fQoUK0a1bN9M5KleujFKK6tWrp1hcmpbS7ty5w9y5c3nk4Ydp26wZqbHun06sbtSiRQvmz5/Pyy+/zIQJE7hy5QqvvPIKW7Zsue9zhYWFMW7cuBSPsWLFiixevJjatWuzZMkSBg4cCMCwYcPo1asXcXFxjB8/nr59+/Lbb79Zn3fy5En279+Pj48Pu3bt4syZMykem6Y9qF9++YUSxYrRqVMnzv/7LzsOHGD50qUp/0Iiku62SpUqia0jR46YHoN5c2bWLPNxnTs7PzYpW7duFUCqV69ubVu/fr0AUqtWLRERCQkJEUD69+8vhQsXliJFisj27dtFRGTYsGECyLJlyxLiR8qXL2/a17NnTyldurQEBwfL0qVLRUQkOjpa3n33XSlUqJDkzJlTWrVqJRcvXrSL79SpUwJIo0aNrD8vy2vcunVLAgICxM/PT/7++2/rc+Li4qxfh4aGCiD9+vUTQMaNG5e8H1QaZ/s+07zH1atXpVi+fAKYtpLBwRIZGenSOYB94kIO0lesbvL7778DUK1aNWub5SPzoUOHTMfu2LGDbt26ER4ezvDhw11+jU2bNtGzZ0+uX7/O4MGDARg9ejTjx4+nSZMm9OnTh3Xr1tG9e3en54iJieHSpUusWrUKgGLFinH48GGioqIoW7YsxYoVsx7r43P37bN06VJy587NyJEjyZ07N0uWLHE5bk1zhw/79uXMxYumNj/glYgI4q5dS9HXyjDjWNOS4cOHU79+fUaOHMnp06ddfl6/fv3o0qULM2bM4K+//gJgzZo1AMyaNct63IYNG5yeY8OGDeTLlw+AwoULExoayp07dwDnd73/+usvDhw4QLNmzbhw4QI1a9bkm2++4dSpU3pMp+YVdu/ezdQFC0xtJYAtuXNTcuRICA5O0dfTidVNHn/8cQB27dplbfv5558B7G5C5cmTBwBfX1/i4uJcfo3Ez4uPjweMrh5fX1/WrFlDpkyZAKz7HKlatSqjRo0iT548lCtXDn9/f27fvk1AQABHjhwhPDycIkWKWM/j4+NjvTpdtWqV9UoXjKvYQYMGuRy/pqWG//77j+efe87UVgI4WqsW/t98AzlzpvhrZtiuANteVme6dDEfN3t28l7v+eefp2bNmuzYsYPOnTszffp0OnfuDMCHH36YvJO6oEmTJsTGxrJgwQLOnDnD+vXrTVevtoKDg6lbty4VK1bE398fgMDAQAYNGsSdO3eoV68e06ZNY8KECdSsWROAJUuWEBgYyMqVK61bYGAgS1PjpoCm3Yf4+Hi6tmhBVGysqf1/pUrhv2VLqiRV0FesbrVixQr69+/PypUrWbhwIWXKlGHcuHHUqVMn1V7zvffe4/bt24SFhbFq1SpKlixJ165d7/s8w4cPJygoiOnTp9O3b19y587Nyy+/zNGjRzl06BCNGzemWbNm1uNr1qzJunXrOH78OA8//HBKfkua5rJP+vVj6datpraHgNe//BJ8Uu+6UkkqjOHytMqVK4ttoes//viDsmXLeigiLaPQ7zPv8dvOnVR+9llibNovdu1K3pkzk3VOpdR+Eal8r+MybFeApmnpl4jQs3lzU1INAv56661kJ9X7oROrpmnpzuZJk/jJZmjVzDp1eHjOHLe8vk6smqalK7eOH6dn//6mthezZaPVmjXgpkI5OrFqmpZ+xMbybs2aHLMZpjhi7FjIksVtYejEqmlaurEpJITZ//xjamterBhVbAoGpTadWDVNSxcujR7Nm4sXm9oezpKFubt3uz0WtyVWpVRppdRWpdQVpdRNpdRGpVSphH3NlFLHlVJRSqltSqmSiZ7XVSkVrpSKVEp9o5QKclfMKc2bywZaTJ061Vo68M8//3R6XOPGjVFK3deUW1sREREMHz6c+fPnW9vmz5+PUirFKnfpcogZQ9yWLTQbMoSzidoUsHDZMnIXKOD+gFyp1JISG/A88APQC5iMUVlmK1AAiAJ+AXoDN4HtCc+pmHDcRmAgEAssvNdruVLdyhNq1qwpgLz55psydepUKVKkiACyefPmez7XtrpVWFiYrFu3LsVjrFGjhvj4+Agg//vf/5we16hRIwHk1KlTyX6tS5cumap7iYicPHlSwsLC5OjRo8k+b2KWql0VK1aUxYsXS4sWLaz/ByIiH374oQBSunRpmTJlinzyySdSrVo16/NPnDghgPj4+IhSylTdyxFveJ9lOP/9J5Nz5bKrWtWvdesUfylcrG7lzsTqZ/P4CnAR6Jvwg2id0L4w4XEpYFLC108n7NsOxAABSb2WS2UDh5Os7alZTyXrP8TbywaKiJw7d058fHykTZs2UqhQIXnssces+6KioqRdu3aSM2dOadKkidSoUcOaWC9evCgVKlSQwMBACQwMlOeee04OHTokIiKff/65NZFVrFhRgoKCZOzYsSIiUrx4cdMvwrBhw6zHjx07VqZOnSqArFixQkRE1qxZI4BMmDDBa8sh6sTqfhfatpXsNkm1YvHicvv27RR/LVcTq9u6AkTkjuVrpVRlIA9GorR87D+X8G94wr8POdnnCxRN1WBTQVooG7hs2TLi4+Np3bo1LVq04NChQxw5cgSAmTNn8uWXX1KnTh1q1KjBzp07rc/z8fGhRYsWTJo0icGDB3Pw4EH69OljOvf69evp2rUrBQoUYMCAARw8eJDQ0FAAypYtS1hYGK1atTI959VXX8XX19faXfL111+TKVMmXnvtNV0OUTN88gmTFy3iZqKmHH5+fLtjB1mzZvVYWG6vFaCUehT4BjiN8dH/PdtDEv51NNfW6T6lVBegC2D6JUmLPFU2cMmSJfj5+VGmTBkiIiKYOnUqS5cuZfjw4Wzbtg2AsWPHUqpUKVavXs1PP/0EQHR0NOvXr+fnn3+2fBqx/iGxePPNN+natSu+vr506tSJH374gddffx2AfPny0aZNGwAST0UODg6mQYMGrFmzhlu3brF69Wrq1q1LgQIFdDlEDSZOZNegQYyyaf5wxAgKFy7skZAs3JpYlVLlgC1ANFBHRP5RSp1K2F0k4V/LT+RUwmbZdz5hXyx3r2qtRGQ2MBuMWgGp8g08AG8vG3j27Fl27dqFiFC+fHlr+5IlSxxeNVsSKMDkyZPZuXMnvXr1okmTJrz11lvcvHnT7jm2z3NlVdO2bduydu1aBg0axNWrV2nXrt19fV8WuhxiOvPzz8T27UsLm+bgHDno3KOHR0JKzG2JVSlVFNiG0QXwAVBVKVUV+AoYAwxSSuUHmgM/icgJpdRC4G1glFJqI1AdCBORB168XYa5N/daygZu376dzp07U7FiRT7++GMg9csG7t+/nwULFlCvXj2OHDnCqVOnqF+/vum4pUuXIiK89957VKlSBYC5c+eyZs0afv/9d2rXrs2qVasYMGAA1apVM/2BsCTLW7du8eOPPxIeHk5Om3Js8+bNo2jRokyePBmlFLVq1SJHjhz4+Phw/PhxFi1axHM2NTMBmjZtSrZs2ZgxYwZZs2alefPm9/V9WVjKISZmKYc4YsQI6tWrR+/evYmOjmbFihX89NNP1nKIX375pfU57dq104nV0yIikJAQagL/2Ox6f8QIcuTI4YmozFzpiE2JDWNUgN2du4R9LYATGFey24FSiZ7XA6OPNQr4Fgi+12t566iAy5cvS4cOHSQoKEj8/PzkiSeesN5kErl782rv3r0iIhIYGCjFixcXEdduXln2lS9f3vKzlejoaBkwYIAUKVJEsmTJIuXKlZNJkybZxValShVRSpluAC1ZskQAef/99603r3LkyCGNGjWS5557znrzKjw8XJ5++mnJlSuXdOrUSR577DHJmTOniNy9edWpUye7m1ciIoMGDZLAwEAB5IsvvjDdvLJ44403BJDXX3/d2ubq92V788qRyZMnS5kyZSRz5sySL18+6dSpk/zxxx8CSOPGjU3HvvTSSwLIX3/95fBc3vA+S/e6dZMVDnIJ4PLaVcmFt40KcOfmrYk1I3KUKNMz/T5LZT16SDxIBQdJ9fjx46n+8q4mVj3zStM077dmDTz5JEyfznLggM3uMWPGUKpUKU9E5pBeQUBLVR06dKBDhw6eDkNLyxYuhJAQwOgPHGyzu1mzZl7X562vWDVN817nzsHbb1sf1sS4GWORKVMmRo8e7faw7kUnVk3TvFNcHLRoAdev8y/QDNhrc0iXLl0oU6aMB4JLmu4K0DTNO82fD3v2cBJ4DvuhVfnz57fO3vM2+opV0zTvc+0avP8+AnTCPqkCjB8/nly5crk5MNfoxOomicvX+fj4kDdvXl5//XVu3bqV7HOGhoYyceJEp/sPHz7M448/jq+vL0WKFGHbtm0opejVq1eyX9PWY4895nQG1fPPP2/9ngMCAnjqqafYsmVLsl7Hcq7Lly8n6/nfffcdw4cPN00RLlGiBNmyZUvW+bRUJAIdOsCFC2zEKIGXWEBAADNmzKBt27YeCM5FrozJSmubN45jtS1f17hxYwHkk08+SfY5E08gcKR///4CyNtvvy3ffPONXLhwQcLCwmT//v3Jfk1biScj2KpVq5YAMnv2bBk7dqz4+vpKzpw55caNG3bHxsTEJPk6mzdvlrCwMImOjk5WnD179hRAtm7dam1bvXq1fP3118k6nzOefp+lCzNnioDEgpSzGatavnx5uXnzpsdCQ08QMLO+4Y2/h6m3OWFJrPXr15eLFy9aZ0uFhoaKiEh4eLi0aNFCcuXKJQULFpRBgwZJXFycxMXFSdeuXSV37twSEBAgZcuWlc2bN1uTlmULCQkxvZ5lYL5lq1WrlrV0Yc+ePSU8PFxy5col1apVk/j4eOnevbsopWTbtm0iYpTLK1GihGTLlk3q168vJ06cEBGRa9euScOGDSVXrlwSEhIijzzyyD0T66VLl0TEmN1Fwswyy/ffsWNHKVmypPTp00euXr0qISEhEhwcLMHBwfLGG2/I1atXHZ7LWXzXr1+Xrl27SsGCBSVLlizStm1bu5+FJd7ixYtLYGCgiBhlEfv06SMFCxaUnDlzyssvvyxnzpwRkaTLOTp9n2nJc+CASECACMi3DiYBOPu5u4tOrDa8JbEm3goVKiQXLlwQEZHnn39esmfPLsOGDZOuXbsKIFOmTJFffvlFAGnRooXMmzdPBg8eLN9//71s3rxZ/P39JTg4WMLCwmTXrl2m1zt58qTUr19fABk6dKhs3rzZlFhFRBYsWGBNbkopefvtt0VEZP78+QLIq6++KiNHjpQCBQqI5Wfap08fAaRz584yZMgQU6KyZUmGf/75p/zyyy+SM2dO8fHxkX/++ceaWAsXLiwzZ86UTZs2Waeu9u/fXwYMGCCAvPHGG6ZzXbp0Kcn4OnbsKIC0a9dOPvvsMxk6dKjdzyIsLExEzIl16NCh1j9QoaGh4uPjIzVq1BCRu4m1WrVq8tFHHwkgderUcfg968T6ALZvF8mVy/q71MDm96VVq1aejlAnVlveklirVq0qGzdulPfff1/AKO588+ZNUUrZJd4mTZrIv//+K1mzZpXSpUtLz549ZdGiRXLnzh0RuXdXgO3HX9vEKiLy8ssvCyAPP/ywtTBwq1at7GIB5MqVK1KhQgXx8fGxfiS3rILgiO1Vta+vr/UK3ZJYE8/vDwoKksKFC1sfFy5cWIKDg03nunTpUpLxBQcHS758+UzFqh39LETMibVSpUri4+MjUVFRIiJSvXp1AeTmzZvWxPr999+LiIi/v7889NBDDr9nnViT6exZkcBA6+/R7w7+f/ft2+fpKF1OrHq4lZsFBwdTr149nn32WUaNGsV3331H/4Q10J988knTWk85c+Ykf/78HD58mFWrVrF7927atm3LkSNHGDlypEtl95ISFxfHhQsXALhx4wYRERFkzZrV+IsLLFq0yFrDND4+3lo4OPHrWo5NyvLly8mXLx+PPPKI9XwWhQoVuu+47xWfIw/6s4Lkl3PU7uH2bXjmGeNfjCza2+aQZ555hkqVKrk9tOTKeKMCUvua9R7Onz/PV199ZU2mljvTtWrV4vfff+fHH3/k5MmTLFy4kI0bN3Ls2DE++eQTsmfPTtWqVa3nAMidOzeXLl1iwYIF1kr/92PcuHHs3r2boUOHcu3aNXok1LFs0qQJAAsWLODs2bP88MMPfPTRRwQEBFC7dm3i4uLo1asX77//PufOnUvqJQCoWbMmzz33nF1StdWoUSPOnTvHoEGDGDRoEOfOnaNhw4Z2xyUVX5MmTbh48SIhISHMnTvXWpIxd+7cgLEKwdq1ax2+dnx8PN27d+fjjz9m165d1KxZU48aSG1RUVC+vDHDKsEajPqiiXnblNV7cuWyNq1t3jwqwLLlyJFDGjRoIKdPnxYR4+ZVq1atJE+ePJI9e3apXr26bNy4Uf7++2+pUqWKZMuWTbJkySLVq1e3LrQ3ZcoUyZWwiNpHH31k95pJdQUcOXJE/P39pUmTJiJy96P5kiVLRERkzJgxUrJkSfH395dSpUrJkCFDRETk6tWr8tJLL0mOHDkkJCRESpcufc+uAMsNp8RsSx1azt2+fXsJCgqSoKAg080ry0KMlsfO4rt+/bp06dJFChQoIFmyZLGWGjx69KiUKVNGlFJSqlQpEbG/efXOO+9IgQIFnN68clTO0Zan32dpTv/+pkuTGJDyAQGm35X69etLfHy8pyMVEde7ApS4cJWV1lSuXFkSL/EB8Mcff1C2bFkPRaQ9iK1bt9KmTRvu3LnDtWvXPB1OkvT77D7MmgU2S43PCgigW9TdOvZKKQ4ePGhdgcPTlFL7RaTyvY7LeF0BWprTokULlFJMmzbN06FoKeHECahVyy6p7gb62FzotW/f3muS6v3QN680r+ftV6nafYiMhBdegFOnTM03lKJp7txEXb1qbcuaNSujRtkuFZg26CtWTdPcZ8wYu6QK8EWrVlxIlFQBPvjgA4+vtppcOrFqmuYehw+Dg9qpV6ZModeyZaa2Ro0apb2RAInorgBN01JfbCy8+SbExNxty5qV2JMneTTRcuuAtT/dxyftXvel3cg1TUsbYmPh9ddhzx5z+4IFfLVxI1euXDE1d+zYkeLFi7sxwJSnE6ubJC4baNlcqSXZoUMHlFLs27fPeo7GjRuneHwREREMHz6c+fPnW9vmz5+PUso0Gyw5lFI89thjDxihY+PGjUMpZYpb8zL/+x/YfNSneXNiXn6ZESNGmJqzZcvGp59+6sbgUofuCnCzihUrMnDgQAD8/Pw8HM1dERERjBgxglq1amWIxf9iY2Px9dVv/1S3cyfY3tkvUgRmzGDkqFEcP37ctGvFihXkyJHDjQGmjgx3xWp71ZjS273kzZuXevXqUa9ePerWrQtAdHQ0ffv2pVChQuTKlYumTZty9uzZe57r7NmzNGvWjNy5c1OoUCH69OlDdHQ0GzZsQCnFrFmziIuLI1u2bLz00ksAtGzZEj8/PyIiIkznqlzZGPP8ww8/oJRi+PDh1n1Hjx6lcuXK5M6dm0mTJlnbR48eTcmSJcmePTsNGjTg5MmT94w5sXPnztGyZUtr/IMHDyY+Ph6A1q1bkzt3bgICAihXrhwrV660Pm/cuHEEBwdTqVIlfv/9d9M5f/75Z6pVq0a2bNl45JFHCAsLA+5+YqhevTr16tVLs3eb05T9+6F5c0j4PwWgQAHkxx/5YMoU/ve//5kO79KlCy+88IKbg0wlrkzPSmtbUlNacVA1JyU3ZxyVDaxVq5aIuFaybu/evdZzNGrUSEREatSoIT4+PhIaGiodOnQQMMri3bx5U3x9fSUkJMRadjBXrlwSHx8vBQsWlGeeecYuvkWLFgkgZcuWlbCwMPn999+tdUzz588vU6dOtdaEjY6OTrJ0ny0wChTbclYqUUTk448/ltmzZ8uECROkbNmyki1bNomMjJQDBw5Yzzdr1iwpVKiQAPL555/LlStXJE+ePFKmTBkJDQ2VF154QTJlyiS//vqr6ec/cOBAmTVrltP/qwehp7SKSHy8SN++pqmq1m39eunVq5fd70KWLFnk3Llzno78ntBlA828JbFaygZu3LjRWgbNlZJ1ton15s2bAkj16tVFxJjr7uPjI5UrVxYRkcqVK8sjjzwi06ZNk1q1aklAQICsXbtWABkwYIBdfJcuXTIle5G7xbLfe+89ERFp0KCBAHLmzJkkS/fZcpRYkyqVGBsbKyEhIeLn52fa98cff8jEiRMFkM8++0xERD744ANrYl2zZo3DmMaPH29awSE1ZfjEeueOSL16jpNq795y5MgRh/9H06dP93TkLnE1sepOJjezlA18UMb/sfNyeDVq1ODTTz/l22+/pWbNmsTFxTFhwgTAqDZlK6lujMTl8sAoN2h5/fsp3eeIo1KJGzduZMGCBdStW5c+ffowc+ZM1q5dS1SiOeQWljgSf92+fXveeOMNa3uJEiWsXyenTKHmorg4aNoUNm2y39e4MXz8MTMHD7bb1bNnT7rZTG9N6zJcH6srf20eZEuO5JSsy549OzVr1mTXrl2MGTOGHj16EB8fby2zZ0me33//PdWqVaNatWps3rwZpRTPPvus3fly5MiBj48Px48fZ9GiRfz9999JxpxU6T5HLly4wODBg63bzZs3nZZKtPwcIyIiOH36NDt27LCe5/nnnwfg008/Zfbs2Xz++efWfdWrVydPnjysX7+eo0ePcujQIcaMGeNSaUMtBUyYAOvWmduyZYMVK2D1am7Hx9uN3ujVqxdTp05NkXq5XiW1E40nNm8uG2jpH03M1ZJ1tuc4c+aMNG3aVHLmzCkFChSQt99+29qdcOnSJVFKiVJKrl69KitWrBBAHn/8cacxDho0SAIDAwWQL774wtoVMHbsWBERadSokQBy6tQpEXFeus8WDj76/frrr05LJcbGxkqbNm0kW7ZsUrt2bXn11VetzxERGTt2rAQFBcnjjz9u/fl8/vnnIiKyc+dOqV69ugQGBkpwcLA0adJETp8+neTPPyV5+n3mMYcOifj52X/837TJesj48eNN74F8+fJZ369pBbpsoC4bqLlfhnyfxcQYKwD88ou5/ddfoUIFACIjIylZsqR1xQqAIUOGpLkiK7psoKZpqe/UKXjoIfukGhZmTaoAoaGhpqQaGBhI37593RWl2+nEqmla8nz1lZFUw8PN7a+8Am3aWB9u376dkSNHmg7p0aMHwcHB7ojSIzJUYk2P3R6a98hQ76+VK+G11+zbCxWC6dOtD0WEDz74wHRIcHAwAwYMSO0IPSrDJNaAgACuXLmSsd78mtuICFeuXHE6KiJdGT8eWrRwvC8sDIKCrA+3bt3Kjz/+aDpkyZIl5M2bNzUj9LgMM461SJEihIeHc+nSJU+HoqVTAQEBFClSxNNhpK4VKyBhhWETf39jyZVEU4VFxK7IyosvvkidOnVSO0qPyzCJNXPmzJQsWdLTYWha2vX999CypX17x44wb56pSUSYMGEC27dvN7UPGzYsNSP0GhkmsWqalkyxscaKqr162e+rUAHmzDE1xcfH0759exYtWmRqb9CgAc8880xqRuo1Mkwfq6ZpySAC7ds7TqodOxrDrDJlMjWHhobaJVVfX18++uij1IzUq+jEqmmac7NnGzekbLVrB3Pngs1U1N27dzN06FBTW44cOVi2bBlPP/10akbqVXRXgKZpjh0/Do6Ko/TqBRMn2iXV6OhoWtiMFvD19WXLli1UqlQpNSP1Ojqxappmb9s2qF3bvn3FCqN4tQNTpkzh/PnzprZOnTpluKQKZJxaAZqmuejYMXj0Ufv22bOhc2eHT4mKiqJkyZL8+++/1rY8efJw6dKlNL3aqi1dK0DTtPsXF2csU22rShV46y2nT5s/f74pqQL88ssv6Sqp3o+M+V1rmubYiBGQqP6t1TffgJMkGRsby8cff2xqe/fdd9P8EtYPQidWTdMMO3aATbEUgoIgMhIKFHD6tMWLF3P69GnrYz8/P/r165dKQaYNOrFqmmYkzzp1jHGrFvnzw6FDkET9g6+++oqQkBBTW8eOHTP8Ejg6sWpaRnfnjjGD6s4dc/uXXyZ5pTp+/Hhes6lw5evry6BBg1IjyjRFJ1ZNy8hu34ZKlYyRAIm1bg1JLHp5+vRp3nvvPbv2jh076poc6MSqaRlXdDQ0aGB83E8sRw5jAoATIsJTTz1FTEyMqT0oKIhPPvkkNSJNc9yWWJVSk5VSFxIWuFuTqH2bZdG7hO2/RPvKKqV2KqWilVJ/KqXquyteTUv3Bg1yPAJg/XqjYLUTvXv35tq1a6a2rFmzcunSJXLlypXSUaZJ7r5i/cpJ+x/Aawlb4kF0YUAZoB8QAyxTSuVM1Qg1LSP47DOYNMmSTEpBAAAgAElEQVS+fd8+qFbN6dN++uknpk2bZmorVKgQ165dS39LWD8AtyVWEXkb+NTJ7ovAWhH5SkRWACilKgJPAmEiMg2YAOQAWrkjXk1Ltz7+2H4GVcGCcOaM0d/qhIjQp08fu/Y5c+bg5+eX0lGmad7Sx1oTuKGUuqGUej+hzdIDfi7hX8uKZQ85OoFSqotSap9Sap9eJUDTnDhwAIYMMbdlymTUAChaNMmnTpo0if3795vaWrZsScOGDVM6yjTPGxLrcqAd0Bo4C4xUStVwcJzlc4bD4gYiMltEKotI5fS+no6mJYsI9OkD8fF323x8YOZMuEcB6iNHjtgtV/3QQw+xbNmy1Ig0zfN4dSsRmWL5WilVEJgMlAP2JDRbFhGyLKZzyn3RaVo6smgR/PCDuW3JEmiVdO/azp07qedg6NXq1at1v6oTbkusSqlGwGMJD4sqpToBezES6UogEugDxAN7ReRXpdRvQBul1GGgO3AT4wpX07T78fPP0KmTue3FF++ZVCMjI2nVqhWRkZGm9sGDB1O+fPmUjjLdcGdXwABgTMLXTwBzgGrAJWAwMBGIAtqLyC8Jx70O/Ilx48oPeEVE/kPTNNft2wf16xvjVi0yZTKWsb6HuXPn8s8//5jaqlSpwvDhw1M4yPRF12PVtPTs+nXjjr/NFSdjxzpexjqRO3fuUKpUKcLDw61tuXPnJjw8nKxZs6ZGtF5P12PVtIzu+nUoWdI+qQ4ceM+kCrBixQpTUg0ICODIkSMZNqneD51YNS09EoGXXgKbGVK8/jqMGeP4OYlcuXLFrsDKm2++SYEkirJod+nEqmnpzbVrUKuWccMqseLFYdYsu0UAbcXExBAcHGzX3qNHj5SMMl3TiVXT0hMReO01+PFH+30bNkC2bPc8xZw5c+zaqlatqkcB3AedWDUtPVm8GL7/3tzm5wd79sAjj9zz6XFxcUyYMMGufenSpSkVYYagE6umpRcXLxozqxIrWBC2boWnn3bpFCtXruTEiROmtg0bNlCsWLGUijJDcGmCgFIqK5ANuC4i0QmD/esAB0VkYWoGqGmaC0SgSxe4fPlum7+/MdOqdGkXTyGMHTvW1NaxY0deeOGFlIw0Q3B15tVnwKtAtYSyfd+SMGdfKRUkIs6qVmma5g59+hgrqSb24YcuJ1WAn3/+mT179pja+rswLEuz52pXQCXghojsAV5JaDuDURilY2oEpmmaixYuhMmTzW1VqxrjVe/DJJv6rA0bNqRcuXIPGl2G5GpiLQz8nfD1E8AfIlISOA2USPmwNE1zyeHD0LOnuc3fHxYsAF/XS4GcOHHC7gaVo9qrmmtcTaxxgGUN3EeA3xO+vs7dcn6aprnTlSvQqBHcumVuHzUKHn3U5dPExsbaDaUqW7asw4pWmmtcTawngNJKqUMYVfwt1W4LAedTIzBN05IQHw9vvAF//21ub9QI3n33vk41ffp0ohMXaAF69uypSwI+AFcT60SMK9NywH/AF0qp8kBejNJ/mqa5U2gorFtnbuvYEVavvq/TnD59mnfeeceuvX379g8SXYbnUieMiCxUSv0OlAZ+FJF/lVI+wEvA8dQMUNM0G5s2GXf8E6tWzVgJwMf1oenx8fG88cYbdu1btmwhe/bsDxplhuZy77aI/Ar8CqCU8geCgJ9E5HYqxaZpmq2LF41CKonLfQYHw9Klxgyr+7Bs2TJ++uknU1vTpk2pXbt2SkSaobn0500p9bZS6jul1FNKqaLAMeAAcE4plfRiOZqmpYzff4cyZSDxYplKGdNYixRx/jwHRITRo0eb2nLlyqWnrqYQVz83vAbUxUionYGiGH2uOYBhqROapmkAxMbCK6/AE0/YlwEcNgySMTNq/fr1HDx40PpYKcWmTZv0MtYpxNXE+jBwRkRuAc8AVzCmtN4Ankql2DRNA5g0CRythlq+vP1S1i64evWq3RjVli1bUqlSpeRGqNlwNbFmxxizCsY41n0isg3jxlXOVIhL0zSAzz5zPIMqONjoAsic+b5P2bt3b44dO2ZqGzx4cHIj1BxwNbFeAMoppQZhdANYJgjkwbh61TQtJcXHw4AB0Lmz8XVio0bB2bNG18B9On/+PEuWLDG16avVlOfqqIC1QDcgNOHxGqVULqAI8IPTZ2madv+io6FDB/jqK/t9X3wB7dol+9Rz584lLi7O+rhYsWIsWLAg2efTHHM1sQ4AYjDGsX4rItuVUk8Dy4H7G5GsaZpjIjB9OgwaBLdtRjFmyWIk1ZYtk336uLg4u9UBBgwYQGBgYLLPqTnm6gSB28A7Nm17MUYLaJr2oKKi4M03ISzMfl+BArBmDTzgx/X169dz9uxZ6+MsWbLQ7gGufjXnXJ6moZQqr5QKU0r9lbAtTpjWqmnag7hxAxo2dJxUS5eGnTsfOKnGx8czxmZ11ldffZVcuXI90Hk1x1xdQeBJ4CcgK3erWT0ENFFKPSsiv6VSfJqWvm3ebHy8v37dft9rrxlDrfLmfeCXWbx4sd0sq65duz7weTXHXL1iHQEEAjcxVg/4NuHrwIR9mqbdj7g4Y77/Cy/YJ9UyZeDPP43hVCmQVK9du8ZAmyFbzz77LFWrVn3gc2uOuXrzqgbGZICyIvIPgFKqEPBHwj5N01x1+zY0a2YUU7FVtSqsXQtBQSn2ckOGDOGff/6xPvb19WX+/Pm6LGAqcvWKNRtwzpJUAUTkPBCOMXlA0zRX3L4NL7/sOKk2aGC0p2BSvXLlit1wqm7duvHwww+n2Gto9lxNrGeBR5VSvZVSwUqpIKXU20CZhH2apiUlPh6mTIHixWHLFvO+gAAYOtS4Us2WLUVf9oUXXiAyMtL6OH/+/HY3sbSU52pXwDJgEEbB64k2+5bYH65pmtWNG/Dqq7B+vf2+8uXh22+hZMkUf9nNmzfz66+/mtr69u2rx626gatXrP8DtmGMCEi8bQM+So3ANC1dOHwYqld3nFSfeAK2bUuVpAowYcIE0+MsWbLQqVOnVHktzczVCQKRQB2lVD3g6YTmvSLioKNI0zQAvvzSWEH1xg1ze+bM8M47xqiAVKrU/9lnn/Hdd9+Z2gYNGkRQCvbfas4pSVyJ/H6eaCzNMhhARELvcbhbVa5cWfbt2+fpMLSM6vBhYxbVnj32+ypVMoZRPfJIqr38b7/9RqVKlYiNjbW2Va5cmT179uiRAA9IKbVfRCrf6zjXF8ixlxkYidFNoGnalSswYwY884zjpPraa7B9e6omVTAG/idOqgB9+vTRSdWNXF7zKgn6f0vL2Hbtgk6djCtVRzJlMkr9DRxoLKWSivbu3cuuXbtMbS1atOC113RZD3dKicSqaRnXRx/Zr5iaWKVKRrHqChXcEs6MGTNMj59//nmWL1/ultfW7nqQrgBNy5hEYOFCePRR50nV1xemTYO9e92WVC9fvszixYtNbX379nXLa2tmSV6xKqWOJbU7hWPRNO93+zb07w8zZzo/ZuBAo/K/m2c3ffbZZ0RHR1sflyhRgkaNGrk1Bs1wr66AhwFBJ1EtoxOBlSvh3Xfh9Gn7/Zkzw+jR0Lcv+Lj/g2BsbCzTp083tfXo0YNMmTK5PRbt3ol1J0Zi1bSMSwT69IHJkx3vr1ULJkyApzy3YPHq1avtili/9dZbHosno0sysYrIc+4KRNO8UlQU9OgBn39uv69IEWNZ6meecX9cicTGxtLSZsmWtm3bkidPHg9FpOmbV5rmTHQ0NG/uOKn26gVHjng8qQIMGzbMrq1Hjx4eiESz0MOtNM2RDRuM2VPnzpnbg4ON1VPr1vVMXDYiIyPthlhVrVqVihUreigiDXRi1TR706YZV6S2SpY0FvUrV879MTkxe/Zsrl27ZmobP368h6LRLHRi1bTEZs92nFQffdSYjpovn/tjcuLUqVMMHjzY1NatWzeeffZZD0WkWeg+Vk0DoxD1Bx+AowX2KlQwFv3zoqQKxpVpVFSU9XH27Nl5//33PRiRZqGvWDXtn3+gSRPYv99+X8+exkqpXjYe9Nq1a3xuc1Nt/PjxFClSxEMRaYk5TaxKqTsunkNExD+F4tE09zp/HmrXhmM2kwz9/WH5cvDSmUuzZ88mIiLC+rhw4cKEhIR4MCItsaS6AnzvY9O0tGfVKihWzD6pBgfD9997bVKNiYlhypQpprZevXrh5+fnoYg0W0klxVEp/WJKqcnAq0A+YK2INE5oLwvMBSoBp4HeIrIhYd+zwAzgUeAw0ElEfknp2LQMxtmd/zJlYPVqKF3a/TG5aNmyZZxLNAwsa9asdOnSxYMRabacJlYRGZpKr/kV8LZNWxhQDOgHdAeWKaWKAdHAciAS6Au8D3ytlCotInGpFJ+W3k2aZExRtfXyy7B0qdEN4MU+/fRT0+MOHTroWVZe5r4+xiuligKFAFNPvojsdOX5IvK2UqoEiRKrUqoi8CQwXUSmKaUiMa5eWwFXgfzAQBGZrpQqAAwFngc230/smsZ//xmVqebOtd/XurWxRpWXf5zet28fiZcdUkrxzjvveDAizRGXEqtSKh+wEnA0f09cPY8TliUqLZ9twhP+fQjImcQ+U2JVSnUBugAUK1bsAcLR0qUjR4zZUv/+a79v2DBjSwNLl3xoU/+1QYMGPJLKS71o98/VhDgGqJaagSRieXc7qqrldJ+IzAZmg7GYYOqEpqVJe/YYN6IuXza3BwTAt99CvXqeies+xMXFMWbMGNatW2dq79atm4ci0pLi6gSBBhjJrHvC4yMYH8mvAq8/YAynEv61DMArnKg9qX2adm9r1xrDqWyTao4cxr40kFRjYmKoW7cuH3zwgam9cOHCupC1l3I1seYF/hSRWQmPb4nIKOASRl+oS5RSjTBGBQAUVUp1Am4BvwFtlFI9MW5g3cS4abUOuAh0V0p1B97CGDWwzdXX1DKwuXOhaVNINN4TgMaN4a+/oE4dz8R1n5YsWcIPP/xg1z58+HB8ffVoR2/kamKNAO4k+rqkUioYCAJevI/XG4DRrQDwBDAHeBbjqvdPYALgB7wiIv+JSBTQGiP5TsJIsq31iAAtSSLwv/8ZK6fG2bxVhgwxhlN52fRUZ+Li4uxWBgAYM2YMnTp18kBEmitc/XN3DmM4FBgJsAJwIeHxP66+mIg8n8Ruh324IrIdeNzV19AyuOvXjYT69dfmdqVg6lSjaHUaERERQdu2bfn5559N7StWrKB58+YeikpzhatXrN8DF5RSjwNTMG4iWbYpST1R09xm714oW9Y+qQYEGNNT01BSFRFCQkJYtWqVqb1OnTo6qaYBLl2xikg/jL5PgN+VUicxhl79JiLrUys4TXNJZCSMGmVstnLnNu78p6FSepGRkfTs2ZOvbf5AFC5cmFmzZjl5luZNXB3HOgQIF5GFYP14vl0pVUUpVd8y/VTT3O7mTWP5lM0O5osULw7r1hlXsWnE9u3bqVu3LrGxsab2cuXK8f333+vqVWmEq10BIwFHA+YmAd+lXDiadh8uXTLu7DtKqs2awaFDaSapigijR4+mdu3adkk1e/bsrFq1SifVNCTJK1alVKFED/2UUgW5O0g/EOOGlvdPV9HSn23boH17SLTks9XMmY4LVnux2bNnM2TIELv23Llzs27dOkp7cVEYzd69ugIs71oBKnJ3SmliFxy0aVrqEIERI4zhVGIzwa5iRWNNqkKFHD/XSy1fvpy337atSwRPPPEEq1evpnjx4h6ISnsQ9+oKsNz5T/y17eagooWmpYJjx6B+fSOx2ibVevXghx/SVFKNj4+nc+fOtGrVijt37taVz5IlC6tWreLAgQM6qaZR97pi7Zzw7xzgBHcH94MxUeCoiPyaGoFpmskPPxj9pv/9Z25XCvr1g9BQr69MlVh0dDSvvfYaK1eutNs3ffp0mjZt6oGotJSSZGIVkbkASqnawF+Wx5rmVl98AW+9BTEx5vZcuYzxqWlkaqpFTEwMbdu2dZhUZ82aRYcOHdwflJaiXB3H2k4plUUp1RWonNC8F/hCRCJTLTotY9u3z6if6mCePM8/D9Onp5m7/okNHDiQ5cuX27WPGzdOrwSQTrg6jjUI+BFjeRSLN4F3lFI1ReRKagSnZWCLFkHHjvZXqWB87B88OE3UT03s0KFDjBgxwm7gP8Dq1atp0qSJB6LSUoOrtQJGAWUSvr6Z8G/2hLaPgLQzV1Dzbr/+alTzP3HCfp+/P8yfD23auD2sByEiDB06lNDQUMTmpluhQoXYvHkzZcqUcfJsLS1ydYJAIyAGaCIiOUUkJ/AyEAfoP7Naypg7F6pWdZxUixWDTZvSXFIFGDt2LKNGjbJLqn5+fixZskQn1XTI1cSaHzgmImstDSKyBqPSVdqov6Z5r2PHjGVTOnWy/+jv4wMffggnT8Jzz3kmvgdw+PBhhg0bZtdeqlQptm7dynNp8HvS7s3VxHoFKKWUetLSoJSqADyMsYqApt2/mBh4911jyektW+z3FytmVPkfMQIyZbLf7+Vu3LhB06ZNiYqKsrblzp2b+fPnc+zYMapXr+7B6LTU5Gpi3QQEAPuUUoeUUr9jjArwAzamVnBaOnbihDGof8IE+8H+AQFG7dTTp+HF+6mj7j327NlD0aJFOWHTrTF+/HhCQkLw8XH1V09Li1y9efU+UA+jS6BcovZ/E/ZpmmtEYMwY4yo0Otp+f968xrTUKlXcH1sK2bFjh8OP+I0bNyYkJMQDEWnu5uo41jMJH/17A08nNO8FpoiIrhWguUbEmCU1caL9vsKF4e23ISQE8ud3f2wpZP/+/bRs2dKuvUyZMixYsEBfqWYQ96pu9SFGHdZ5CQn0g6SO1zSn7tyBXr1gzhz7fXXrQliYcbWaRp06dYrQ0FDmz59vV/avYcOGLF68mJw5c3ooOs3d7nXFOhzYBcxL/VC0dCsy0pjnv8GmHnqePDBpErRtm+YG+1vExMTw7rvvMmWK4xWKcuTIwcqVK/FLQ3UMtAenP5doqefff6FzZwgOtk+qQUHGSIB27dJsUgUYNmyY06Tap08frl69qpNqBuRKH6u/UqooSRS0FpEzKReSluZdv24sM+1g2WbAWHp6wwZ48knH+9OILVu2MG7cOLv2/PnzM2PGDL3oXwbmSmKtAJxOYr+4eB4tI/jxR3j5ZfvyfhaPPQbffQdFi7o3rhQUERFB9+7dWbhwod2+8uXLs3PnTnLkyOGByDRv4WpXgLMi14kLYWsZ3fbtxsJ+zpJqo0ZGpao0nFRv3LhBy5YtHSbVCRMmcOjQIZ1UNZeuNM+hVwnQ7mXmTOje3fG+1q1h2DAoX969MaWwW7duUbNmTQ4ePGi3r3PnzvTt29cDUWneyJXEGi4iI1I9Ei1tskxLdXQDp0YNGDvWKKySxkVGRlK7dm27pJo7d24+/PBDevfu7aHING+k+0a15IuJMaacOprn36cPjB9vFFFJB7p27cq+fftMbU888QRr167Vy1Jrdu71rj8D/OOOQLQ0RMSYPZUzp+Ok2qaNUQMgnSTVsLAwvvjiC1Nb4cKF2bRpk06qmkP3WvOqhJvi0NKSKVPAWX9i584wa1aaHpua2N9//013m77jokWLsnv3bvKm4ZliWurSXQHa/dmyxZjv78hXX8Grr7o3nlQUFxdHu3btuH79urXN39+ftWvXUrBgQQ9Gpnm79PFZTXOPXbuMqalxcXfbfHygfXvYsSNdJVWAoUOH8tNPP5naPvnkEx5//HEPRaSlFfqKVbu3CxfglVeMcaq2vvgCXn/d/TGloqioKBo0aMB2m+/3xRdf1Hf/NZfoxKolLT7eGPT/88/2+7p2TXdJFeDdd9+1S6p58+bl888/R6WTvmMtdemuAC1pCxc6Tqr16zuuq5rGnT17ljk2pQ19fX2ZP38+BQoU8FBUWlqjE6vm3PXrMGiQffu8ebBunbGESjozbtw4YmwWNPz6669p2LChhyLS0iLdFaA5FhcH77wDFy+a20+ehJIlPRNTKrtw4YLd1erEiRNp2rSphyLS0iqdWDV7V64Yfae2NVQ//DDdJlURoXfv3kRGRlrbihQpYjeGVdNcoROrZvbrr0YVqn9sJtyVKOG4WyAdEBHat2/PsmXLTO0DBw7URaq1ZNF9rNpdCxfCU0/ZJ9VChWD1asia1TNxpbIJEybw5ZdfmtqeeuopunXr5qGItLROJ1bN8NFHxgqptp58EnbvhnQ6KP7o0aO89957prYsWbKwaNEiMmfO7KGotLROdwVoxsqpH35o316lijGFNTDQ/TG5yccff2w3CuCzzz6jTJkyHopISw/0FWtGt2aNMdDfVrNmsGlTuk6qe/bsYf78+aa2KVOm8Ho6nPSguZdOrBlVXBy8/z40aWKUAUzsu+9g5UrInt0zsbnByZMnqWpTgLt06dJ6FICWInRXQEYUH2+U/XNU9T80FF56yf0xuVFMTAxt2rSxax81ahSZMmXyQERaeqMTa0YTH28UVFm+3Nzu4wNDh8LgwZ6Jy42++OIL9u7da2qrW7curVu39lBEWnqjE2tGEhdn3Pm3TaoAa9cay6ykc3FxcXzyySd27V999ZUHotHSK51YM5KBA2HRInNbpkzGGNUMkFQBZsyYwZ9//ml97Ovry/HjxwkODvZgVGnH2bPGLOeYGIiNhYceMoY52xIxKkrGxt49tnt3x6v1HDgA335r/N0XgQoVjIJqjsyYAX//ffe2QI8eULy4/XH//GOsYylibAULunl+i4iku61SpUqi2Rg50vIeM28bNng6Mrc5duyYAKYtJCTE02F53J9/irz3nkhIiEi9eiKdOzs/9q23zG+f6dOdH6uU+dg7dxwfN2+e+bj27Z2f85lnzMfu2OH4uN9+Mx9Xvrzzc94PYJ+4kIO85opVKXUaSPy356CIVFBKPQvMAB4FDgOdROQXD4SYNh06ZNz9X73a3B4QADt3QsWKnonLza5cuULdunVNbdmzZyc0NNRDEXmPc+dg9Oi7jytVcn6sbUGz2Fjnx/r6GlerFvHxjo+7nzUnXS2H6+myuV6TWBNsx0iiANeUUgHAciAS6Au8D3ytlCotInFOzqFZhIdDrVpw9aq53dfXGPifQZIqwOTJkzl79qyprWfPnhRy9Dk2HbJ8dHaUcAoXNj++ds35ee4nsWbOnLzEajv6LymuHns/50wJ3pZYTwFrReQmgFKqOZAfGCgi05VSBYChwPPAZo9FmRacPWv0m9omVYBRo6BaNffH5CGRkZFMnz7d1FagQAEGDhzooYjc69w5o1jZgAHQuLH9ftvE+t9/zs9VvLjx9zhzZuPvc1K1v9u1MxKrr6+xObsyffxx+OADY79SxixqZ7p1M4ZeK2VsjvpXwYhr/Hjja6UgKMj5OVOFK/0F7tiA00A8Rt/XReAtoF/C49cTjumS8LhzUufK8H2sFy6IFCtm35+aKVOG6lO1mDBhgqlf1dfXV/777z9Ph+UWBw+KZM9u/Pe3bev8uI8+Epk5U+Tbb0V++cV98aU1pLU+VmAO8CcQAIwBZgG29/EsH2TsLuyVUl0wEi/FihVLvSi93Z070LIlnDljbi9fHvbuhSxZPBOXh5w6dYr333/f1Na3b19y5szpoYjc59AhePZZuHXLeLxqFUREOC5S9sEH7o0tvfOaKa0iMkpEvhaRL4ElQCYgPGF3kYR/LR9aTjl4/mwRqSwilfPmzZv6AXsjEXjtNbBZspnmzTNkUgUYMWKEqXh1jhw56NOnjwcjcp8jR8x9nLdvgx6u6x5eccWqlHocCAXWYcTUHuOG1Y8Y3QLdlVI3MboHTgPbPBKot5s1C1asMLfVqAFLlxqdXBnMjh07WLBggalt/PjxGeaG1SuvQJkyxr9//glVq4KDmbxaKvCWK9bLGFeo/8PoBvgbaC4i54HWwC1gEkaSbS16RIC93buhVy9zW758sHhxhkyq169fp0OHDqa2Rx991K4tvXviCdi3zxicv3Vruq1V7nW84jdORP4BHC6DKSLbgfRZZTml3Lpl3IKNS/T3xs8P1q+HIkWcPy8dGzx4MMePHze1zZ49G98M+EcmWzawGRShpTJvuWLVkkvEmNdnk0SYNy9DjVNN7LfffmP27NmmtjZt2lCzZk0PRaRlNDqxpnVz5hiTshPr3h3atvVMPB4mIvTr14/4RKPRS5YsaZdo06Njx9w/EF5zTCfWtGz3bnj7bXPbI48Y1ScyqEmTJrF5s3nuyMSJE8mejot2gzEf5IknoE4d422heZZOrGnVwoXGb1F09N22rFnhm2/S9XIqSTl+/Dj9+vUztdWrV48mTZp4KCL3eeEF462wbRs88wz07u3piDK2jNeTnx58843jFVXnzDHG12RAt2/fpk6dOpZZfIBREnDChAkoT1fkSGXz5hnDqRJLalqolvr0FWtas3GjMQLA1pAhxoTwDCg+Pp727dvbFVkZOnQoj6fTZbsTe+kleO89yJXLeFy4MLRv79mYMjqdWNOSRYugYcO7cxTBKFQ9Zw6MHOm5uDxs9OjRrLCZGFGmTBmGDBnioYjcq2BBY6myM2eMLoHQUGO0neY5StLhbcTKlSvLvn37PB1Gyho/Hvr3t2+fOBHeecf98XiJy5cvU6xYMdO0VTBWYS1ZsqSHovKMEyfgwgWoXt3TkaRfSqn9IlL5XsfpPta0YPFix0n1008zdFIFmD59ul1SXblyZYZLqgClShmb5nk6sXq7W7fA5k43mTMbowIy+MTvM2fOMNZmaNmwYcNo1qyZhyLSNINOrN4sKgqaNjU+31lkymSsqPrCC56Ly0t0796dW4n6m3PmzEnfvn09GJH73Lmj+1G9mb555a1EoFUrYwmVxN55RydVYO/evXz33XemttGjR2eIOqs3b0KxYvDWW/Drr56ORnNEJ1Zv9O+/UK+ecWWaWLFi8OGHnonJi0RGRtKjRw9T27PPPkvXrl09FJF7ffml8SFm3jx46ilo0cLTEWm2dH2P0x4AABYjSURBVFeAtzl6FGrXNpKrrc2bIQNckd3LhAkTsB31MWTIEHzuZ7nPNEoEpk41tz32mGdi0ZzTidWbnDgBlSsbpd5tzZkDDz/s/pi8zI0bNxhvWSUuQaNGjXjppZc8FJF7nTtnP4w5g1yopyk6sXoDEaMKcYsW9kk1f35YsAAaNPBMbF5m+PDhXEu0PnOuXLlYsGBBup+2alGkCJw8CWvWGFeuefLYr7KqeZ5OrJ62bZsx/9BmOiYAFSoY+/XHfwA2bNjAp59+amrr27cvQW5f29izMmUyBos0bWpe00rzHjqxetLRo1C/vuPfjurVYd06yJHD/XF5oRMnTvDKK6+Y2goWLMjbtmUTM5jMmT0dgeaITqye8tVX9supWHTqBNOm6YGKCWJjY2ndujXXr183tS9atIhclsojqfK6Rim+6Ghj3GjOnI4Xuo2PN9aVEjE2pYyF+xwJD4e//rp7bNGiRgldR7Zvhxs3jPMXL64rVqUpIpLutkqVKolX27RJRCnL75Z5+/JLkfh4T0foVZYuXSqAaRszZozDY//7T2TfPpHYWMfn2r1bpHZtkerVRSpVEunUyfnrtm5t/q9ZssTxcRER5uP8/Z2fc+JE87G9ezs/tlw545gqVURu3nR+nOY+wD5xIQfpK1Z3u34dOna0X0PjxRdh1Srw9/dMXF5sxowZpsfNmzdn4MCB1sd37kCTJnD4sHHXHOD8eaPqk61bt4z7hBaOrkAtbD8wJK4pnpjtfbOk6hol59g9e4xVzLt1MypD3msxhPM3z3Po4iEu3LrAxdsXuXnnJrfv3CYiJoLbMbeJio0iJj6GmDijC2r1a6sdnmfsjrF88+c3CJIQqyCI9d/EbbbmN51P+Xzl7dpX/rGSUT+OSvobcOK9596jZbmWdu1HLx+l7Yr7X4rozQpv0rNKz2TFci86sbrTqVPQsqX9jarhw416qrrDzM7GjRvZmjgTAiNGjDCNAvDzg99+Mw/9vXTJcWK1XVzBWbIE+79xzo61HT6b0okV4MABI7H+/rv9OFaA2PhYRv84muV/LOfghYPOT2r7GjgfTXH86nF2nN3h8rkSi4iJcNh+OeIy+//Zn6xzXoq45LA9MiaSX/755b7P92KpF5MVhyt0YnWXAwegUiWjwyyx/v1h2DDPxOTl/v33X+rXr29qe+655xwWry5f3pxYL192fE7bxHrnjvPXz5IFAgKMBOvv7/zvno+PMfxYKWNLaoXtIkWgVq27xzrrXwWoWdOYbOfjY2yVKsGAAfbHnbp2iocmP+T8REkQhLj4ODL5ZLLbl1GGsKUGnVjd4eZNaN3aPqlWqgQffeSZmLzctWvXqOhg+e5x48Y5PL58eWNimq+vMY/C2SSskiWNRRgsCTOpQRdTpzq+OrTl6wt79977OIDmzY3NFdOmuXbcjH0z7n1QEmLiYxwn1iSuZrWk6ULXqe2PP4yOsQMHzO0NG0JYmB5O5cQ777zD5MmTTW1K1eXy5U3kyWN//MmTEBkJpUtnvMEUIsLav9bS7/t+/HX1LwCqFK7CQ7kfIl/WfOQKyEWgXyBZM2clMHMgAb4BZM6Umcw+mfH18aXBww3w9bG/xvrryl/8e+tflFLWJGv52lFbYuXyliPQz35Ry0u3L/H39b+T9X0Wy1mMfIH57NojYiL449If932+AtkKUDjH/c2ucLXQtU6sqWnePKMEka1HH4WDB/WNKif++usvHnH4GfkXxoypSP/+xiB5zezPy3/SZnkbFrVYRLm85TwdTrqkVxDwpPh4mDDBcYdY5cqwaZNOqk78/PPPVHe4tsgZAgKK8thjOqk682jwo/zaVdcR9AY6saa08+eNTrQ9e+z31agB336rp6gmoXfv3nZtPj7TyZy5KN9+a1RT1DRvl/7rrLnTnTtGyT9nSXXLFp1Uk7Bhwwb27zcPxSlZsiTz5r3JH3/opGqx99xe0mMXXnqiE2tKOXjQmHd47Jj9vpAQo5hKUuNwMriYmDjeeaePXfumTZsICfEnA64N6NDec3up8lkVKs6qyPwD84mOTWIgruYxOrGmhB9/NCpR2RanLl4cdu6E+fOdj//RiIqCEiW+4OhR853db775hoceSt74zPSqx3fGygkHLxyk4zcdeWPlGx6OSHNE/7Y/qLVrjZHctgoWNBJutWrujykNuXgRypW7wPnzHU3tTz75HE2aNPFQVN7nWuQ1np7zNPvOm0e7dKjQwTMBaUnSifVBHDtmjFG1VakS7NpllC7SkhQRcYZTpwrZtPozbtxMPfMnwc3om1SeU9kuqT4a9CgvPpx60zK15NOdfskVFWXc/b9xw9xeo4YxnCqjjVJPhmvXrlGnzvOAeUbam2/2oV49o4BHvMQTHRtNVGyUaYuNjyVe4omTOOIlnniJp1zecmTzy2b3OuE3wjl88XCyYiybtyzFchaza78RfYMdZ5I3j75wjsI8kf8Jh/vWH19vd2Nq0u5JnLx20u7YyS9NxkfpayNvpBNrcg0cCEeOmNveegtmzNDFVFw0duxYTp06ZWoLDMrBoSe2UmpyKS5HXOZG9A0nz7a3u9NuqhSuYte+/vh6On/bOVkxTnlpCr2q9LJrP/3faRoubpisc3ao0IHPm37ucN/LYS8TE3/vZQGWv7Kc+qXq3/M4zTN0Yk2Ob7+1n0TeuLGx4J/++OqS//77f3vnHmVVdd/xz3d4zMgwPGcsLxUFqpLoSgRaklilJErSRiMGtGlaVNSkCCiyFMiyGhddSqSlrphgCFqNQ5L6AFeiGNOo4RUSRYowkIgoAQQEHBFBhpfD7P6x9+WemXvnzqP3OfP7rHXWuWfve87+3nNmfnff39779/uI+Q0Www8aNIiBUwfyykevtOqada4uaXlb69VtmbKFIb2H5FqGkYK29ReXDebOhSuvrB/vrU8fePxxM6rNxDnHww8/zKGIG6W8vJz169cz8oKRrb5uezCsEz8z0YxqAWA91ubiHNx1F8yZU7+8qAgWLYLy8tzoKjCWrF3FL363kCe/+2S98mnTptG1a9ekvsfiDsWUdCw5tRV3LKZzh84UqejU1kEd6NKpS9I2+5f157JzLmvVYNgZ3ZIPQJZ1Lmv1wNEFpyeGPYwxZvAYautqE8qFGN5vOLMuntWqNo3sYkFYmsuCBTBpUv0yybsEbrklvW21UZ5ZWcU1L3wOFhyBiOu0tLSUnTt30rNnT3Yf2s3y7csZ3m84FaUVdC/unjSknWHkAgvCkk5+/WuYkjiAwZw5ZlSbSXXNB0x48Uqoqm9UAe655x569uwJ+BHzb17Y8jQbhpFPmGFtio0bfQL3htlUH3wQpiUuwTQScc7x5R9N5FjRDlhVv+66667jzmRRwAyjgGk7Xv1M8PLLMHp0Yv6OmTPNqLaApVuWsq7meXgdqImXl5aWMnfuXFsIYLQ5zLA2xoIFMGZMYvKkcePgvtZlmWyPOOeYvXI2HAcazKefOnUqp5+eGBHeMAodM6zJ+NnP/EBVwxxVV10FTz9tkZabyeYPNjPooUF+KeYrQCRxZ1lZGXfccUfOtBlGJjEfa0M2bYLJDXKNFxX5n/73329zVZvJgaMHOH/++f7gENBgksa0adPo3bt31nUZRjYwwxpl3z6/gurgwXhZcTEsXuzLjWZz87O3xg/WUC8cwMCBA5kxY0bWNRlGtjDDGsHdeCPaUT+D5NszbmLXp0ph27KE959Xfh59y/omlB88dpB1e9a1SkP/bv35y96JifTqXB0rtq9o1TW7FXdjWL9hSevW7F5DzYmapHWpkMSogaMSyp1zTPjFBJa881NfUA38rv57Zs6cSdeuicFSDKOtYIYVbwx+/r3r+eYLL9Qrf/wzMLFoPlQmT/A+67zHmHPtDQnlf6r+E6MrR7dKy9+VT+aFyYnJ7Gvralt9zcElI9l42x8oKUmsu27JTWw+sLHF1+xAJ969+QT9GkT8q67Zz/K1e6Azvqf6q/r1vXv3ZsKECS1uzzAKiXY/eHWy7iSXzPwSFz9QWa98bV/4l68CKVyqv/lN+vWsXJn+a77zNnz8cfK6I0dbd82TJ6GqKrG844lyhg0+A3aSYFQBJk2aRJcuyZeeGkZboSAMq6QvSKqSdFzSOkkXpevac36yiB9U/pazIm7VOuBbV8AJ68+3mF69oGMnB88l1vXv35/bbrst+6IMI8vkvemQVAIsAY4CtwN3AYslDXHOnUx5chMc2LePbbNup7Iaov3VdX3hrc1d6PPnbnTsEL9FR4/A/v3x9+1xzzF9euLP6P1H9lP8+gCOH4uXlZdDyWmJGmprYe+e+PEnNZuYPn16wvvq6uooe3VAvZ5nt+7QrVvyz/be7shssSNHufvu6STrKLo36uD9AaeOS05rPJ5MdTXxz+SKWLhpekKv/fDhwyx55ImEc3v16kVVVRW9evVKfnHDaEPkfRAWSWOBZ4EZzrl/lzQbuBv4knMuaeDO5gZheXf8eM5avDiteo3knDhxgk4WANwocJobhKUQXAGxxMe7w35X2NdL3ynpW5LWSlpbXV3dvCt//etpEWikZvny5WZUjXZFIRjWhsSGk+p1tZ1zC51zw51zwysqKpp3pc9/Ps3SjIY88cQTXHrppbmWYRhZJe99rEAsKVLMEdi/QXmr6dGjB/PmzYMDB7xD0lYCpY2KigquvvpqSktLcy3FMLJOIfhYS4Ad+JXmc4F/BU4AgxsbvMpIoGvDMNo9bcbH6pw7BowHDgPfB94Hxv9/ZwQYhmFkikJwBeCcWwk0nijIMAwjj8j7HqthGEahYYbVMAwjzZhhNQzDSDNmWA3DMNKMGVbDMIw0Y4bVMAwjzZhhNQzDSDN5v/KqNUiqxq/WagnlwAdNviv3FIpOKBytpjP9FIrWluo8yznXZDCSNmlYW4Oktc1ZqpZrCkUnFI5W05l+CkVrpnSaK8AwDCPNmGE1DMNIM2ZY4yzMtYBmUig6oXC0ms70UyhaM6LTfKyGYRhpxnqshmEYacYMq2EYRppp94ZV0hckVUk6LmmdpItyqOUhSfskOUlLI+XnS/p90PiWpMsjdVnXL2mIpGWS9kv6WNJLkgaFuqskvSPpmKTlks6OnPdtSbskHZX0S0kZz4Uj6bWg8UhINnlJPuoM7ZaE5+sk/TCU5dWzD+1uDxpj2/qm9KS63xnU2UNSpaSPJB2WtDJrOp1z7XYDSoC9+PxZt+Azwf4Z6JAjPQ/hsyQ4YGmkfD3wITAZ2AQcBLrnSj8wClgBTAmaHbAM6AMcA9YBU4GPgZXhnM+G970EzABqgcos3NMHgRuA74Q2t+SjztD2/UBNaP+H+fjsg6bt4fn/Q9jGpNKT6n5nWOezwElgHnAT8Fi2dGb8jyWfN2Bs+CO+MxzPDsdfzKGmgVHDGvlHnx+OJ4bjG3OlH+jc4Hg/PmXO7aH98aG8MhwPIv6FMSLUrQQ+AUoyrFX41TV/FYzW5jzVeSFwFLgzZljz8dmHtrYDPwHKImWN6kl1vzOo8ZzQxk+BzoQvnGzpbO+ugFg3f3fY7wr7c3KgpTFSacyJfufcidhrScOBXngD1FKtHYEzMqkV37urBl7DJ6G8Kd90SioCHgXmA69HqvLu2UeYAByS9L6kG5vQkwutQ8N+BP4LtUbSA9nS2d4Na0MU9vk8By2Vxqzql3Qu8Et8D2ZqC/VkS+th4HLgVvzPwNkt1JINnTfgf6lUEk/v3h3o1AIt2Xz2jwDXAP+M/7L6caT95ujJhtbisC8FrgVW4107DfP8ZURnQSQTzCDbwn5A2PdvUJ4PpNL4YYq6jCJpKPBb4Dgw2jm3R1IqrdG690JdLfFeQUZwztXi/aUvSRoH/C3wYp7pPAOoADZEyv4J6JdCZ86evXPuvthrSZ8FphO/P8n0lKWoyxTbw36Vc+5ZSRXAaOLGMrM6M+2PyecN34PZF27cJPxPgG3kbvDq74GZ+G/IDfifrUPC6+gAxiGgR6704w3B+3iDM4v4IEZfvKH9X+LO/1XhnGEkDgotyrDOMcB/4X2S9+J9pXvzUOdQYFzYvhvafzFoybdnfwHwPH7g51a8m+UI/ksgqZ5U9zuDOgVUhb/Tm4FXw7P8dDZ0ZuyDFcoGXAJsxP+keQMYnkMty8M/VXS7HvgU8Ifw0LcAX86lfvysgIY6Xai7GtgatK4k4vgnPgp7LPxzlmdY54hgjI4CH+FnLozIN52N3NvYrIB8e/Z9gV/hQ+0dAdYCY5rSk+p+Z1Br7N4dC/fuH7Ol05a0GoZhpBkbvDIMw0gzZlgNwzDSjBlWwzCMNGOG1TAMI82YYTUMw0gzZliNvEHSqEjEpFG51mMYrcUMq5FxkoSZS7bdi5/8/lrYDuVY87gQru/DED5wl6RXJE2IvMe+CIyktPclrUZ2eAO/4gn8csHYUsH1+InYALucc+uAkVnWloCkrwHPhMO9+OWR/fBLInfi1/QbRqNYj9XIOM65sc65kc65kfgoTjFOlTvnHk3WA5R0b6TscklvhqDVT0sqk3S7pPckVUv6gaRTnQVJnSXdHQJEH5cPzP1zSQNIzTfCfgXQzzl3kXOuDzAYHy6P0MNeFjlnWdAYq5ekyZI2hB7vQUnPhRgLMX3XRz7b1yStDgGW35Y0thW32sgTzLAahcRi/Brw04Dx+PXf9+GXrJbjA2/fEHn/Enw0q8HAW+HcbwCrJfVM0U7s/+JcYIKkwZLknNvqnFse6nYBb0bOeRPvwtgajh/Cx1S9EL8W/ShwBfB7ScnC0D0F9MYvsxwMPCPpghQajTzGDKtRSExyzp2HDwEHPnjJZfhANTtC2RcB5FOwfDWUfcU5dyE+ruYHwJn4eACNsRAfeb4Pvof6NrBX0mOSzgRwzj3a4Bq3hJ73v0kaiA+aAvBt59xQ4Czgj/hwgN9J0uZ/hM82FB8/tANwR6qbYeQvZliNQuL5sN8e9gecc6udc3XEDetfhP1fR877H0kOOIDv2UIKX65z7mV8AJfHiAc9Ph3fG14lqUsTOkcQD0/349D2MXxQkMbafiq0vYv4F8enm2jHyFNs8MooGJxzsZkCtWEfnTkQiyakBnuANSQGK363ibbewIcbJPRAZ+MDO5+JT5myurFzG7S9AW9Uo7yXqm0Sg0YbBYYZVqOtsiby+j+dc0+BH1QC/gYfRjApkqbge6pLnXOfOOe2S1qFN6wQN+hHIqeVRl6/jjfkAv7bOfdA5NrDiEe3jzIe2CipH/C5ULYp9Uc08hVzBRhtkjDIFMsU8KSkLZI24rOcrgBSpYq+GJ/hsybMKPgjPv0I+ClisUGrrfjg2QCVkl6VNM45tw1YEMq/J2lHmB3wIT5+6akU1hGmS3ozXLsrUIfPLmoUIGZYjbbMWHxE/s34waMB+FTH8/BBxRvjEfyg1Va8z/Zc/HzWRcAVzqd7wTm3Hx9FfyfQE+/X7ROuMSXUbcD7Z88G9gA/ws9WaMg1+Gj8xaHda51zVS3+xEZeYIGuDSNHSLoeeDwcnu2c2547NUY6sR6rYRhGmjHDahiGkWbMFWAYhpFmrMdqGIaRZsywGoZhpBkzrIZhGGnGDKthGEaaMcNqGIaRZv4PjPCDP24vx8cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T_list = range(T)\n",
    "plot_list = range(num_comp)\n",
    "plot_list = range(4)\n",
    "Loss_list = cumsum_cost_pca_mean\n",
    "plot_total_loss(plot_list,T_list,Loss_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}