tabalho_final_apresentacao.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>SER-347 - Introdução à Programação para Sensoriamento Remoto </h1>\n",
    "<br>\n",
    "\n",
    "\n",
    "<h1>Classificação simultânea de diferentes imagens - trabalho final </h1>\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "### Maria Antônia Falcão de Oliveira"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O presente trabalho teve por objetivo elaborar um programa capaz de realizar cálculos a partir da discriminação de classes por meio do algoritmo de Máxima Verossimilhança e classificar diferentes imagens simultaneamente."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"_trabalho_final_intro/figuras_apresentacao/diagrama.jpg\", width=450, height=450>\n",
    "\n",
    "**Figura 1** - Diagrama da proposta de trabalho "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importando as bibliotecas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################### IMPORTANDO AS BIBLIOTECAS ######################################\n",
    "\n",
    "from osgeo import gdal, ogr # leitura de imagem e shp\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "import geopandas as gpd \n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# importar constantes\n",
    "from gdalconst import *\n",
    "\n",
    "# informar o uso de exceções\n",
    "gdal.UseExceptions()\n",
    "\n",
    "# biblioteca de funções relacionadas ao sistema\n",
    "# sys: System-specific parameters and functions\n",
    "import sys\n",
    "\n",
    "import glob\n",
    "\n",
    "from numpy import *\n",
    "from numpy.linalg import inv\n",
    "\n",
    "import os\n",
    "\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- O MaxVer é um classificador estatístico baseado em Bayes.\n",
    "- Por ser um classificador supervisionado, requer um conhecimento prévio das classes de interesse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Acessando os dados no formato shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################### ACESSANDO OS DADOS EM FORMATO SHP ######################################\n",
    "\n",
    "# informando diretório onde estão armazenadas as amostras\n",
    "\n",
    "dir_shp_tr = 'C:/Users/Maria/_trabalho_final_intro/amostras_shp/samples_tr.shp' # amostras de treinamento\n",
    "dir_shp_tst = 'C:/Users/Maria/_trabalho_final_intro/amostras_shp/samples_tst.shp' # amostras de teste\n",
    "\n",
    "#acessando as amostras treinamento e teste\n",
    "shp_tr = ogr.Open(dir_shp_tr)\n",
    "shp_tst = ogr.Open(dir_shp_tst)\n",
    "\n",
    "layer_tr = shp_tr.GetLayer()\n",
    "layer_tst = shp_tst.GetLayer()\n",
    "\n",
    "# identificando o número total de classes \n",
    "tam_tr = layer_tr.GetFeatureCount()\n",
    "tam_tst = layer_tst.GetFeatureCount()\n",
    "\n",
    "# checando uma determinada posição de um objeto \n",
    "atr_tr = layer_tr.GetFeature(5)\n",
    "atr_tst = layer_tst.GetFeature(5)\n",
    "\n",
    "# verificando as coordendas do objeto selecionado \n",
    "nx = atr_tr.GetField(\"X\")\n",
    "ny = atr_tr.GetField(\"Y\")\n",
    "\n",
    "# identificando o ID\n",
    "fid = atr_tr.GetFID()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verificando a estrutura da tabela de atributos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>classe</th>\n",
       "      <th>num_classe</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CAF</td>\n",
       "      <td>1</td>\n",
       "      <td>324022.500002</td>\n",
       "      <td>9638402.5</td>\n",
       "      <td>POINT (324022.500001629 9638402.500000015)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>CAF</td>\n",
       "      <td>1</td>\n",
       "      <td>323527.500002</td>\n",
       "      <td>9638192.5</td>\n",
       "      <td>POINT (323527.5000016287 9638192.500000015)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CAF</td>\n",
       "      <td>1</td>\n",
       "      <td>327247.500002</td>\n",
       "      <td>9635897.5</td>\n",
       "      <td>POINT (327247.5000016286 9635897.500000015)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CAF</td>\n",
       "      <td>1</td>\n",
       "      <td>326272.500002</td>\n",
       "      <td>9635807.5</td>\n",
       "      <td>POINT (326272.500001629 9635807.500000015)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CAF</td>\n",
       "      <td>1</td>\n",
       "      <td>324577.500002</td>\n",
       "      <td>9636182.5</td>\n",
       "      <td>POINT (324577.5000016291 9636182.500000015)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  classe  num_classe              X          Y  \\\n",
       "0    CAF           1  324022.500002  9638402.5   \n",
       "1    CAF           1  323527.500002  9638192.5   \n",
       "2    CAF           1  327247.500002  9635897.5   \n",
       "3    CAF           1  326272.500002  9635807.5   \n",
       "4    CAF           1  324577.500002  9636182.5   \n",
       "\n",
       "                                      geometry  \n",
       "0   POINT (324022.500001629 9638402.500000015)  \n",
       "1  POINT (323527.5000016287 9638192.500000015)  \n",
       "2  POINT (327247.5000016286 9635897.500000015)  \n",
       "3   POINT (326272.500001629 9635807.500000015)  \n",
       "4  POINT (324577.5000016291 9636182.500000015)  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir_shp = gpd.read_file(\"C:/Users/Maria/_trabalho_final_intro/amostras_shp/samples_tr.shp\")\n",
    "dir_shp.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verificando o número total de amostras "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "O número total de pontos de amostra é: 500\n",
      "Número de classes: 5\n"
     ]
    }
   ],
   "source": [
    "num_classe = (len(dir_shp))\n",
    "print(\"O número total de pontos de amostra é:\", num_classe)\n",
    "\n",
    "# verificando o número de amostras por classe\n",
    "for i in range(tam_tr):\n",
    "    aux = layer_tr.GetFeature(i)\n",
    "    auxclass = aux.GetFieldAsInteger64(\"num_classe\")\n",
    "print(\"Número de classes:\", auxclass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verificando os tipos de dados das colunas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "classe         object\n",
       "num_classe      int64\n",
       "X             float64\n",
       "Y             float64\n",
       "geometry       object\n",
       "dtype: object"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir_shp.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transformando o dado shapefile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################### TRATANDO OS DADOS EM FORMATO SHP ######################################\n",
    "\n",
    "# parâmetros para construir uma matriz para as amostras\n",
    "# número de classes\n",
    "nclass = 1\n",
    "for i in range(tam_tr):\n",
    "    aux = layer_tr.GetFeature(i)\n",
    "    auxclass = aux.GetFieldAsInteger64(\"num_classe\")\n",
    "    if auxclass > nclass:\n",
    "        nclass = auxclass\n",
    "    \n",
    "# número de amostras total\n",
    "nlin_tr = layer_tr.GetFeatureCount()\n",
    "nlin_tst = layer_tst.GetFeatureCount()\n",
    "\n",
    "# construindo a matriz a partir do numpy para amostras de treinamento\n",
    "matriz_shape_tr = np.zeros((nlin_tr, 4))\n",
    "\n",
    "for i in range(nlin_tr):\n",
    "    aux = layer_tr.GetFeature(i)\n",
    "    matriz_shape_tr[i,0]= aux.GetFieldAsInteger(\"FID\")\n",
    "    matriz_shape_tr[i,1]= aux.GetFieldAsInteger64(\"num_classe\")\n",
    "    matriz_shape_tr[i,2]= aux.GetField(\"X\")\n",
    "    matriz_shape_tr[i,3]= aux.GetField(\"Y\")\n",
    "\n",
    "# construindo a matriz a partir do numpy para amostras de teste\n",
    "matriz_shape_tst = np.zeros((nlin_tst, 4))\n",
    "\n",
    "for i in range(nlin_tst):\n",
    "    aux = layer_tst.GetFeature(i)  \n",
    "    matriz_shape_tst[i,0]= aux.GetFieldAsInteger(\"FID\")\n",
    "    matriz_shape_tst[i,1]= aux.GetFieldAsInteger64(\"num_classe\")\n",
    "    matriz_shape_tst[i,2]= aux.GetField(\"X\")\n",
    "    matriz_shape_tst[i,3]= aux.GetField(\"Y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Acessando as imagens utilizadas na classificação"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Arquivo aberto com sucesso!\n"
     ]
    }
   ],
   "source": [
    "###################################### ACESSANDO A IMAGEM ######################################\n",
    "\n",
    "for img_file in glob.glob('C:/Users/Maria/_trabalho_final_intro/imagens/*.tif'):\n",
    "    \n",
    "    # lendo a imagem\n",
    "    img = gdal.Open(img_file, GA_ReadOnly)\n",
    "\n",
    "    # número de linhas e colunas\n",
    "    linhas = img.RasterYSize\n",
    "    colunas = img.RasterXSize\n",
    "\n",
    "    # quantidade de bandas\n",
    "    nband = img.RasterCount\n",
    "\n",
    "    # acessando as bandas\n",
    "    band = img.GetRasterBand(1)\n",
    "    img_asarray = band.ReadAsArray()\n",
    "\n",
    "    # checando abertura do arquivo\n",
    "    try:\n",
    "        img = gdal.Open(img_file, GA_ReadOnly)\n",
    "        print(\"Arquivo aberto com sucesso!\")\n",
    "    except:\n",
    "        print(\"Erro na abertura do arquivo!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checando as imagens no diretório"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CBERS_4_MUX_20160703_3b_567e8_20m.tif\n",
      "LS_2016_5b_B2457e8_15m.tif\n",
      "LS_2016_6b_B14567e8_30m.tif\n"
     ]
    }
   ],
   "source": [
    "# Acessando o diretório contendo as imagens que serão classificadas\n",
    "\n",
    "os.chdir(os.path.join(os.getcwd(), 'C:/Users/Maria/_trabalho_final_intro/imagens'))\n",
    "os.getcwd()\n",
    "\n",
    "imagens = os.listdir(os.getcwd())\n",
    "\n",
    "for dados in imagens:\n",
    "    if dados[-3:] == 'tif':\n",
    "        print(dados)  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checando a extensão espacial das imagens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################# FUNÇÃO PARA CHECAR A EXTENSÃO ESPACIAL ################################## \n",
    "\n",
    "# transformando coordenadas geográficas nas coordenadas da imagem\n",
    "trans = img.GetGeoTransform()\n",
    "\n",
    "def world2Pixel(trans, nx, ny):\n",
    "    ulX = trans[0]\n",
    "    ulY = trans[3]\n",
    "    xDist = trans[1]\n",
    "    yDist = trans[5]\n",
    "    rtnX = trans[2]\n",
    "    rtnY = trans[4]\n",
    "    column = int((nx - ulX) / xDist)\n",
    "    row = int((ulY - ny) / xDist)\n",
    "\n",
    "    return (row, column)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Armazenando o número de pixel por classe e criando a matriz com valores espectrais"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################### ARMAZENANDO O NÚMERO DE PIXEL POR CLASSE ###################################\n",
    "# associando o pixel da imagem a classe - utilizado para o caso das amostras variarem de tamanho \n",
    "# funciona tanto para tamanhos != quanto para tamanhos =s de amostras.\n",
    "\n",
    "# criação de vetor que guarda o numero de pixel por classe\n",
    "pixel_classe = np.zeros((1, nclass))\n",
    "cont1 = 0\n",
    "cont2 = 0\n",
    "\n",
    "for i in range(nclass):\n",
    "    while cont2 < nlin_tr and i + 1 == matriz_shape_tr[cont2, 1]:\n",
    "        cont1 += 1\n",
    "        cont2 += 1\n",
    "    pixel_classe[0, i] = cont1\n",
    "    cont1 = 0\n",
    "\n",
    "################################## CÁLCULO DO NÚMERO TOTAL DE PIXEL POR CLASSE ##################################\n",
    "\n",
    "# calculando do valor máximo de pixel armazenado por classe (corresponde a tamanho total de cada classe)\n",
    "# se utilizado números de amostras != para cada classe, é necessário saber qual classe possui tam >.\n",
    "max_pix = int(pixel_classe.max())\n",
    "\n",
    "\n",
    "################################## CRIAÇÃO DA MATRIZ COM VALORES ESPECTRAIS #################################\n",
    "# considerando a localização das coord X e Y do Ponto com as da imagem\n",
    "# matriz construída a partir do maior valor de amostras, para deixar as demais amostras de mesmo tamanho, \n",
    "# a matriz foi preenchida por zero.\n",
    "\n",
    "# criando matriz com valores espectrais da imagem\n",
    "matriz_amostras = np.zeros((nband, nclass, max_pix))\n",
    "\n",
    "cont2 = 0  # percorre as nlinhas da matriz_shape_tr\n",
    "\n",
    "for i in range(img.RasterCount):\n",
    "    band = img.GetRasterBand(i + 1)\n",
    "    img_asarray = band.ReadAsArray()\n",
    "    for j in range(nclass):\n",
    "        lim_pxcls = int(pixel_classe[0, j])\n",
    "        for k in range(lim_pxcls):\n",
    "            # utilizando a função que percorre as coord da imagem e associando a localização das amostras\n",
    "            pix_xy = world2Pixel(trans, matriz_shape_tr[cont2, 2], matriz_shape_tr[cont2, 3])  #chama as coord. xy da imagem\n",
    "            valor_espec = img_asarray[pix_xy[0], pix_xy[1]]  #coordenada pix_xy[0] = x e pix_xy[1] = y\n",
    "            matriz_amostras[i][j][k] = valor_espec\n",
    "            if cont2 < nlin_tr - 1:\n",
    "                cont2 += 1\n",
    "    cont2 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"_trabalho_final_intro/figuras_apresentacao/cubo.jpg\", width=450, height=450>\n",
    "\n",
    "**Figura 2** - Representação da matriz tridimensional que armazena os valores de pixels para cada banda por classe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cálculo da Máxima Verossimilhança\n",
    "<br>\n",
    "$g_{i}(x) - ln \\: p(\\omega_{i} ) - \\frac{1}{2} \\: ln \\left | \\sum_{i} \\right  | - \\frac{1}{2}(x - m_{i})^{T}\\sum_{i}^{-1}(x - m_{i})$\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "Onde:\n",
    "\n",
    "$i =$ classe\n",
    "\n",
    "$x =$ dados $n$-dimensionais (onde $n$ é o número de bandas)\n",
    "\n",
    "$p(\\omega_{i}) =$ probabilidade que a classe $\\omega_{i}$ ocorra na imagem e seja assumida a mesma para todas as classes\n",
    "\n",
    "$\\left | \\sum_{i} \\right  | =$ determinante da matriz de covariância dos dados da classe $\\omega_{i}$\n",
    "\n",
    "$\\sum_{i}^{-1} =$ sua matriz inversa\n",
    "\n",
    "$m_{i} =$ vetor médio \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculando a matriz de média e a matriz de covariância\n",
    "#### Parâmetros utilizados para o cálculo do MaxVer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################## MATRIZ DE MÉDIA #############################################\n",
    "\n",
    "# calculo da matriz de média por classe - matriz bidimensional\n",
    "# as colunas possuem os valores da média de cada classe com cada banda\n",
    "# cada classe vai ter uma média por banda \n",
    "matriz_media = np.zeros((nband, nclass))\n",
    "for i in range(nband):\n",
    "    for j in range(nclass):\n",
    "        # soma de todos os pixels / total de pixels\n",
    "        soma = sum(matriz_amostras[i][j][0:max_pix])\n",
    "        matriz_media[i][j] = soma / pixel_classe[0, j]\n",
    "        \n",
    "\n",
    "###################################### CÁLCULO DA MATRIZ DE COVARIÂNCIA ######################################\n",
    "# cada classe vai possuir uma matriz de covariância\n",
    "# calcula o valor banda x banda\n",
    "# cálculo da matriz de covariancia por classe\n",
    "matriz_cov = np.zeros((nclass, nband, nband))\n",
    "# somatório de todos os pixels - média \n",
    "for j in range(nclass):\n",
    "    lim_pxcls = int(pixel_classe[0, j])\n",
    "    for i in range(nband):\n",
    "        for k in range(nband):\n",
    "            matriz_cov[j][k][i] = sum((matriz_amostras[i][j][0:lim_pxcls] - matriz_media[i][j]) * (\n",
    "                        matriz_amostras[k][j][0:lim_pxcls] - matriz_media[k][j])) / (lim_pxcls - 1)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cálculo da Matriz de Covariância**\n",
    "\n",
    "$covariance = \\frac{\\sum(x_i - x_{\\mu})(y_i - y_\\mu)}{n - 1}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Armazenando os valores espectrais por banda"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################## ARMAZENANDO OS VALORES ESPECTRAIS POR BANDA ##################################\n",
    "# criando matrizes para cada banda para guardar os valores espctrais da imagem por banda\n",
    "# lista de matrizes, cada banda/matriz vai ter um plano !=\n",
    "\n",
    "matriz_img = np.zeros((nband, linhas, colunas))\n",
    "\n",
    "for i in range(nband):\n",
    "    band = img.GetRasterBand(i + 1)\n",
    "    img_asarray = band.ReadAsArray()\n",
    "    for j in range(linhas):\n",
    "        for k in range(colunas):\n",
    "            matriz_img[i][j][k] = img_asarray[j][k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementando o algoritmo de Máxima Verossimilhança"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################### CÁLCULO DA MÁXIMA VEROSSIMILHANÇA #####################################\n",
    "\n",
    "from numpy import *\n",
    "from numpy.linalg import inv\n",
    "\n",
    "import math\n",
    "\n",
    "# calculo do maxver\n",
    "matriz_mxv = np.zeros((nclass, linhas, colunas))\n",
    "aux1 = np.zeros((nband, 1)) # matriz onde estão arazenados os valores espectrais da imagem por banda\n",
    "aux2 = np.zeros((nband, 1)) # matriz da média\n",
    "band1 = 0\n",
    "\n",
    "for i in range(nclass):\n",
    "    for j in range(linhas):\n",
    "        for k in range(colunas):\n",
    "            while (band1 < nband):\n",
    "                aux1[band1][0] = matriz_img[band1][j][k]\n",
    "                aux2[band1][0] = matriz_media[band1][i]\n",
    "                band1 += 1\n",
    "\n",
    "            band1 = 0\n",
    "\n",
    "            a = (log(1 / nclass) - 0.5 * log(np.linalg.det(matriz_cov[i][:][:]))) #\"np.linalg.det\" calcula o determinante de uma matriz            \n",
    "            b = (-0.5 * transpose(aux1 - aux2)) #\"transpose\" permuta as dimensões de uma matriz\n",
    "            c = (inv(matriz_cov[i][:][:])) # \"inv\" matriz inversa\n",
    "            d = (aux1 - aux2)\n",
    "            e = dot(b, c) # função \"dot\" realiza a multiplicação de matrizes\n",
    "            f = dot(e, d) \n",
    "            matriz_mxv[i][j][k] = a + f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Criando a matriz de classificação\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################### CRIANDO A MATRIZ DE CLASSIFICAÇÃO #####################################\n",
    "\n",
    "# criando matriz de classificação\n",
    "map_class = np.zeros((linhas, colunas))\n",
    "aux = np.zeros((nclass,1))\n",
    "class1 = 0\n",
    "max1 = 0\n",
    "cont = 0\n",
    "\n",
    "for j in range(linhas):\n",
    "    for k in range(colunas):\n",
    "        while(class1 < nclass):\n",
    "            aux[class1][0] = matriz_mxv[class1][j][k]\n",
    "            class1 += 1         \n",
    "        max1 = np.max(aux)           \n",
    "        while (aux[cont][0]!= max1):\n",
    "            cont += 1           \n",
    "        map_class[j][k] = cont + 1           \n",
    "        class1 = 0\n",
    "        cont = 0\n",
    "        max1 = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizando a Classificação"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAEICAYAAADBWUaVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvX+UZddVHvgdK0KipMaSjezoya046eexUWPUBA39eg31xrGTafDKBHpmoIFJMD/WqHtNw6DCOG0gK4gMjC0mQ+EsN3Q7Y4Z4IHZ7jCs4DkkDIU41a7paCE8H07KcqQYhiefBNpbsttqRhXXmj/f2fd/dd+9zzr3vve5Sqb61alXV/XHuufees+/+8e19QowRO9jBDnbwQsWLrnUHdrCDHezgWmJHCO5gBzt4QWNHCO5gBzt4QWNHCO5gBzt4QWNHCO5gBzt4QWNHCO5gBzt4QWNHCO6gGCGE+0MIv7LA9i+GEF43+TuEEP6PEMKTIYQHQwjLIYRPLOrak2v+wxDCYyGE14YQ/t0ir7WDrYO/dK07sIOthxDCdwP4EQCvAXAZwAUAP7Po68YY99K/3wTgbwF4RYzx6cm2Vy+4C18H4PUAfhHA2QVfawdbBDtCcAc1hBB+BMBbARwFcAbAlwB8M4BvBfB04tR5468AeJQE4MIRY/zvJn/+rat1zR1ce+yYwzuoEEJ4MYB/BOBYjPGDMcanY4zPxhj/ZYzxLcbx/1cI4f8LIXwuhLAeQthL+94YQng4hHA5hPCnIYQfnWz/6hDCh0MIT4UQPhtCOBtCeNFk36MhhL8ZQvgBAP87gAMhhC+EEH4qhPC6EMIT1P7uEMIHQwifDiH8eQjhnZPte0IIvzPZ9pkQwq+GEG6Zw3lfE0L4yKTfF0MIf2fez38H1wY7QnAHjAMAbgSwVnj8vwbwKgAvA/BRAL9K+94N4EiMcReArwXwO5PtbwbwBIDbALwcwI8DqOVuxhjfjbEmei7GeHOM8Sd5fwjhOgAfBvAnAF4J4A4A75PdAN4GoAfgawDsBnD/jOddD+BfAvjNyb3+EIBfDSEs2jzfwVXAjhDcAeOlAD4TY/yLkoNjjL8UY7wcY3wGY4Fx90SbBIBnAdwVQviqGOOTMcaP0vbbAfyViZZ5NrZPYP9GjIXVWyba6n+KMf7upE+bMcbfijE+E2P8NICfA/BfznjeAMDNAN4eY/xSjPF3MBam39Wy3zvYgtgRgjtg/DmArw4hZH3FIYTrQghvDyFcCiF8HsCjk11fPfn93wJ4I4A/CSH8+xDCgcn2/xXAJoDfDCH8UQjhrR36uRvAn1jCOoTwshDC+yYm+OcB/Ar1qet5PQCPxxifo1P+BGNNcgfPc+wIwR0wzgH4TwC+reDY78Y4WPI3AbwYY/MSGJuViDH+XozxWzE2H/8FgPdPtl+OMb45xvjXAPzXAH4khPCGlv18HMCdjrB+G8bm9dfFGL8KwN+VPs1w3gjAbvFdTnAngD9t2e8dbEHsCMEdVIgxfg7APwRwIoTwbSGEpRDC9SGEbwkh/Kw6fBeAZzDWHpcA/C+yI4TwFSGE/z6E8OIY47MAPg/gy5N9fzuE0A8hBNr+5ZZdfRDAJwG8PYRwUwjhxhDCf0H9+gKAp0IIdwB4yxzOO49xZPzvT57H6zAW4O/DDp732BGCO6ghxvhzGHME/wGAT2OsPf0gxtoc4z0Ym4R/CuBhABtq/98D8OjEtDyKsWYFjAMpv42xwDkH4BdijB9p2ccvYyyE+gAewzjQcniy+6cA/HUAnwPwrwB80Dnv8xhzIEvO+xKAvwPgWwB8BsAvAPieGOMjbfq9g62JsFNUdQcvRIQQ7gTw0zHG77nWfdnBtcWOJriDFxxCCDdjrNHtv9Z92cG1x8KEYAjhm0MInwghbHaMAO5gB4vC92MsBH/7WndkB9ceCzGHJ6TU/4hx+tETAH4PwHfFGB+e+8V2sIMd7GAGLEoT/EYAmzHGP5o4ld+HMZ1iBzvYwQ62FBZVQOEOjKOKgieQ8L989Uuui1968sV47pab8Jo7P42PPXkbbnjczpt/7pabAADP7gJee+un59bh//gHS+6+5265Cc/uQtWnZ3aP+/DaW8d99XDD4083jpU+f+zJ23D95frxL3pqerxuR/dH8Jo768/gkcfG/Xl2F7Drpi/W9l1++isb23if9Jf7DdTb+dIjz9WuL32Q62rofc/umrbJ/fH69qVHpvxk6ZuG9FXa1H0uaVPft77OM7tvarT5pUeeM7fre/nSI8/hK15Trm/o/vH9yD7p71e85kXVPQv4fnLXlXuwcMPjT1djPwd+j1Y70g95Npef/sq5zN9HHrsNTz/1xGdijP5EzGBRQjAY22p2dwjhXgD3AsB1t96K3T//E7h0+CTGpH5gz+n70F8Zsy42VwfVeb31iLMnTmHP6aN48PDJxkWWjx3B2ROniju65/RR9Fc28NKfH6C/soEzowu17Zur4+2b9w+m/XnzuD8PHj6JPaePum33VzaAJ4Arh/ZjNAzoYUyqA8YpCHI/ADAahvF13jywmqo9i0uT+9b3unzsCO7uA0tr57F5/6A6r7dRMIoVRoPL2FwdYHjgYvX/lUP7cfbfncLysSMAUF17+dhP4+6+39bd/fH9Aajaa9OPCk+gMbJ4bABjsl/uOtLm5uoAvfWI0ZunjernJcf2PtB8hpsP3IX+8XIPz/q5vcl+8X5pm/vK544Gl4EAbB4dYBem4wPA+BlNnlXvn9X7zddYP7cXGr31OB6rk+ciYx+Ta8jz5r+B8fPmMWO1zc92/dxe3EDXajNnNa67HX/S+WQszhx+AiLNxngFxqz7CjHGd8UY74kx3rNv95VqYgsuHT6JM6ML1aSXH5lMlw6frCYjYzQMONjbh+VjR5ICiq8DjF/QlUP7zTavHNpfH2QYD5aDvX3Z9qVPGv2VjarNpbXz6K3HxoTmYzdXB7hyaH8lNDV0v4cHLjYGKjAefKPB5eqHB6ts4+uOBpex+cBduHJof+06PGi9ATwahvFEVgJQrmlNFAtyD5urA/MZDQ9cbPzw/XrtAeNnz+9Ct9/b2OV+RNoIQOln6f7RMFRCRQtA7if3Xe/3npWF3npEbz1iae18dX2rzep4eiaV4J7s997r5uoA6+f2Vvtl3FvX8mDNz1mxqMDIX8I4MPIGjMm0vwfgu2OM5hu45+4b44Nndlu7sOf00YaA5O16UmoNDqgPGNH0lo8dqV64HCsTXdo62NtXncsvy5swbcHtsEYIoHHPLNBZE/SgB5Z8nRvCnAazaAmi9QHjZyEaZ8kAXFo73xAa0q4M/rbaIDDVTOU9La2dr55/V+3SEi6l56e065zGVwqrHbm2vEse50DzntbP7UVvPZoCWwur1Li2ND+rjRQ8DRFojvcU9Jy/7vbN348x3lPcgMJCNMFJgvoPYlyU8+MA3u8JwBzk4ew5fbT2IwJwae18bcKLBglMv4aiRcr2g719NQEov0fDUNNq+KV72lxqm/eVBsZCVQQNH3Pp8En01qMrcFL7SmFpCSIYWAAC0wHnXVOev/xmbUDaTA3+FPic3sYuDA9cHE/64w9X/3vHe22tn9tbaXf9lQ33nFRbOffCLAKQr6uFmSV8ZezwWNJ9Hw1DTQPj9sViAHztUfbpfnlCtOTeSo4XLEL7Y2yJjJGUJghMtT5t3vJDFAGnwX4z1v7ETwekv0KsDQrYbADaaYQs1HU7Vw7tr2lclW9GCWDtCrC0tFKhbUH6AaChBS6tna80Me07Y3+PJQhKNCiG1UYbLUsfm/NTtW2/tF+ijY2GoXXblg/P0gBFCPJ4snyblvbsPRc+NqXF8z7dR6td7kcbDdDCntNH8eh9P7r1NMF5w3pQpRNawNofMB4soiGmoIMygC1gSiGarNW+9C/lJC7VBHvrsaaJtn1eDBaAojFrWL5H+S1/pz4WPJG0X4+Rm7zesdw2Q0zLWUz1HIYHLrb2H/K5QJlvk8d3rj2GHhuicVvPj8HPTGug3Ddtrp8ZXSiae9Y4l20yjzwfeRtsGU3wz7//Z2qaEWt2oo1YgQiOlmqIf1BDtMCSrxD7GC2/oAX2Qep+Wj5Ovi+thQGo+b/0S0/56ljr5f+t/mrtgs/RAzqnzVkfC44yA7ZJaWlTixRO8/LdLaIfVnTe83MzdPRWfzS8cZszd0ueE/t+22h8nt8/dTxHlGf1CW4JIXjzrbvj3a//4WrCyQ1qBygLC0/4LR87UtFNGHois1DJhedZEMo5gC1UtLbD53n91W2xT07MTwuaHqP38fNL7efnxRNNtCSgXAiJEJT7kTb0ZG7TZs6MzvWHJ2bXYEiu/dL9pYJXP7MSWGMPsM1ifZ5n8pb0VwtXDvSVUF/Yx9+FKrMthOBXhZfEl/78T9Q4SsB0opZobVqYsM8vBU8702DfIGuGOsCSuk5pvwXa7yiRUbmvlCmQi+ZqP58WTLnJm/OBCWRyecERrXW0ESgl6CKE5q19prTfHPiZlbg0tPYoYH8tYxYBzR8WAQszcfuUzN2U8MvNz20lBJmmIkgFL3iSa5ORhYpoUyWBEA9sWmvz0aIpeALLetk5c1baFIFtPadScGBDwBof0Jz83iS2hEVK09POeQuzUmlK0EYQzqMfXbVBQUoTTAUvUubvLAI5ZepafjoRiCm3lUV3K8WWpMi0haTtyEM9e+JU9SMO/oO9fRUJWnD2xKkxqZT8dP2VjYqqsbk6SDrz2+DS4ZMuW17+Z4oBC+RLh0+2frFW9JkDKnw/JffGVCJ+PpurA3eSCZlaE4YlkJCi/5jbJ21UZppxHAufLrQahndfJQJoXjQX3dZocLl12ymBxe9AAj9Cek9pjhZRnv+W/5lcL7hyaH8V3OitR+w5fbRKTuCEBgazIhjLx47UxsEsmSNdsSWEIOBrStrvNhoGlyrDQsj6WpVEpFJgQejtl/vQ9BfpO5Dn3qVgCZBchEw+BHpSpHiM6+f2JjNYANv0kg+PB34/s36YNLTwmZXM3hUpIdfVt1l6L6X+Q8A3r63ILvNteY6y+0p+yzHyA6DG72WUZB6lUJIVlsOicodb4frLzQdQOe1X7YmiH2gjUHIYAGzu4LygBa34L4Hx4Lh04kLVh+VjR9BDxPJ6M9+35Dpaw+Tr59Db2AUM6iTuFEaDy+jDDyJY5m5vY9d4Aq5OHePraJqBnPHB0D7DNgTrnEO/axR4K0SOAcrkcfZzP62PXQqV+2IV6GHcjpy/fm5vNQ68rK0SLB87gksn0spH2wixoLceq2UOu2JLCEGgTh7WaVEW+isbWF6fPNzDwJ7TAxzs7atRTOYFdtxeOnwSB1fqVB3PR6L9IzrAY1FhUtAJ7ClUQgkTIbFq99USNl3T0RiaGMxOdE9z5SBK6fV10MUSiNdCmJUK3pKgCefmWlF8uR7QjQ8q52ghm2qrJrQOT7e3yS3n/ZZALQmYYAhaDaYbtoQ5/OyuaW4u+/F0Cpr8lr/5AVUvUpmc84B+EVxQgU1wmdzyQoV+w2aBhxIBWALpiy6EwBOJPy4inHLE2BJYSfX8v6S66Q8cH9uV1J3SHGf1LXrQ7WrSdcl1UwUaNNj3zOOvixvgyqH9tR/ZxufyuNbzSWtt2lffFpwPz9QtDd5v9aMLtkR0+IY7d8cvPngjgGaamg5C1Er9KM0qFZIvDdeXgMnbPFCEbyeDSvfP+rJxKp+HtiaO14Y4zS14EWIr5SsVvd184C4zPUxrRm1pI7MQmxcdbb5a0AJWiOse2Zqhi4NosODRKaqlSQUpzm0bGgwrEdKnVF+2BUVGcod1zixQr2IC1NXtNiF1fqml9cs8ocV9yeXoMvmb+w/ktT+m4FhtzxM5mgxDCxUputA//nCDO6YFX5c8XT355y3MZiFjpzDvvnIWiVe7Tweq2rqHRBgd7O3DmdGFIl+dVjC6kJ6t3HeuHbpInuCWMIdLoNXgrpwifU6OUJzbZgklLbSW1s5X/ZfzS/I8rSR0y9SZNQra5Xye3OLDLeEBliBFMWmTuQKURUstk3IemLew9gjP+nr8HrpGXDdXB66fTn4LNUZqEQq6zkn+kQiznp/ziAZrbDkhyFogkP6qyFfCezBapdcUmXmX6NEOf53gLsgVIEhpfpZgzGmH2g+okTKTU362zQfuwvq5vTh74lRVAFe4apbTvurPA3dV17Vg3U8bgcLXK9HwrMIK1xJWwVuGLoqh+YEAXDK9fIz5JwUhMcuPfNBFUC2tnW/4vHMFjUuuK8fJXOH25i0It4w5fEP/p6vIsC7kqcs6CfTDYBO6a8i9LVKcRf0/9ydXlTrlB7QS6flaJfCI37zN+s2TS78LMaMkMs2/SzGr/26r+P9K7ts7hrNugKZLAahH/wGbmcD/MywBlJpX3Adxd/AyFJZiIb5Fy2LzFBurLauMnG5nW/gEvyq8JO4Pb6gERUr78x4qb+uaiN0VqWiWCCauiuNVtymF9hW2FX66Dp3OfxY/U2kSvPVxEqe9THRZM6ONn+xqVHlZ1DU4cOHB85Fa1Vg0OCDXZaxbygXPIx6j7C/WVYb0tSVAoi06htdfnT5npblays2sQnBL8ASfu+Um4HPjF7oHRytipSXcLMgLrKJKuLaCXRdCHQ/YC5X218fswY0uAtAC05Fk4gnpubSkvh7Usl7K5jlguDGZ1MfH5tpwY7YCCB647Hyb6jRYXYz2ODxwEeurg4qAbEELSLl+f2UDI4zHCQubPaeP1gIhlhZfAkvo8W9gWvlcfL3DjanLZPxBtxMRLh0+OeUNHp5erwSllWTmbeVtCU1QosO5m7PKQV0L6EEksNYg0ZFiwTwoL4yUX0+vF2JVp5Hjuj5Tiyf2fY8tN7SZtlFYXVfPQyqYkmu79PhFQq+hItAato4GL8Llw2NE+iRavZc9kmtPUCLggLL7knb/7w++5fkfHZZ1aYWfxD6JEgfuvAMcOXgvUqv/TG4t2d4Gmjyei4KmtLrN1YEpAEsd2Po6wHgwy2StTLeONBTNT7SQyp/2ns1WEICV5o2p9pUi2Iubor+ysTCf99kTp2qC+Mqh/a1ykhmsbZa4VRi5AMi8lKCZhGAI4dEQwsdCCBdCCA9Ntr0khPBbIYT/d/L71lw7vIC4Zo6z1pKKEgOLCZ+nYPWnJHo7i/ATgWUtA5pq1yveaiXFC9oMMnZZ8HuQyGVbAcgCVKLQwDRyqiclL+uphZpXxfpaQd/DmdEFnBldqARF6sMjGTfemjoptJkfIgj549Jl3OZcWdaxAi+owzSdeQjCmczhEMKjAO6JMX6Gtv0sgM/GGN8eQngrgFtjjMdT7fBCS1oVB2CSjdvc/CycwlybgF32qoS20vUYqVOoAxwMXTfQokwsIsfaCw6VmLQabEpz//neNXcuR/LmAg5C8G6DXNTX2m/RVto+ezGFpb0uZmkOlq/QKry7qOsCfuaJ5Qqblzm8CCH4CQCvizF+MoRwO4CPxBhfnWpHyusLeJDrqGqbst2CeUaL+WV4Zemt/60orKCECsPw1gzxfHqzfgRKn598wCyNoUsRUf3Mcgt18/q6nlCcZ7HU1HX00qVA++fPz11SNK15UYqUzz0XdBRwkeCu98MfSw6GpDK0UprktRaCfwzgSQARwKkY47tCCE/FGG+hY56MMTZM4hDCvQDuBYCv+MpbvqH3tn9gDnKdP8iO46sVGPGqvOS0uVRgRMAcPKA+8b02vGosi9Ds9AAs8etYC7/nwBHPkqUROG/WuoZQcrzrzNsPyFomYx4fa/64yDW6aJJAWovztC0WxED35TJLx5GlDXpCcGntPH47fuCaCsFejHEUQngZgN8C8EMAPlQiBBkcHWZovpG1Rq9glrC5da42d7tQUjxCa05wesd0jeRqTbD0WckEyCXH8/H8vthPl4O1WE8JrHzilKAr4e+VQnMp2xCD24BdRG01QD5XKw9t+6aVgUWZ5NaaQtazFe34mq47HGMcTX5/CsAagG8E8GcTMxiT358qbU9XOhZz01tTtUv0UpBKt+OXnaprmJuoXJ3Zy/nlnxRmCaaUOJxzKMl1luvIkgciILjEVApt79GKiotAtMzgWUr3c71CiejyWBWBsgjr5OyJU5XQvtpMiKsFmYu54Cff/6w56oLOQjCEcFMIYZf8DeC/AvCHAD4E4E2Tw94E4NdL29RakKcZbK6O1wyRvEUA1VoHJdCmgVfuSgsnL4/XqwlnHe/dK0NH5PTLbmsO6YnTdiIJ+bnkPM7RttYd9iD5xm0Gtqa4pIjSXasv83oqmw/cVa0ZI7nSKToLMD8XRRUdN5aXKAXPly594/uUZz8vocysEP2R1mvrzJsFMkvGyMsBrIUQpJ1/HmP8NyGE3wPw/hDCDwB4DMC35xp65LHbsHzsCPprzZXbLNR5duNyP22yRHLOYfZr5Zzxwu5niL9uae28mX8r5+l2rO2633twtNbXLs5p+Z07V2cQtF0TRNchzKGNcOoCJhpXGS1O3yo+3tpGRUnS2kiuZHwXaL9Z9f/QPj7n2hgNA/prsxHhNcQM3nP66NzXiWE/pNx7DxE47Fs019032zU7a4Ixxj+KMd49+dkbY/yZyfY/jzG+Icb4qsnvz+baenbXdAU0LSxymgGr0TrHsOSLwUJBIMKLC6W28Vd5pGk+N0cgZg2Qq98wry81qFPZNaNhKCaucjuibZdUCBFfVNsABNOi2hzvFW1NwfMZsqlr+dGuBvj98MdUkAt0sKtongJw+diRmtXFub6yfxa00fLm9S62RNrcDXfujsMn6ivFlRCOdZDAc6qmkOJEpbQSb5F0gRcU4fuxKs1YmMX5rDW+Ns5wKyrOffTKklmLuWvohZCsdLGSCeFVvWa6DF/T65NeUlL6IPc272CHhhU5tQIRDF0BWsAfbzlvngEMjw8KoPEegbJUObY02tLgtk0VmZf+/E+4a5MCZRHVkgim9wUTeAEAnX+rj83xBLmf1v3ktMvSZQRyyJnReoLr/F9PU+NJUUKO9sq/N44roMoA+ZQ8DoTkPm4eV02jS3TVIyMLtLBdWjtf3ZsusWVBu18EMjfmJcA9KhR/3HR1a7lXXR1GTHqpXNNWaG+rytIlJiybiXpAdH25YlpaX1xr0ElQhpEKgggsE680h7iLM7gk71pDP0NLm00JQD4OaEZVeRswXfCnK0qeHU9M6x3oPpQ8L+8Yy40A1N8fp8dJJJ2foX5nqeKq+ke2e8+lrRXgoaSd3sauKoipBSCDXT3iT2+znOes2FJCsAvkq2FpShbk2D2nj7p+HhZgWpiJT8SD1khkMDIBWKKtwvnS5cn5unx9rQGWJJhbvkOZgFbUWG+THGUZnExZWlo7X1Fh+B70hLQmsQgeaaN//OHOTnZJ8JccY77eaHB5vAYybUt9eLRAAuw1sb1VDeVYXlJBxozVLkfT5Vy+HueI9zZ2ZesL6r/bwvLveYKG55z0i90S4lvlPlsV3oHxWBYydqoWoe7PPAIzW6KeIFBWG42jUrqmWUrNFzWbq+EK+JycHzCXviPQRSit++qtx6reWsq/Yt1LF1Sk2UPj/3MmX46PJQJM+4BEwAuHz0uhqxZmP74rWXMvh9oH6/jDWJpcs0pbOzQ97szowmSN6qON+5J+6zWjrWdSZVCs7HO1WB6fVTS5g/tCorsllapT8ycnLCzhZwlljcr3qorDcn9KfHw8bgDfPSLtyPvpYRstvu7BEnqpiNjZE6dqAkawuTrAntOTVLNh8zwAje06NS0XKLAmvpXmVmJKaDqNnqBtIJpLf22iocD3cQmsIhajYQCG9kDlZ9JfG//ubexK0lAE4nPkIqKMHD2net6DjaryCT+v5WNHJu1OP5o6PU+0j+X1I7XJX3vnh8a/9wwnz3+1LIpd8r4b73QyFvVHRfreRfMr9Qe2cSuNhgF9+dv56JUEK+U9lC4NeunEmC526fBJXDfj4utbJjAi5fX5BZd+OdsUVqgEAjn6OZqmBZ8FL3iiK5roKJ3ldPfWKGHo52GlF+aEGlAv+CD9saCPs77QOhCi+1rMDZysU6zvx+qTtT+VSglMn4FO5+NxpseOZgnkuISldCMOTKQ+ZNxPThn1LJVSIZjra6mQ5OO4ug0wFYReHQBrLW7AHpte+pzu47YJjLAAtHwGgO+bsLQ1+ftgbx8O9vY10nLYL8SLppfy03TfN1cHtewCaVf2lzr/LX9Vf2WjIQD5+aSyFSxhlXKaV5o09Uefx77NlPugJDVNaCy5Fd+031EgvsXU8xVByf6+1Gps/eMPo3/8YfQ2dtUI1lwAVfosvr8UJJB2sLevFiRIHW/9bTEK5P1Yz0ajJPXRAvsGNUH60uGT5nuWMXnp8Mlav3rrsfKTHuztq3zLfG8ps38RaYlbQgg+s/smAHk+U26y6yCAHkC56HPOb2LlMOf4jHzuLNDLh8rvnLYs6Wt8falgnAMHdfgDwRNO2uHnoEvAexCBci0qO/Mk5XHjvaeuRViXjx1JBto85JZl1W3psamvI2Z1avxbWrSVRQX4LhkRypz+JhotB9TazIcUy2Ee0eEt4RN87a2fxoOZEj+5L4AVbdLZJ3wcO75lsWlZJCZFXD0zulANgP7KdHGiql0crUXERK3nAVxyP57p04UILimFOhCQu74+jmkdvI/9cJJGKAssCSyisl5pTRzd8wYL/FqO6ommqbUEf3IOD1ysBQCAid/zgfpxnm9V+yFr7ajKO9aKbcMDFzFCPe9c+wc9En7tnSl/udVv3XeGNQatvG1pr79Wd5doHqHGaHB5/B6GzXx1njvbLmPkiw/e2Pi6dImm8fq3QN3M1m1Wq9ORMBRoP5G1JiqXKGKhpgWsFji8vjL3Q0O+nqksDSuTITUJeYKktEGZgF5EXPs6+RwvY8TKEhFoMnppJF6Dn6nlJ+a2rcgvP6+SwgxW8VQhOHvLVFomrq6cLZDnCdSzkORjLO/RI+5brg8rUcB6Z17Gkx4XutSc7k8XeAGSRfgEt4wQ7L25ngVdWu+OJyvDCjhY5rY1+L3Jw5BULV3q3HKCcxteZNnqKzuXuX3LsZ4zs7wMFi2IREDre5VrcFVr1m68SQE0F0riduQ8/jDM4wuvhbj1zDwKjBZiJZkvlqkqbWiLJLUJGkkGAAAgAElEQVQ8gn52nhC0so88jVAfZ0X1vSUBBKmxOm/oD0JOFmybwAgwdaJapq0OeFSmnhGFsjTKVDHKsyfGi8qcGV1o+IjkWkLqZee6t9ZDSngC9oSxIqyW2ZSKyqZIwNZg1efrrzf7ErX/ScAaLd83H88BI/0uqg/Zxi70jz88l0wFfawnAOX6/Jv77BVYqLUzKaMmpcAEwlW0sjiW1s7X/LJ6ZbcU5Hq6Xf3etUDU1+fgjwhXEXyWmZobP22zfkoi2osIgljYMprgFx+8sfpfkzUt86girPb2VRMxVVaI8xMtSoRHrdAaD9MrSrW/lNmrv+KeGaNpPHx+aXDGgyw4xJWSgSZtARibYTlBpYUkg7VKwPYLptDWH8Tjgp9/ig7F1BivTL8FK3hivS9vjRhgqo1bmh4vLm/lc1vXFuQWh/KO4fWQGTxe5Tlqa8C6zxJTuc2CVMvHjmyPdYdfe+una19kiR5pDVBraEBzAHgRJPGDsAmkhS0w1Tp1VenN1UGV5mQNOotg62kfDO8LXhvAVJvQO1++9rOYJ177pTw0De0r5XL7klcqKI2+WlFCK5dbQwtAWW7UQ5eItab5sLZmvZsUtSdFeUn5KPl6pZB3UfoOrA+0x8rQaCMAgfnkMOewJTRBXmPE+uKk1hvgKDDQzARI+QCtAehRYEqDNEz+BOoBFvZhpsr2C9gfxH3RSPmEcpNB5zQLrHVV5IvvPTsBt5cb9JbfTGsC1kSwnp8V7NGBEsui0L7A2r0UpKtpMFle2rXAmpLnj/V8dnKd1Pvl989E/lwpMUsL1dt5vy44a1kCuWeQQk7QbZvAyD3/+XjJTR3NlL9TQshbKEnMOe0PZIqLvpZGKkolsLIRrDY0xOnLUWoGm+I5gclCpIt5rE1Dbbbq46wotPiWtOPdi0jrMlECMbmta1TkYGVGi4mYMqc9U7pN9WsNvjZQ50hqIegJOC0MpU8Cfb8smEuKxwL5UmNa2OfMbq8EvvXx1YKwjf8w5+YCtokQlHqCrFbrhwqktTEWJDyB+dxcRFUXQgWawgGop12lXracX7XvRJ1TVTO0acH3xBNO98cThLrP4vuyJpPWzHLas6U5erQcjgqLNuj554CyRda9tX4tOgj719qCBYYlLARWFDZnfegosY4OM1KaILMLLN+g+Pr0u9cRaB5XXnBR5hU/B29utIXlZ9dYuBAMIfwSgL8N4FMxxq+dbHsJgNMAXgngUQDfEWN8MowXHHkHgDcCuALge2OMH811QigyrGGIYNApOgJtEuu0Msuhz0ERr4qt9dWxzHQNPeDZwc2D0ZoYFjHZ2q/vy9ue28+ml6VtajOWJ3zKsa+FcKlG7EF/TLTGJRAB6gka0fRKltu0BGPKJPaOt/rB92RpR15AqatJbJnYXl64Zw3lrDCBN++kHznqk7f/amiCJRkjvwzgnQDeQ9veCuDfxhjfHkJ46+T/4wC+BcCrJj/7Afzi5HcRKirDsM7rEyG0vE7ZCphqjbr6hzcozATyw81jzG0Ea8DwYLYmYWqgWfwt+d8qWMoOdw8pM1omhlSWsXxp3n2ONahJgON40yTz7rFEAGoBwRokCxotdIYHLgIHgNFa/f6kPenn8MBFrK/6Cyx5mqE2PVkYmabw6vTd6L54z0J/UCwfmrmOyioapjJfQ19LvxvLjVKfJz61TPu/LXgReIEIP09ALh87AgwXs6SBoMgcDiG8EsCHSRP8BIDXxRg/OVlb+CMxxleHEE5N/n6vPi7VPpfX118U/op4/LuSiFNOrU7t05qgZTYDZYRkDe1HS5GrpT0grfl5tBsPbO7wNS1BarXFwQ3v2Xh+wtQEYQ0O8AWX7MuZtxz9zJnBlvZnmb7ajwb4z9vzCZZCtGHdT4t24n10PFdN2+ystmNTTGlP2/OEIFuGXh+vFVn65SLYJr9fNtl+B4DH6bgnJtsaCCHcG0J4KITw0LN4BoAdsLCIuIKSyKlVXikFqZShU+hK6RxtghJ8rNyLJu/mvuReu7lAypVD+xvHcYUQT5PM/e9lUAhSSf58jCcgRoPLjfVC2tBBctVq9LU8yLOzxoK+967+MI15FpqQPpUU4ZgHstpe4jxgtsXGcuiqCT4VY7yF9j8ZY7w1hPCvALwtxvi7k+3/FsDfjzH+fqp9WW0OqEcfPXCaFfu0BLZan9b2UvskD5mR89OlNEJ9rmhPHhE5Vd8v1aecQPb6UdLnkvu16v7pAIfVJvvLdAAiReydB7Q2qbU/S9NNaYOspeXoIgLfDVGPRLfxUwvk2qkMqlJwOmKOkpX6EFjzzssft469Gj5BC38WQridzOFPTbY/AWA3HfcKoLE2eQM3PP40Nn9eHLAXzCoXTHlZAhUgOAHsOT2doOxH1A9sNAxJOo0nCK0oJ2Dz4NhMKqUvyHHeF9ETeCkOV8qE9bQ0zg7I+S8tWP5K8T0Oj08iwGj69qzndPbEKSzjCEaDifA8FLG5fhdwaLJG9bl6G6WwTGbepvf1VzaASZqaJ3jZdE6VESs1f62xNhpcRh/jCj192ONBj8eUidzFx2ZF2asq25hWufHAQq0k6quTGbouLZFDV3P4QwDeNPn7TQB+nbZ/TxhjAOBzOX8gMK4naAkm+dlz+qgZeNCqsq55p5EynXQuqQU5l3OMvS+cXthHfqe+ltYAtiaOJ5xqfktVGNX6m8GmcUqgpsxsS9gurU0XYOf9WthIu/rd6kCTl5+cguQuW8JJNDxrH+fSlvDsdFFdoMxfrcdQzuWgz5fxmENKQAramqzVKn5OcQeNkrRHFn46cDKP4hoaJRSZ9wJ4HYCvBvBnAH4SwL8A8H4AdwJ4DMC3xxg/O6HIvBPAN2NMkfm+GONDuU5IxojASnHy6A+50kkMMWu1KaC/ToBdzl1XOtEv1HMUC6xqKm2QMoNz5nLOR6jBZrDFI7Sc8ZZp3NCeHZqHQH+8rEmjzWSg+Ww1ctFgLej0hPZoMLPAo1VpcPDMKl7RJkjoXZ+RC0S2uZYX8EtdSyOnPW4LsvTNt+6On//4jZVZCtipWsxZsvJ+S1R8b+Fy7YOwktgt0iifZwVqNIGU29NoG1TREWuthbAAYr6ihuc/bBNdtsA+NE+YaGjStDd5OMvDix7nosH6Y6SFIEdj2RfpRb/nIYBSH5YSLmFbWP1IZUel+uq1WWJ2e3OX/Y6A7effFgUUgLwA06RNLqYwzz4AqCpBW5wqPi4FOZaXH7TMQgFfTyKYOe1NJu3m6qCmuXjC1uNztdUSNcQU97RPS5haVac9WNWDGG38giVFAkrqCDLaRn9LBI+FFJewTR/EBPfOSUVyBVZ/uc2rVQZrHtfZEkLwNXd+ulHz7czoAs6MLtS0goO9fbUKLfrH8x/oF2oJIXmYLCg8X5gWJtaAsQSRxe73BISYalqoeVxAvhZrcXw8By08wdXWRJdraH8nE5xFqFtJ/DpQIYsXCTwBqPvJGpslxHg9E8sdofuma+zpdoB0FZhSaGqUhrwj6zpiDsv6Orm+dO2vp3AwIdsT7LlMEQ2es3ye1BrlOT4v/+CWMIe5snQuMduDRbEQCNGSzW0rT9lqy6OdpPJcxWzRFVJ0+ykfmkBHG7uYpqnrp4RpznxP9WEWTpuXR2z5bNnEt67JS3pafj2vr7kKLYISkzRX7KCERJ3LOmkrbCykSszJdr2N5ydTgLpQcLSf3RK8vKSA7N9WlaWB8U3KEpmcqeFpPJY2s7R23qzAzHw7rlfoRZS9KF1/ZaMi9FoapAxW1nIsDVCqEW+uDpLRx+GBi6ZQ4HNSkWe+j956rLRs69lKf9tEiKX/0gcR3F1XaPOED09yqbeX8oWNBpdr70X3yQqEyN+lQl/Ggc61lp/SNkqQM2OtxALRxnJmY04AetusSLjwXS1tjf/Xf2sKjQWvnucs2DKaoJClBTmNpg1RlB9oSrBynrLs56+Ojlp78KKzbSK7FgFXX7ftYNAr7Xk+wtRztojMuuoIb7fKM+VgLeCtide6KACQjgB7qW66X7pCjNwPm9CeZmfddy6KbBU54PtMVdixkKqXyNsA39WQE0QaB3v7zKBhCvoapde0eL7bIjosucMMXRUGsGuVMUqcslZ7Al16iE1fXftPzk+hizmZo9VYbXmmbM6FYMEz1dkc44+DRNGtKi1dS1VZYEEg5e6t0ln8/rQw06avt/qdJQSr9g0BV8Ih9JAzc60qM57/rU0wIxdsagMRgkA6vc2i3Xh0nDaZXdvCHJbF160H2FuPDZPXmtibq4PKQVwCTyjptXT5C9yWsW71uUQoaRNaawJ60nFUWW/rr2w0Fq/SboQUr09gFbuVKDoXL9hcHcxkCltorH1Cwkf6xRWFREhYgQ1Nm8lFgjdXB9h84C5zn5i83n5GKSWmSxuCqxWR1bhyaH/RR760fyWE6nliS2iCeqElhqV58URsWyCBTUHLLC6pTKzrElqFXFmrKHWyCzS/T8MiL6fcBEIvsmq+lXITmTStgwFaU9IL1M+S52uZghZH0CJZWyR2KxCSW2jIgqcdWuDyVqUE41q5s5XmsqYaHscud8y8wGZqm9Q4C56GyPsY28Ic9oQgR3OB5gT1Mj8EKZXaE4JAXhDq7BRdYIEHr4dS3yaAmk8qte5ILiii0dZEZyHI70WblPJeuOJwDqksGtY2ue/sN/Si9ezD1AEgWUlOE6G9gAkfwwKtjTncppxWV/rNLP69ttBzITUH52WC83oxwDYRgl7anK4CzZqRVWVagzU3Pl6EluUP1LDW3dDbtU/HG9geJSYnLHmRd32+9l/pcy1tMqXtWT7BnGYK1Hl62lfI8AoYAHZl7ZSAYfoL35el2VvLFAD2kp+WdsvHpqrIaFjphoJ5msKMq2UWW8ENq4K6Z9527aeez9vCJyjQlBV5SFLzzBMYqcq2y8eO1B7YntNHa4JU9rn162ShdXUN+d9yansDtpScbAkaWeBbg/1vOa5fyXWsfqV4jfL/+rm904IDw1D5iaR/4p/rrcex8CJ/nV64XHDl0P6kb9Gqus3BrJLV7kaDy1V/hgcuTt/3pH85zTCHxup1k+VBSwWg+B2ZipM6f5ECUFNeWLjxdj2uLOGXy0vO+QPn6S/sWkprrnjksduwfOwI+mt1MvLyup3poQMolw6fBA6Xqf7iGxtf58J0G8qDHqKK78FRYNjUOkrIrhZyWiEfJ7Byk3ORYdlvRZFZ29HCrzJ7SRhKFRPgAg6u7KtV8uEo8tR1MTGVh0eLzOT+ygbgaFpa45AF5PnvnrO6Az/rEkEpKFnwic/TQnB44CLQssyap2l6VKF5ImXC7jl9FD00i5aIpVVC6Lb4gTnS9KXDJ7G8foRM8bfMcIdbRBN8dpdfcl1ebKnkL/kSzlJNV4jS3BYw9dVZE2pWcmdpZRQd7S3ROjXaaDq6Pp01IXXgSrQG7gub+m2wuTpA//jDY7+e4aNkTdDTuJmozi6FVH/kA1HKeeT2rIh1SRtef2ok7fU4t5p7JWlp3hxKFY/15mebiLAupDErtoQmCNQjmDWuF/ZPF1E6ITnDNieplOVuIUUerq1vQAVfpb85s6s0AOH577wJwP61nBDx6DAWMXr93F70MY3qSgFTgBb3Icd0auDqZyrnHFybBpNYY2RYlaUb7VPqm+b9YTheZGnzAdtFUT2LSbYLQObrwCd9V8cirQVyez3YBUdLtX/v2HlmT7DWlYpEayHUW48NK+zM6AIO9vbVtNUSPmAbGo3MvUfLbs/FlgiMyOLrVgBEoE1NEZrzgvcFLS3132YwW2AtWEd/U8GJVB09r1iAJjtbfUndi14GwFu+NIWDvX1m5oa+rxz9pLHS2wRs6lsC0Irea3I8t9GF5mP1fR51CBk6Ut9lwaQc0ZpNVO+9Wz4+JnenPpizBEh663F7lNK6/nL9gVoCUBzKs8ATdDp4IpDq0Snkyn+l/meIj84r188LBLEpxhVnNMTZr6+vfWEp3yH7Gc+MLmBzdWD6XvlZcd6oNeg5OMVCQe6xzVoi8iz0R0C0Q7mOfk9s9nKeLwdHGtphS1jnSY512zZzGl+Ji0fGf46qMhqG6hl4AQzRALkie0kwg4MiqZzm1BhaPnakupd5mMNbQgi+6Kmnq7+tm+JBKcn/YopKwYUSeINEX3NeK3CV+uS8DBgGR195wsqP5zcsKQYgGo/8sOnMf8vHQiYJL3sANL/y1rucNaqn6xAyhxKoZ9ukIsuy2p4F/nBYmSqzgkvSM7RAlm0sjHtkugNNv1vKJyhjOuc2krJVgB9s1GPKKpull8KQMSPH57TPnAJSmk+dw5YwhyV32DKDU8JI+EK5FepLkHMoe+3rNYlz0KYmmzNidllmHZBO82Kz1/MdecEByy9obedCElZlboEmL2s3R4472DbnWEjPDA5IeNkkHpeR+2AVUZjn0pceuP+WZqxXBfRWVmxDUE4JJa/IQwraHBZBaLmYpL9eBWkNHoML5wmGEH4phPCpEMIf0rb7Qwh/GkK4MPl5I+37sRDCZgjhEyGEgyWdeGb3TQ3fTak2tkgBWCLcul6bv/ByLdYA9ETrMvFKneZMmZE+WPe+tHa+0ga5mG0JpFBpbi0QOTZ3DEMEYG6VN0u78N6x5Ubw6hXOCxYVx+oHgMZcyUVdgcWt1lYCEbDcV+4Pc4L5f0ab8dYGJebwL2O8cJLGaoxx3+TnNwAghHAXgO8EsHdyzi+EEK4r6cjm6nih9TamqBzb1sQqfZA86LxzuhZVkN8cebY0DdZUmMIhvkA9CdmHJ791Jog2qXifmJh6v5zHJg0/A2tRdfYtaQHFNQdnCSgJpP3quSh/puXbunJof0Pz9gJM1t8AGhpoFzTe4eR//ijqzBbRzHK1BdmcnUVZEG1O2m1TA5HBc9Xrz9VYDJ7RdfH1+wF8Icb4j9VxPwYAMca3Tf4/A+D+GOM5JKDT5hYJS3O0BJmOpgJ+eS8NzxzNpXBZbVjanLVd18zTbaS0K8t01UUaPKGp+6P9k54Jz+bmPMxLfV86wsul2RlCqvaez6JN31wASL9TLhfXpTiBhS5+2hLzmI/RhGi+dptcZz5e5vK1WnwdAH4whPA9AB4C8OYY45MA7gBqK0M/MdnWQAjhXgD3AsCdd1xdumLXwVMx1Ifj/z3HrLfmMVBetYX3W6ugachEHqlzOQrrCU4PvLKa1iRrwmtjl7moOvdLoya0W2Rh5FDLCRaNyvEvLa2dx2gNGG40+15VrznQvIbu5yyC3BN+mw/cNeZJHoo1YVMzJYfNYgIC/WFP+XCBKWXIE2yaUpTTQAFMeL0Dk0doHo+0gGz0f2juao2u0eFfBLAHwD4AnwTwv022W/FqU1LEGN8VY7wnxnjPbS8tspjnAvFJ8CARs1oGiiWYUilxbaAjrto32KAH0YI/Yqpyri3TZaR4AUc+9fX43lgAsunKucBsGjMpWc63hKimvViwtK+UINHX0f/LfXEtQU3UFcf8aBjQP/5wg6qiNVgu3mr127rHVARZPy/9t7xnMYctQcPCq7+yUfFKU24ZFkIelYXrMLaFmMmi7elrsOWTok5ZOcZevrL0eR7opILFGP9M/g4h/FMAH578+wQAtmtfAZhE+asKL/+wDWpZLGJaUW5zW7K08CE9U7i67kSDqwQYTdpKiCpng+a58fHeQK/VR9zYhc1zMgkfrp0/AtDHVNvafOAurKMuGKRPKaHGmS5cJUf+NnNuQebhqooAH6oXTOB753vzzHPJhLHA9yP966OZ07x+bm9Nq9So+idZJOsRo5WJOXwoLYQ40spaEH+UmS3B0FkdXWAJSU9rswIfVrkzzzzW23i7xvgaP9rmVhropAmGEG6nfw8BkMjxhwB8ZwjhhhDCXwXwKgAPztTDGSFLcXr7uoC/bDIwWJiVaItLa+dbpUBZGqLlp1taO580x4H6F1Qfy8flnP7r5/aagY+ugY6aKetcj8n0OstGP3epSs7te+a/V6UHsBe6svqYM4k5P1760XYZzKqGoaMFWZriPDQmTZ7W2p3uQ8P3mhgTcnxKiC4SJRSZ9wI4B+DVIYQnQgg/AOBnQwgfCyH8AYC/AWAFAGKMFwG8H2P14d8AOBZj/PLCel+AkgEgLyxFldFCpJSg3QZyHfa/ccBBNC/tp2Nsrg5qZq3li/RY9rxKnqCq9WeQePn6NYKx4rJpc5L/Zo1Q+z4l+l25ApxVAfvHH66tdldplJPy/wIp79XFf8cfhC4ZJOvn9lbX33zgLrNadgpWHUKpzuO9WxkLs2Zb6XM50qwrvniQPqb6wdxGbQLrH2B+VXO2BFm6bXTYevBtAx2c32jljAL16Kouod+1OKaGR1aWfYyaIE6UY7LO87RMfa+pdrQpaWVqAPWFjfQ1vOhxLfCSqBKtOYEeaVxDm8pWUIOJ0R7aBkGsxaDkfrgf3nH6HIY8K6/SuIDfRW495xLkhJ7lktA+Wz5OtqdKagHNEv7i23/0vh99/leWvvnW3fHzH7/R3MfamafVdaUJWBWmATQmZFX5hOqksbD0SgeV+gktASjbtXCR47QQ1Hw3fR9At+COlVEhPjtvqU0+jvdzX/X96DJcso1zyi1StPbrNdwGqjqQ7mdOqHWJ/OYEn+536hygnuEC2B/n1MdToJ9zrrCpPk5rZ/x/quiHXBtIC0Fu1zKJvWrysxZQ2DJC8O7X/3D1vxXut0L8XYSffsAW543/B+olsxg8Qa3V2EoX1hFoYWENJl3+3dJ8OM2rVAibEXFqw+IalnAQ5Th9PY/XODxwsSqbzwKMNSZLQ5V90m9GRR9S+0s0Pw1LKOp1knmcyBom3sdE+pMTgkD9oy0CjQWrJQhlvFgadi744KW28bFesMnimMo70ILQ0wKtPmgsHzuyPYSglNLiB2yp6KUcohJ46w9bA4m1EvmbaTYeediqM5gjQafA52oBBDQXK2rLSbTAJqhcw8pgsb76uh6f1b+UqQ40Sc8CL983JdS05lxaqUZfo/JjOilz1toolvC0iNApAcAuHK88mLUuj7TNKI2+pkjO3n3qcWa9e6Du+mg7t/mYbbHGyGtv/XRDAAL1iFSp6l4K/ipurtZrEwrPTgILm6uDmqYnObQSeZMIn14LRdbaYOgBooMhHiwh01/ZqIIGgJ9nmmrTOq4WoCBTUk9kLlWl27WKkEq0sjqGIr1c/YV/SxTd8xW2RdtoLDDl90lwJlUdWoI0jWo3mb7zGEtVUDl74lTNdPR4c10DQPpauvwVMJ1/1rvnMaEj4ZZA5srxXtks/f+8MmUEW0IIAvWwOidKWxGgeT0AEWJWWXItlGTtDBFy0jdZ8N3q0yyTVgSJxXeT/fwb8AnEHizhK8LLov5YSJGhAVSLK0k+72gYGkul6rZcTTlD2al4k46Glir7rvvOz5LHibSjTTqvxL/8bfVd03cqDiBBk45zJmNDMCa0XJ5zVhtWShtro0yqB5puEqA+Ps0PrmGiWxFiuX7qA9EVW8Yc7r35Pnd/KbG5yxdC/I3seNa+Fi9aKRqkDt6kJpk2mb3orLe8I0Onr7E52LX2ndY4rbV4NfQym1X/6LmVcjXZd2cFEVJ5yLwmsmcui3mdC4bofGo9BhqcxJa0Ge0bK/GDlZbFYkElsK7VlVEhbef4nOwvTLEQ9Do0uT7oZ7Ut1h3OCUFgvpUlrCiX5+PL+fe4fx53sDRKzCihwLBDWvs1PXqLxy2UY9jJzx8Fb0ALdOQRaJZxst5hafEKHRjh4guSZ2sdz9DUHU2T4XvX4CUUrPFgfZAs4e09K0GJcCsRktxvwF7wyruut8/KN9ZFMCzfn1yftWEZ23qpBp6PHgdRu85e0EKQtbg2XxIe0DLB5ctlpR1poXhmdKEmMD2qgocSoegd40VG9STTgR6vGKd1DkcUuQ0PekJ7a08Atraun6UW7tJ/FsBaMGvNIyUEc0EFva3kveaqwLCfkIMBKe1OC4E2mluJ1piivGhY+7Xl5IHdEHoc5oJBur9W8HTW6PCWWW0uB0uTyGmHWnVeWjuPJUwnt6U58XV01WhvIlTO/DW7H100wRxxmZGKCleCDdOS7pyPzMeJ0NqDo+ivNctp6Tb5fBEYS2vnq+fAwkYTXfk87x6HBy42cqUZ/ZVxLjNWKbdbpb/V6DGrzZJmnIVikXz1gksWNlcHyZXnuB+pqj0WuvjBPAFo0bzaQGtl/H6tlQKt7ZUb4BCqVMdlNNkgljaaO6YrtkxgxEOlCRgvj7lKFjjKnHOGy3WsmoECmRzLx44UCzXLQZyL2mkhxguae9fl9LPSPlmQe6uEpyrsyTm2El2W52+9By6XroWf90602S6Cg6POoqleObQ/W6ihosKoPFoWhpJix6wADhyl3rd1fbmmFUHn+9apYIwS09iC12ZKaWibMZITzN4aLnp1PK7yI79Tpjpfe14BkueNJri0dh44Mf5bvo6jYRj74ZzqHwd7+8wvuNaY9MuSL6ZV/LQ6djhotGNNFEs7Wz+3t6rEIuB+6nZEY+H9Yl5xWas+6iuvac3DMplTHxnuGwBg0LzfqnRVwpzcg6PZ2m/aj6m1K23iyrFLa+eB47tq/EA2tUbDgN4h38yVa88Cy/9n1SjsghRV5mohda0SvmF1zKS+YPWuhwPXNaDbLOU1dsGW8AnKQks5pPyCAo7ylkAmU0oDnHWSAGhEWT0zM+UH9DJTmDRdmlOrJy0LMPHdaJqI9VxZCHaFR9Ph56Wd76zFecJY+zX1olbzgiat58BCWnC1hFpOeMyjH8vHjtQ+NJWbhdxLFlfQIk6nriHYFusOd4WYbgd7+6qfLgIQaOYoWuanp96XmKBcUaUtfUXMMqn1p/slfC0RFNaE1JPe8k3xs9OVaiwhlyMdlwgaFuRi3go2Vwe1+7La1XXuWACK6S4VeKzyZfNAl8i/kKPlmafMwXlqPVcDZ0+cMi0LIcfzx0iOLfVVLuJZbBlNcH94Q1aYaMf6rKic1CAAACAASURBVNVbvPblGiXRwVxkrEtQxDpf02J0JNWLhloaIjvp20apPZ4X4GuDqYCAdhfkyNc6uqvzqK0qz1r7TmncXd9Vl/Ja+mPYlS9XEoFPtQPUF/yaJS9fwGNBa+apSjZexFyUFCu4tm1yh4dPTP1zKZoKY5bB7Alc9pVpkzPXnzYk5XkIR90Piz+o4e1L9UVfS/fdMm24Xa+EvWXO52DxBEthXVuu3/Ujx210qTajryVI0VrmbTozPYWFjS420oaN4QlBQY4oXpo/DGwzc9gSgByh05iHaeNRYKyJkVPZ5ykAS/NbxbSyzGQxBeVH9vHfJZB3UMob9GAJX+9ZWMd6xVm94zU4uMTP1wrKyDGz+DpzfeKPQ+n7TglHy1RMmY8507LtvZcKZ2sZTh3tbUMInxVbQgje8PjTAJoCyeLKlfjfSn1ROhor/gouhCBgs5iLAJT2qS1SA9DLBhE6hsdJY7OQBZoWAtY1vHtk31ZqKVEvrzlF+WGhx6asJ8RzQkd8qpLLa5mwLJTb+BCtMlklmiH7QedB/UjRSzRSvDtORADsoGRKIOkq4F3vqc39dMWWMId1dDiVgmbB0ia0Sasjg56vS5vAbTSgkuPaalMcGQZQrZObQs4EtSLLVt1Aq3af95y5r4xUqS8AbvRXm7tM8LaKt1o5zilzWdcY1PdSAs6X7kKDkSi3JiFbE700b7gNLB8cUPcPal9cri3dV90u0LS4WNh2IYUv3CcYQtgN4D0A/jKA5wC8K8b4jhDCSwCcBvBKAI8C+I4Y45MhhADgHQDeCOAKgO+NMX40dQ0JjCwCPLB1ipwnBAE79zFFr2ABkBOCVuHJVN/bappe8U6+vpfjafVB4PUltV0LOU3haePXkxxh71lrYZTL5dV+zVkCbXp9lFygxKLJWHUCdU7voqg0ll9QwCml1jkpiJDVedZAM42xNGWOt42GYfHl9Scry90eY/xoCGEXgN8H8G0AvhfAZ2OMbw8hvBXArTHG4yGENwL4IYyF4H4A74gxJh0eWgh6mh0wux/Qa0eueWZ0oaaFcjQrN0m6BDtS52jBogm+loCyBLkloJjDxSlv3vHWNT3oitAcpfVMdS/Q4RUx1U52q4yVLozgRcSt5RJK88Cre+4YHeacaBbEJYJh3tAZWOyf9zTBEkFYrZDnLL/AUfE2Gu9V0wQbJ4Tw6wDeOfl5XYzxkxNB+ZEY46tDCKcmf793cvwn5DivzUVqgqWwCMlAnf5RGom2NBVtbsuktPJ4vb5pSDtW9gOfowWMlK+X43N5n1xwItUfeXYWZaZLFFdrzKkPEgsQq9qxR/ZOMQFyZdG4jbbmsK6CbVkFbSgzFrqamFogWdlAFqx3wtA1FS0hWNI/oO7PvKrR4RDCKwF8PYDzAF4ugm3y+2WTw+4A8Did9sRkm27r3hDCQyGEh57FM9X2RQQZLMh1WGuyvnbyYksi1PzSc5PCIv+W3ruuIKzXGtbt6Amm70kCG/IjwR7JqZXc5Vz/5Fw5rw0JFqgLP6s6jhaAXn+42rbumxaQVnCMz+mCbIBmcLm6XsolIhXMZ0XbvGANnU2U+9tD2+IRGmL+Sh+veu5wCOFmAL8G4L4Y4+fHrj/7UGNbYybEGN8F4F3AWBOU7V3NXWsVqxSXzruO9TJHw+BWiNFtikYxAjDcuIgR7Lxga/B7fdLrbOgobNX+MJ3DPMK0CEIP0wotVuSWhYHnC+XUNan4Ivm+6+f2on/84cqE9Xx/uoCqfi69Q7HS8vgZpASstFWl11He8vq5vcBq2q/bdgz2VzZq1W7kPkoL4/LfXYvhCliT40WMSs9l6OVnAQBDcius+eudWNZFfw3Vc2IXyeYDvgYKTE3xRbkGioRgCOF6jAXgr8YYPzjZ/GchhNvJHP7UZPsTAHgR4VcA44pHHp675SZcef1sjmkNmSi6QrQe4OwHWj62H6NVu3RSqZ/SE0Bd4AVjPM0xFbmtcoEHzaKiWuvqbeyqHcft1a43MW/6tI0F6mhwGf2Nh7Ek5dfVdWv9onsW1Gg3w2lVmysPxNrxetywEJfzS/y1OZ9gqZAs0fI8dPErM1hQjAuWNgWUV6C0Mf/omXMxi9HgMvqoz5FUYMlzQ4g7aAmT/Sem+61ADPd/ngIxKwQn0d53A/h4jPHnaNeHALwJwNsnv3+dtv9gCOF9GAdGPpfyB2pYA600cscTpnL4Y3/W39Zf2TArtQhS0U/pnw4uANNqIqXpaUwVsJbz9K6fa7e/slGrpwfYlZpT0KuYWVFXFmbDAxexvjrA5jmgj6nfkTVSpuBYJpf1zMcaiO/H5PEj46YN6Tln8nvPvKsGJxqjaITDAxfTGkMLpGrvyfM6uDYOAi7hvE0vcpZX0HxKgTzzVJk0oKk9j9/neO2ZXER83hphSXT4mwCcBfAxjCkyAPDjGPsF3w/gTgCPAfj2GONnJ0LznQC+GWOKzPfFGB9KXYPXHbbSbXQFCgs8eEWYeAUxLaFaugiPdV0OLnjnt3GAc/VrxqyRcd1nTZfxruVxDvk4q6+yoJJX0kyQSrsTNLRVwz9oLc/pVd6+1vA4hvOuMOMJPwsl68kwrL561+D7bVgnaFb5Acq4icA2yR3Wi697mKWQQQ5dzpeJphdq4v2lEUCLuOpx4rQPkvdZ/fPM+pQQ9Ph12tdltcEaMt8LX1/3zVrcquT+U/BcGF7FbOlHikM6DwHKrhprbFg+traCUC8joe/bo6sIcoKQg1gedcmCxzvlcluldBy5v1l5glsibW4WcCRTb/fg5Wm2WY+WJ+Ke00cb1+e2uPip9FfAka7aMgBMLqU0vdL+pu6f+5DSAsWM99Yl0TQc+bHuL9cnqeysTaz+ilqrWJXP123r52tFhC3IOdYk1tQmfY3SthmSwrd+bm9WAHbB2ROnKsGnn5mXLpjar9Ffsdc1LgmSaHaDwBKAkhOtlQN5T7MUrBU8bypLAxPyJsbkTV352Qt4WLAGevXVXCvTNMyvGI7WJ/FwUGuLTbbpwJiuwetV17VycuV8L0LKGpanwaTS4wSVcDlH/SFt0MoCmfqSxutHcJ8sIaX9sFIj0OMZimM+h0qYrw6qSs+5tUAsaDoVt53SEr1nyv3S9RA1vOhrDsvHjqC/Vve7SWRcnoE2fXPpiCys+dj+8YcbS4jqe6iOlXdNbIvhgYsYrVxGfw1YXq+XyauCJmhGqC0qVBc8r8xhwKdr5JAzG63IYM780ax+LiHEgkr3XS/0462+po+12s0RVHO+RR2dTcHK7eU2rbY4QAE0faj63uQYTQDXqVYaOb+vF9ji/vCzYo0kVdZNuxd0u6lnr1daS6E0e4KP/77HlpPneKR0gWcyW5qilfWikaMM6f3aZ+i9w9+OH3hhmMNCwAXK07f0l9sTaLxdiMfel/3Kof2u01r+rigSkz4z30mEnlTUtQa3HKcXBdICEcibTxaBW7Sx2te8xUfFMkGsicHPQZ6d9Enft9asvH555o9MQMs9Im3KucICqDT5SdVpvt6e00drgsUio+sq2Lo/lubYBl5VGV2FWlNeSv2Hugp5xSGd/FhC0RKMnP3CWSGWK4IFnNZCxbfIwpDnrfd7Vmx5TTBl8i4aekLW/CaJEkGsEXpVdEtgraZnBR0E3jKHKcqJFQDhRciZEuOtX+JFFHVKoNXnFHWlQcZGWdEJwB4rMrEtyoyYV6lIp04HzFXbSZnP0p5XSVmQW1VOQxYJO3viFA729nXKaRZ4lblLgibafcGBGT5fimKk1uBhgWyxMGbVBLeUT9CLrlp+ly6Rui6J8XLOOJl8MC3vvT79GudSjXKlkaxjtBbCpF+LU5cykVOwTOHexi5snvMFVKoCTa0M1sbFypeohaqG9V6r6Dr2Yx2Taxq+QB4fmlfJJpVwFpfWzuPSiQsV7UoLsBzpV/xZFfl742K1Ep/0p8RHuLR2HgfX9mG0On0mqdXavFzg2lhBnUgO1PN1u0DeqbUAlgV3jRulYbIGLv9bHznexr57YPJM7/tAp/sSbCkhmItoWttyNJG219eTgpnr7CfytC5GicbnMeMtTVIXKJDrpz4eAjErvGds1f2z4E3szXN1rbK/MiZFi9CqzhvWtekciVkEm+eH033jQBBX2xZNRMi/y8eOAMN6ul6XD2t/ZQPrqwM3UMPBM/6t30Plqli/CxrWwkteGXrRBCvK1vHumqAIJRFgJUiVRksJzzb8xHlbhFtCCL7oqafNyWBNWu185mNzsPhffK7naxCtT6ckyYDrAqtumwZz5eQ6Xv05SxDmiND62XrlrWSf/uBIYQXmu1WC0AlUWEEI60OW+1+3Kcfw+5AJKc+Iq2/r+wSl45V+VK0PAo8hHQSq7nNoZ+hINowFfs/jFE87G0Q0waoQxrnuNBK3sveG317uWil3Rpesm/7KBh5tfVYdW8IneMOdu2PvzfcVHZuKxM0Kz4QRH4tljnTJYyw9R2da6Ikwa31DL3J55dD+ivagj7UIt1bB1BT3TUe6rQ+gR8C22uS1hzl90YrcWumCevJ1HU85LVuv+GftT13X8o1p6KIHpQIwtegTgMrHKP0oxawFITzwh3ThRVWvBqS8fmnUF5iPSpyrFi3HeJhl3YSSKrqlC8mnJm2OqqGFA1BW8sj7+Gi/EWcTCHR9OqtaiecLFng+UH0/LFh4H9dU1PeVgiVYNUoqaLdxOfB1vWeRKlIA1IM31kfKW75TwBkoJYJV55TzB1Y+4vpDWwLW9mUcbYu0uRJNcB6UAwuaw6bb1xqHfgGCLkUhU/vbBnBS0JOV14/Q1ZYtjh6b1qnAlC4xD6Q1FiDtF/SEAfdLw6pqzQETuRcRzt4SoZbw9II3AkvIWVpfyRoxOZSY7FrTtbS7NhBB2DbHGLDXdOkCVlxk/s0qBLeET7AElm9wHkLCS5MS8AOXjBWOELdh9HsL28g+7hMz6ud1rwLPF1kFCRTHTgYvO/+t/mgOWFXPTx+nBI9eaGf52JEx30wJI6AeWPCCYj2MJ0pVoQXjjARPaAoqP+Vgo3avevLydUu0In1MTtMrhVScsahbMq7Xz+3FpZEE3i5gXthcHbTKvqk4mooCo7VYz8fNc0zT0677IGbC84Ys7aGKTjqDxwuspM4XisOVQ9PlN+WB84TtUrDBivR5OcOCrv6pkgnFVY4ZMmj7KxtJzUWuU0VBVaHWkmekzbDUOVxj0NIKN1cHjXJmAus+LW1oc3VQrefMbXjP1Hs+XvDDasMLylnQH0X90RZyOpCvxNIVXoHcNmtva4oMMCV766QEgRw/z3JaW0oICo2Df1KwonE8wOQh6onK58txfE2rki1rTyW8vxyWjx3Bwd6+SvCxSa772SV3tL+ykS2J702Q3nqsTXwL3E9tBor25r0/yX7hLzpnQQD1EvmeP5P/TkX4xSSW8bJ+bm/SBJdzdSRZnqlExUvgCcg2Gp/VR2677fhIjVXJlEkdw4wBL5+Yj7PgadDaItJ/y73OcxH2LeETLFloyRs0KUGpI4+pyK8Fz3z16qXlgh3euYtG7tnpAqe8Lxd0SbXH55dOVB3k6Br95+syzYijnLnMGvbf8fPQa/LKsQKmE2lfqheM02Z9rj5liRC1PnKpUlUlkGwoK4OnTfktfa4gR/q3xtG2CIy0WXLTQ8qJ7h0jA9qjvjBSg1HalSKiVhulUV6rn971Svdp4Za7jhWhLem3nKMDC0A3bdbiQeqsHy9YZt1jSTXq1H3pNlIphRaNKNc/r5+pPntYlBkMYOqzRf2+2qy93CVAwpFlgHyCt29uPyFooZQ6obd5k9jjmek1F3QkMOXUTvkwStCWB5jqE+/3oJ+l3ldybTlWV/zwjtdmZKrCd6r6jC4YyvdjRbf1++lSoJd5iNKGwPOR5Z4LX6ukFJXgWgpByUnmwgmphaVS5bi6RotZEM4qBLeMT5D9dqlIsNz8mdEF80WzwLIExJnRhZrGxgNsNAw1rVB8hRpe9RCu4DEahtamLzu1S9fE0H0o2d62zdQ74UnQX9mockxT19Z+UGufhrwH2S+mpvxmYSgCTvuLdVUfL6hmlR+T7VqIyz0z2tS40/dq3b8lHHPLmeb8el3Bbcq7FkHmBUZ4v7U9h5RyMY97LFloaTeA9wD4yxivMfKuGOM7Qgj3A/gfAHx6cuiPxxh/Y3LOjwH4AQBfBvA/xRjPpK7xzO6bGtt00EN+ywPJ+TZymgx/0UupCbXjhk2hIBSENgIsd63cfaT2l5jNufalmIJ3PlAfzLr+m6a3zOoT9YIxjRxl4xzpg8UN1PCeC5tjMn68fFkREJJmZi2gVGL68rH6o+2hdMH0tjh74hT2DI+OaUcdtbg2Pl4uv9bGQmqDkoWWbgdwe4zxoyGEXQB+H8C3AfgOAF+IMf5jdfxdAN4L4BsB9AD8NoD/LMb4Ze8aKbK09bDGSwm287cJdGqVdwzg89m8YIuVAtbmZadygT3kggdtfXpAU8DoklZeepflC2wDa4Bb5q/uNxeWSJmxcqygS4Fe3Qe5ZklQoOQ6qcCA3qffG3ME2TTNBe1y6DrXOODmuQus44H0x4EF/FUxh2OMn4wxfnTy92UAHwdwR+KUbwXwvhjjMzHGPwawibFAnBsO9va1WlOEIQn3OgrJJp9w0YQgfenwSfRXxmuJyPkWBYTP0W2fGV1omJXyv0xSHVgRcLFQjcopPaHDpJ6Dtc87nqkuzM3jIgRWkc2212GkBCC3Iz8snDdXB+YE9a7Lpq12b1gTULZ7gixnAncpYtA2mCSmvnyMmKtZajZKxFsoXDLXqjlRmCVSfaCcNWHM/hua5dkTpzAahpqbiTEPc7iVTzCE8EoAXw9Uhf9/MITwByGEXwoh3DrZdgeAx+m0J5AWmhU84WJBBE1V462j9gHUyxxpYci+lZR5aU02EaBA3XwXv6T4deQFM2/QQuoedTWVlMBkQcJtyt+jweVaZoY2QUU4ap6YrgSsr+354FL3xM9Ct6ufkzVBq/WnVTCDPzhc+dqzEPjDyT5PEQyWX5D5dB5KBJ0Qh3W/cgJZOI0lPuY9p49WAlAEn3zoxPfXX9morZNsPm8qX8Y+Qvlb+w11dWluV3z0UhoMmFZdF8zD9VQcHQ4h3Azg3wP4mRjjB0MILwfwGQARwP+Mscn8/SGEEwDOxRh/ZXLeuwH8Rozx11R79wK4FwCuu/XWb/gbT9WjwzriJ3+XoNTnoM2pVDsl5iZHJsdFWMeDSmeZ1Au15tdT5vV0Z7lfvmfAj7RbkWLrfcj/qT61fXe5+7GiqTqC7JmsXhky3pa6D6aCcGQUSFNCuMpOG5OYTXjpm7YWvBRIPmdzddDwn7MGpStlV9uJO1rz7Ro54l1gLd0KoMrr1uXmeFlOFoRXhScYQrgewIcBnIkx/pyx/5UAPhxj/NpJUAQxxrdN9p0BcH+M8Zw+T5CiyFh+oNRA0pNkFieqHgDetXnQ5oIsenICfgl39vFYjuHUs0gJ7dw9WcEHoF5TjukoV4MAzjQX7fS3qshY96T7qtvS98Hvkfl/uoSYhrc9JzD4HfO6NIIuZds0NCdWC8PSZQzkQ1BaBizVjoXR4HIl8HM83oUXUAghBADvBvBxFoAhhNtjjJ+c/HsIwB9O/v4QgH8eQvg5jAMjrwLwYJtOMSG27YMdDUOt+EAb5IRKClZmiTURWYilhIcV2ZOJnJrobbXl1PENM24V1XVrdJSOz7sNWBvrrzXNSC5c6sHKRumhXqWbg2Y5X5a3v43/z9JaPRMvJQBzbIk2vrMS7W40uFyNB2BazKEErAF6GvSVQ/uxfKz+jufxEbBQEh3+JgBnAXwMY4oMAPw4gO8CsA9jc/hRAEdEKIYQfgLA9wP4CwD3xRj/deoaWhP0HOLzMgc9bS7lX/AisDxxRBCmosdAvaS+qPba3LVMvZT2Y0VN9fW9fdp0tTCPMkht4d2HZa5zTbycpszbtTC1FmCq7c9EO3UBiVLNSvpiMRIssKvFg7eGtS4BJ8JF0gk9yH17Cy55S3ZqszrlQrCK8WrXhcwfub9tkTZ3862744HPNQsflJh984C0L+o3DzD5O2UGA/ZiON6gYrZ7rsS+oMRks3x6OX+jJVx0zTuPbyb9956NR1fJmc9WP7TfFahnfeh355n7MpFLzTdL4GkzkPvK19UrrpVcy0Jp9kdO22NTW1dLl/2CK4emRU+18BNwNkxJLnFJnjELWPENekIQGL/vbVFZWgtBoBkYKSW5Wkhpf5bzWH9lrYkF2MJBBpgVyODzLK0ztV3alj5J6f0U2nw4vMT+VPSSNd+2gZnU8blUM625egvYWwGfnA/LGit6hTSN1OJCAq9YQEl0s0TbYzPe0m51dFxbFZbiof2gwjnkKLl1nyUCTmPzgbvMStNWUEtrs7NqglsmbU7D8lfxYGgzwb1jue5af2UDB3v7aprZntNHcbC3r+pLSugIR6tEgDSoHYkUKCaFav4hg+lFkkZWCubYic9RaCCe1ib0BS+h33tWnq+z1rZR41Dz99inKX1hzVuXEZPjmbrC+0qEs3dPpQJQriOUly70DqFsydgUGpFQeKRdaTsnAAF7nArFRd9DqoCDFoD6fG/xdmlXVzKS7bp/qTW/u2BLVJZ+0VNPA5nxsLQ2XqPVW9pQo8Ts8vhsJcfrgdRY71bKNq2NJ6alvWjoSWEl6pfch1cWK3e++NaAaaaIJ8ilyrZXATv14Un13xOmlXa9lqbpsHlnmcRcZh/wP1QuL48mcYkGyO0B0+fWBqxtL00ousPjuxqBCJ3O56XOtUmpswofALbg31wdAINpCTEWeOO+TSvPVOu8rAb0cLF6R7LyH4BqYXar//PEltAEn9l9U+0rqdFWq+lStmkW6GIJ3ksSMq7Vv1RalzVgOe95XthcrZNqSwZbziSfFZZG6F1T99/UVtQ9pVLnvEIZgpL1QuTcrgn/3nu2VsjzPqKs5ctPSkFIvVPR5qxiCdJHeS7yMxpcRv/4w9XfUrS3f/zhShjqeTEaXM4qMfOKFm8Zn+Ddr//h2jadK1qSU5uKoGrMyiH0KltoPxWbaLp/QN2HyLBUfn7pbLZ7vq1UZNLTpOS5MC/PQ+kX2RJkum852o8VwdYrqVlR4lS7ejtv0z6oNkVEtYbImjkH2zREOOkit/q+JTLKgTXLamByce4Y0TalbTZDLZ8vB/cEwnTQ26o1eYyF4lP9ElicRt62LaLDFlnaS6bPCcGSLIwS5Cg5lnamB46AB4Ymp3JaF2slKX+iFgheP636bpYvEfCd9Cmzqa1ZogMUXuQ3Jbz0PenjvI9l6TOzoCPKHL22IK4ES6CxX9capxwU0sKBaSFehLdL9RjrnKW180myclvklvP0lqzwFiaT/mxrIShoQ49JfbnaoOSanknLaXMWDvb2mZEuj7/lTeaSii2W9uRpyrrSiO6nhZSTnCF9EJqFri5imXZyHsMTgHyszrix+uhpmynoPulzNKXEOi6nXWvT3NMEhc6lhWibCkn6uhxZzi0HO6+AhF5hUMAKgy5ybCkZ24IiU1pZupTSIEhx1QTzMon5pYg2apm0MnAtQWgNLotr2IYHKNCD3DPJLOGd+yKXaoRa4LHWA9iTn+/FMos9Erf1AWB4ZjBgr0XNWhhgpzoCdjS25INcyh2tAjzORzBlLXnkY8EisjHmgRzpe9tSZDR4gvALTA0wbzAIrSBFPM4hpyXJMp3ewJIJItSXNgOwRPCJaSnHSqku6zpSmUOqdUgJJaFjsJCznik77y1qk+zTUWsWXjqa6GlDWrDpwIROn+MKQSlIu2dGFyrqVKV1T56LN26kzJQ8J135RigxXokuDmoIJUsqEHElIn6uYnKfGV3IFobg5yJ/W+9xEZWoZ21fn8NlteaF540QFPBg1jwwC10jxSmB2OC1GZpQyvnNv3MZIylNJuWzBPyPAJcl0pof++OW1s43fJne80zRVvRxJbC0wdS5noBJReo50iy/pXSTRirqLOWmuH9cr7DE2hCen9U+MH5Xco+y37u3kuvNyqDoItC6CC9RJFKLoM2K5405vEiUmsfW4LQiu2IO6wog4qts63D2TOJSPxb3ywJnV2iw6Z4yiXV5K6uiDmtzpbCW77T8czqgAjQpIiUmsiXI2wRTtNmeq2qk/XhWhgTQzM5hFoK+R6+vui/SLjB/f9+8IfduEb+3fWBkUfDoELOYyBrW+rQygDm1Lpea1jXarWkVDE1dAJoBpYbG45BwvYikVeIKsOsY5gRNToB5NCO+p1SgoSS4I23oSLEXLClBTnvWY8MLDgDNZ1CiOafoJ1sBLJy11ST3O2tgZEtkjLRBaVAgB226VO0O7cmQ4nd5xzPkhWlNhek2Hlfw0uGTOLhiF2PQE08/Hxnsy+t1jqHun1w7JWzkOVlRbI+y4VFvqsyAQdOxnyObW/uX1s7Xsld6G7swWrmM/potnFKCMDWu+isbGAHoY6PqD+/LQVN4+N2nPgj8QRFBqD+ylvbN19Uf+K0k7FKQfjLvsS0NKIUXnCZomQQAagICaC/0BGwKs4al/weaPDCe4CnStzWpOYqpSahyf3o7awGe+WpdQ/axKcfQApGFZ4p4y/frmXsWSl0DXhEOLTRy2ij3mZ+bR71JCdfctbXQ1FVggDp9RBZF11QdOVeb2paLo1Q4Lqq+X+p61tx5wVBk5gXLlNAEVr29VBhq81LD82lYWSNW+hyQnqBeFRyBV6FXV+/I0ZG8Sa21RLkWV7xJlSerorHOIt6l5mZKm/Q0yRL6laxyqDVind1knevRcDzTlrV46wMg79viklrXZ8xDCPKH/moIQitDRubNC4YiMy9oikNJrbbcMV6UlQf22ROnshPYU/GtybS5OqhRSvj6Vn9kEOkBOxqGBs3EM6+tPFYdZfWEj1BNdFt8btWnwWWsn9ub1JKs82S7J+T6KxtJTqPXnoAj5JZ5nUIqmu+hhPDsHWNpuJ7AiSIUdAAAIABJREFU1SiJvsr4WpQA9PogdC5v9bku2HJCcJZAREnbstJbiZYn+/Qxcq4nbPSLkS8YlzkSXh5rgd4L1Zw/+cqn1mjVfU4NVksj08UDRLgI303O8ap8SFEJzjwZDcclry4dPlmk3aYCBSmhUiJs5VxuRwRH6RhMmehtNFbLF8w0JeYZajeB7m9OkOv9nrCxtgtvlPfp/+cBaU8Ha/ha4hech29wSwZGcqp8aZYHm2AAsOf0dADouoEeRO3XwQRPO/RS9liL0uaaJfwaQuswsGc47aekTXl97q3HWjBEo7QYqhe91u+BMzeqiapKMAkNZBlH0EM0I6spAVTCbatydrG/EsQS7CrVxGS9DH289bzaZhzJfVaBEdSL2crHhgWhhdEwmO/f4lXy+h9MJcq5b6xrAs3sGWlj3lohl/7n3/x+AQAfnO0629InyC/Ygvb7pYRgaWlzhq4cA/i+KIHlE+SXr/ueqzbMfj4r8KDh8dqAut/I4rwxbQSor0/M0ALVEiolEVotgDX1xRKsukCFhZJrWx/o1IdZ98dr32ujRCtN+VH1NTRS3E/ZD9QtHE94Lso0zvnkF06RCSHcCGAdwA2T4z8QY/zJEMJfBfA+AC8B8FEAfy/G+KUQwg0A3gPgGwD8OYDDMcZHu3awFHXBd8E9zgp8MLoIPasvbGIBU3oH0BwslgmdQqqPB3v70N+oBzpKJqhoC146lY4q6tSwdWeySvu9jWkRUJl4OW3HQ6kZqt9BblW8Uk0xF9G1nmuu/VKN0vqYbq4O0MP4w5MqrMFVbXShBC/dDoCpjYl/ed6mMAc+qnECmwPaZo2eFEpWmwsAbooxfmGy/vDvAvhhAD8C4IMxxveFEE4C+A8xxl8MIfyPAL4uxng0hPCdAA7FGA+nrjGLJijCL7fkoMDz780TTG7WPhwPHodP02z4Gvq4WWgkXRfUtvxsWkiYpnIi00ZrpG37wdCUltJnlNLa2MzPRXBnIVJb7aVgrZ9SUkQEsDmgAiurxCqqkdIEu5jL/PG1SN3AfCgy2cBIHOMLk3+vn/xEAK8H8IHJ9n8G4Nsmf3/r5H9M9r9hIkhdPHfLTcUdFkewDnCUCrNUUKMr2JHLL0kGVYnzlvN5PTNDBLjlj2Eitgg0jiRraIe6VP1dP7e3dr4VKMht17AqxWiKkhVdt/psQQus1DOXya8DP6Kp8z2wAONray07J6yt9yBruPD/coy5EBF9PKz7669sNNwRbaCJ3wwr+GG5auYdJDl74lSj6IcWpPOYw0XR4RDCdSGECwA+BeC3AFwC8FSM8S8mhzwB4I7J33cAeBwAJvs/B+ClRpv3hhAeCiE89OwzX3CrbDB0oENjz+mjpno8D5U5BQ7XS+K8/vK3+QpaL50HPi+65A06S0PT0Kaj+PZSx+XgRW1FkPQ2djX2y3vLafOp+9D+QKvSi7Wd2z174pQrCLoIFt3f4YGLtfei2+RnoxclYoqRV53IM4Et6AwUIE830de1+nC1s1DmIXSLhGCM8csxxn0AXgHgGwF8jXXY5Lf1JBsjK8b4rhjjPTHGe66/4eak30o0v5z2xvtTvr95fD2sh8+VPpgUvbR2vvbVTFV34a+p/q1NAl5tDKDI80RTsfxj+n9rUnrmsKelcJu6PWmTq6zwZBcBxu+ONSdP69T9YlOOnwe3x8/dmvCWb7biP3ZYeJ41Sr5fHcRIPVMtDC2zsCQTRUOX1gJ8l02qJFwptnKKXiueYIzxKQAfATAAcEsIQQIrrwAq/+8TAHYDwGT/iwF8Nte2Z6Z0NVsX4esD0tU2eBLx/s3V6UI7ucHgJbQLJUIEn2RhsICQ64hAKI0sCmTxHGmTJ29qovJ98m9pU1/DWvJSIB8v1vq1cNTn6swNT3im/JfyMZF7Tn1ASmBp1SL02J+on1VqwXJg+qEUH/GlwydNU71Ecy0di23asTh+88AiCzyUBEZuA/BsjPGpEMJXAvhNAA8AeBOAX6PAyB/EGH8hhHAMwGspMPLfxBi/I3WNm2/dHW+7/z4AixNeVwtCj9E0HU4z8iJx7ATWaUJAfeKWpHmljvVgCZrUWholEO2HHfbWdUSTtaJ+Ot2ulHrE9+CttZKi6ejnJ8LNaouDEUz/sZaszC1czs9KlxMTsNXDz5Wr3DCs8lxAuzL3+ngGn+ulpc4LHDBceO5wCOHrMA50XIex5vj+GOM/CiH8NUwpMv8PgL8bY3xmQqn5PwF8PcYa4HfGGP8odY0b7twdv/jgjV3v4aqgJLpl+R69geqRkIHmYNTmszVRukJPUD3JUyuqWWBzj4/liWmlb2lBo9vUQkvX4WNYgktfQ/fBEnhee7wsgFfe3zuXoSPyOW2bBaFmRPAHwuIMeksAlNJcvKKmliDUH7FZhaA197gP26Ke4D133xgfPLP7WndjoeAiAimyLNMBWKu0tBEv59e7jrWmhxzLiyAB9mJG+lq5fVqDtASn3IdcPyU0LLpRipqiNTr+26q9V6I1e8EjNnFZsFpaIFDn85W4GvR1pDAvgz+0zAf04NFOSo7lAiApjfFq5BbvFFC4SpjVx+EJpsrxPqG4LK2dx57TR5OOb6FD6MlY4gOU83VgItVn6/xcoEKDNSddip7vw2qbn5123nt+P84x1vuZSpPqswXLz8e0IhHKlqDUws577kxX4m05sH9QjhdfZ8plUBL46DL+5yUA5+1f1HjBC8HSBzyPF5oSNloYCbu/S1sCzUWzJpKe/F4FF+sc3WeOBFtCgLfJuhx8XdmuTU+mllgEa+8aQNMPxjxFCYRY1XHkXKZuWe9D34duR7TbVMBDvxdNpbHuS7teJGBiCXaLA9h1PIvAFMqWpswsQmAtmorzghOCbaqraKQqbmjtzYKn4fCg598cqdUCzaNfSHtA02fIfDreVpmZ6hoW9LUEltBbP7e3+hGzl38YQgURAexFggWsvfAEl3sWYcmmqdwrt6PpIbyPV5fztFsd8WVNTO6pv7JRRd4teB+0NtxEfh48rjT4I8IVWUq5tDl/oD5mkRDWxKx4QfgES1LquvovSh3AupBozjy2fIACESCezy/lL/TgaR7eusCp/mlYprsFvaCS1S+rAG4JT86qst1fqS96laqqU+K3S/ErrffchYANTAW95Rfk9EkrACfUGmuNGSmAmwrY5ZCikM0Cr909p48uPm1uO6AkOuW9NK3hyf8l1Vw0UoLI0uxYa5EJJhPA8p15AlBfh7VLwDbpckEXbs8ya6WNkj7pvrHGxs+MJ20qU0aI9Xp5U9aWmNgOpLMlUhFyvi8mNnManPXeRTucBdYz0MEeLdBkDWX5rZEKpFwtDc+CF8SZRz3BF4QQ9OBNJs7c6K9sVBkfMnis/MVSYeiZm55WwBOIJ5w2n/UxqaCJ9rnptvkYizisOWnedXTGSA5yjJiQHLGVwc4VTKQ/lrtgz+mjNeJ6iuem+Yf8bFPuAcv3qc1j/ax1Bogg95GwFqq3gmjynK4c2t/Iu7WQKodlFXxtS6eZFxalYQIvQCFo+S20pse+IIHHk2oLz3RjX5hMahYCHDzxAhNWoKIkSMHny/acT86j6cwCj8wM2P4nFtRyrnyMeKU+Das+o3VNoJ1fzgsoWYGRNrSY1PHWR1mXyUpBU4y0QG2raS0yOLKoAMnzRgimghKlEAHn+f94YgnEvAKaL6LUmawnmo7gykQTIaYFkea9ybHarGXhaEUcOdLKTn3uC8MiUcs1LQ3PMu/YjNd90kESj6O4uTpofLRkVTXZL8eKxs5FJmS7Bx0d5veRCmjoffyRYo3PE6T8vNsKW2ASPDOWcrAyQSSIYC3+DkyXPpgXrpbpPI/rPG8DI11yCblWnXcel9MHZmO7s/DT5Fwr0AHYKWYMLyjhrcxW48U5ZGkmLFv90deywOdwezpbhAVdKgAi24VEzgEMq6o0+wu9D5zAe+/W9Stz3ElzA5rCu62W58HLIpFnJ/20AiQa/Nwsv5+11GsKJWb21SqacN3tm9s/Y4S/aHpNX4ZV4+xqwyvdr8v4s7DjbZYwYIGWEnjs3+MMDC8NjieTB84k4fYsaEGn/ZN68gJTIVfTjlUkV2BVeRHklhtgWKb1aBjcjwA/Q6v47LyEHqOkmIJgczW9nITcv2Qk5WglJYKwja9x0fNxViG4JRda4vxDzg1dxhFcOuEXHb2aArB0QtZoCMP6cZrqwgLC4g32VyZl7Fdtao2UWO+vbADihMfF2vkWLN+eJjavn9uL4cZFYDBuW6MSvKhTYrgd7oNFYtY8QCtwpSd7ziXBAkBrTNW5Q31WswzZqHnITDSXHFKVvvU46a9sYA98Glg1J04Ay8fSNTstH6OleCwqRa5NO/O65pbyCXKZoIO9fVVktqKJGH6LXJHHRYAjj7njhE4j2QlAfYLpSSS+O+2fAqZCxMoosCK+er9oc1ZKlrTLgRSP1Gz59riPeptczxKAkhkj/kq5tjXArfebisyzAGwT5U1RmbTZb+1P+Ubbwgo86fEj92Z9EFJ0Ig9W0FDMaXZF5NrtUoewzfHzCsJsCSH4yGO31QasHoSVs9lQ0682d6ltlWqp2MG147RJKELLmnysLVrmsm6zTb6pjvbyeZrgzCa1Jbi1NulFpTW0Rmj12ZpMpe+Blx1gdBk3Fp/SQs7FUIJUQAbwP6SzzgerIAXTbtoSqRcxP+fd5pYxh6v0JmNFMBlUukTPPIiSKaRMXg/m2sSTZabEvOeJuXlu+jWvRXVXx78q36ExscTc1UEIgZhxm6uDsSmLqcAVM1mgtT4xvU2z2CBW91c2gNVpAKjWFpp50bwA0PqwaVbmNIKSDCBen1Yfzz5AjZTrIAcdLfeqyOSgg0Zef3IBJu0qEu1JL+8q/6fI0pJRkwosMhaloLAGePbEKVy3HdYdvvnW3fHA59JfTysNqjSNp61KrssFdUFqHRSGjhJrbU8763PR2bZUC25fC9JUIMSL+nqk65KiqqxNdonK87v2ovvCEPAIy3xPjEUEP3Kw0ghL4KXTCbj6eGqMc61CXi4iVT7Luh4wm6sqN4e3RXS4ZMlNHdUqnSSlS3HKsRrWmiVtYF3b0jB58ll8QEvD6AKLyqKhCzBoAef1KYVcpFOTr6281rbgc/nDVvqhsHiSVwtefngpNGGaaw3K/mSq4ERAalaDF7n3UPL+LCHXRnnZdkJQ+8vki6ULjs4LOZNXJiMXRS1FqTYo0FqIXE9ra4JSYWidBzQrJOtzuF+pNnnCspnb1hTk5zvPkuwiANoKspKS+IsGvyvrWVrzRXMH+f51cEo0O2sOWNt1ReqccGMsKuo7qxDcEoERoF4EU+dwli4JyDjY21dVZrawfOwIDvb2ZQWgHCd9y6H0y235MzkTQAcqAHuJRq8Nb7vOxuCgBvsFPcEn4EAJp+ulVg1MgQXUvAWgbv/5Bg5+6QXJcrncwPR5ar8skPar82JOvI0FYJfoswdua57t5pAVgiGEG0MID4YQ/kMI4WII4acm2385hPDHIYQLk599k+0hhPBPQgibIYQ/CCH89dw1ntl9U5bGALQzDazkd5kQekUzC2IG6EnNyfx6e2n/hPqjwRFVTqGbV44uC9kUZ5D76UH8hZZZq53rbTXoea4Tzc+wLXR6YwlmrQxjwaJW8cdba4MWcgGPtmhLTWtLfWlT+n9WlGiCzwB4fYzxbgD7AHxzCEFGx1tijPsmP6J/fwuAV01+7gXwi7kL3PD40wDs0j9AfTCWPhT25emikd7XrzTX1KLxaKS+sCKgPYEqx2hOXk4Qeql2mmvIv/kYryKMjhxbFVKA8URjAa+jzfrdMuVCX5vLlZXCKnk2C7iqDeBzAK8G5KNjramsLSegbJ54WjuQ5l/ye5kXQZq5t1oL1BrivJEVgnGML0z+vX7yk3IkfiuA90zO28B4feLbc9dhQrSnqUjCeO5BcNkrgbxUq+yQ3nfp8MkqAd/K3MhBBIE3iVPr6ALlHDsGFyGY90TVEWxL27M+XnrtYy726WkmHIXU4BJnFmRsyP7+8Ydb+/FY4Mu9pp6n5vPljm8LXRxD/uYPqSA1ZjxhpQWhJfxkLvF8Sn2gOBiiBVtq7uYyXhalFRb5BEMI14UQLvz/7Z1/qBzXdcc/h7ixo9itHFtp/SIFJ36m8Q8apXWjJ1IJx6S1Y9JSQUEJpWlDwDb4D/+IU8UU+uOPUAylLxRELYNLUyiJmzaixrS4xk54glpKnUZJ/ENunkA4qtwqaR2nVNDE8e0fe8++s2fPnbm7b19232q+sOzuzJ2Ze2fuPXN+H+As8ERKSe/cp7PIuywiF+ZtbwO+bQ4/nbf5c94mIs+IyDPnLu1t04daIoCK0hssskRF54qKz+ixEWfmQ9eazm2PadJt6T7rwxZxWzb8rGmS2/0RBxcdH0WE+HPa3/psSpEfem3LlSyspAFOxep5bVYT5bL0mfjFqJzC4oHnWTzwPDcv7AxrE08i3XpJJRA9b5stxhaa3yiO0UbWDPmEZoIcvWAGXv5BTsMmzk/FU2UMFlZSGEan3+Nwh/58TdlwJo0qIphS+lFKaSewHXiviFwP3A+8C/hFerWHD+Tm0SwcoiwppYdSSjeklG7Y8kq7vs8vQCvmWnHXikFtSnr74G2qIfswIqWyPY8VS0pJAEo4uf/BXtynSxO/8vR1RQ7GLi6bpkphuRi7rcaPsEQUS9ELvoi8z2ps63N4Xacdr/6OFPsqKvlFr6GFXgT2/RgHeh/sM7TERe9HZD3eKAuy6jZ9P2FtPmsi1SgZqrXY2/7WIlIveL2gelpEJSeaiJfu87HIlqBuZEjsyC4yIvIHwP+mlP7EbLsRuC+l9CEROQR8OaX0ubzvReDGlNLLpXNefs3l6T0vvn+c/jeKp6qH8m+Ukue7rfMLwz6C3go3FFebXQq0XU2KIz2/GmLseZsIVuTIHDkql5ySo2OiLCp6nhpECSCiRA+lbdHxbURbn2PkdzkOojRittA7DL5cLRfoicrqA9cO5D2cNKLi6tZZ3Do1W08JO2/0HrbN1SYOzxJHy8nXwidMUfjIkFJfNtxFRkS2icjW/PtNwAeAE6rnExEBfh14Nh/yKPDRbCVeAl5tIoAAPzjxepWurRQ3Wwv/EEuFuxX6NlIRwHI9pRhXq0eshZ0wkXGjxjJsRbH1wF4/EnsjRONtU2n4baU45qZ+llAigDUiqj/WV6NrwkYaTWrObee37699hlbiiJJxtJ074gj105aNupQCrw0bmZqrJnb4CuCzIvIGekTzb1JKj4nIUyKyjZ74exxQHvgfgFuBVeAc8LHazvjklOBy5S2VF7n3htdjo5hRe+4IfQ6QNbbevpVsei915rYP5uT+B9mzcjt77hzduduKhTpJQzFo2fyG/n+LtsntI0AiLqjP3R4uu9ZoCi1/Xf9MSkTaOuWWQtqiawJw9BJWH1g7RxMHuLq8xLkHUiN3FmXE8RxgaXxD7lRjcIElB/jV5SVYGrw39nlYMVQNexZHDh6Cg4McoVVRNKXisoiIWJvOrrba48mDvXUziiGk1+aTre2aMHMRIx5NPm2RtzysKVlLyTM9rNirBEEJnY2xLJV79OLzeqCRKX4xRwk9R+U422CjCnRRR8YPhY0HtiJtFPgf/W7q/zjW8dpjorb2fpZCGPfufm5AxPVx1jBeVEn0bO1/y53780f9bJuHOnd9FFGt+kZRG9LYJkpbEdpzfDVi+D9/8ZPzEzbnDQtNdWChnZPzaEuZr5lH2q7hw4l0Eo5CBG1fvE6lRhc4KdQ42Db5k0HMGbURRd8GBp/5pA0MtURSDSKWi/Nhhm3nsUSsKbQxIpw+Xlu3KaJQRYhDDiMOLMqxqES/NimJwuodm4iV109G+22/mzjA6BxzETb3+tY3h2FAipKPVAlXPXIHNy/s7FsWvbO0/w+DBNCe3/enVJSnCdFDVT2jNdaUUrzDWmLSjSSAXj+n1t2SA7vnkvy5Sy5FJaOJfelNMkJmVHgx1hdNspbhyFpfS8C9NVkJ0bl9uwYS36pKxN+vNj1lU6U9VU9YF6xRrOqRXrBEvCLxvAmltjZMb5LO0zOXTxDWdHk6QNU5RXrBElGoIRYDhHAvsHcwe0ukTFbovpWnr2NhX1YG74+vE5UKANjCoF5x8cDzfU5gkaNg6npEOr9xYImsFWf9vbVO102coPVJXGTwHtUYaqzhoY3jHBXezzFClDxWt5c40pKbSVtbf11/PWDACj3UfrmXB/LM0v+wsC/19XwANy/0vleXl4Y4QMuleaLhOdWS/s5zelsOH4O9Sz+W0hbeSjzpa80EJ1jCkYOHBjzwozf7OBZRy93ZN6F9M64uLw0oejXEy4Z6WUtYTRJQT+ibdIlqqBiV8ystdr+9FEpn9ZE+WL+EkoN7U3uLSRLAtozMHk2c/DSyxlhEXHRJHaLPTh2ZLaK5eeTgoYEXOYwecTIJb4QabCSBhRkigkpMVISyLG8k/nqR1SMiLLpY7bH6xrfiTU3RaY14qBUh7BvYjuXxM8eHRPP1iLyjuo54xb49h9b/8OFnkeEmcnVpg/rfle7huOLsJFQGG51xphQB5F8IPulFkx68VqdmCUqT1T6CnrvtxViLUvhd07qy/Z9EhNBMiMM/zOvLPgxbB3XgLdBPVb/mJe9vhIrSCwzGoTaJzoPGjZ6VrO9OsG+w/Zm90hedvdjsobrGM8tSNKCc3P8gV3EHZ4ylzyezbEKbpbUG9vjFe46uqR8MfEKHPgElznQSGUcsVKSKMO6YbIW9WUIbt9WkB7dQI8aWw8f6ojDECWSj7X6/uj9FiDgw/b+mKz4+dP5JoMYvc1LZhmbCOnzh23ekhU/c3V/4NVbWqF1kZbJ1SfzC8n5eUWRJKUqg5GlvJ5hyUrYORI0bgi8Ar2NqeuiRlVDRFnXho19qrLNRGn7bj1FRcneKfCQ3Y17AGqOJdXfxvqgeUcRTSRcYET8YthBrHzzjUYra8G5jtUTQri9FW8VAD3vsqbvv2/wuMm3F1y28ktf6GHkiOEpqfX8NnYA2A6939VCC5gmfx6juM4r1vulKvmURSsW+mxT73s8s4vhscfUa3Z93qSnFy44ytlmB7Xub20yJ8Pl7HGVbHzWNlr6o7Qs7Sivn0/XDMPFqI4RtmZXaRPno+PUSwZnRCY6KyPG5lAmjJo2PbRst1FI4XF/cNRlRFG2uPKMiCqIvtbN9GBelzNWldP3+eqvLa/nvSmUbff9KPoUK1aU1FSefNtr646NzFD4hgkf0PCK3kbaEpzaZgiWAem1LkGpdUGrbqWuYT3On59A1Ha1Bf41JVZucCSJ44qVt1W2jh9rEZY3CTZV8BZvaD5TPLBxj/QFrEfW7tHjGRVSEva19G7xF0+tL16MzKtX8mHTkzHrQdI+iKBXFKM/Vjrc0r6wBRL8jI0Kpv6PO1yY/wVLfrK+s11P6ZxoZ0CZhFIEZIYI1KKXxqX0btOUis2+hc/t29V0NfEywfTiaYbqULMDuq9WVNOU1tPC+bSVro1q+ozRZNmOytZZHC8O2iWqdaGqrCJYQ+ufYZFiy97opUmZaOsKIcC0cvaTxZdWW0zHigHUe9X1IjW5b56zPv2edifXlvnjP0TXVhDNyeVgCGuX20/kZcYujElBL4AbSl+XtpWJPk8JM6AQvvnRH+v4LF7W2i/KUNen9otC0Nn2FLsqoYLdCJ5VOWK839It1I/WBVqFdqkVikzJYRAkTYFgXZ4liJIbaEDFr4Ghzo7B61+jeaXRGifjOKqyBI+pzdN9LaoFon98fGfVgWIdn9ag6X6Jn6HXd0Zqxur8m5+zouKh/FhpYUFMBb2ElrTt2eCY4wXe9/TtD27yZX1NaedS6ktRwYjpxS8p7TSxp968uL/XfVgsricfPHB/Qe1h9RxSuZ8e4XlhOQeFT3LfB1xhWDlC5B4UPF/N1R2rEYO94PtCPlTRgSCgR5Y1OXxU5X5e2W1jOuoS2uHivXvGEUT+aCNgSLV+Gom91zs/Q+sZ6H1A9V5P0ZJMPR4kPSjhy8FA/HLO0FppKfW4EZoIInnhpW1HctQ8jski1ERXbtjZFTySS7Lnz9r44AWsJU5XQFdMXZUSp/KOJY8dTU2skEqdKIrtfkP78JfeaJiLkoS+BNlg9lS4Khd5L69heIjw/bn3g6gPXtrbxUUjRfYdBQhZF9ERp/u13TaSOxtDrMW26R6u+Ka2TPXfe3n/5K0bR9a4ul2vM6Ll0DikzERHFSRHGmRCHrYtMxEp7fz8YLsjj2fEmcday25Ffmq/T4FOT6/6ma3zspT19n6vIiu3Hp9DJZ4P4RzWGlLgx7UdJLK5NC1Vyhym5VkTQZ3zzws6+hdK7e9iUT6XC4xuRcQbK6atqC7J7V56obYk4lu5vDWzwgOXQ1Lk5knL8+CK3GwsVrUfNRu39b2tcaWoI3VyIw4oSAYxgWfVSiFwJ6q8WKaD9JC/5/dkqWh5XPXJHmCFar+2P88pnRUT8avRiStD9OScZ6xkt7pqYVYWKXGpEiBZeG0e5kbGrJQLoXXza0MSp+v5bo9QkooCueuSOgdj3mj4patQzUd2RpnKcba47vi54SS+ox48SttqEmSGCTZkurP7AB4iXFkrTW6otZ6DqWYZ0f9nKavugGX1tn3VftFi824KHVlWLOJ82cWbxnqP9vtnx25RYduH5RafwKZ7auB9rvCipNSLYl1ibOKXn1761cVXrRXQN5aT0Xo5CDM/t2zXUx8gYMk4Mtr3GkYODKe49Y3Fu367hbOtOzVCTzMIaAyP1Tgl67ii1ner3tf86Fr2OfltL8iRE4pkQhy++dEd69013DaXMgTJLbduMgsiwUQqnK4lnHlFMsB1HpMssWdBUlI9C2ZqshD6Eylvv1Fprq79FccC1sJxuTURMjXXetvV9jwivoklsn5So7KvPQXsCXD2uJBJHGaWb0CYiR+oPb7Tw+2GQ67VzoqS4fzYGAAAG3ElEQVR7924sPkxVfyvsM9fM6RbK8UUETbfbb9u3PXfevm5xuDqBQq4x8gzw77mq3DuAz9Mrt/mvwG+llH6Q6w//FfALwH8B+1NKp2qu0aSIhfXXHbWuLVAu/6j7a8Vs9Su0+QQtQW8KLbIcrh2fVp5rCqOy19drKgasg3vXchZ6JXx/PJXhcf64iABGBN5GNvi+th3bhkkSu+L5960tUH2J7j3Qfs1GtcwYfW4ihD4BbhMBjM67wHMDoXmluiNK+DxT4JkY7+EBvRyaQ8Rsf6wXt9fzhNLWIlkvRskicxfwAvCT+f8DwHJK6fMi8iDwceDP8/crKaVFEflwbldINzqMaLHowhhngcCgQtZnzIgJXrPC1xLRJi6o1NdoLEq0BpzAdw9fK+pLpGTWyWUn6pbDx3rjN1mSV1TMLGReWbznaEggm14UTeNuQ9RGF2gJJWKyurzE6tP1xZtK0Hts/3vVRIkj9LHQdhzeAdwnpfWo0YEqIYzuY8nhuX9u4z7T379/2JCo6d/0eCWWbUyMZl/ypWx1X5Sq3wYx2Hk+CZcyRZU4LCLbgc8CnwbuBX4V+A7wMyml10RkN/CHKaWbReTx/PtpEbkA+A9gW2q4kGaRiSyzHm2Om/rfwou/ilHEOVjjzhQrT18XFnXyD9FzfXpMZPaPRIfICuvR5NhtHZHHcTguOUXXVifzE3fUF5m/775f2jebBMBivVyijteLgf3zr5Qr2EUqjZK1O8I4VuIaXaXntkp6wOj6SgR1Hqt0VUqCYKUvz0HqnPc1h0vwluVJiMO1RPBvgT8GLgHuA34HOJpSWsz7dwD/mFK6XkSeBW5JKZ3O+04Cu1JK33XnvA24Lf+9nrW6xfOGy4HvtrbafJjXccH8jm1ex/WzKaWx33St4rCIfAg4m1L6qojcqJuDpqli39qGlB4CHsrXeGY9lHyWMa9jm9dxwfyObZ7HtZ7ja3SC7wN+TURuBS6ipxP8DLBVRC5IKb0GbId+/e3TwA7gdBaHfwr47/V0skOHDh02Cq1+giml+1NK21NKVwIfBp5KKf0m8CXgN3Kz3wb+Pv9+NP8n73+qSR/YoUOHDtPEepylDwD3isgqcBnwcN7+MHBZ3n4v8KmKcz20jn7MOuZ1bPM6LpjfsXXjCjATztIdOnToMC3MTNhchw4dOkwDHRHs0KHDeY2pE0ERuUVEXhSRVRGp0R/OFETkL0TkbPaP1G1vEZEnRORb+fvSvF1E5M/yWL8hIj8/vZ43Q0R2iMiXROQFEXlORO7K2zf12ETkIhH5ioh8PY/rj/L2d4jIsTyuR0TkjXn7hfn/at5/5TT73wYReYOIfE1EHsv/52Vcp0TkmyJyXF1iJjUXp0oEczzyQeCDwLXAR0SkPWvlbOEvgVvctk8BT6aUrgaeZM049EHg6vy5jV6Y4aziNeATKaVrgCXgzvxsNvvY/g+4KaX0bmAncIuILLEWBno18Aq98E8wYaDAcm43y9DwVsW8jAvg/SmlncbXcTJzMaU0tQ+wG3jc/L8fuH+afRpzHFcCz5r/LwJX5N9XAC/m34eAj0TtZv1DzwXql+dpbMAWesk/dtGLpLggb+/PS+BxYHf+fUFuJ9Pue2E82zMxuAl4jF7gwqYfV+7jKeByt20ic3Ha4vDbgG+b/6fzts2On04pvQyQv9+at2/K8WZR6T3AMeZgbFlkPA6cBZ4ATgLfSz3Hfxjse39cef+r9FzCZhGfAX4XeD3/v4z5GBf0os7+SUS+mkNuYUJzcZQsMhuBqhC7OcKmG6+IXAz8HXB3Sun7IsVMvptmbCmlHwE7RWQrcBi4JmqWvzfFuDYqvHWG8L6U0hkReSvwhIicaGg70timzQlqiJ3Cht9tZvyniFwBkL/P5u2barwi8hP0COBfp5S+mDfPxdgAUkrfA75MT+e5NYd5QhwGyoyHgWp46yl6eT5vwoS35jabcVwApJTO5O+z9F5c72VCc3HaRPBfgKuzBeuN9MLyHp1ynyYBGzroQwo/mq1XS8Crys7PGqTH8j0MvJBS+lOza1OPTUS2ZQ4QEXkT8AF6hoRNHQaa5ji8VUTeLCKX6G/gV+hlnZrMXJwBheetwL/R08v83rT7M0b/Pwe8DPyQ3hvo4/R0K08C38rfb8lthZ41/CTwTeCGafe/YVy/RE+E+Aa9LLPH87Pa1GMDfg74Wh7Xs8Dv5+3vBL4CrAJfAC7M2y/K/1fz/ndOewwVY7wReGxexpXH8PX8eU7pxKTmYhc216FDh/Ma0xaHO3To0GGq6Ihghw4dzmt0RLBDhw7nNToi2KFDh/MaHRHs0KHDeY2OCHbo0OG8RkcEO3TocF7j/wH1YwCYJQStpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.title('Classificação')\n",
    "plt.imshow(map_class)\n",
    "plt.colorbar\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Criando a matriz de confusão\n",
    "\n",
    "**Tabela 1** - Matriz de Confusão\n",
    "<img src=\"_trabalho_final_intro/figuras_apresentacao/matriz_confusao.png\", width=450, height=450>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################### CRIANDO A MATRIZ DE CONFUSÃO #####################################\n",
    "\n",
    "# construção da matriz de confusão\n",
    "\n",
    "matrix_confu = np.zeros((nclass, nclass))\n",
    "cont = 0\n",
    "\n",
    "# matriz transposta da do envi, saída do programa\n",
    "\n",
    "for j in range(nclass):\n",
    "    lim_pxcls = int(pixel_classe[0, j])\n",
    "    for k in range(lim_pxcls):\n",
    "        pix_xy = world2Pixel(trans, matriz_shape_tst[cont, 2], matriz_shape_tst[cont, 3])  #chama as coord. xy da imagem\n",
    "        value_class = int(map_class[pix_xy[0], pix_xy[1]])  #coordenada pix_xy[0] = x e pix_xy[1] = y        \n",
    "        matrix_confu[j][value_class - 1] = matrix_confu[j][value_class - 1] + 1\n",
    "        if cont < nlin_tr - 1:\n",
    "            cont += 1\n",
    "\n",
    "# matriz de acordo com o envi.\n",
    "\n",
    "matrix_confu_inv = np.zeros((nclass, nclass))\n",
    "\n",
    "for j in range(nclass):\n",
    "    for k in range(nclass):\n",
    "        matrix_confu_inv[j][k] = matrix_confu[k][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculando o Índice de kappa e acurácias\n",
    "\n",
    "Índice de kappa é dado por: $\\kappa = \\frac{P_{o} - P_{c}}{1 - P_{c}}$\n",
    "<br>\n",
    "<br>\n",
    "Em que: $P_{o} = \\frac{\\sum_{i=1}^{\\kappa}n_{ii}}{N}$ e $P_{c} = \\frac{\\sum_{i=1}^{\\kappa}n_{i+}n_{+i}}{N^2}$. \n",
    "<br>\n",
    "<br>\n",
    "$P_{o}$ é conhecido como Acurácia global.\n",
    "<br>\n",
    "<br>\n",
    "Exatidão Global: $P_{o} = \\frac{\\sum_{i=1}^{\\kappa}n_{ii}}{N} * 100$\n",
    "<br>\n",
    "<br>\n",
    "Acurácia do usuário: $au_{i} = \\frac{n_{ii}}{n_{i+}}$\n",
    "<br>\n",
    "<br>\n",
    "Acurácia do produtor: $ap_{i} = \\frac{n_{ii}}{n_{+i}}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################### CALCULANDO O ÍNDICE DE KAPPA E ACURÁCIAS #####################################\n",
    "\n",
    "# construção do índice kappa\n",
    "soma = 0  # valor total de pixel da amostras de teste\n",
    "cont1 = 0\n",
    "\n",
    "vet_col = np.zeros((1, nclass))\n",
    "vet_lin = np.zeros((1, nclass))\n",
    "sum_diag = 0\n",
    "diag = np.zeros((1, nclass))\n",
    "\n",
    "for k in range(nclass):\n",
    "    soma = soma + pixel_classe[0, k]\n",
    "\n",
    "# somando as linhas da matriz\n",
    "for i in range(nclass):\n",
    "    for j in range(nclass):\n",
    "        cont1 = cont1 + matrix_confu[i][j]\n",
    "        if i == j:\n",
    "            sum_diag = sum_diag + matrix_confu[i][j]\n",
    "            diag[0][i] = matrix_confu[i][j]\n",
    "\n",
    "    vet_lin[0, i] = cont1\n",
    "    cont1 = 0\n",
    "\n",
    "# somando as colunas da matriz\n",
    "for i in range(nclass):\n",
    "    for j in range(nclass):\n",
    "        cont1 = cont1 + matrix_confu[j][i]\n",
    "    vet_col[0, i] = cont1\n",
    "    cont1 = 0\n",
    "\n",
    "# calculando a acurácia global\n",
    "po = sum_diag / soma\n",
    "ac_global = po * 100  # acurácia global\n",
    "\n",
    "pc = (sum(vet_col * vet_lin)) / (soma ** 2)\n",
    "\n",
    "kappa = (po - pc) / (1 - pc)\n",
    "\n",
    "# calculando a acurácia do usuário\n",
    "ac_usuario = (diag / vet_col) * 100\n",
    "\n",
    "# calculando a acurácia do produtor\n",
    "ac_produtor = (diag / vet_lin) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Salvando o resultado em txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################### SALVANDO O RESULTADO EM TXT ######################################\n",
    "\n",
    "num_arq = 0 # variável que modifica o nome do arquivo de saída.\n",
    "\n",
    "# Acessando o diretório onde estão salvas as imagens\n",
    "os.chdir(os.path.join(os.getcwd(), 'C:/Users/Maria/_trabalho_final_intro/imagens'))\n",
    "os.getcwd()\n",
    "\n",
    "# Criando lista com os nomes das imagens acessadas no diretório\n",
    "lista_nomes = os.listdir(os.getcwd())\n",
    "nome_arq = []\n",
    "\n",
    "for dados in lista_nomes:\n",
    "    if dados[-3:] == 'tif':\n",
    "        nome_arq2 = nome_arq.append(dados[:-4])\n",
    "\n",
    "#-----------------        \n",
    "\n",
    "# salvando a matriz em arquivo txt\n",
    "vet_ = np.zeros((1, nclass))\n",
    "for i in range(nclass):\n",
    "    vet_[0][i] = i + 1\n",
    "\n",
    "with open('C:/Users/Maria/_trabalho_final_intro/saida_txt/result.txt', \"w\", encoding='utf-8') as out_file:    \n",
    "    for i in range(len(vet_)):\n",
    "        out_string = \"\"\n",
    "        out_file.write(out_string)\n",
    "        out_string += str(vet_[i])\n",
    "        out_file.write(out_string)\n",
    "    out_string = (\"\\n\")\n",
    "    out_file.write(out_string)\n",
    "    for i in range(len(matrix_confu_inv)):\n",
    "        out_string = \"\"\n",
    "        out_string += str(i + 1)\n",
    "        out_string += str(matrix_confu_inv[i])   \n",
    "        out_string += \"\\n\"\n",
    "        out_file.write(out_string)\n",
    "    out_string = (\"\\n\")\n",
    "    out_file.write(out_string)\n",
    "    out_string = (\"Kappa: \")\n",
    "    out_file.write(out_string)\n",
    "    out_string = str(kappa)   \n",
    "    out_file.write(out_string)\n",
    "    out_string = (\"\\n\")\n",
    "    out_file.write(out_string)\n",
    "    out_string = (\"Acurácia Global: \")\n",
    "    out_file.write(out_string)\n",
    "    out_string = str(ac_global)\n",
    "    out_file.write(out_string)\n",
    "    out_string = (\"\\n\")\n",
    "    out_file.write(out_string)\n",
    "    out_string = (\"Acurácia do usuário: \")\n",
    "    out_file.write(out_string)    \n",
    "    for i in range(len(ac_usuario)):\n",
    "        out_string = \"\"\n",
    "        out_string += str(ac_usuario[i])\n",
    "        out_string += \"\\n\"\n",
    "        out_file.write(out_string)\n",
    "    out_string = (\"Acurácia do produtor: \")\n",
    "    out_file.write(out_string)\n",
    "    for i in range(len(ac_produtor)):\n",
    "        out_string = \"\"\n",
    "        out_string += str(ac_produtor[i])\n",
    "        out_string += \"\\n\"\n",
    "        out_file.write(out_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Salvando a classificação em formato tiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "###################################### SALVANDO A CLASSIFICAÇÃO EM TIFF ######################################\n",
    "# obter metadados\n",
    "linhas = img.RasterYSize\n",
    "colunas = img.RasterXSize\n",
    "bandas = 1\n",
    "\n",
    "#salvar a banda em GeoTIFF\n",
    "# definir o nome aquivo\n",
    "filename_output = \"C:/Users/Maria/_trabalho_final_intro/saida_classificacao/classif_mxv4_5b_B2457e8_500x400.tif\"\n",
    "# definir driver\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "# copiar tipo de dados da banda já existente\n",
    "data_type = gdal.GDT_Int16\n",
    "# criar novo dataset\n",
    "dataset_output = driver.Create(filename_output, colunas, linhas, bandas, data_type)\n",
    "# copiar informações espaciais da banda já existente\n",
    "dataset_output.SetGeoTransform(img.GetGeoTransform())\n",
    "# copiar informações de projeção\n",
    "dataset_output.SetProjection(img.GetProjectionRef())\n",
    "# escrever dados da matriz mxv na banda\n",
    "dataset_output.GetRasterBand(1).WriteArray(map_class)\n",
    "# salvar valores\n",
    "dataset_output.FlushCache()\n",
    "# fechar dataset\n",
    "dataset_output = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################# LIBERANDO MEMÓRIA #############################################\n",
    "    \n",
    "#fenchando as imagens para liberar memória\n",
    "img = None  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (geospatial)",
   "language": "python",
   "name": "copia_geospacial"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}