class16_bayesnet_dataprep_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 16:  Bayesian Networks\n",
    "## This notebook prepares the discretized data matrix from the phospho-flow cytometry data from Sachs et al. (Science, 2005).\n",
    "\n",
    "In this notebook, nine Excel spreadsheets from the supplementary data files from the Sachs et al. article are loaded and combined in a single 11-column data frame.  Each column is discretized into three classes (low, medium, high) based on the values in the column. The breakpoints for the discretization are assigned on a column-by-column basis, selected (in a single global optimization) to maximize the sum of the mutual information of all pairs of columns of discretized data. The initial breakpoints for the discretization are the column-wise quantiles (0, 1/3, 2/3, 1). The resulting discretized data matrix is saved as a tab-delimited text file \"sachs_data_discretized.txt\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read all of the 9 Excel spreadsheets in the subdiretory \"shared/sachs_datafiles\" into a single data frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir_name <- \"shared/sachs_datafiles\"\n",
    "data_file_names <- dir(data_dir_name)\n",
    "\n",
    "names(data_file_names) <- data_file_names\n",
    "library(gdata)\n",
    "suppressWarnings(data_matrix_orig <- do.call(rbind,\n",
    "                     lapply(data_file_names,\n",
    "                    function(data_file_name) {\n",
    "                        data.frame(read.xls(paste(c(data_dir_name,\n",
    "                                     \"/\",\n",
    "                                     data_file_name), collapse=\"\")),\n",
    "                                   expt=strsplit(data_file_name,\".xls\")[[1]])\n",
    "                    })))\n",
    "rownames(data_matrix_orig) <- NULL\n",
    "data_matrix <- data_matrix_orig[,1:11]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "define a function for computing the mutual information of a Mx2 matrix of factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## df should be a M x n data frame of factors\n",
    "joint_entropy_df_cols <- function(df) {\n",
    "    stopifnot(is.data.frame(df))\n",
    "    stopifnot(is.factor(df[,1]))\n",
    "    \n",
    "    n <- ncol(df)\n",
    "    M <- nrow(df)\n",
    "    \n",
    "    ## use \"xtabs\" to tabulate the number of occurrences of each combination of logical\n",
    "    ## values for the n inputs; divide by M and call the resulting matrix \"probmat\"\n",
    "    probmat <- xtabs(~., df)/M\n",
    "\n",
    "    ## define a matrix \"hmat\" that is the negative of probmat element-wise-times log2 of probmat\n",
    "    hmat <- -probmat*log2(probmat)\n",
    "    \n",
    "    ## some entries of hmat will be NaN (do you know why?).  Set them to zero\n",
    "    hmat[is.nan(hmat)] <- 0\n",
    "    \n",
    "    ## return the sum of entries of hmat\n",
    "    sum(hmat)\n",
    "}\n",
    "\n",
    "mutual_information <- function(df_twocols) {\n",
    "    stopifnot(is.data.frame(df_twocols))\n",
    "    stopifnot(ncol(df_twocols) == 2)\n",
    "    stopifnot(is.factor(df_twocols[,1]))\n",
    "    \n",
    "    joint_entropy_df_cols(df_twocols[,1, drop=FALSE]) +\n",
    "        joint_entropy_df_cols(df_twocols[,2, drop=FALSE]) -\n",
    "        joint_entropy_df_cols(df_twocols)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "run a couple of simple tests on your mutual information code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"independent ternary RVs:  this should give a value close to zero: 0.000052\"\n",
      "[1] \"identical (as in perfectly correlated) ternary RVs:  this should give a value close to 1.584963: 1.584962\"\n"
     ]
    }
   ],
   "source": [
    "## this should give a value close to zero\n",
    "dftest1 <- data.frame(factor(sample(1:3, replace=TRUE, size=100000)),factor(sample(1:3, replace=TRUE, size=100000)))\n",
    "print(sprintf(\"independent ternary RVs:  this should give a value close to zero: %f\", mutual_information(dftest1)))\n",
    "\n",
    "## this should give a value close to log2(3) = 1.584963\n",
    "dftest2 <- data.frame(dftest1[,1],dftest1[,1])\n",
    "print(sprintf(\"identical (as in perfectly correlated) ternary RVs:  this should give a value close to 1.584963: %f\", mutual_information(dftest2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "obtain an initial ternary discretization of the columns of the matrix, using 1/3 and 2/3 quantiles as breakpoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "discret_breaks_initial <- apply(data_matrix, 2, function(data_matrix_col) {\n",
    "    quantile(data_matrix_col, c(0, 1/3, 2/3, 1))\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "optimize the ternary discretization breakpoints for each of the columns so that the summed (total) mutaul information for all pairs of columns is maximal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "final  value -6.834433 \n",
      "converged\n"
     ]
    }
   ],
   "source": [
    "compute_sum_mut_inf_given_discret_breaks <- function(data_matrix, discret_breaks) {\n",
    "    ## discretize the data matrix\n",
    "    discret_data_frame <- data.frame(lapply(1:ncol(data_matrix), function(col_num) {\n",
    "        cut(data_matrix[,col_num],\n",
    "            breaks=discret_breaks[,col_num])\n",
    "    }))\n",
    "\n",
    "    # compute the mutual information for all possible pairs of variables, and sum\n",
    "    sum(apply(combn(1:ncol(data_matrix), 2), 2, function(pair_col_nums) {\n",
    "        mutual_information(discret_data_frame[, pair_col_nums])\n",
    "    }))\n",
    "}\n",
    "\n",
    "opt_res <- optim(as.vector(discret_breaks_initial[2:3, ]),\n",
    "                 function(discret_breaks_internal_vec) {\n",
    "                     discret_breaks_matrix <- discret_breaks_initial\n",
    "                     discret_breaks_matrix[2:3, ] <- discret_breaks_internal_vec\n",
    "                     \n",
    "                     compute_sum_mut_inf_given_discret_breaks(data_matrix, discret_breaks_matrix)\n",
    "                 },\n",
    "                 lower=discret_breaks_initial[1,],\n",
    "                 upper=discret_breaks_initial[4,],\n",
    "                 control=list(fnscale=-1, trace=1),\n",
    "                 method=\"L-BFGS-B\")\n",
    "\n",
    "discret_breaks_final <- discret_breaks_initial\n",
    "discret_breaks_final[2:3, ] <- opt_res$par"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "make a discretized data matrix based on the optimal discretization breakpoints selected by `optim`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_matrix_discrete <- data.frame(lapply(1:ncol(data_matrix), function(col_num) {\n",
    "    myfac <- cut(data_matrix[,col_num],\n",
    "        breaks=discret_breaks_final[,col_num])\n",
    "    levels(myfac) <- c(\"low\", \"medium\", \"high\")\n",
    "    myfac\n",
    "}))\n",
    "names(data_matrix_discrete) <- colnames(data_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as a check on the discretization, for each column in the original data matrix, plot the original expression levels vs. the discretized categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2dC3zVRN73057SewuFlkopcnG7\n4lrUXWC5qICAouKigqCLjxfQIoK4ykVWUXS9sKCw7IKr8ujCvi4rgvqKj1AvFcFXfBAEXW9Q\nC+VUrlpKaYH2QG/znpxkknOZJJOcycnl/L+fT8/JSSbJ5Nf/L8lMJjMcAgAgajirMwAAbgCM\nBAAMACMBAAPASADAADASADAAjAQADAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIA\nMACMBAAMACMBAAPASADAADASADAAjAQADAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAA\nMBIAMACMBAAMACMBAAPASADAADCSCpyE1TkB7A6EiApeCatzAtgdMJI6Z2YVJFaUnLQ6G4Dd\nASOpM2PgZu7kqElWZwOwO2AkdTrs80t0INvqbNgLKDtGEpAChFEk45hfoup0q7NhL6DsGEnA\nOiCMIvfcVsPVT5podTZsB5Qdw8DXIBCGzNn5PT2dptRbnQ3bAWXHMLCRQBhAD1B2DAMbCYQh\ns7XN6hzYEyg7hoGNBMKQSStacMjqPNgRKDuGgY0EwpA59dqYtKvXnbE6G7YDyo5hYCOBMIqc\nWFmUY3UeALsDT440aPvq0d7tJ1udC9sBZccwsJFAGDJzeqVNeBvu7CKAsmMY2EggDJnRq09Z\nnQVbAmXHMLCRQBgy8KBaESg7BhNURgJhCMCDagWg7BiKZCQQhgg8qCYDZccwsJFAGDLwoJoM\nlB3DwEYCYcjAg2oyUHYMA1p/qwMPqslA2TEMaP0NGAHKjmFA62/ACFB2DIOu9XfL8wOyEjIH\nvNgaq2zZA3gFXxEoO4ZB1/p79i9e+XLf168Wz41VtuwBvIKvCJQdw6Br/Z29O/B1GB7YAgAR\nupuWzB8DX3VZZmYFAJwLnZHuHvRxXevJT4dDrR6AoOxIgq5fu8aZuf6F+XN9McwZYFug7BgJ\nbb92bdVVx+LylaUaCBiAArqLc5xWf/t5Gm5hABroQiROq7/9ZJXG5XUY0AudkeK3+juvweoc\nAI4Aqr/Vef3Wz/dDGYkAlB1DkYykKoxK9ferI+3GiEGDNdNcW06rj94y0qYpAhedH4ND1QVT\nXQyXHd0aL1gKdWFUqr/H9p1rLx7swnG95mgkyvqX/hCg413RSBkdYnK09LDVxXDZ0a3xgq2j\nJUxE9ffzOQKe3sYENY1Fub88P/M1jUSFtAHTOLurp2COgednBRfrX8dUFuWe25OdLobLjmPv\nN7iiWSwqOvRDbvS6YCOpC0Oo/q5YJ9DxAq1dxJirEtulJNyrkYg6YKb0/6Dig75amyNgOyON\n8vToxU4Xw2VH2xlp2niEBi3QSERvJHVhVKq/bRcwfTgugbtOIxF1wLTf7//Y10F/NmynS1/+\n1v0mjUTUuuguO67oJZB8If06MeGtlAXzkj7XSERvJHVhVKq/bRcwOfxxaF0mqQMm/aj/44iB\nF9hsp0s+r8tvNBJR66Kbb1cIdLBbUQAt+9XFa7TS0BtJHZXqb9sFTAIfMO00ElEHzPjRFc0V\noyfoz4btdEnkdUnVSBSHZUd05G9/P6GVht5I6sKoVH/bTpjAlTVBIxF1wBwf5+E8447rz4bL\ndXFP2bGyY5fc7jUaieiNpC6MSvW37YRhGzAINRw0VEHlcl3cU3a833+pTvyrRiJ6I2kJo9j6\n23bCsDaSQVyui3vKjnm8LpdoJKI3krowc39GNTel59wfZ1ckbq7+2imMm3VBbio7BsrUSRqJ\n6I2kLgy3B904YscXwwnPAGwnDMOA8dYiwy+wuVkXBGXHyM2I3+rCcHta2x1AqKJj5CK3ChMt\nrtcFyo6hm5Gm1ITh9rSkNiJUkyzP2oFbIf1SM6uxhW3AvJ3nK8sp/Fh/NtwaMNHiVl3o7v65\nGxf8eglqfWSoPGuN2C42pQfVFmIH24DpsbJ14LKXfqU/G24NmMCmoOwYuZnAh5Ywby6489LR\naFnel5GL3CqMSKLvZLqvwUO/+6ZKgXyb6LJcbI0+RdAF/1I4fvPKjriJUNK59OvEBJZGohTm\np0bCTLsYqalWRBAG/2ohp6Y20nmlrwxHZd3o8zFPPCVxXejXMZOHx4sIuuBfz5FTm3drV46b\nCNmtkXNMb+1UsIuRBnJExpJTUwfM2uSkjVszVtHnAxu6i010kejGy3G7RiIoO5KhN5JuYcrE\nU1vaL+jXMZP9O0UCzwUS8K+D5NT0Z976WuQ7bSA/tgsYdCOX8KhWGig7kqE3km5hSvGboOfR\nrxMTfuCFadJIFIe1UwiNm6GZJH7KjhKMjaRbGIwNA6bf7zSTuPwWxldL4ndTiLODRw6On7Kj\nBGMj6RYGY2nA7CgjcfkNxNnfBK3I6BZGpeNMK3WpTiYXGMlkB+U+fsqO/3ehiGAk/OtDcmp6\nI+kWBmNlwBzREy9cOyMBo36ltumbw1XcJ5XUrOWa5TXjp+w4FXcQJBgJ/5pHTq2j1s6JwlRx\n++kTf2ooYNSv1DZ9czgGuhjGLkaSyOaN1F8jkcurv6u4LWafedWv1DZ9czgGRnJm2ZFIIcdd\nppnGvOpvjLUBowv2tzA2fXO4int+HTVPGtLFPdXfjGozDVd/Y1xupDOzChIrSk4qLLXpm8Mx\n0MVNtbxMjeRIYWIQMDMGbuZOjlIeqdCWbw7HQBdH1vLWZOuRpYuuyim66m+bVvNWcb+fQs0Y\nQwHTYZ9fogPZCkttq4vpRnJkLW8V9z87qXlFny501d+Eat7j4tOZTsWGjytqqrhrxxMYeDlp\n7hWGAibjmF+iaqVX8G1b/f3MCmrmmFB2tO0JxrxKGLrqb0I171P4fFZAnzfWnH2M2OP5+eR+\n2p8JWpM6YO65rYarnzRRYWn8Vn+rlx0JJ5gz4ok+7yL6rLEmFkZSFyaG1by3y2NpFPWVp/9b\na70pfSVyOsvThh+wiZyd39PTaUq9wlKCLu/ixry9KPdgAjEwknrZkXCCsUMToVgYSV0Yk6t5\nm//yoMSIyyWyesrTN8kp3iNu442FEtMekKffJ++ROmDWKHkoAEGXTeL1L9PCVvExMJJ62ZFw\ngrFDE6FYGEldGJOreT/XVTbuFP0O6QMm9xO1pbat/jbdSOplR9s+XzPfSOrCEKp5f8RvPBp4\n9hTOZ7qMZKB7zwioA2bjkE2Vam8O27L6+6AuPROD3iJmVHa07QmGWDlFRmflFDaSujCEDiKX\n4nfwC40fGMbGRhL3qX8PVgZM2yd6WsVvC1qTUdmRcIIpF++1sy181bxKV5wZM5K6MPo7iGx9\nXa5gnb1UmvzvfaTEJ+8YJXHRLyVS8+TpgXKKRVrHREEcNs5k+2KfetmRcOJdIdb+tLOw85NY\nGEldGMoOIlu/lh5n/b+LLpBoVyBPPyo/8aom7mq0XON2TpE8vVjrQPRBP3zJQ92SCmaTen3h\n8XoVO41hbaQSuQbz/P7y9P/RWm86uTbzCXJqRmVHm/bMGwsjaQkT0UEkJliY/9GVz99SHbw5\nUAfM5D4bKjYW36WwtB8+lshFrANmkfxArGCoPL1Ra73Vcg1myX3y9HpyakZlR5v2zFvFLXmV\nwIsvk+Y+asxIGsJEdhCJCRZmra43MxnUUhiGvnaKfyLyXabC0rbJ9yityCJg9mZ5SCQkEmdr\njcZHA6OyI+WJN9bUpOiJzxxDbe3UhaHsINKFRso95P84SDizCmww1Ug2roRRh/LEy4KZcj1b\n35Hy9Fqt9R6Rb487d5Onl5BTs32xT7ODSBcaadXo8ubyUQryItR0SmmJy42kXnY0uWfe5idu\nlbjwFxKpneTpy+QUbxO38aJ8ezz+Vnn6NfIedQx9qSqMCiFGyniInitY3IoYRe8tDMdN17kH\npYCpkGtbynZIk7uIAxh484Pu3BIlEhLk6aAEQ3XmkASjsqOAaT3z2vEBPjYSjTBEgoX5sq8e\n/qB/Z8ygDhipXk5rlNFwQgLmgPTC+7cKt+mz5HfiyU2HC8jr3aczWxowKjuq4NYrtdSywUph\nYg/9g8fAp1JjXhWCdXlP1z/+cv07YwarsqMyLOJluy49Y3tFslSY2EMdMJd86/8o665/D/Fd\ndlSGxQP81lWLJK4fJpHdS56+R06h+mCHEnoj6Rbmafx/t/B9JONQB8zLnZ+tLcldTb/l04T3\nblxoJCm3tGXHb3ETod5BM09J3b0evLCHRFKePL1U7hGWPLDIBLnGrddv5OnllNmiRHf1N70w\n1fRvyE7dRbnN2EFfzXuwT/JIchsMMo9iIYPeu3GhkXSXHUlNhD7SIws3Qu/BMITeSGwK1WS6\n6jihxwjqgKkrKZjfeVGzdsJwQm7tOtC/+b3iJivvlWNbdlzbzj0nGKnRauAzSmEUcLKRCm+r\nRVXD+urfQ0jAdNax4mIDO2MGlB3J0BtJtzA6hulwspHe4T9a/6J/Dy43ku6yI8btRtItzGP4\n+LTfwXeykdT7slDB5UbSXXZ8Ew/e3VOet7b9QnpG2/uWV2oipFeY0/iKpN0rTOG/dWw3NjDr\nIFIRlxtJd9lxC+7LokietzaF/oXV8RfbWxdsJCaFagXe1WPQ2MCokw8VQoyUpaeywd4BgxOy\nKDvup+/e089L+nfGDHojMRHGOTDq5EOFYF0+66UHrQGTzSS2ZUfnQG8kEIaMel8WKrhcFyZl\nR+fAqoNIFdwqjIhWJx+KuFwXJmVH58Cqg0gVKIRxdMsGo7hcFyZlRwUcqQtdB5EquPyBrGFc\nrguTsqMCjtSFtoNIRdwqTLS4XBczy46O1IWug0gV3CpMtLhcFzPLjo7Uha6DyOj6b3OkMFSA\nLmTiUBe6zk8I/betxU0+ehDXOD1cfqc8uac02c8mGjEKGNCFTBzqQmckQv9t28QmUD2uIq7R\nskp+Vj99sTz9HdX+TIcuYELaTJISxKkumsShLpTdcSn33zZ2sGZzw1EjrzHEwOHG1rvsbs0s\ntacKmDSvDFmo+NRFm/jThdJIyv23/Vm7u6DElHRDJCYbW8+TrZ2n7TRHvTBo+inQRQfxp4uB\nAUt0Y7TwOGChdhoSFKMuUHOC/2h6lt0Gg3CyLmbiSF3ASOr09ztpa3FP7YQGcKwummXH6HCk\nLgEdQBhFHui3ryT5EWI/qFHjWF20y45R4UhdAjqAMIq0zeQu381ucyE4VhftsmNUOFKXgG9A\nGGXa5oww53rkbF2Es+4RA2+CauNIXaCMpIKZ97vIwbogSZrkG46y2ybGkboEQiSlBnVlsC0l\nHCkMj3zDy2iDoThWFz//uKGiae/Yf1f+/gZ228Q4UpeAkSaYe+ad9h9j6z1eZmy95WuMrUfC\n8AuPFDhZl04/+z9+6oKOG3hhQAtH6hKwTtthb76JZ14nY/iFR5eTVen/2JeBfjLBSI5EugZV\n+220p5uVWbElhl94dDnT+36w98N+9/xwpYE3b1wJNtKTgTu7Ky3Nix0x/MKjyzk7/9zE/Jm+\n9TPNuet1HthInbZsm9jyp8WW5sWOGH7hEYgvsJE8jU3F6EyupXmxI4bfBHU5J+5bj5ZMB10k\nsJF6P9xY/MWu9pbmBXAONw8rR19fcYvV2bAP2EhrPHtXexLMaeXsaI5BbSaJtJ/8H0czrM6G\nfZCMdKAJnai1NCu25I8mPl9zMh28/g9vjtXZsA/SYMwshqx1IWmlrVZnwZbc3/+DvZsGTrM6\nG/YBG2njkE2VcAsTSa7P6hzYk7OPnevpMhvEkQgbjJnlpqN/I8PL0W9DTCmvwOaFkL/d8T3r\nMpIrdDEBZ+tioqoshKmdq3Nv8gpsAmaePU8w1upiTqt4Z+uC1258qFtSwezGqLYVBgthokjL\nxkip7zMvIzlfF3NaxTtbF7z25D4bKjYW3xXVtsLgc3Z8QnqHqb4By9BR7nF0ut1OxWysOidt\n3rLctD/6777vzcq5/0zLEwVJF/2dv1Tzm+H/gpKQ9/Z8bvYzfMqz92XnruS84owoyWdfDHCF\nLiadeJ2ri9SJPv8+9XeZhrdDgD+am649Utlvzp9uQGuzh6EPzlE8u3N3Nr7H3d64gTuGHv+v\n6r0DZj/Vtaz+k+4hwshJyHub1vQeV+NPOb/vvqPX88IEZkR5EP+405QykuN1MenE61xdpOrv\nQ/6Pgx2NboaEP48NngqE3i/8okPLtEdTfQ8rv43AfY9aA39e1Mt/LJsLur3pn/tWiDByEvLe\njqI2Xg507gcI/cBPBWZEeRDmlAWcr4s5J14H64JDZNXo8ubyUUuMboaEP48HuWaEKjyteTuL\nv71gy8A3lLPhlf48fOAm8oryhycIs58LSULeWxu/zJ/asxehJn6qTTGxpbhCF3NOvA7WJaz6\nm+OmG91SBIEzjD+THxag2+d2ap02O61OORte6a9TPUInK7q/5Z+7XhCmDW2jEQaJwnQpQ2gv\n5xVnRHkQ5pQFnK+LOSdeB+uCjSTVw0R78yzD52zcdUcr+81Cr2dej97MHKaSDa/0VzK55tDV\ndz9T+FH9p0W8MMe49xtu1iPMnP6VP49jFTAmlQUcr4s5J14H62Lyc6Sa8entpzai2sTFqIZT\naRIbdNT1t6Rl3VHf8niXxAtLAw/YnsjMfVWPML6StPzV3BE2AWNSJYzjdTHpxOtcXQJGygif\n4Wje+QGhymRGPa6ZUBawCqa6uAg2ugR8Y3KXxTFm1lXVtRNZ9RJlQlnAKtjp4q4TLxtdAjrs\nCSb6nFnMyVszUsccYbQxE8oCVsFOF3edeNno4nwdzMWEsoDzcdeJlw1gJABgABhJBXeVBQAz\ngfhQwV1lAcBMID5UgLIAQAsYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABg\nABgJABgARgIABoCRAIABYCQAYAAYCQAYwNxIezjSeIhtt6SihWlGduZ1vtXJkoiIxxev4gRp\nE3EwjtLEBCOR5h7me2A2dIgUA3UcyDew3RhClkQEixKn4gRpE6mAkzSJkZECkph1rrD7edk8\nI1Fgc3HYG4kCM7ZripG83D87p09qRFxp3prXentyn0b5XFc5/8JAHOsL69DLl5wNTissQOuK\nPL3WI/GLX0kY7ANxazuL4240PZCTfksdQkEbF9YVh+iQ1ijNy52CUE27fawPUw8hkoSPR+I/\nvrYH+p+IV3GCtPEfjJh1nFEnaWKSkUYf/L54NuKGbP7R88fTH3PVIVckYSAOdFvJwbxvgtPW\nCQvqkzc2/iMPiV/8SsJgH4ib2LBBGHdjfv/9R4fOQHVBGxfWFYfokNYYsnl9fitaMZj1Ueoi\nRJLw8Ui8XNvMfifiVpwgbfw5FbOOM+okTUwyUjlC7xYibhNq2tvYtF3sWBnrIgzEgWoL+jwb\nklZcUOtZebLVh8QvrzTYB+K24XE3ur/v39H6kI0L64pDdEhrbEIteVvRsBdYH6UuQiQJH4/E\ny83h+NNfnIoTpI0/p2LWxYw6ShOTjMSPzpGIuEqEVvcpviXMSMJAHP7TRHJtaFpxwYZh6YM+\nQuKXVxrsA3H7cefM/KAFPEEbF0f3EIbokNbwb3TKzMMp5LGmYkWIJN6w8Ui83PDLS1DcihOk\njZfDWRcz6ihNzLsilRbwB1GXtIM//YYYSRiIA1V1HlkSklZc4NuKGv/WEYlfXnmwD7lL9C4f\nI/T14pCNC+uKQ3QEr/Fhz6VjWB+kPsIkCR2PxMs1fJf8WdyKE6SNP6di1sWMOkoTk4w06tBu\n/jbWi2oTvjzzJFceYiRhII7WYYuqc8uC04oLTmd85FvRFYlfQYN9yLrMGlT10/Bp8sbP4HXF\nITqC12jKOWcd64PUR5gkoeOR8Mf3cJ+meBUnSBv/wYhZxxl1kiYmGWlxTsbkxkC2FmXnPjei\nMMRIwkAcS/s2o9XdTwWnFRagVYWJvf03wMJX0GAfsi6+qe3Tb67DG2/onY7XFYfoCF4DTcpm\nP6CyLsIlCRmPhD++hh7Pxqs4QdrwI6wIWccZdZImJhmJuIg4m22dPmmIjplMhwkzgIokMnEq\nDlEbnFEnaeIyI0UO0XH2SOctLPdgALsYyY7iELXBGXWSJiYYiUM1xDFQniJ2+0tOa5TIITo2\nJU1pY7kHAyhLIhOv4hC1ETPqKE3s3YAEABwCGAkAGABGAgAGgJEAgAFgJABgABgJABgARgIA\nBoCRAIABYCQAYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABgABgJABgA\nRgIABoCRAIABYCQAYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABgABgJ\nABgARgIABoCRAIABYCQAYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABg\nABgJABgARgIABoCRAIABYCRAN5yE1TmxDyAFoBuvhNU5sQ9gJMAIZ2YVJFaUnLQ6G/YBjAQY\nYcbAzdzJUZOszoZ9ACMBRuiwzx86B7KtzoZ9CBgJCo9kQBdFMo75Q6c63eps2IdAiEDhkQzo\nosg9t9Vw9ZMmWp0N+4DPtVB4JAO6kDk7v6en05R6q7NhH7CRoPBIBnQBqMBGgsIjGdCFzNY2\nq3NgM7CRoPBIBnQhk1a04JDVebAV2EhQeCQDupA59dqYtKvXnbE6G/YBGwkKj2RAF0VOrCzK\nsToP9gGekACGaPvq0d7tJ1udC/uAjQSFRzKgC5k5vdImvA13djLYSFB4JAO6kBm9+pTVWbAX\n2EhQeCQDupCBB9VhBJWRoPBIBnQhAA+qw5CMBIVHMqALEXhQHQY2EhQeyYAuZOBBdRjYSFB4\nJAO6kIEH1WFA6291QBcy8KA6DGj9rQ7oAlABrb/VAV0AKqD1tzqgCwF4BT8SaP2tDuhCAF7B\njwRaf6sDugBUwMUZABgARgIABkC/diqALgAt0K+dCqCLMjWgSwhwrgWM8DRcqUMBKQAjZJXC\nq8MhgJEAI+Q1WJ0DmwFGAozw+q2f74cyUhCSkaDwSAZ0IQK1mWFgKaDwSAZ0AajAIQKFRzKg\nC5nG2V09BXN8VmfDPmAjGS48vjrSbowYNFgzzbXllIcHupCZ0v+Dig/63gu6YLCRDBcex/ad\nay8e7MJxveZoJMr6F+XhuUiXfI7rwUyX9vv9H/s6uEEXNvGCjaS78Phn3Hqml34xTWVR0aEf\ncl/TSFRIGzCGC9Vj79e/jqnMTUjyJP63RiJqXdKP+j+OGHhPy3a6MIoXw8Xow2UCnS40ugWT\nmDYeoUELNBJRB4xhbBcwPThPInepRiJqXcaPrmiuGD1BfzZsp8u08Ud+ZhAv2EiGC48FF+tf\nx1TeSlkwL+lzjUTUAaNbl+p1Ar+5nX6dmJDJFZ/HnauRiFqX4+M8nGfccf3ZsJ2RXk/kuIT/\n1UhEbyTDhUfbGan1ioTEW7QSUQeMbl3+nCOQWEi/TkzI5W9RtW4fdFypGw4aqoexnZGWdjyv\nKP0NjUT0RjJceLSdkVZ2WPxM8ocaiagDxj26dOWNNFgjURze8t57M0KXPqORiN5IhguPtguY\nu1I4LvVJjUTmF6ptp0sOl5LMnaeRiE4Xbq6LKmH+nf3i4qRPNBLRG8lw4dF2AdOb/weP1khk\nfqHadrpcwOV24q7WSESni7cWGX5Py3ZGQn/qcu6LWmnojWS48Gi7gEnkjZSqkcj8QrXtdFnC\n68Lsllc3TZUCV003aw9moqf622Dh0XYBE7jlSNBIZF6hegd+hPdLPWvFgG5cSjtuoEYial3e\nzvOV5RR+TL/7efi5Y1f6deyDic+RMK43kk5wE5iU7mbtwSDpvC5dNBJR69JjZevAZS/9in73\nTbUCXewWL1RQGimawqObjeQqXQJG0qqTpzZSou9kuq/Boz8bttOFCkojRVN4tJ0wDI3kKl0y\neV2YPZA9r/SV4aism/5s2E4XKuDWTgHzn5fYTpd2vC5aY3hS67I2OWnj1oxV+rNhO12ooDeS\n7sIjxnbCsDWSe3RJ4HVJ1khEf4Kpr0W+0zp2/zxu8WHgKmYqbbekpv1BKxG9kXQXHjG2Cxi2\nRnKZLokaicy7UleKbRA7XmDWHgzyCK+L1mHTG0l34XGpXduUsTWSewrVcKUmk8/rcpFGInoj\n6S482u0M84P4WkeZEDD414/k1HFYqIYrNZlUXpfOGonojeT4wuNAjshYcuo4LFTDlZoMo0oY\nqdZOb+ERYzthAoXqNI1E5hWqMbbTha2R3HOlDjQp02qVHI/V37yTkrTSxGH1N1sjuedKncCm\nEoau+rvl+QFZCZkDXmyNXGSpMD9Xkrj6duLsmqAV47BQzfj5GlypwzYjfqsXHmf/4pUv9339\navHcyEVWCvOTh1wwIpMWdBpgVKgmnGBwK+d8mwTM4Z0iQsDgX8fIqal1OTOrILGi5KT+/Ljd\nSOqFx+zdga/DhDKZlcJUca99TODDj0hzl3HN8pqMCtWEE8xj2Lda7UNjxK/JZ5XryKmpdZkx\ncDN3ctQk/flxu5HUC4+ZQi1yXVbkImuNpAsjRlLXhXCCOY2vSFqPJmJEg9jqujZQFvDgXwrd\nuVDr0mGfP3QOZOvPj9uNpF54vHvQx3WtJz8dTjgDudxI6rrY9ARDJkE7Xuh1yTjmD51qpcou\nu5apgxhEjpFx5NSsqr8bZ/J90OTPJZzHXG4kdV1seoJRYNwMzSTUutxzWw1XP2miwlLCLe+7\n4wXSbNKh6G78yD5wpZYe4FeSU9MbSavw2FZddYzYnby1RhoyisClxLkDDBlJXRebnmDqe+aQ\nSE4hzu4T9H+l1uXs/J6eTlPqFZYSbnk3iW8OZ/7C0CGZx17eSIRLZwj0RlIvPNr0Ul1/US8S\nGe2JswcaCRitQnXECeYwbjplYQ+0VdxdDxK4gzj3ZkMnmDVKHgrgqFteNPhGzST0RlIvPMaw\n+vubMok178rTP2mtt2udxIBr5OnvyakZFarn/oxqbkrPuT/oivQcbsxrYd8EMbjlzVXtwMqm\nt7wNY4lDTXTuRpx9m64TLyWX3ecAACAASURBVDaSeuExhtXfheR/tmYnTsPke5XUNHn69wp7\nYVOo5vagG0fs+GI4IXcuLztuHLKpUvnNYZve8lZxDy0kUHIfae5UfbpgI6kXHgmX6hXiDVOS\n1tvLxhmw0LRNMypUc3ta2x1AqKJj5CKXG0lcVXF5xC1vuRif2Ra+LVDFFfel5nxjRlIvPBIu\n1eUrBDowEObodcRra/si4uwHot8hq0I1t6cltRGhGsKLp1Ya6QCXmkagXTvS3JSEFnlNRm0Q\nCbe8L4jx2c7CN2TNPMFgI6kXHk2+VH/GPUO6uN4+izR3Yn70O2RVqOZuXPDrJaj1kaGRi6w0\nUts7ciHx1eUSlwyVp1fLKYI7jaQfpeOhbkkFsxsVlsbhLS82knrh0eTq78+439BfcnsR7qN0\nw6hQ/eaCOy8djZblfRm5iHXAdCP/szV7G+hLXm8MOTW1LpP7bKjYWHyXwlLb3vKmpRNII85N\nNWYk9cKjCmyMpAcDQ0NEwKhQLfAT6cTMQpe2Uvm68Vf5mtztGnn6ZTnFf4jbqJabv3/znTxd\nS94jfSUMX/30XabCUpve8tZdoOdxSV9DtXZikOrPnMuNZKku27lEEgkJxLm50e9Qx5X6kP/j\noNLNQQxveV+XTyr3zpOnCfcIobwqDxB7ftCwtuvJqVm92Of1KnaU6HIjGcbluqwaXd5cPmqJ\nwtIY3vJeJd/0ZxTK05oVvtPl2qvz+8jT88mpdQx9qVp47If/WZGL2ATMa+sI/G0Vae4DMa1s\nUNdFBZcbSdqryuASpt3y7s3UcaVO/HX0O9RhJPXCY9vke5Q2wESYfGIjsKQ04uyh0e+QVaFa\nBTa3dnqI6a2ddI9So502BLeeYKSWDaqFR7RBt5H275TY8oU0ueusVoaq5PWKZ8jTdVrr6YNR\noVoFJV02B9UfyNfhN4jx2PalLMDyRyQKhsvT/1dOodBwWRf0z9cCnxa9IWtnI6kXHlHTKaUl\nCsKcSiEf1dNaGcojr6doZGMwKlSroKTLOfKlNTFDnlYqbUj8Vq5OysqXpx/TnzM1qHW55Fv/\nR1l3/XtgYaTPuQQS5LkJnaLfoQ4jqRceVQgWZt+E6yVGXyuReok8PUZO8Vf9O2MGo0I1AVIP\ntLVLSE+WF2bfTJz9ka4jYQu1Li93fra2JHe1/j0oGWmH3ED51fekyY+I1fTN6+XL+twpEp0H\nyNNL5BTb9GcyAt3V36qFRyLBwqzV1RmJgW46mcG0UB0CqQfadZ4iEtnnkubm/UbvwTCEvjbz\nYJ/kkdU6tnxavAHNI7+CX59JDhLN2rfB8tU5+xx5Wms0bp3QG4lJ4XFtsuuMBLqQqSspmN95\nUbN2Qsyj+PCCOoVx05VaarQa+Iyy8OjCgAFdFBLeVouqhvXVv4dgXdYlXUKi4y9JcwsM7IwZ\n9EbSXXhchVvzBq2zNnM8PRf30XMkjIltodqFRnqH/2j9i/49uFUXbCTdhcev8HX4fHneWq1O\n/YNZbO8zjIhuXUgdRK5tT7xXITPaynexY9tBpBuNpLvwiAkRxn1G0q3LPPx/DyoLrO1E7EP5\nB+LcRxyhC5MOIt1oJN2FR4zLjaRblyaxB8YuQbr8j5544S7TdSRsMa+DyCWExwJuKgpgIzEp\nPLrQSEx0aSVeeirz/0KcfVz/zpjBqC8LAgfFJ0OdgnpXclO8YCOxKTxmraDnJnsLI8JEF6VM\n/Fv/Zk2GUV8WKgTr8pauK/Vv9e+MGew6iFQkWJgd9O+5+rlP/86YYYtRF941UiY1F0Z9WagQ\nrIuvjMhTbxJnl+vfGTNYdRCpgg07/KMgDkddoMIZ72nFHlYdRKrgVmFEQBfGuFUXug4iVXCr\nMCKgC2PcqgtdB5EqUAgzdZf+zZpMbAvVCjhZF8O4VRe6DiJVoBCmq4HW9iYT20K1Ak7WxTBu\n1SUWo5o7UphoAV3IuFUXMJJJuFmXkMc7OvfgVl1i0R2XI4WhIk51SfPKEAMoDnWhMxKhO661\nYkdgKT2Ia5weLj95Te4pTfaziUaMjETQ5U2sS0/iGq7QJfi11adICeIwXuiMROiOa5vY6r/H\nVcQ1WlbJjYGmL5anv6Pan+kwMhJBl8/ELjsLRxLXcIcuP13brvgzhPYoRU8c6kJ5i6vcHdfY\nwZpv2IwaeY0hBg43tt5ld2tmqT2jMlKc6jLmxm9fyf1R2UhxqEtACu3Co3J3XH/WblSXmELq\n7V+bxGRj63mytfO0XUsY0EWFlKMIPXGHipHiT5eAFJqFx+gwWng0OmLfuBnG1osAdFEg/RBC\np7p+rWykqHCkLgEpNAuP0eFIYXhAF6Utjdx6Fq05fy0YSUKQQqvwGB2OFCbAnt8mXfpZb88Q\nL6sNhuBcXapvSN2L0LJMiBcJQQrNwmNUOFKYAJfdV/5o8hN7pl7JaoMhOFcXhNr4Ubh8X7Hb\nYBCO1EWwjmbhMSocKUyAdvWonvMRR55jgHN1CdA0j+nmZBypi2AdKDyS6bgbtS1HqDKL1QZD\ncK4ugeqXcs7IUKkUOFIXwTpQeCTzQNEC/+f2Kyaw2mAIztUl22BDOzocqYsghbmFx2nkYYI1\nebzM2HrL1xhbL5LmF3gL3fX7Y6w2GIJzdansP6UB+cyJFofqgsUws/AIuI6zs4q/Nc1IjgTE\nAAxReu5fIXaCCIjx87jMlOuOWJ0V+wG6qHB4ZIbVWbATASONGfrFzivNeVbiaEAXgJaAkVK/\n8ZcfPVZnxX6ALuqc8lqdA/sgtP6uRqgR7ngjAF3UWQXaSAhGOobMq8x0MKALQItgpF3Cc2qT\nnlQ7FtBFDdOaCDmSsBf7rM6OvQBdFDC3iZAjgRAB9GNuEyFHYqIU0b9U6uXotyGmlFcw451W\nJrhAF1OaCDlbl8Da8gvVUW0rDBbC1M7VuTd5BQYBA7ooYUYTIWfrYmIZiYUwUaRlEDCgizLs\nmwg5WxeTb+2OT0jvMNU3YBk6yj2OTrfbqZiNVeekzVuWm/ZH/8nu3qyc+8+0PFGQdNHf+Us1\nvxn+LygJeW/P52Y/w6c8e1927kp/QViYYT9cogvzJkLO1sVkI9107ZHKfnP+dANamz0MfXBO\nq2I27mx8j7u9cQN3DD3+X9V7B8x+qmtZ/SfdQ4SRk5D3Nq3pPa7Gn3J+331Hr+eFCcww7fAM\n4xZdWFd/O1sXExut+vPY4KlA6P3CLzq0THs01few8giS3PeoNfDnRb38x7K5oNubiB+sN1gY\nOQl5b0dRGy8HOvcDhH7gpwIzojkC0EXpEMyo/na2LiY2WvXn8SDXjFCFpzVvZ/G3F2wZ+IZy\nNrzSn4cvkyTyivKHJwiznwtJQt5bG7/Mn9qz13+25KfaFBNTArooYEr1t7N1MbHRauAM48/k\nhwXo9rmdWqfNTqtTzoZX+utUj9DJiu5v+eeuF4RpQ9tohEGiMF3KENrLecUZ0RwB6KKAWdXf\nDtbFxEarfM7GXXe0st8s9Hrm9ejNzGEq2fBKfyWTaw5dffczhR/Vf1rEC3OMe7/hZj3CzOlf\n+fM4JgEDuihhVvW3c3UxsdEqn7Oa8entpzai2sTFqIZ7ViUbXumv/pa0rDvqWx7vknhhaeAB\n2xOZua/qEcZXkpa/mjvCwkigiyLmVH87VxfBSK5qnPnOD/57j+RmBlsCXVRwzRuybHQx8YGs\nVcy6qrp24g0stgS6KNL4ULekgtmNjLZmKWx0cUWIhHLy1ozUMdDTQgRMdZncZ0PFxuK7GG3N\nUtjo4kIjATEgY7f/47tMq7NhHwJGygifAQQAXRTJPeT/ONjR6mzYh7AykkuKA4wAXRRZNbq8\nuXzUEquzYR8C8bEnGKtzZCNAF0Wk88t0q3NiE+BECxhBelPLhq2CLQGMBAAMACMBuoFKmEhA\nB0A3UAkTCegA6AYqYSIBIwEAA8BIAMAAMBIAMACMBAAMACMBAAPASADAADASADAAjAQADAAj\nAQAD2BppD5chNCBJvf4nsQ8X/ENgYZqRHXqdbnfQJYKAJHQ4QhzWRuI3uetY9ddDx4oBg39U\n35SResPPxg6PYrCOA/kGthszVHRpebBT8uXl8afLnpADrsz2Bn1Jkw4KGjOM5PV/7MgI6lWM\n/3HjiANHbrjWtPOEvU/OKrq8WPR93fQB8adLiJHO9g/oIX7JcxwUNGYZaVdKUMD4f7SlbUXo\n23ZS3oXBONYX1qGXLznr5f7ZOX1SI+JK89YIC9C6Ik+v9Uj84lcSBvxA3NrO4tgbTQ/kpN9S\nh9BrvT25T6N8riveqDhMh7RGaV7uFIRq2u1jeqi6UNYFTVjgX5wUf7rsCfR1Lxwfmj0loIf4\nhfCkk4LGDCOV+xr3XzVaDBjxR9veswi9cqGkiTAYB7qt5GDeN/6Zow9+XzwbcUM21wkL6pM3\nNv4jD4lf/ErCgB+Im9iwQRh7Y37//UeHzkB1nj+e/pirDmxXWFccpkNaY8jm9fmtaMVgpkeq\nD2VdeE49Njj+dBGMJBxfaXEDr4n4xSNMOilozDASPzrAiCNyoTrwwx8vT+RulTQRBuNAtQV9\nnuUFLUfo3ULEbcILaj0rT7b6kPglD/iBuG147I3u7/v3tx417W1s2i723CysKw7TIa2xCbXk\nbUXDXmB6pPpQ1WU5l/xZ/OkiGClwfIfP/ZbXRPzikScdEzRm3dohFN7z8utdb/9Rvi0VBuPw\nnyKSa/mZ/GgeiYirlBZsGJY+6CMkfskDfiBuP94mP3ABz+o+xbeImogjfAjDdEhr+Dc6Zebh\nFPJ4U7FBTReEfPMviD9dBCMFjm/kUp+PK28Sv/iF0qRzgiZmRlp+/n/4L6yJMBgHquo8skQ8\nuZQWBBILC3xbUePfOiLxK2jAD3mbXT5G6OvFqC5pBz+YVGC7wrriMB3Ba3zYc+kYpgeqExVd\n7t4eUkaKG12kK1JpQUbgTdt54he/EE86KGhiZSRf1uuB1ymxJsJgHK3DFlXnlvlnjjq0m7/d\n9eIFpzM+8q3oisSvoAE/5G3OGlT10/BpqDbhyzNPcuVe7gxeVxymI3iNppxz1jE9UJ2oGGni\n8P3Hpw2JP10EIwnHhyQ9+K9Vp/Ckk4ImVkb6SnzBH2siDMaxtG8zWt39lJdbnJMxuTGQWFiA\nVhUm9vbf/ApfQQN+yNv0TW2ffnMdQouyc58bUdjQOx2vKw7TEbwGmpTtY3qgOlEx0rHxmclX\neuNPF8FIwvGhECPJk04KGhOMpAqx6p7tcwLSMB0zre3uHXSJQDASXVoniONCI0UO03H2SOct\nLPegG9AlApsZKWpxWBtJY3tPERPUMO33NnKYjk1JU9pY7kE3oEsEvCSUx+cIcezaggQAHAUY\nCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABgABgJABgARgIABoCRAIABYCQA\nYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABgABgJABgARgIABoCRAIAB\nYCQAYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABGAgAGgJEAgAFgJABgABgJABgARgIABoCR\nAIABYCQAYEDASJyE1dkBAGcSsI5XwuLcAIBDwdegM7MKEitKTlqaFwBwLNhIMwZu5k6OmmRp\nXgDAsWAjddjnnzyQbWleAMCxYCNlHPNPVqdbmhcAcCzYSPfcVsPVT5poaV4AwLFgI52d39PT\naUq9pXkBAMcCT44AgAHYSFvbLM0GADgbbKS0ogWHLM0IADgZbKRTr41Ju3rdGUvzAgCOJaiM\ndGJlUY51GQGcA7TNjESSou2rR3u3n2xlVuwHBAwZaJsZCQ6ROb3SJrwNd3ahQMAoAm0zw8BG\nGr36lKX5sCsQMGSgbWYY0PpbHQgYMtA2Mwxo/a0OBAwZaJsZBrT+VgcChgy0zQwDWn+rAwFD\nBtpmhgGtv9WBgAGogNbfgBGgbWYY8KhRHQgYMtA2MwwwkjoQMGSgbWYY0K+dOhAwikDbzGCg\nXztNIGBIQNvMUOAapAEEDBFomxkGGEkdCBgy0DYzDDCSOhAwZKBtZhhgJHUgYMhA28wwJCPV\nQGUDCQgYMtA2MwxspKeh+psIBAwZaJsZBrZOVik8wicBAUMG2maGgY2U12BpNmwLBAwZaJsZ\nBjbS67d+vh/KSJFAwABUYCNBEyEAiAKwDqAbaJsZCZaicXZXT8Ecn6V5sR0QMGSgbWYkOESm\n9P+g4oO+91qaF9sBAQPQgo3Ufr//Y18HK7MCAM4FGyn9qP/jCDwvAQBDYCONH13RXDF6gqV5\nAQDHgo10fJyH84w7bmleAAcBbTNDkeujGg5C4wYCEDBEoG1mGCCFOhAwZKBtZhhC5ydzjbds\neHWk3RgxaLBmmmvLKQ/PcMC4XBfDbTPdqovQ+UktMvy8ZGzfufbiwS4c12uORqKsf1Eenv6A\nqRX4ne10ac9x+cx0Mdw2063xEvVNy9j7o90CYxYVHfoh9zWNRIWmBcyjuDFEIf06MWFSQu8i\n7hmNRNS6GL6DsV+8pCZ6UqKPFyzF23m+spzCj/Xnw3bCTBuP0KAFGonMC5jTlQL5F9GvExN6\n8sdxuUYial0MY7t4OY9LSOBGaCSiN1KPla0Dl730K/35sJ0wb6UsmJf0uUYi8wOm4GKz96CT\nZN5IeRqJqHUx3DbTdvGSxOvSXiMRvZESfSfTfQ0e/fmwnTBo2a8uXqOVxvyAsZ2RAlfWRI1E\n1LoYbptpu3gJ6JKgkYjeSOeVvjIclXXTnw/bCYOeLDj/Fa005geMPY0UfcCI6G6b+fVCgfOu\no18nJjA20trkpI1bM1bpz4ftjPRaUlJS4icaicwLGIzLjaS7beYqsSI5pTv9OjGBsZFQfS3y\nnTaQD9sZ6WJemHEaicwLGIzLjWS4baZbdXFf9XciL0yaRiIIGDLUuhhum+lWXdxX/c22UO34\ngPl0hYgQMPjXl+TUOmozDbbNtIsuEoyNpL/6Gz/Bt6WRmJ15HR8wY3qJCLrgXyXk1HH4WICx\nkXRXf9v2CT5rIxnEdgETuOXN0EhEp4uRtpkviA5OMlAzbCpWV3/b9gk+QyNF05jXdkbinZSq\nlYZOFyNtMyvWCXS8gH6dmGCX6m/bBQxDI0XTmNdKXXw3EZswd+5GnH1n0JpxeKW2S/W3W4WJ\nFit1qeKGjCLQdwBp7gCuWV6TWhfDlVNujRdsJP3jAImVDV1cKoyIIwOmitOFESMZbpvpdiPp\nHgfoKfxvKKBfJyawNZLugHnaBrrEwEiG22a63Ui6xwE6vlMgr5h+HTOZ0ldEEAb/mkdObV7A\nVJcJdLJQlxgYyXDbTLcbyfA4QHYR5g2xUeRCQRj8631yancHTLVHj4/SW+U1qXWByqnwzYjf\nhscBsp0wj/HC7NBI5LyA+V1fiS5F8vRSYuL6Wonre0ikZ8vTJXKK4NGm6SthoHIqbDPit+Fx\ngGwnDPpPuzTNJj3OCJjqXu0kPDIJCUE/5BQ3EbfxyTqJJcvlaYUzjbtrM4kwNtIao2Np2U4Y\nhMbN0EzCKGBanh+QlZA54MXWyEUsdPlMV1GHRcft7q7NJMLYSLlab/AoYaUwp67oSyKnM3H2\ndUEdazEKmNm/eOXLfV+/Wjw3cpHLjeT46u/9YqVQmWAk/OswOTW9kTYO2VTpuCf4VdwzKwg8\n/Dhp7hwTnpdk7w58Hc6JXKSgy9mS8RIXXSNNTviQlLh+6DkSaSkSCR55ur2cYjblIanBqDbT\n5Cs1C/qSz0ZjyKnpjSRuR3+GrDXSfvrEn5rwvCTzx8BXXVbkomBdapfgSsSFTw/+rUTyefL0\nXVKKhR8Rd7VI7mRtzGR5eiPlcVDCqDaTcKV+Pkcg0W6NVhP4uGfX+YlhrDXSv8qoWWrC85K7\nB31c13ry0+GEB9nBuqxLIt5rduxNmlvYj/bwTYBRbSbhSl1p10arvJPaaaWhN1LjQ92SCmY3\nKqSy6aW6SlcZwpCR1AOmcWauf7v5cwm9DAXrsrb9QnpGW3nvw6g2k/JKbROYVE5hI03us6Fi\nY/FdCqkIl+pqfIZx9xN8rervtuqqY8TewUOMlKwnmwZK8Myg1kW9bSbllTrmVOwkMfxm4uyq\noBV1tGzgL8bfZSqkIlyq/4zvebsaOiQmxMBI6gFDeaV2oZHU22ZSXqljzU+Jev4PqbpafEjV\n34f8Hwc7KqSy6aW6inuI/pZpqiEjqQcMZfW3C42k1TYz4kptjzuYr2qpKdUXL9hIq0aXN5eP\nWqKQyqaX6oNc+xwCqemkuVmJLfKajAKGsvp7bUf6OpGye35NmTUzoNZFvW3m3J9RzU3pOfcH\nXZHscQdjXi1vWPU3x00npbLppbqtdB2JAdcQZwc/VWUUMIQr9VIcMEF9WaztrOOgFvfVkZg1\n1Lqot83k9qAbR+z4YjihYxxr72CuHU/NFcaM5MXUkNNFXKqPm/S6wMheEh26yNPPaa03Rr72\nJKfI0/eQUzMKGMKVmlTN60IjqbfN5Pa0tjvgL9sTigpureWVGq0GPpUK1YRL9eN4d10MHlUw\nhysllj8j0W2UPL1eTlFH3MSn8qVn6d/l6Z3kHTIKGNrqb/cZSb1tJrenJbURoZrkyEVuN9Il\n3/o/yrorpTL3Uv2NruNj8eopo4Chrf52n5HU22ZyNy749RLU+sjQyEVuN9LLnZ+tLcldrZTK\n3Eu1jRtnMmnM60IjqbfNfHPBnZeORsvyCP25WmmkIwl64qydoepvdLBP8shqxVTmXqptbCQm\njXldaCTxf6Ga5idSQxnWRnpJbnc4Yao8vZWY+IBcQHj6Foluv5Snl8kpjgStSG+kupKC+Z0X\nNSul0n+p/lnO0Xf7pMn9TaS0Ry/oKNEhW8KTKk/nyCkM9GgfAaOAkapoIhe53EiGYW2kcXIH\nfe2L5Om/aq33sFw/1/9SeVqhUoveSIW31aKqYUr/Q9pL9SnpcdaP7ciXk4flJ14tkRvzk09e\nz8BQX2owCph+OHuRi1xuJPW2mSqwMNKxay8lkdWdOHtW9DvUYaR3+I/Wv6im1bxUb9Z1DzqE\nuJNG2Wk/Bz1nVrpWGoRRwLRNVqhed72R1NtmquDWooDxDiJFQgJG4TJExhFNYTQCZkO8Gkm9\nbaYKbjeS7g4iMS5vU6YRME2nlJa43EjqbTNVcLuRdHcQiQl5gS0xnURSKmluch/9O2NGbAPm\nLV3/+AH6d8YMal3U22YSOIxbfFxoKGMhfM4lk0hMIs1Nyot+h7HuILKG1FfCihUdJhNnb9K/\nM2aYFzCYYF2aNulotFpWqX9nzNBbm6nQNpMAy0arzeuJjSmLbiXO3hb9Ds3sIHKHWGGf9UvN\npF2VnvNah3kBg7Hhm6AUUOui1TZTEQVdWhbLz4BunCFN/vEbre09OUWi82/l6Q06s6WBiR1E\nvko/3HuvN+g3GyMsCxibQ98GMfAZZeVUEA2j5WdAKcXy9Bqt7d0jPwO6+Ap5eqX+nKlhj85P\nqslPjKzEvID5Bx7isbtm0qm76DcbI6h1UW+bqQJFvDjyDsbRvQgZx7yA+RYXDHtrJnVkwIio\nt81UITheDk2fQiJjOHH2q/p3xgwwkgKxDRgFnGwk9baZBN7E92095XlrU4iv1F0wijT3Yns/\nFjBspG/FrhCytc+8NsS8gMG43EjqbTMJbBGrDzKL5Hluer5m2Ej/EPszbNddM+mFjLsDZYB5\nAYNxuZHU22aq4NYH+LG4tYv7gFHa9r/1b9ZkqHWhaZtJxNVGCsmvzl243EhMAkaBd43cMpoL\no/7+VHC1kdK8MuyNtDPzDkIHitYS24BxDoz6+1MhxEgZI+kpsnKwYkojLQya8ZTOXSgFjA+/\nAPGs/2Ry7nHxxwmdWzeL2AaMc2DWQaQibm3MK1yA9vw26dLPenuGePXvQiFg6lPI12e9PjWJ\n2AaMc2DU358KIUZKIz4wIjPQCUa67L7yR5Of2DP1Sv27UAoY6fX4dlwCxz2GfzXo34MZxDZg\nFHByywYmg3fvu5n4HOm8kcTZ5HGnYwOtkdrVo3rOR+zcRAvNgEngEjnOQPNpU4ltwCjg5EoY\nMwfvdqQugpE67kZtyxGqJHSSr4WmMP77uQTuFv0bNhUIGDK26PzEkboIRnqgaIH/c/sVBjro\noTESx/1B/4ZNBQKGDOhChtZIzS/wFrrr98cUUlF2O0XeA88DWvmIMYwCJhpdHBowVMShLnSP\njQjdTq3FjRB7aO2By0rirqbaTewwrzsuDV1OD5dHi03uKU32s0nsmKcLqdFqEI7XRThU7nZy\nz/QYQrdT28RGqz2u0thDIl9GskmttwSjgNGvS8sq+V376Yvl6e+Y5CdqzNPlM7HRauFI4hqO\n10U00kO91xC7gscodzs1drDGSHm/9Z+YEkddY4CBw42sdc01l92tOXxfe0ZlgSh0Wbhw1Ehj\nxwe62FAX0Uh7Doy9+J8qYw4rdzv1Z9IQ9yEUJXiInQtpkphsbD1Ptmae+m7XEoaOaHTpm5ji\nUF2022bGny7YSAhtvTbzd4toZNSN0cLjgIXaaUhQDPdORzSNeSlwrC7RtM2kwJG6yEZC6Mel\n1zLYYiSOFIYHAoZMNG0zKXCkLsFGMgtHCsMDAaPGWaZbC8KRuphwoo3AkcIEU/Us4w0KOFeX\nD685iNCg/u+w2l4ojtQlyEhN8xhsj4QjhZE4tHQgl89ygxKO1eXN5PtPIrT9wSRzOix0pC6C\nkQKlgHLOyMh0FDhSGIGjyy/juKmbzemYz7G6XLxM+P7LrxltMBRH6iIYKdu8qinkUGF4Xroi\nodsfPuEIY5YzwbG6JFUJ3z8mMdpgKI7URbBOZf8pDchnVnlp2n+Mrfd4mbH1lmv2c0sJd25p\nm//TLCM5VpeMH4Tv3Qbe06LAkbqI5jk7q/hb04zkWB7p5bn6n3WmGcmxXDdb+J5htyaUFiKZ\np/Tcv4KRwmnbMbNrMvfScavzYTN2dyzZ1Xj2y0ntNceKiB9k8xwemWFhPmxL6/+bltdutNW5\nsBlVd2RxXPrNe63OTQ0ACgAAHlZJREFUh42Aq5A2zR9OtjoLtqPl8CHG42M7HGykt/N8ZTmF\nH1uaF9shjZ33mZW5sB+gSyTYSD1Wtg5c9pKVnVnakK7lga/GOR6LM2IzQJdIsJESfSfTfQ0g\nTAgP5vOjI33e+7wtVufEXoAukWAjnVf6ynBU1s3SvNiOtvmddvnmJv1B5UWtuAR0iQQbaW1y\n0satGauszIodea59719utToTNgR0CUeqtauvRT44w0TwUhKjd2ldBugShonV39G/C+fl6Lch\nppRXYPAqHt+U9+lOW4jdSkWxVdCFvFVH64LX7sp/nGH6SJaFMLVzde5NXoFBwJjznjnoQsbZ\nughr54u6XBDVtsJgIUwUac14OZwJoAsZZ+sirH3Amx94JYnpWBF8zo5PSO8w1TdgGTrKPY5O\nt9upmI1V56TNW5ab9keEzt6blXP/mZYnCpIu+jt/qeY3w/8FJSHv7fnc7Gf4lGfvy85dyXnF\nGfYDdCHjbF2wDQMjgDUxfaWaP5qbrj1S2W/On25Aa7OHoQ/OURy5j7uz8T3u9sYN3DH0+H9V\n7x0w+6muZfWfdA8RRk5C3tu0pve4Gn/K+X33Hb2eFyYwI5oj+HlcZsp1R6LZAgnQhYyzdcFG\n6u930tZicneyBvHnscFTgdD7hV90aJn2aKrvYeWB77jvUWvgz4t6+Y9lc0G3N/1z3woRRk5C\n3ttR1MbLgc79AKEf+KnAjGiOYMzQL3ZeaWDMKHVAFzLO1gUb6YF++0qSH2F9a3eQa0aowtOa\nt7P42wu2DFR+xZ8/APHPwxfWEnlF+cMThNnPhSQh741/B48XxrPXf23lp9oUE1OS+g1Clcxb\ne4AuZJytCzZS20zu8t1GN0ImcIbxZ/LDAnT73E6t02anKXcwHnTUneoROlnR/S3/3PWCMG1o\nG40wSBSmSxlCezmvOCOaI+CqEWpkXjYHXcg4WxdJjrY5IxgPS8nnbNx1Ryv7zUKvZ16P3swc\nppINr/RXMrnm0NV3P1P4Uf2nRbwwx7j3G27WI8yc/pU/j2MTMP7ba/bvDYMuZJytS0AO854L\n1IxPbz+1EdUmLkY1nEpVRtBR19+SlnVHfcvjXRIvLA08YHsiM/dVPcL4StLyV3NHWATMLqFv\nJfYPHkGXSJytS8A6cse8RjdjJ975wX8Hn8zitTMzO/6OPaALGTa6SFJU+220xxWtv2ddVV07\n8Qarc2E/QBcybHTBRnoycIJhXqVpBSdvzUgdw/wph/MBXciw0QUbqdOWbRNb/rQ42s25i+Cm\nh664i2EE6BIJ1sHT2FSMzuRamhfbYfL4SI4FdIkE69D74cbiL3a1tzQvtmNPMFZnxkaALpFg\nI63x7F3tSTBn+BIAcD3SlbmuCZ2otTInAOBgsJEaH+qWVDC70dK8AIBjwUaa3GdDxcbiuyzN\nCwA4FmykDL7F6neZVmYFAJwLNlLuIf/HwY5WZgUAnAs20qrR5c3lo5ZYmhcAcCzYSNLztemW\nZgcAnAk2ktT+O6qX+QEgToEWHgDAgICRoBEiAERH2Buy0AgRAIwQ8A00QgSA6IALEAAwAIwE\nAAyI3kh7uPAxLLx0G12YZmTnlBu3GQGR+BJo6vU/CX3VtD3YlR/TeB73EoovKYJRkmX377KS\nLn7DUbqwMFL4HPFADuQHz6zM9uJaDdTyYKfky8uNHTDFwB2hO7YFAZG4Xceqvx46VoiY+fn8\niMZt5517Bb88jqQIRkGWHzvOKvf+NflDJ+liopFCRDjbn5eJr8yYext6sej7uukDTDtz2PBM\nvUfqNG1HRmBicUd+PGP0RVJp4lFkXpZtKEUwCrLccTe/cMbVTtKFjZG83D87p09qRK/19uQ+\nzeey7YH+J/K5rgitK/L0Wu9PNXuK2Pfegd51aMIC/2pJ0tEIA3OsL6xDL19yFm+KK81bIyzA\n2xC++JWEwT8Qt7azOA5H0wM56bfUIXH/gR0L64pDdkhrlOblTkGopt2+qA9bv0hCxOxK4SdW\ntBdGLJl1dVu35QjFkxTBkGVpyQqIc/Kgk3RhZaTRB78vnl3n+ePpj7lqL9c2s9+JwPHUJ29s\n/EceQqXFDaKRxv478HXqscGSSsLAHOi2koN53+BNIW7I5jphgbgN8YtfSRj8A3ETGzYI43DM\n77//6NAZSNq/tFFxyA5pjSGb1+e3ohWDoz5qAyIhrtzXuP+q0f6Jp5N6+/jZrYX/RH+4DKF4\nkiIYsiyHOekFU+fowspI/jvbdwub9jY2bee7gZ3D7RNEqPWsPNnqQ4fP/VbsDXbrRYExb5Zz\nyZ9JKgkDc6Dagj7PSptC3Ca8QNyG+OWVBv9A3DY8Dkf39/35WI+k/UsbFYfskNbYhFrytqJh\nL0R91AZECpQQE0cc8U902l4UuI//NPkE+oQ7hOJJimDIshzifDiBc3RhZSR+OI5EtLpP8S18\nLodfXiKKsGFY+qCP0MilPh9X3uT/fe0/hbV88y+QVBIG5vCfNJJr5U1xldICYRvilzz4B+L2\n486a+UEMePD+pY2KQ3ZIa/g3OmXm4RTy2FMmIt/D8HBb0KdJn/snpgcqX5aieJIiGLIszWmB\nW7vtg846SBeGV6TSgrqkHfzATl6u4TvxguPbihr/1hFlBAJmnr+ElHnav8rd20PKSMLAHKiq\n88gSaVOBwxcWiNsQv4IG/5C7SO/yMUJfL0bS/qWNikN2BK/xYc+lY6I+aCMiBUWMf+IP/ruY\n5s4v+Xy+qQNRPEkRDFkWNOEO/ueD/ZykCysjjTq0u3h2bcKXZ57kyvlcPtynycudQaczPvKt\n6CrJhP56HT85cfj+49OGSCoJA3O0DltUnVuGNxVILiwQtyF+BQ3+Ias0a1DVT8OnIWn/Z/C6\n4pAdwWs05ZyzLuqDNiJSaMSc7jUHlbU77v+1mauKJymCIcuCKtrf9/2PS5PeRg7ShZWRFudk\nTG5Ei7JznxtRyB9IQ49nG3qnI7SqMLH3JkkmdO0ifvLY+MzkK72SSsLAHEv7NqPV3U/hTfHJ\nhQV4G8JX0OAfskq+qe3Tb65D4v4DOxbWFYfsCF4DTcr2hR+D6URGDNri+d+7RvO/WvKfjScp\ngiHL4r90jEpLuoSPZefowspIRtYkrsW2hp80ZMdMC/pKinzYFkocSRGMliwO0sXlRoocsuPs\nkc5bWO6BDhsYyS5SBGMHIzHShYWROFRjpKeHp4jvPhnalCKRQ3ZsSprSxnIPdOxRf88rnqQI\nRkMWJ+li7yYkAOAQwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIAMACMBAAMACMBAAPASADA\nADASADAAjAQADAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIAMACMBAAMACMBAAPA\nSADAADASADAAjAQADAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIAMACMBAAMACMB\nAAPASADAADASADAAjAQADAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIAMACMBAAM\nACMBAAPASADAADASADAAjAQADAAjAQADwEgAwAAwEgAwAIykAidhdU4AuxMIEQgYMl4Jq3MC\n2J2AdSBgFDkzqyCxouSk1dkA7A6+BkHAkJkxcDN3ctQkq7MB2B1sJAgYMh32+SU6kG11NgC7\ng40EAUMm45hfl+p0q7MB2B1sJAgYMvfcVsPVT5podTYAu4ONBAFD5uz8np5OU+qtzgZgd7CR\nIGAAIArgyZE6W9uszgHgCLCRIGDIpBUtOGR1HgAHgI0EAUPm1Gtj0q5ed8bqbAB2BxsJAkaR\nEyuLcqzOA2B3gspIEDAk2r56tHf7yVbnArA7kpEgYIjM6ZU24W24UANaYCNBwJAZvfqU1VkA\nnAA2EgQMGWjMC1ABrb/Vgca8ABXQ+lsdaMwLUAGtv9WBxrwAFdD6Wx31xrwtzw/ISsgc8GJr\nTPME2BBo/a2OemPe2b945ct9X79aPDemeQJsCLT+jobs3YGvw/AgO+6B1t/RkPlj4Ksuy+J8\nAJYDRlJBs5uyuwd9XNd68tPhUNsZ90C/dipodlPWODPXL1r+XF8scwXYEejXLjraqquOwatc\nANzaRQVUfwMiYKRogOpvQASMFA1Q/Q2IgJE0qFErO0L1NyAiGUk1YAgcLxOYcdlImzFi0GDN\nNNeWUx7m00arv1+NwZHqg6kuQBg4RNQDhsDjuMY8c669eDCH4wrmaCTK+hflYWaVqtXJqVR/\nj+0bk6Ol58EuHNeLmS5AGNg66gGjQsHFzPLChvG8u/+kkaiQNmDyGtSXR1R/vyOe27MHUu4h\nViwqOvRD7msaiah1AcLARtIKGEVsZ6QE3khpGomoA+b1Wz/fr3zLS6j+3iKe2zOLKPcQK6Zd\ndOFF5y3QSARGMgo2knrAqGA7IwXuNxM0ElEHjHqLD5Xqb9vpModLSuJe0UgERjIKDhHDTYRs\nFzBsjaSOSvW37XS5hNflJo1EYCSjRF39bbuAYWukxtldPQVzlNrSqVR/206XFF6XThqJwEhG\nwUZSDxgVbBcwbI00pf8HFR/0vVdhqUr1t8t1AcLARlIPGBVcHjDt9/s/9nVQWKpS/e1yXYAw\nsJHUA0YFlwdM+lH/xxHlviwUW3+7XBcgDGwkrYBRxOUBM350RXPF6AkKS+f+jGpuSs+5H65I\ncQ82knrAqODygDk+zsN5xh1X2tcedOOIHV8Mvz9ykct1AcLARlIPGBVcHzANB5WfVXN7Wtsd\nQKiiY+Qi1+sChCBXf6sFjApxHTDcnpbURoRqkuVZP64Q6PArJntgBxjJVOA5ktqm5mo8qOZu\nXPDrJaj1kaHyrKW9BJIKqXIbO8BIpiJ0fqIVMCq42UjeWqTel8WbC+68dDRalvdl5CI36wJE\nInR+ohUwBF7AZ95uJuXMKLEPmJ8aCTPBSPGF4Vu7clwWuIBldhjANmDezvOV5RR+rD8bYKT4\nAhvJ8QGzfIqIEDD4l0JgUAdMj5WtA5e9pKPi4Cv8itz59OvEBDCSqWAj6Q4YjF2M9Mh4ESFg\n8K/nyKmpAybRdzLd1+Chz8cq8cW+lO7068QEMJKpYCPpDhiMXYwkcR4fMFM1ElEHzHmlrwxH\nZQbKgbbTBYxkKthI7gkYNI5L0HrTnD5g1iYnbdyasUp/Luyiy1COiMIAPmAko2AjqQeMSo+i\ndgmYIMbN0ExCHzD1tch32kAm7KLLV2JvT2WBV/AT8K/d5NRgJKNItXaqAWPXV6ordpIYfjNx\ntjdoRfMDxi5GktjMG0mr8QoYySh01d82faX6SAL5voVMu6DraRxWfyN0wRDNJGAko9BVf9v0\nleoqbj994k+5ZvmHedXfGBsaiektLxAKXfW3TV+pjoGRHFmbWZOt50rdxciVGgiDrvrbpq9U\nx8BIjqzNrOL+h1RI/PBj0txXDOkChEFb/W3LV6pjYCRH1mbGQBcgDEdXf8ciYPTWZp4RT/R5\nF9FnjTVV3Fe11JSCkVjg6OrvGBjpzKyCxIqSkwpLCbWZ86SiB33WWFOlp4jEgZFYgI2kO2Ds\n0DgzBmfeGQM3cydHKY1aTqjNbBJ318XSEwwYKdZgI+kOmFdt0DgzBgHTYZ9fogPZCkttW5sJ\nRoo12EgQMGQyjvl1qVbqpsy2tZnXjqfmCjASC7CRHBowphvpnttquPpJCk08EaE2s0ZszNap\n2MgRseFHrrgvgW49SXPP51rkNcFIRsFG0h0wGGuNdClp/Mbhw0lz+xky0tn5PT2dptQrLCV0\nEPl0jkBigdGjip6W5xeSKB5MnL0yaE0wklGwkdQDRgUrjVSdpOeClGnkCf4aVUli2EHk8UqJ\n3RXy9Bmt9Q6TG/MeI6cGIxkFG0k9YFRgHTB3yteQon7yNHmErHq5Uu6OSyQ65MnTs+QUp4JW\npA6Y3E/Ulsawg8jO5LOD5rAHfcnrjSGnBiMZBRtJPWC8XsVehlgEzOmre0h0bC/hSZWn8+UU\nTxC38YZ8rzLtAXn6ffIeqQNm45BNlcq9KxE6iMSwNlK1fBW65CF5WrNTzzPymaS6Rp4+S04N\nRjIKNpJ6wPTDJ7LIRSwC5n/13KFxBobMiIDR0JeEDiIxLHRp/CuxUNPtGuLst6LfIRjJMHRD\nX7ZNvkdpAywC5jP7GkkdkzuI/F/ukt8QOOcXpLk9O0e/QzCSYSj7tdug20g18u3HD/LkfkJr\nPYT262r2348uy6rQD335ULekgtmkHiBlTOsg0rEnmDhEGvpSPWCaTiktCQmYPV9iPksm/7en\nSym+JFdvGC1U64M6YCb32VCxsfgu/XsAI8UX2EhMAmaDrn98X+L2TskF4qPHNAvHRqFv2cC3\nMvwuU/8ewEjxhdSygUXArFW4DJGxctwT+urvQ/6Pg4TqbS1YGGkbl0YiKZk0NxnKSFYiVX+z\nCBgXGmnV6PLm8lFL9O+BhZGa3lxHouhW4uxPo98hGMkw2Ei6A6ZJrD7IDzZSBqltjgJFFjZG\n01397S/cUa6wBDcRIo+P5Pu9rECvwdLklRu1tnuT3Douo6s8/QJltigBIxklrPqbPmAewysE\nvcC2LonUKrJv3q9IcwtZ1L4ZhTpgpEfRNZQrHMSNVi8kLm6aP1ci63fytObl5EX5gdHts+Tp\nTZTZogSMZBRsJN0BcxpfkYJeqa4mPihcmH0zcfa7rI9FB/SNVgOfSi88qhB8a7fr/F4kkjoT\nZ9+nf2fMACMZRWq0GviMMmAU6Lpa/2ZNhjpgLvnW/1HWXf8eQsqOWStIPLiUNPcmcm1mbAAj\nGQUbiUnAKOBkI73c+dnaklw9B0B41XxtSg49aeR7wtgARjIKNpL+gBFxuZHQwT7JI6t1bPlR\nQtnRhbWZQBhSEyG9AYOhMFLhvw1s11yoA6aupGB+50XN2gkxpLIjGMn9YCPpDhgMhZHeNWJQ\nc6EOmMLbalHVMAPFlpAyUhqxNpNMoSMeCwBhYCMxCRjnQB0w7/AfrX/Rv4cQI+lpdLAYKhuc\nCDYSk4BxDow6iFQBjBRf0HUQqYLLjaTe358KIUbqVEnPI2AkJ0LXQaQKFEaaukv/Zk2GUQeR\nKhhvFa89HJh5gJGMQtdBpAour/5W7+9PhRBdvMRrT/5fiLONdkPDAjCSUeg6iFTB5UbS6u9P\nEZfrAoRB20GkIi4PGDP7+3OyLkAYdB1ERtcdl3sDBnQBROg6PyF0x7UWj0bRg7jG6eHyI8Zk\nuc/pfjaJHUYBA7oAInRGInTHtU18GaLHVcQ1WlbJ7ZmnL5anv4sqt8xgFDCgCyASdXdcY8k9\nswczauQ1hhg43Nh6l92tmaX2jAIGdAEEKI2k3B3Xn7VbjyWmpBsiMdnYep5s7TxtZyMf6AII\nUBopKowWqgcsNLbeuBnG1os1oIuLACOpENLggNE2g3CsLkAkgfiAgCGT5pUBXQA1AvEBAUMm\neP9PMdpmEI7VBYgk4BsIGEVO8B9Nz7LbYBBO1gUIQ7wAtRxqRs2HWtTTGsXJAdPf76StxT3Z\nbTAIJ+sChCEYaVsB1/O9Qq6rOZWfTg6YB/rtK0l+RHNcPEM4WRcgDMFIF/3F9zy33LfoElP2\nMe0/xtZ7vMzYesvXGFuPRNtM7vLd7DYXgpN1AcIQjOSpRz7uDDrhsTg3NqRtzghzrkeAqxCM\n1Gfpmee5v59Z7MjXxs3DzGcCgLsQQuR/u/jLSF25gs8tzo3N8Kq8JgEAwYjn2uaDzajpoIFu\n7VxPDRgJoEAw0p7fJl36WW/PEK+lebEjT8OtHUCDECKX3Vf+aPITe6ZeaXFu7EdWaZvVWQCc\ngGCkdvWonvOhmmSLc2M/8qDKDqBBMFLH3ahtOUKVWRbnxn68fuvn+6GMBGgiGOmBogX+z+1X\nTLA2MzYEqr8BKoQQaX6Bt9Bdvz9mbWYAwKmYeK6N/o0ML0e/DTGlvIIZL4QwAXRxIzY3Uu1c\nnXuTV2AQMNzLUuOGxJlRbw3jfF2ASAKq/jwuM+W6I6w3zSJgokjLIGC8J6XGDVtSo96atFXH\n6wJEElB1zNAvdl7J/BkS/x87PiG9w1TfgGXoKPc4Ot1up2I2Vp2TNm9ZbtofETp7b1bO/Wda\nnihIuujv/C0Mvxn+LygJeW/P52Y/w6c8e1927krOK86IksaHuiUVzG5Ebe9HuyUJV+gChBEw\nUuo3CFUyb/nN/5dvuvZIZb85f7oBrc0ehj44p1UxG3c2vsfd3riBO4Ye/6/qvQNmP9W1rP6T\n7iEBIych721a03tcjT/l/L77jl7PB0xgRpQHMbnPhoqNxXdFuZUQXKELEIbQ+Um1/9TL/Irv\n/981eCoQer/wiw4t0x5N9T2sPPwS9z1qDfx5US///3hzQbc3/XPfCgkYOQl5b0dRGx8m6NwP\nEPqBnwrMiPIgMviXkb7LjHIrIbhCFyAMwUj+k5nPDCMd5JoRqvC05u0s/vaCLQPfUM6GV/rz\nBMr2fKTx/3YhYPZzIUnIe2vjl/lTe/Yi1MRPtSkmpif3kP/jYMcotxKCK3QBwhCMtMvrLedY\nt3IOnHn9/7wPC9Dtczu1TpudVqecDa/016keoZMV3d/yz10vBEwb2kYTMEgMmC5lCO3lvOKM\nKA9i1ejy5vJRS6LcSgiu0AUII6xfO5ab5v9j4647WtlvFno983r0ZuYwlWx4pb+SyTWHrr77\nmcKP6j8t4gPmGPd+w816AmZO/8qfx7EKGEmZ6VFuSMYVugBhmPwcqWZ8evupjag2cTGq4VQ6\ntQqKhvpb0rLuqG95vEvihaWBB49PZOa+qidgfCVp+au5I2wCRnq1j13p3BW6AGG48KHCOz8g\nVJkMLymGA7qYScBIGeEzHM2sq6prJ97AYkugC0BLLPr+jjEnb81IHcOknQboAtASiI89wVid\nIxsBugC0OP9ECwA2AIwEAAwAIwEAA8BIAMAAMBIAMACMBAAMACMBAAPASADAADASADAAjAQA\nDAAjAQADwEgAwAAwEgAwAIwEAAwAIwEAA8BIAMAAMBIAMIC5kfZwGcI72qnX/xS2aGGakb25\no9P3gCwCESOsxLMursEEI/Fb3XWs+uuhYxFqebBT8uXlaNvF7S753OD/nmIEkwP5BrYbW/bI\nh37yZTyF9YhjXVyDSUby+j92+M/BLxZ9Xzd9QEvBkhNLuraYdhJ1wMl5DymLURmJAgfo4hpM\nNNKuFIQmLPDPSTrInUE+7pD0jxVGKFlfWIdevuSsl/tn5/RJjYgrzVsjLEDrijy91iPxSx4F\nBXFrO4sDkjQ9kJN+Sx1Cr/X25D6N8rmueKPi2CXSGqV5uVMQqmm3j/Vx6mRPoMv7wJEGer4P\nHIl/qu2B/ifiWRfXYJKRyn2N+68aHZhx6rHBzT3/fub5Hs1SwAgjlKDbSg7mfeOfOfrg98Wz\nETdkc52woD55Y+M/8pD4JY+CgriJDRuEAUnm999/dOgMVOf54+mPuerAdoV1xbFLpDWGbF6f\n34pWDGZ9mHoRjBQ40oCRAkfi5dpm9juB4lkX12CSkfiBE0YE+lBbziV/hkr9v9+TA0YYoQTV\nFvR5lp9ZjtC7hYjbhBfUelaebPUh8UseBQVx2/CAJN3f9+9pPWra29i0XezOWlhXHLtEWmMT\nasnbioa9wPow9SIYKXCkASNtE8ZamcPxl4Q41sU1mHhrJ+Kbf0FT1+W+vxXKVyRhhBL/+TO5\nlg8YfoiTRMRVSgs2DEsf9BESv+RRUBC3H2+aH82BZ3Wf4lvEgBGHPRHGLpHW8G90yszDKZYP\n1y4YKXCkASPtFzrkHn55CULxrItrMNlId28nlZGEEUpQVeeRJeKZt7QgsI6wwLcVNf6tIxK/\ngkZBkfuK7/IxQl8vRnVJO/gRtgLbFdYVxy4JXuPDnkvHsD5K3UhXpNICr9TzvZdr+M5/vY5n\nXVyDyUaaOHz/8WlDWrouqXuuUK61E0YoaR22qDq3zD9z1KHdfFnAixeczvjIt6IrEr+CRkGR\nA2bWoKqfhk9DtQlfnnmSK/f6jSquK45dErxGU84561gfpW4EIwWONNhICD3cpymedXENJhvp\n2PjM5Cu96POLPBdvl8+8wgglS/s2o9XdT3m5xTkZkxsD6wgL0KrCxN7+koHwFTQKihwwvqnt\n02+uQ2hRdu5zIwobeqfjdcWxS4LXQJOyfayPUjeCkQJHGmqkhh7PxrMursEcIylAfK7B9mEH\naeySmUzHUjaGYKTAZOQBx7EursFtRoocu+Tskc5bWO7BGEFG+iQxfGEc6+IaTDCS4iafIi6q\nYTeoJCKNXbIpaUobyz0Yg5dFONL3U2eGLYtnXVwDNCIBAAaAkQCAAWAkAGAAGAkAGABGAgAG\ngJEAgAFgJABgABgJABgARgIABoCRAIABYCQAYAAYCQAYAEYCAAaAkQCAAWAkAGAAGAkAGABG\nAgAGgJEAgAH/H8V+QBbcSozFAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "par(mfrow=c(4,3))\n",
    "invisible(sapply(1:11,\n",
    "       function(col_num) {\n",
    "           protein_name <- names(data_matrix_discrete)[col_num]\n",
    "           plot(data_matrix_discrete[,col_num],\n",
    "                data_matrix[,col_num], log=\"y\", xlab=sprintf(\"[%s] expression category\", protein_name),\n",
    "                ylab=sprintf(\"[%s] expression level\", protein_name))}))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "add a column indicating the experiment name for each row to the discretized data matrix; then save the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_matrix_final <- data.frame(data_matrix_discrete, expt=data_matrix_orig[,12])\n",
    "\n",
    "write.table(data_matrix_final,\n",
    "            file=\"sachs_data_discretized.txt\",\n",
    "            sep=\"\\t\",\n",
    "            row.names=FALSE,\n",
    "            col.names=TRUE,\n",
    "            quote=FALSE)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}