aula-11-cholesky.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fatoração de Cholesky"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrizes positivas definidas\n",
    "\n",
    "**Definição (baseada em autovalores)** uma matriz ${\\bf A}$ é _positiva definida_ se todos os seus autovalores são positivos ($\\lambda > 0$). \n",
    "\n",
    "Entretanto, não é conveniente computar todos os autovalores de uma matriz para saber se ela é ou não positiva definida. Há meios mais rápidos de fazer este teste como o da \"energia\".  \n",
    "\n",
    "**Definição (baseada em energia)**: uma matriz ${\\bf A}$ é _positiva definida_ se ${\\bf x}^T {\\bf A} {\\bf x} > 0$ para todo vetor não-nulo ${\\bf x}$.\n",
    "\n",
    "Para a ordem $n=2$, temos os eguinte resultado: \n",
    "\n",
    "$${\\bf x}^T {\\bf A} {\\bf x} = \n",
    "\\left[\\begin{matrix}x_{1} & x_{2}\\end{matrix}\\right]\n",
    "\\left[\\begin{matrix}a_{11} & a_{12}\\\\a_{12} & a_{22}\\end{matrix}\\right]\n",
    "\\left[\\begin{matrix}x_{1}\\\\x_{2}\\end{matrix}\\right]\n",
    "= a_{11} x_{1}^{2} + 2 a_{12} x_{1} x_{2} + a_{22} x_{2}^{2} > 0$$\n",
    "Em muitas aplicações, este número é a \"energia\" no sistema.\n",
    "\n",
    "O código abaixo mostra que esta multiplicação produz uma forma quadrática para a ordem $n$. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAAZCAYAAACl4o35AAAACXBIWXMAAA7EAAAOxAGVKw4bAAALJklEQVR4Ae1di3EVNxTFHhdAkgrAHZihAqADSCoI7gCGCjykA0IFhHQAVOCYDiAVQNwBOWetu6OnXf2lt7vPVzNrafW59+jcq8/qrd87+vHjx52thKOjo9fAehfXfVxfcb0E/mvEGpQBZUAZ6MqAzj9d6VXhCzHg8+ujrWwOTAfeAC83BXdw/x7Rfdw/WIhTVasMKAO3hAGdf26JoW9ZN0N+fbwhLp46WC9wf4bO8RRBgzKgDCgDPRnQ+acnuyp7KQa8fr2lzQHJszcC8nGCnbcUwapXGVAGDp8Be67R+efw7X1bejjr1ydb6T0+Pjh1sJ6Z+3+cfL1VBpQBZaApAzr/NKVTha2EgZBfH68EYwkMvpyoLySWMKdtlAFloJYBnX9qGdT2a2Rg9OvNvJBos2heoriLXc+5nb+vtHnP4aXRJ0cy3Kh83hcG1ZPGgNoqjSetlc6Azj/pXGnN7TDg+vVmPlYQitGB50gvuTHgv1K+xkbgmYWJu60rYHuC/I+Sr/GyDMAeaqtlTXBw2nX+OTiTaofAwJxfb+rkAB14jH5wAR6e2s3999wndrTjG5r8N8g/cj0DbfkvlL+jrbyQNIhA/n9MIP+nIUP/NGFAbdWExr0IqbHVXgBWKkH/msw/hFHKFdrp/FNpx5zmpXbK0bF0XZ9fb+adA3SALyDyCf2SBjNG4yZh+N6DPRLMCeJqRh9PDO4anDPFmrUAA2qrBUg/RJU6/xyiVbVPIb8+2RA9n4CVx8TcOY/BfYIfC/oluBnh9yvwo42d04N+KlVyIQNqq0LitNmEAZ1/JpRoxgEw4PXrnc0BFjy+XMencS56XIi/4HqF616PhZALLGTzHQKGh7iomxh4SvALrm/QOxz9I17FcT1w+L6RUf61cuckY82cgt+kkGOnJIF7qqS2Co+pPZkhWc2ax4rOP2lm7DVX9JKb1qu6Wpv1azg9Pydn73kEy8/N+Vm85PH4HLc393aMfC6GV3aenU4o58aAL/aJLm4SuBl5wTyT9sqXdiUxZPOdg0FPSXu3DWSRO4Ie+8M6Jr8Jp5SN6w2uDybmycXAncTIL+IU7by2LJGJNlGsgjkWQ9ZWbfUC2GkvnnSt3VZBrDEbSXmNrdA2ef5B3aB/xcqJF6ForEhfa2Pob+bXkFU9/0BGCqdBPynhNKa31FZGbhBvig0hp8pOaJ/r117MKX1CnWZ+LQszBXIRe24ThvthMXLyuJDIpPfFLmMaIVgu9VGPzjhOmuaeAoY8xDTKuFGRdi1iI7vl5oCbmvc2Nty35JRcjVxQF67Jxgl5WZyiftRWBTKTsNpchdLQXzU4XdmQ19tWHBvjOEJ6MoaICWENtkrC6nLou0efimyFdjljJehfhtfmY8XX59L8Uq7m9EFWrU/HOE3yE8N98pzey1aQm4R3jks3D7KKfJpyEHL8Oog5tU+olzWvuP2172VzQGDjwiwVkDc+yUuexCgjaZPNQUb56ERsg8BTig/S3hejjvcp122DuiSKct2L/eJmyM3n/c4i78p071Gf3PEHodyn+GacQj7x2hMeOSBpYx71I5Ry6rVlrkzUT8I6w9eh2Ip2H0+QkOaJGLo78Y812CoJ6wz2prYCP8ljBXWD/hUrl76gXjb/aMN+E6v39E7kS2zazM0zTeYgg6dq/oGMGKdJfgI5WZzG9Foc5spNwivyGSM09WkjM8evg5iBL1gufUG9LK6svk9OLGRzwIVyZ2GmElxk7UwU2zHyvQuKURgst2WZ+tTlfZonDlwkiE/N3k2JK3fuHu2Jzatrrs1cHmXgmgxM1kVoxilkcfA+FgxI870MKhnzpMyOTZ1oP1Ev2VYxmSgvwmrjttMGW7QPdpu5NOTsxVaubuilz+6MLbcO7xGC/i9tUK+ZrUSmxJCdhFXqu7HBlm0rtEseK6gb9K9YuYtZ7tEuyD/KuXiMm3GkZ0/vRF4sRvvqOQgymvg05AQ5dfuC+kl+gnoxTrP0Co6YXKkncSpeqW/HaFtsJ7RN9mtbJ9MxzLFykYd6MRvQlrMnnSfmRQ9uBNxv9+NnJdw8uPnMbhqAgQs/g/cLhAyOc9SlsaT+0GiJPwbHKXCN39KIPC7aDN9xNeMUOrJ/VyKF0wFpxp8UmSVYMyAUVd2nrWyA0Esf+BXXIzvfTafw6raJ3efKTMUa05tbbvQmj5WYf8XK5/AlcsV5hxOpBP4qLL/4jBuGnZeQpULPGHqJp8n8k8OZsVcTn87RK1wm2kqq8/skksbg2KBRwuhN9mtbbQxzrFxkoZ6sk9511dS11xeecg3/JHAsghBfWmkmuegNGwMo4fcKyMLnVCu7NTJJHsNkI4LyF7ik/KbWSv4CF0nnlzGNGwMDjQPWDr045VPM5HclenDaQOYsVpuknumlbAW9fKr7hOsCfjLZYDfgdUJbqcwY1omiPhmlYyXmX7PlhVzZc+C1ocHO68OMIxXYe88/Ps56+7RPL9cfWQuS14ot+nUMc0J5Fldcw3ANXypo3Iz/jTdsJo5RQCfn5DU6OQBw58B72RE/RD1JI7suQD4NzGO5wdCIn9gSUU7d1wabXbR42mAj9mukX5vrDWLmkeivPTmlPuj5CB073+6I/Oac1sr0Yd2XEaGffrSIrWgfXBxop8DBz6jHUMvrKMhK1MgMYbVUdElCd/H8gz7PjgUB6isv4Qo4+ZRuP4HJU9lefxUW2Lv6tI8zchrykxJOxU6MfXpr5Ibw2rp7pKG7yK9jmEPlNVwZG3ADxlOhYbMwfH0yhDLzLa5LXPx+gXe42DkOPuZ9Bih7YMjXf/LFK/tIAlVvAmTyKXq23NHHBn/i4oaEsnisQUdk3iSE5E4qezKMDB4H7iywnuo72WhLfDIx7JThhhuDgQ+nj604JUcPoMM9sZDjM7EhcVVz6vQhV6YXKwWlhq3aSvoH/JzM+dnqM9jtb+ZX8jo7rmpkEhPDHNabkrS/aE9s2ePKwZ46VoL+BZneckcfO5c8VoQJyKBN+b5R9hxCGWhfylXP+cfLmfTbYG/t0169jWw1wWv3J5QutZPhKXtdFSzQG8Q8V17DFdryVOg3XO9Gn0ZieCkxN4YQOrf3xcBYea4+qd9CrpGR/eKUYOgVh/qGMp4M2G/B8372ZdFcfCG9ubJYH6EZVoNtE7YCVk4GfAlptAvSHOQkZXzpp4RTaQM5wXEn9WIx5DTHarB1txX0BP0rVh7jJlYO+Xxomn0ROdZWytGeduzOleiLxcAyyynym/uJjcWn166Tk26NF/K62ymGOVaew0+sLnSNL5ue4KYm/BxpHCuPNO9WzFOQtWKb4MKu7gx4OSFdIE1nZeDpwfjLkENO3Z+J3hJxHbBuxlYYePyoiS+jXlvc0V68/8vKq01W26oT1u62ivlXrLyWeMjnUy7/XWxyepcpuztXqXhCnHXykwFaSG8qdrdeB7zd7RTDHCt3Oai85zrzBbZ5mv2rjGjEJyF+JsGdJtM8Kr1EB4bjtVg56lYFgoaA2Y8rqgQv2DjGGcr5NMod/E4A50c7GZk3Mb2Z4obqvbCWYOnRJsYZyrmR4/HcN1w8Ih/GC2xV9c5OTC/0ZIdeWLOBZDQA5uBYiJVnqJpUhWzOedW/CjsRvHBGjDOU9/LpoC1LaemFtxRPSrsY5lh5ig63DmRyTfkX1yPMT/LPB5yv+JHZefbmAI0WDegQNwdv0ZlV/NbComSocmVAGdgLA2Zy5js9F5bC4fQOc5F9UmQVa1IZWDcD8GtuBLjhHR5ecM93D17hureZzQFAB08s1m0CRacMKANbZgDzT5en3C1zoti3zwD82nsq9D+5zvG5iRKZogAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle a_{11} x_{1}^{2} + 2 a_{12} x_{1} x_{2} + 2 a_{13} x_{1} x_{3} + a_{22} x_{2}^{2} + 2 a_{23} x_{2} x_{3} + a_{33} x_{3}^{2}$"
      ],
      "text/plain": [
       "      2                                     2                       2\n",
       "a₁₁⋅x₁  + 2⋅a₁₂⋅x₁⋅x₂ + 2⋅a₁₃⋅x₁⋅x₃ + a₂₂⋅x₂  + 2⋅a₂₃⋅x₂⋅x₃ + a₃₃⋅x₃ "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing(use_unicode=True)\n",
    "\n",
    "x1,x2,x3 = sp.symbols('x1,x2,x3')\n",
    "\n",
    "# ordem do sistema\n",
    "n = 3\n",
    "\n",
    "x = sp.zeros(1,n)\n",
    "A = sp.zeros(n,n)\n",
    "aux = 0*A;\n",
    "for i in range(n):\n",
    "    x[i] = sp.symbols('x' + str(i+1))\n",
    "    for j in range(n):    \n",
    "        if i == j: # diagonal\n",
    "            A[i,i] = sp.symbols('a' + str(i+1) + str(i+1))        \n",
    "        elif j > i: # triang. superior\n",
    "            aux[i,j] = sp.symbols('a' + str(i+1) + str(j+1))        \n",
    "\n",
    "# matriz\n",
    "A = A + aux.T + aux # compõe D + L + U\n",
    "A = sp.Matrix(A)\n",
    "\n",
    "# vetor\n",
    "x = sp.Matrix(x)\n",
    "x = x.T\n",
    "\n",
    "# matriz positiva definida\n",
    "c = (x.T)*A*x\n",
    "\n",
    "# expressãoo quadrática\n",
    "sp.expand(c[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ainda falando sobre o caso $n=2$, observamos que os autovalores da matriz ${\\bf A}$ são positivos se, e somente se \n",
    "\n",
    "$$a_{11}>0 \\qquad  \\text{e} \\qquad  a_{11}a_{22} - a_{12}^2 > 0.$$ \n",
    "\n",
    "Na verdade, esta regra vale para todos os **pivôs**. Estes dois últimos valores são os pivôs de uma matriz simétrica 2x2 (verifique a eliminação de Gauss quando aplicada à segunda equação). \n",
    "\n",
    "A teoria da Álgebra Linear permite-nos elencar as seguintes declarações, **todas equivalentes**, acerca da determinação de uma matriz positiva definida ${\\bf A}$: \n",
    "\n",
    "1. Todos os seus $n$ pivôs são positivos. \n",
    "\n",
    "2. Todos os determinantes menores superiores esquerdos (ou principais) são positivos (veja Critério de Sylvester). \n",
    "\n",
    "3. Todos os seus $n$ autovalores são positivos. \n",
    "\n",
    "4. ${\\bf x}^T {\\bf A}{\\bf x} > 0$ para todo vetor não-nulo ${\\bf x}$. (Definição baseada na \"energia\").\n",
    "\n",
    "5. ${\\bf A} = {\\bf G}{\\bf G}^T$ para uma matriz ${\\bf G}$ com colunas independentes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretação geométrica \n",
    "\n",
    "Matrizes positivas definidas realizam transformações \"limitadas\" no sentido de \"semiplano\" do vetor transformado. Por exemplo, se tomarmos um vetor ${\\bf x} \\in \\mathbb{R}^2$ não-nulo e usarmos o fato de que para uma matriz positiva definida, a inequação ${\\bf x}^T{\\bf A}{\\bf x} > 0$ deve valer, ao chamarmos ${\\bf y} = {\\bf A}{\\bf x}$, a expressão anterior é o produto interno entre ${\\bf x}$ e ${\\bf y}$, a saber ${\\bf x}^T{\\bf y}$. Se o produto interno é positivo, já sabemos que os vetores não são ortogonais. Agora, para verificar que eles realmente pertencem a um mesmo semiplano, usaremos a seguinte expressão para o ângulo entre dois vetores: \n",
    "\n",
    "$$\\cos(\\theta) = \\dfrac{{\\bf x}^T {\\bf y}}{||{\\bf x}||{\\bf y}||}$$\n",
    "\n",
    "Uma vez que a norma (comprimento) de um vetor é sempre um número real positivo, o produto $||{\\bf x}||{\\bf y}||$ no denominador acima é um número real positivo. Se ${\\bf x}^T {\\bf y}$ > 0, então $\\cos(\\theta) > 0$, e o efeito geométrico da transformação é   \n",
    "\n",
    "$$|\\theta| < \\frac{\\pi}{2},$$\n",
    "\n",
    "ou seja, o ângulo entre ${\\bf x}$ e ${\\bf y}$ é sempre menor do que 90 graus (ou $\\pi/2$ radianos). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Método da Fatoração de Cholesky\n",
    "\n",
    "Trata-se de um algoritmo para resolução de sistemas lineares ${\\bf A}{\\bf x} = {\\bf b}$ através da decomposição da matriz ${\\bf A}$ em dois fatores simétricos, ${\\bf G}$ e ${\\bf G}^T$. O método é aplicável apenas a sistemas cuja matriz associada é simétrica e positiva definida."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Passos\n",
    "\n",
    "- Primeiramente, é necessário checar se a matriz associada ao sistema cumpre os requisitos da fatoração de Cholesky:\n",
    "    - **Simetria**: a matriz é simétrica quando sua transposta é igual a ela própria:\n",
    "$$A=A^T$$\n",
    "    - **Definição positiva**: averiguamos se o Critério de Sylvester é satisfeito. Ou seja, verificamos se todos os determinantes menores principais da matriz, constituídos pelas $k$ primeiras linhas e $k$ primeiras colunas dela, são maiores do que zero:\n",
    "\n",
    "$$\\det({\\bf A}_k) > 0\\text{, onde } k=1,2,\\ldots,n \\text{ para matrizes }A_{nxn}$$\n",
    "\n",
    "Concluídas as verificações anteriores, decompomos a matriz ${\\bf A}_{nxn}$ em uma triangular inferior ${\\bf A}$ e sua transposta ${\\bf A}^T$, a qual é triangular superior.  O processo é descrito abaixo:\n",
    "\n",
    "- Para uma matriz $A_{4x4}$ obtém-se o fator de Cholesky da seguinte forma:\n",
    "\n",
    "$${\\bf A} = {\\bf G}{\\bf G}^T \\Rightarrow\n",
    " \\begin{bmatrix}\n",
    "a_{11} & a_{12} & a_{13} & a_{14} \\\\\n",
    "a_{21} & a_{22} & a_{23} & a_{24} \\\\\n",
    "a_{31} & a_{32} & a_{33} & a_{34} \\\\\n",
    "a_{41} & a_{42} & a_{43} & a_{44}\n",
    " \\end{bmatrix} =\n",
    " \\begin{bmatrix}\n",
    "g_{11} &             &            &              \\\\\n",
    "g_{21} & g_{22} &            &              \\\\\n",
    "g_{31} & g_{32} & g_{33} &             \\\\\n",
    "g_{41} & g_{42} & g_{43} & g_{44} \\\\\n",
    " \\end{bmatrix}\n",
    " \\begin{bmatrix}\n",
    "g_{11} & g_{21} & g_{31} &  g_{41} \\\\\n",
    "            & g_{22} & g_{32} & g_{42} \\\\\n",
    "            &             & g_{33} & g_{43} \\\\\n",
    "           &              &             & g_{44} \\\\\n",
    " \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "- Este sistema pode ser resolvido comparando cada coluna de ${\\bf A}$ com a multiplicação de ${\\bf G}$ por cada coluna de ${\\bf G}^T$ de maneira que ${\\bf x}_{kxn} = {\\bf G}{\\bf G}^T_{kxn}, \\, k = 1,2,\\ldots, n$. Através dessa sequência são obtidos sistemas simples, em que cada coluna $k$ terá os valores de $g_{kn}$e, após a atualização do iterador $k \\rightarrow k + 1$, não será necessário resolver novamente as $k − 1$ linhas do próximo sistema gerado.\n",
    "\n",
    "De posse dos valores que geram as matrizes triangulares, é possível seguir para a última etapa, a qual consiste em obter o vetor ${\\bf x}$ fazendo\n",
    "\n",
    "$${\\bf A}{\\bf x}={\\bf b} \\implies ({\\bf G}{\\bf G}^T){\\bf x}={\\bf b} \\implies {\\bf G}({\\bf G}^T{\\bf x})={\\bf b} \\implies {\\bf G}{\\bf y}=B.$$\n",
    "\n",
    "Aqui, usamos o mesmo processo utilizado na Fatoração LU. Primeiramente, obtemos o vetor ${\\bf y}$ após a resolução do sistema ${\\bf G}{\\bf y} = {\\bf b}$. Enfim, obtemos ${\\bf x}$ através da relação ${\\bf G}^T {\\bf x} = {\\bf y}$. Fazendo com que, assim, seja obtido o vetor ${\\bf x}$ que resolve o sistema linear proposto no inicio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algoritmo para a fatoração de Cholesky\n",
    "\n",
    "O código abaixo é uma implementação de um algoritmo para a fatoração de Cholesky por computação simbólica."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAACDCAYAAAA5xZrQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dWdIVN7aFDUHUowP/jqh34xnYZgTgGbgZAfYM7GAEhJkB9ggwngH2CCiYAXXfKwL8Rz3eh8tdn5ASZZ5sTzYnm6UIncxUr6Vma29t6dx49+7dR3OaGzdu3FH6r+fMoy5t1etGnbvdjIARMAJGwAhMicBNEbovZN812GcTZPaL0vgawraknaDcTsIIGAEjYASMwEeij68baGTgRG9lGP2q9yrn+O/Mf/CrMoYbvSMC+ufgyI5gBIyAETACRmAdCMAQ3q4U5a6+v8EN8ecXsi9lvxTBe6XnZEaE9IkSe650/5gsUSdkBIyAETACRuDCCIi+/aAiPEHSenOusigTqPd9E9G5EHa6RsAIGAEjsAYEctHu1OWBFf65mmgU9+J+LQuxRZz8UPYzEV3cbM5EwNieCZyjHQoBj5NDNfcilZ2FkGbc6I95LeR+X98oMCFGDvuvckOsfNtENEdq+LuxHY6ZYxwPAY+T47X5EjWeS7QLhwlHWhh1YLhPiOjPiYhGz3/pWauMpDhoFENoa02Xf22kHTpOja3S+0X2iezz+KTtbIzAphFQX550DlJ6jJOf4hh5FtPfNEYu/HkITM6Rxs70jYjl55UiJa2n3yvucKkoJRVGaaAABTd7JXsyiXf5Fwkd52VKbEmLDfQkMWDx85fsl8eB0zXdKQJTjhPmrJcaJ5x24HgE34yVr/m2ORYCc3CkaDLRYavmOzn8qY53nTwi0eWITIkjVZhXshDSpyls/uzyz8Me5H0ybIVXUOfOcHukdyQDtJONEdgyAlOOE3DImQUkZzAFNgdEYBBHqsk0dBQRshLhq+D2o/zzDsZqDa4SWz1ek9KruleS9GcTAjNhC9FMZ4jTwid3ayqO3Y3AKhGYepzEhX5eVyQ2bfNiHtbvO0PgZt/6qCP+pLCILuq4zZBMDFMS01bSf1H5husMRFRxv5E111MBaMDnJNiyCKoslBCzY9jLtjECW0dgknGSg6B5CyYBbvfklEIezu/7RaAXIY0Ejr2AJOZrEmHAjT6uwiU3uBoIZkEolSYiYL4T53NX4dK7nG36ILAAtiycUBBLnGmfYjmMEVgVAnONE81jMBjoEDxSHpasrarVlytML0KqDvLv2BHDxrqKd7LyioSx7Qaje4p3l44ny+QMhxM25mNnfK5vm/MQmAXb2E7sa58sjs4rpmMZgYsiMPk4YWzIItb9XOPFc9hFm/dymQ/aI4WYqrNATH/Qkzt0cw4S4tqo2Ulc+X9bU9U6t5pgdmpCYA5s1b5IDDjfi/jdxghsHoE5xkkGCswBF5tzYqGNocii+HUvCPTiSCuVTXukBVcaJ90/YketBPfn1hBQeyK6Z680EFG+ZdNe6daq4/IagUkR0Fi4Lft3w5jgyJ7NwRAYTEg1ucKFsuKCK2WTHQNRZf90DtPVMbv85yjTXtI8wS5ODiyWXugdBTCOw9C+ufRhL/V3PYxAHwRK4yQyDG8VMdcbYJzwXT0n3yd9h9k4AoNEu1ldIZp0nIeaaNGCK50PzcKd9ao0UUJi8oYzYvX3TM8X6sBhr67LX2FtGhDogR2KEyyQwLwwcfIovv1iBPaMQI9xwpbUjwr3Rs9PZZmzuPo0J65ysjkCAmf/jZo6EAeQ6TxwKvfcgY7QXVxHI2AEjIARAAHRwEn+Rg2uFM7lXyaiwGpjBIyAETACR0Rg8B5pAknEk31Szk0l5aPk5acRMAJGwAgYgcMgcO4eaQBIxLTxuMthEHRFjYARMAJG4NAInM2RHho1V94IGAEjYASMQERgFEdqFM9HQBvViMS5XsxmIgQkIUF5zsYIGAEjsCgCJqSLwl3KjBtQPPGXIPGHETACRmB7CJiQXqDNxI1yBtd/uZRhL0y+0Of30Yl3/gDBl0BkGPnVCBiBdSJgQnqZdoFgtP3d3GVKddlcfxHhTH9icF9F4QLw0v/aXrZ4zt0IGAEjUI+AlY3qcZnb9QsRDXOkZZTza9j4Z6A74lI5p2xjBIyAEVg1AiakCzePiAPcloloBXctLPKjVGCU/rqvEtKfRsAIGIF1IWBCunx7cEdn6R7b5Yuw3hwjF/pQJQxi3vWW1CUzAkbACLxHwIR0+Z5w32LdetAjEf1NvtzdbEWjepjsagSMwMoQMCFdsEFEKNBG5VrFwxmIpOwT2Z9kn8m+lAWPYPTOfihna/nXHy6E9tWTAGFjBIzA6hGw1u6yTYS27tNls7x8bpFI8m9B34rTfBUJKN8515n+TYh/VMDgF4hq+PKPETACRmClCJgjXbZhjnp+FO6S/6zNufGSMpH8Ppe9kVkffVm2bzo3I2AEzkRgFkIKByLbKMY7s6ybjiY87qgCJeKx6QoNK/x3Cs650GTQys2JanL30wgYASOwOQQmF+1CRIVClxgvAKWwaK8+EBdynZCTG/tmiEAxvO/lhhu40cNp68b+QJ/ICScaueyT4v5Rpf2P1Ceovo0RMAIbR2ByQio8OsV4mkDhSJhMeV7JFoRU73u94YbFwT3ZQxmIpNqb9g1EU+8sjmh3+sl38v9VT5SLjtgnqLqNETACG0dgDtFupxhPkyf7ZSiSvK3BD8KazC5uuBGROOG8UgUP8uTs7EPhwL/dIOKm7XErzNH6RFFxvxgBI7B5BG6pBmhHMqnlGpRnVSwSDIhGLzFeXSaaUPd4ww2Li8Np66b2hUjq/ezbnHbaJxI8fhoBI7BNBJjTAkMARwqHwL4Uz1FGEx4iPGzgwERYkxgPIo0YD79eJhLlvdxwA9h/9Kq4AzUisLM+0VhPexgBI7AJBNiOCnovcKRTG4gGYrwXeqazgLihgNTLxAkz3XDTm/j2SnzhQLEuV1pEtHL8CpdEnmkh8lpFZSHx2ZAFyMLVWyy7PfWJxUBzRrUIeKzVwmLHEQhMTkg16Y8S48UJE0UU9tHCDTdKM7yPqOclo7JqaRVrqs5pZfNlIrhyY+Fxu0pE5Q42cPwQ3rBQycN0+SvO5ozqRH331Cc21wZ7KbD60tRj7Y2w4cwzuh2lEwh7wcz16EbgZneQ6UOoM38hWyie6J2jIclAQLjdBo7sb9ncT5/rMCrzD5VyNxXse3k07o8qDYgE4oGfExGNCaFoVSLACgsxeaJwHAlC65nB+5dsMF3+Kdwanyr75vvEGnF1mT4gMPFY4/+EX2scPmY86h3FySDm+5Cj346EAPuY72T5j0zO9LXaGJbwq7Jd5Z7SX3VnVUv9n3WlqzAMtkZM5c+AJC24zyIc8WR/qnG7k9zkn9ouuMU4jf4pXvWpeBBzFjDk2ccSdjX9pVoff3/oR8biAxbqs1OONdLiqF4Ys3pn8a/PD/n5fd9Y5G1+Sx+DjDrHjUERdhZYq1pEqlg0kyGojUZh4aa7lIzQ6OU40HVKKK6cyaPEkUZ/3NN+a4qTu+Xvdf4pm+IZ8861pQu/sS9H7y9j8XP8SRGYbKypX8OF5obxUzde8zB+3ykCgwnpTnEYUi04x19F7NgX4Z9M4Mzy4z55Woh1H+UO+XskmHCD1fiBQFfT1Xf1/lk4UgxiYJbDrf6EsTECR0Rg6rGWYxjThkjfy939fhwEbi5RVXU0CM67ue0SdcmIW9r3bONK24hsXlw0nHPDajcQV2H2jSxcZp0JCjgqU+I8q2G6/Kvh/W0E9o7ApGNNYxNdD/QUHmVzw94xdP0qCMxOSNXRIDSfqpPl/+wxy3ulbrN+xkEDAYPrPDGx3q2inkgAIZgFoVQ89lr4TuLbuwqX3ot8FA4iiUj4ceGYvXT5Z0H9agR2j8BcY43xJ4tY93ONufyPGXaPqSv4AYHZCamy4uhKo3jzQ1E2+QahRNsU8WzVcHa2jxYf4qC7SgOuHeKImBaNXI7+sNo9GZxyh9giYq7u0xCNeK3+IZB/jMDxEJh8rGUQMnbva+yt8pRBVk6/zoDArHuksVM9j6vBUHy5QXTodGiH3pWF++L8VXWfUM6rN09VQgYOXHdVqeh+E6HLaxWxgehWTZ0bRJK8+O/OQETj99uEX5d/NZOlvnfW7kvB5nwmRGCqsaYiISH6H9l7adxlxeRIms3BEJibI+W8YyF6jJMpRyc4C4k7nCoKM3Wiy+JvtuRfMkqn0a8UcP6PJLotiXdVPuo0+cIgpssi5IXe2TuFiMPxB/y6/BXuIkblYvG0p3a/CI7OdDkE2sZSJMhvVRq2dpJhLPL9e3Lw8zgIwJF+HKubnpPUXh0R8WJVtAkRYF8vJzKlP7tWPDiu2r9Ya/ObpNADE2FAqUzUhTLnBsIKtzq1QakBolTCNQ5s8uryJ8wlzK7a/RIAOs/FEegaS0iMftT4f6Pnp7JI1riZLCeucrLZMQL/THWDkP43fqRn8hv7/FadKuz1ZQmhIv4g+4YA5USVIxxweX+qg57sNbT5ZWku/Up5q8dgvlFZ4RQnNUrzk7YEu/zb4s7st8d2nxkyJ39JBLrGkvyZt0pz1yXL67wvgsB/Uq5ni3ZF6O7UETsSljtKMtz8URi5wUlh884HoX2OX/Qvwm/oJXGeQbyrerAyPRFVb6g+kxZ1x+0+KU5OzAgYge0icBYh1eQIJ4myy8leZZw4v9aKraR8o29EHliIKcSWfUTSgegM+os1hV+NiStT6kVdMHDSJdFrcD3oz17b/aDN6WobASNQg8BZhFSTI/uciC7hLh9W0mU/DFtn2FfgL9bgWOHcSKNWO1XuWzL5MRg4UysclFtvr+1erqW/jIAROCQCt0bWGkKIaJZbPVC8gTheQWjr0o3utX514TfkhngXThQFK/Z54VBtIgI7bne3sREwAkbgo7M40oRbnCDhShMHyhPierZB5JtzrHovlI7a/M7OcJqIaXFA/SGqNgMRaGvbNr+B2Ti4ETACRmByBPgnF/YqOeOH6nauCNQrM01y7A1y+w6KQ5wb3YOotlfd80DCAQzBkssSrGyUg+N3I2AEjMDOENCcjwSSOxFujBXthuMqShACDDGt/vvIzqBrrQ6cKNf2mYi2wmRPI2AEjMC+EBhNSCMcnA3lSrwjExG0lJOId1+9xLUxAkbACBiBRgQmIaQioHCkg8XCjaXaoMfBFxEbbLHliyzJDXvoaKzbGAEjsDEENMezFVprJiGktSnb0QgYgSoC3HjVOBirgf1tBIzANhCAkCKORdP2yGLZbbSWS7l1BHwsaustuNPyS1qCoiRn4DG8ozhqmhDgaPxhKy+cUuH4C2c/ETnxtDECRmAGBOJExX/N2hiBNSLwiwjnz1gVDnpw8j/Iayz0hct0P2I17hzphSvh7I3AlhD4SoXliJSNEVgjAldZoVjwcZd6uM41c/drAwKjLmRoSNPORsAInCLwpZwOz5GicOUJ+rRzzOUirEt/HtKUjzhR+mcycFqlv7dMHn7WI2BCWo+LXY3A1Ah8pcnq0JrtcVJ/Khy8Vzx172pO74lw7y2mjYsc7k/ngh2bngiYkPYEysGMwEgEDi0mi0T05dEXEyP70ODoEW+IaSdnGonob8rknuJZ0WgA2iakA8ByUCNwDgKaoO4o3mG5UdUfLdA7mpx/PQc/xxmHgHDnspjG/48m9UhEUTIKWqj65t2mJwI+R9oTKAczAiMQgJC8GBF/61Hhcrj9bHNGBIX7VJM0AXHntyJMJdG0wvD/ww9yd7nR5q3HSUamjfIa5eJaVu577VqoQSBpB4hqnUERjgVf+AcrPdOxyLqwdqsgYEJaAWTvnyMH75omhi011V0Vtvc+1ZYq1lVW9Tf+veltj4m+K6nF/WPZ+ROKxKU9rBBLlHIgrjzRes0JLMdJwj6j0sGf9i/uIh+bttIi/T9jOhDIXFlIn2WjsK8UFs6Tq1xPrjKVW1G2ckx/9UWAldM72S8EJpfQ2+4UA7UxkxoDMLSx3v9O7zxlGPCIdP6WRRRX9AV9P0/feifc6/Qd445KmzSzdNhLK/Le+jvYyfKHBrupU9+6qN5waz/0Db+mcCp3+EcnyiTDPPmsrnxyfy1bHS9FH5YfnCOJFH1A72PSztNhvBZjuq58yU3hGKNFuZK7n+eNS2EJ9y743vkcqYA4kkEb7ykV1sqUiaG0MlWH+FOW1fdbwlQMK+5k6s6ZjUkbcVkqC9xbek/5bf3JJJtzK1uvz5DyM3n/PiTCisJCAJPSDdKYF5ED7Cyi2jvnEFl4Vo+TjEn7mvErCxGFAQocc2eh3o8r4pG3zYQI3JowLSe1fgRqB68GYtO+SVGjcyeGnmmHiUGZMVkxMQSRWJH5hl/ipJUm4w3XZHjRVXcIyLXa83p47FXE+FGleKh6wHG+kUX8Oagusf3rjpOMSluYsieKuPa+LFxmTrjldWpoB4WlL34lu7fF6mmFo4vqDNPA3ILhffLrD01IA7aH+Rk1eEFJnRJivIqJYSOtxqTVpQiykaoMLiYLIqQXmzQiPBCbswlOHCvsX3KcpESAx6StdH9S/MeASjoQClmkHn0WbPRF2uXsepHvxkzrfvUUdYGQfhwTSs8p0nUaK0RgzOClOiudGFaIdKlIrICPqrHLoqt1ERH7FCJKuD7E+miOogHbGk9hVm2yegWxq76ZzPuKYBvrFtNlIRsIqb7pX3D9fYgo6RKOOJsxYKfC0pfoG5SfO4FLCxO5tZmrzLPYlhqYRpZE8frP9AYh/W/8SM/k5+cZCMSO7omhJ3YTTAw9c7pYMIhDp+j8YqWbN2MmPhStak1se5Ru2CNHTMkEz3dBFKJbrVgu+sHxM8n2PQaioNOYmD/ia+r5s75RyEttTT1wP+s4SVPaTP7yIy/mGMTN9K97sn0NcXLC0jfeRcLFenK8J/QJfaO89pdspyg7FVhx87C0V3W/OgUd+vxPHoHOa63dCTQqhSMDmpV10IDmKVvV1sONQYBlkim0/fSOHwOPP39+wrc6wWKanuQX86bM5M//Z4b89U69cE+2pLWbwjU9Fa8tbeoMHtSbgbJovZvKPIW76jIIpynyXEsaqjvEpFFjV370MSbJ1MfoIyW89N2oLY6f7Cza3ko39fPNPROedU/Vi7FWwrgu3FrcKKtsdY6kTQq3vmVVHObnsMDpG6ctnNICSwV59xEcqc10CEAM0HzNxVLV1U+bvD4QWMXvfT5suqKHHkG5sUFslKetMo06ZxYxaUp7zzfeDBFB5ZDv4Z2Jq818J88HWQC4hXzs4HWV+VfFcvnlCCfa3uJe0BhGL6DL5OmEsOqve/0D9reqYI5pFzZr8E8iXcqSxlPu1llG9QX6Yu1+dWfkHgFMSHuANCDIaicGdSRWcZs0W53UhDkcFpP/UQ2TXi0xjRMbfjnhRAnmWfRjqc/eX6NYDv+IMaLfE21v+SNmTaJWvdoIAYhoITpfOyJqw+oCnjGFKcZV7C8wIXCvJ/vsmf+k+9UUIhkT0oTEyGdsrNVODOqQe11hj2y5WaOzf4co6agG7ufTuspHIlgQ2kgQ4UiZEL+Tf0lKEccXSjYQ28IoHIR40DGQIvIxX5ijaJetGvpHoWwU+wVjrHGfXX5BnKsnolhMUFh6/zrN781pknEqTAxCoW5ioNGYGPArTNvEoLCsnPjfxiNPwgVWG36BmypWzhuux7lFp+8zcTeZb+XBOU32xhHV0e9xK5k4VpJYruCmYrwQVmOG4xzhGEgpsj+qCMDhFRhWPdf8rfaGiLJ1lm89JbdcslHaTlN4rnm8kdkqlzu62rMQUjq+LH/d85MsopqXsoklH13oFSfgiWHFjXOBoh39P0iDMlAT7prYmBThJB7L/hGfHJYvuFHNGxDiwIWQjr55T0ex4FCDifPLkGMgKepkT5XhB1nmPOxz2ZNFhNwK0XU14za/atgR30hJEIFuyggbuMnb6hvVPW+20+hnydTtsye/WZ8QOPbPJtGUVDp0HhoqpMdTlvQBIWjnpafc0NAsueub8IDGKhWtvknKlfK89FP1AR/qxQo8TBKUKb4Xd9/qGxyK70uUW/mndqAt6KyltorlPmnDVFaFb/RLYfb8VP03ox05RzvQX2SVdHncD/kGQ9LIbIGp3OifEFb650W1vZU/ik3Fnbd6L41dfSexdd091o1+Q7DqE1blAMvBGq990p4rTMIupR+/mR9D/8rro3fmqUB08U9xeMqczEdyO5veKC75KOmQdkHoJiFYShgi0aXS3thxFDesYmPh6Jy7umRZ9fHEMGJipV9swaqdWSjVXnK+hfJPVUZhAOEIR1SmSnON6aiObMPkzENt2ysM47+WkLX5TVFnykceU6S1VBqxzGALLUi2WNTLjf6V4w5Vg74Ux67iNwuuukXM2fSGPGQFxTzHX2C1W1XalXG4ekvsOsBUTa6KjgbWrq6yUt0b5fPyK0RaVVAu8I3YLLSj2okBWGqHtjZs87tAPS6VJZi9uFTmK8r3kcqCOK7Uf1ZUvqmKAneU9h7RIg4X3GssrElrmHZYU3n6YP+XAoEt3GRhhOt1/EjbaYw18E/77IV+Scd8NAm9uVWUbIIXTbhUGJtv/J6otLdlBUBx4qYzstIoaem1xbXfpAjQjmufGCat8MSJsQhktXtoo/H7WOP5b+YGxvaOwYBIjbrgfgFsYHI+WyCfybJQn/mkLTH5B6asLUybH31yCnoDIf1YFnXo/2vLsI9fLBSDhUkYhQBW5UmMe6LSTpg6o3QgxFZprwNnObctTAzLoTE8J/r+o+HRdhkDHJBwwC3s0oyd0OcGRXMx+8iPmKPnzmtr6Y+gN/9QXaGd4Waj/+p5JTuVBm8nq03GTYYGV8WCejOdE2Isy55C4o6aotp9QgTAXsntXRw3IWInSdFnWyct+rViQVzS4pP9MwjOZ11xFWYzRnWBK0Vz3+P4Aq0m3GFs0IZu3Fa6QLFWkaWwGUNv/leVgHZOf0WgGuvsCTg2OBNJIKT6ZlV/rTRNRGktm00gEPtxa59VGCQ17Pt8mfq33NjX2asIlAV2qK+ehzNxLqPNw+JJ38UF921+EwH1m9LxFlkFTOHOAmMSeoNod3HT1HE0oSCvHvPPBovX5agZNrUheLT5HQSvr1RPtidqjfBhAENUuKElJ7hc3hBERXlEhUfj8I0sHAUrYP5mDC42mC7/FO6ST+qpcj6QRZsVonooozqH7SpVOr9MIGDQ5jcWJOH9RGkg0s372dhkNxVfGMCQnSxihMlk9IZr48iElTAr48bBL38bI3B4BDQoO8WTCsN+FLer1GpIxskN1flPGMwJVLkj2uXoWDHZxrAcAQsa3fGbMgQOo8s/pV19Kh7EPGlEVr3rvinn6P8JJd+8znUZ2W06BJbA+1J9aTqUzktJ9WYMM15vDOZIFfn96dPz8p4tFpWZLfEeCa8Vlx5F30yQFbRxIJACrGt1j8ZuLRGNYKM9ya0+ORGFsCH2q9uXzve2WPQygHPT5Z+HDe8x7/xC+JMw5zj0GQcKc07SjnMmAmPx7hp3l+xLZ0LSO1pX3VNCgwlp34RTBkd5GpdDtDT7KXCGbUQSIDi2VUtsNalBMLFV6Q+iJ053l9z1jfZ0biB+BbHt8s8jLvGu8phK9gRafYF90iRZoP3hbvJFUc+U9hlsS31pMCE9p8nUYVjJhzsyz4nfN86WgO+qkzBbJeffVe6p/dfSpmoPuEAIIFshjUbhCHPdGOCDBwfIcwPBDERUaXBRySvVvUSMY9pws/fyiOm9yz+F83M1CLDfnQz743doQ7V7n/6T4vm5AgRmJ6TqGKy0Pl3LhLgCzDuLEDErXd7dGckBZkMgEijSh7ChSNRmEM8yKdYaJkmlB8EkXDD6hkjznbjRuwpX4noVhsUol5SgOJLC6fO96fJP4fxcDwJqx1y0zjxZ+teS9ZTUJelC4GZXgAn8OSfng+nDgESrsRDfDYvq0DMgkC4TgYDdFtGC62wyEFr2MdsMHOVdiJ8skhoIbxLxQTBPbkTSpMs/pDDxfq44g/0Vz2alCMT+xLaBj6istI26ijUrR6oOgoiKfQCLKrpaouzfqRlaDr6+L7U9IlA4KAzvcNglUWXw2cCPyp3uQEYcS5/OucdqDZgMWxeOcTzUHQGpc6umD+F9zdhSOiWuNQbs8q+m5+8LIhCJKOc878V+ccHSOOtzEZibI2XyLFT5KSQdR/aI/1Xaq43iwNrDwoO/lOKcJBIJJvcTLqoXIOsKlBYCbeJdFkEnotdzqhHHCvfU1u3LXnX5n5On4yyHQBzrjA3GCHMj7zYdCMR+vyoaMhshVWXZ9+HQeWFix0Hslc7KsXJnkkgTVB720n+AW5Rl4Rf2SvZAdK4y3BBdBkWKzG2Lr4lA5ntb1Xq0iX2rYVu/RZBZUL2VzRdWcMR8/97lrzA2F0aAeVC26c++01Emzg/zF1+0rU0LAmulIXOKdvl7mqrMnxUX5+fShARkpQ12AQUhIR5PJuNiEmnzU7i9GOq++ZWp2jgnNrRlqZ232FiqE7fzUHREuydGfiwK8759EuYMB8S9PyrtN3p+KkveXJ6SxkWXv4LbXAIBtRmE8XO1VeI4H2btxlEnH3UZ3jCrpCFnE1J1EgY05+VO9mnkh8IEV1NVDar7DzJHJtjSxKP0wl29sRNmQcMZu0a/UsBtf2x+fzSHX+14W997UqSgvzaJdnGvHmvJ4Rj8rvFAfqUxkifS5Z+H9fviCNDvw3ynccAiq1aBUH5I7qrXPhJ+FzoGqseUZpU05CxCqoaHAMI5IbKoXnPGxPm1Bjj+hVE43LH5pECYQoSrONdFhAO+CCMWJydi7q1CEdt8b4oUtA//SFR33g8uvLSdsdW2c7knQYD5Lo1niGLpz77Vh9I8yrMkfdM3OgZhDo3h2O45NAcb55NV0pCbapzBRg2MeBZxBUSRVVduYL1PRJORSEIoASK/2JyOxvEC/I5uGFB72B8NSmWqC/0gibVO+sRGGztxnHVcKW7sB9sYARDgkg3+7Bt9EUTzEEIIZjDZPPo2uWXPIpzc9qJjkFVv+Ouaacit4dUpxWCSfK6OwiFxDprDUW7RqlkAAAciSURBVF3RQUqhPnywn0PHYjKCgBIfNzbdbd5z+Q+2AkScIMLCSGVm9cy+eFoQ0ab0ByYRTGrv91/b/U0SFepWNXVcajWMvw+CQJwHm+bCVhQUd3c6Bq0V7u+5ShoyipDSUTSZMrHAbbD6KjiQOlzGdKy69HbotpmJWO3+jfA/oiJF4jjziQ4OHMKaxHg77Jqu0iUQUL9iobonHYNRMK6VhtwcVav3keEqUfFGLIkm2qjJROmw/4SyEhMT/03KhB1Mm18Ks9Wn6kZ9E7ezhWowuJ9SUNpFj9qVt/yKPXDCYmI7/qInFokGdd+EUf9OHHdVtAsGW2q/TeB95EJqXEBEk47BqHn1SDgKt8VpCLr8TACI4c7+P1IVnPikA4fiBhcQQ40wZPHAReW1BGloenOHV3k5+xaOYegdScQb2eI/OOXGwgpxL6JdwhX9Qn6b/tcLlT/0d9Wp+KeTiMELuZ1osav+NjMhINwhNvQ/+uNdWRZlo/83VWksZlQHyo6CZj5GUr2o21tZjs4EfYPFCuaMWhFQuzG3cSfC8P8jbUiZfb37eUdoCGfnZgS4qLx0C1Rz0FX4JEUKJgGIKIoUiVtDMsGCANF/IVHQdzJX6UXPQpFCcYr4mf8aXykzq978qBILySdrLOxeyyT8ITYsatibf6XvxBTkBAm37yMGvBdXVcbwtX4x/KyPmD8LziR9Y4GZFmLUC/e96RjMiumlEh+1R5oKTSfWu8VaCZDznkwKmzFq80Aozymw4ub7i0wkW7usgUkOw0SXJu2cqAZP/8yOANxa6wUv8m87RtLmN6rwIpJD/wYRgskWWVO+3AxWpKkx1BiwKQG7z4fAJIR0vuIdI+W4Mj3cQkT1ZvGwRUWKRDzhcOC6qcf1MXrrqmrZeThfpb3KSlyVfrT5sZePNAXJS5fJtdVDWBO6Lsj25W9Cuo72hCt7WleUOElveg+opV5JkWJrRIgJGcOeHAbFo+QWHPwzLwJxXLCAyReg7MkXym0iZtdt0o82P0ovf8SsWBsj0IrAzVZfe06CACtb2bTXUZcmigb5hBDCxMkCMWLjJf8Kw15drQZsm19dIZZyi/VicRCUJyj/UnlPkY/aCsKPhSPF8Ezi3uDgn3kRyNoAYpo0x8M2gT5PLniJfa5W+tHmN28tnPpeEDBHOnNLZkQDQpn+17JvrhCYVe4BdVVA9Ya4MLHtVZECDpT6YeBMH4U3/yyJQK/D+XEM1ko/2vyWrIjz2jYCJqTzt19QkVY2rIZPjAYyk3HTtYCr3ANSmQulh5MK1TvsUZEi7JMKCxYK/HnDiUShHgq7ToWAMO9UeIuEsiT9ULwkCYGbrfWbqozKn74fuGY9q7d/hWwU5uTS+pR/m18K4+flETAhnbkNNGgfazDAnSF+rbu5iMH1tFoMwsoNm0/Qq9gDUp2sMfhBlJvEu9Um9Pc6EGg7RtLmN7r0GsMoKzXe/iV/FtGMaZ5XsteywbT5pTB+rgcB75Eu0BYiPIkYMmCqppabURwGFRZi6j0gQFiXScpFaHWm9l1XCV2a8J+fLPwyW/yDitwgcrV+E0GHFCoskkUYWXCVLltR3unPP95W82vzq4b19+URMCFdrg0YRKw+C6PBBZGEWDaZtAdUXJmogLiVTEwn7QGloxkhTJtfKRF/DEJAE10iniyO0j/CDErDgXePAOM7jUcufgh/o7b7Wh+wgrcOWOdLVZl90OqZNCbhpv3R4nagtgJHQlm7z9Pm15ZmXz+l37n/0zetjYZjkmSPNBHVjVbDxZ4JAcY7/3ZVe/vXTHk62QsgYEK6HOhMttxOku+TwqFCBMeYtn2eNr8xeabD6o37Pylx1XfPihSJkCYxb6q2n0ag10LYMO0DARPShdqRPQ8RFXLj8H7aKxl9rZzSLfZ8qlVp86uGPeOb/Z8HxFO9TvZ/5Aa3vXdFChZHXwnnNvE8ENkYASOwYwS8R7ps4zLxhn1SEZqu/dFlSzY8t9b9HxYOshwzeFtNus2vGnbl3+yNmhtdeSOttXgsQGUL/Qe9o+UbTJtfCuPnehAwR7psW8CJwqlheJ4cewk+2/jx/s97yQKLIxsjMBgBLSjpO9iTf31q8xuckSPMjkBOSF9qFVTN8A816ImWaDWQv3sjAAfDChQDZxoOhoevjf3AVarISUS9sdJPU1xhgEjXYt1p4HQqRmC1CIg2ojCGYmGtgZCiMNE0oXu1XQvb2Y6B8KhR4EbZH/UkfDaUjmgEjIARWAwBlELZzqo1t+JkfiJaqA1tx1EIgLWIKAsXxKI8bYyAETACRmDlCGjubr0n3cpGyzcgXClKBY3nR5cv0vQ5asFgRYrpYXWKRsAIrBABrsdaYbH2WyQRGIgoZys/idKA/VbWNTMCRsAIHAABc6TLN3LQ1DMRXR5452gEjIARmAMBc6RzoOo0jYARMAJG4DAI/D+S4T+Tk6F+qgAAAABJRU5ErkJggg==",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\sqrt{a_{11}} & a_{12} & a_{13}\\\\\\frac{a_{12}}{\\sqrt{a_{11}}} & \\sqrt{a_{22} - \\frac{a_{12}^{2}}{a_{11}}} & a_{23}\\\\\\frac{a_{13}}{\\sqrt{a_{11}}} & \\frac{a_{23} - \\frac{a_{12} a_{13}}{a_{11}}}{\\sqrt{a_{22} - \\frac{a_{12}^{2}}{a_{11}}}} & \\sqrt{a_{33} - \\frac{\\left(a_{23} - \\frac{a_{12} a_{13}}{a_{11}}\\right)^{2}}{a_{22} - \\frac{a_{12}^{2}}{a_{11}}} - \\frac{2 a_{13}^{2}}{a_{11}}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡  _____                                                                ⎤\n",
       "⎢╲╱ a₁₁          a₁₂                             a₁₃                    ⎥\n",
       "⎢                                                                       ⎥\n",
       "⎢              ____________                                             ⎥\n",
       "⎢             ╱          2                                              ⎥\n",
       "⎢  a₁₂       ╱        a₁₂                                               ⎥\n",
       "⎢───────    ╱   a₂₂ - ────                       a₂₃                    ⎥\n",
       "⎢  _____  ╲╱          a₁₁                                               ⎥\n",
       "⎢╲╱ a₁₁                                                                 ⎥\n",
       "⎢                                                                       ⎥\n",
       "⎢                                      _________________________________⎥\n",
       "⎢                                     ╱                      2          ⎥\n",
       "⎢                 a₁₂⋅a₁₃            ╱        ⎛      a₁₂⋅a₁₃⎞           ⎥\n",
       "⎢           a₂₃ - ───────           ╱         ⎜a₂₃ - ───────⎟         2 ⎥\n",
       "⎢  a₁₃              a₁₁            ╱          ⎝        a₁₁  ⎠    2⋅a₁₃  ⎥\n",
       "⎢───────  ─────────────────       ╱     a₃₃ - ──────────────── - ────── ⎥\n",
       "⎢  _____       ____________      ╱                        2       a₁₁   ⎥\n",
       "⎢╲╱ a₁₁       ╱          2      ╱                      a₁₂              ⎥\n",
       "⎢            ╱        a₁₂      ╱                 a₂₂ - ────             ⎥\n",
       "⎢           ╱   a₂₂ - ────   ╲╱                        a₁₁              ⎥\n",
       "⎣         ╲╱          a₁₁                                               ⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# implementação de algoritmo simbólico \n",
    "# para a decomposição de Cholesky \n",
    "\n",
    "B = A[:,:] # faz cópia da matriz A\n",
    "\n",
    "for k in range(0,n):\n",
    "    for i in range(0,k):\n",
    "        s = 0.\n",
    "        for j in range(0,i):\n",
    "            s += B[i,j]*B[k,j]\n",
    "        B[k,i] = (B[k,i] - s)/B[i,i]\n",
    "    s = 0.\n",
    "    for j in range(0,k):\n",
    "        s += s + B[k,j]*B[k,j]\n",
    "    B[k,k] = sp.sqrt(B[k,k] - s)\n",
    "\n",
    "# saída \n",
    "B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tarefa \n",
    "\n",
    "Converta o código simbólico acima para uma versão numérica (ou implemente a sua própria versão) e aplique-o na matriz abaixo para encontrar o fator de Cholesky: \n",
    "\n",
    "$$\\textbf{A} = \n",
    "\\begin{bmatrix}\n",
    "6 & 15 & 55 \\\\\n",
    "15 & 55 & 225 \\\\\n",
    "55 & 225 & 979 \n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cálculo do fator de Cholesky com Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matriz A = \n",
      " [[16 -4 12 -4]\n",
      " [-4  2 -1  1]\n",
      " [12 -1 14 -2]\n",
      " [-4  1 -2 83]]\n",
      "Matriz L = \n",
      " [[ 4.  0.  0.  0.]\n",
      " [-1.  1.  0.  0.]\n",
      " [ 3.  2.  1.  0.]\n",
      " [-1.  0.  1.  9.]]\n",
      "Matriz L^T = \n",
      " [[ 4. -1.  3. -1.]\n",
      " [ 0.  1.  2.  0.]\n",
      " [ 0.  0.  1.  1.]\n",
      " [ 0.  0.  0.  9.]]\n",
      "Matriz LL^T = \n",
      " [[16. -4. 12. -4.]\n",
      " [-4.  2. -1.  1.]\n",
      " [12. -1. 14. -2.]\n",
      " [-4.  1. -2. 83.]]\n",
      "Norma || A - LL^T || =  0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAECCAYAAADXWsr9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAIoklEQVR4nO3dwYuc9R3H8c+na0SpaYU0SEhC9SAW8WDokIvQQ7CQetAezUEsCDkJEbx49R/w1kvA0BRFEeJBJEX2EJGAxmxCtMbVEqQli5ZEgzVBqUS+HnYOaXd1npjnN8/OfN4vGNiZDM/zfbJ555lndpKfq0oA5tvPhh4AQHuEDgQgdCAAoQMBCB0IQOhAgJkO3fZe2x/bPmf7maHn6ZPtQ7Yv2P5g6FlasL3T9jHby7bP2j4w9Ex9sX2L7Xdtvzc+tmcHn2lWf45ue0HSPyT9XtKKpJOS9lXVh4MO1hPbv5N0RdJfq+q+oefpm+1tkrZV1WnbmyWdkvTHefj+2bakn1fVFdubJB2XdKCq3hlqplk+o++WdK6qPqmqbyW9LOmRgWfqTVW9JenS0HO0UlWfVdXp8deXJS1L2j7sVP2oVVfGdzeNb4OeUWc59O2Szl9zf0Vz8gclje07Je2SdGLYSfpje8H2GUkXJC1W1aDHNsuhe53HZvM6JJjt2yQdkfRUVX019Dx9qarvqup+STsk7bY96OXXLIe+ImnnNfd3SPp0oFnwE4yvX49IerGqXh16nhaq6ktJb0raO+Qcsxz6SUl3277L9s2SHpX02sAzoaPxG1bPS1ququeGnqdPtrfavn389a2SHpT00ZAzzWzoVXVV0pOS3tDqGzmvVNXZYafqj+2XJL0t6R7bK7afGHqmnj0g6TFJe2yfGd8eGnqonmyTdMz2+1o9IS1W1etDDjSzP14D0N3MntEBdEfoQABCBwIQOhCA0IEAMx+67f1Dz9ASxzfbNsrxzXzokjbEb2RDHN9s2xDHNw+hA5igyQdm7F+UdEfv213ffyT9ckr7WvVbnZvavi5K2jq1va36Yor7uixp8xT3J0lbprivaX///inp86o1/+Drpja7u0PSXH18+X8s6eGhR2jq8NADNPb40AM0NPqBx3npDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAnQK3fZe2x/bPmf7mdZDAejXxNBtL0j6s6Q/SLpX0j7b97YeDEB/upzRd0s6V1WfVNW3kl6W9EjbsQD0qUvo2yWdv+b+yvgxADOiy5JMa9ZxkrRmwbbx8rDjlSOnvVoYgB/T5Yy+ImnnNfd3SPr0/59UVQeralRVo2kvegjgx3UJ/aSku23fZftmSY9Keq3tWAD6NPGle1Vdtf2kpDckLUg6VFVnm08GoDedlk2uqqOSjjaeBUAjfDIOCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EcNWa1ZVu2Miupd63uoE8/JuhJ2jKi8tDj9DU09+st8rYfHhB0r+r1hwgZ3QgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EmBi67UO2L9j+YBoDAehflzP6XyTtbTwHgIYmhl5Vb0m6NIVZADTCNToQoLfQbe+3vWR76WJfGwXQi95Cr6qDVTWqqtHWvjYKoBe8dAcCdPnx2kuS3pZ0j+0V20+0HwtAn26a9ISq2jeNQQC0w0t3IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQYOJ/9/xTfCHpcIsNbxB/WlweeoSm6tL5oUdo6tTXNfQIzRzdM1r3cc7oQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCDAxdNs7bR+zvWz7rO0D0xgMQH+6rNRyVdLTVXXa9mZJp2wvVtWHjWcD0JOJZ/Sq+qyqTo+/vixpWdL21oMB6M91XaPbvlPSLkknWgwDoI3Oiyzavk3SEUlPVdVX6/z6fkn7JWlLb+MB6EOnM7rtTVqN/MWqenW951TVwaoaVdVoc58TArhhXd51t6TnJS1X1XPtRwLQty5n9AckPSZpj+0z49tDjecC0KOJ1+hVdVySpzALgEb4ZBwQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQjQeUmm67FF0uMtNrxB/P2b+f7fr099XUOP0NRoy3x//9bDGR0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBJoZu+xbb79p+z/ZZ289OYzAA/emyUst/Je2pqiu2N0k6bvtvVfVO49kA9GRi6FVVkq6M724a3+Z7zR5gznS6Rre9YPuMpAuSFqvqxDrP2W97yfbSxb6nBHBDOoVeVd9V1f2Sdkjabfu+dZ5zsKpGVTXa2veUAG7Idb3rXlVfSnpT0t4m0wBoosu77ltt3z7++lZJD0r6qPVgAPrT5V33bZIO217Q6l8Mr1TV623HAtCnLu+6vy9p1xRmAdAIn4wDAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBvLqGYs8btS9K+lfvG17fryR9PqV9DYHjm23TPr5fV9WaVdGahD5NtpeqajT0HK1wfLNtoxwfL92BAIQOBJiH0A8OPUBjHN9s2xDHN/PX6AAmm4czOoAJCB0IQOhAAEIHAhA6EOB74AmxxeqVwyoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAECCAYAAADXWsr9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAIdUlEQVR4nO3dQYhd5R2G8fdtjESqQayzCEmoQkUqLhI7ZCNUCBZSN7pUrCshKyFCN26le3duAoa2IIoQKSIWySJGAhozSaMYR0uwFKcKiQSJWbQS+XcxtxCaae+J+b575s77/GBg7jjceY/65Nx7ZzLHVSUAG9uPxh4AoD9CBwIQOhCA0IEAhA4EIHQgwFyHbnuf7c9sn7P93Nh7WrJ9yPZ52x+PvaUH2zttH7W9bPus7QNjb2rF9hbbH9j+cHJsz4++aV6/j257k6S/SvqVpBVJJyU9UVWfjDqsEdu/lHRZ0h+r6v6x97Rme5ukbVV12vZtkk5Jemwj/PezbUk/rqrLtjdLOi7pQFW9P9ameT6j75F0rqo+r6rvJL0q6dGRNzVTVe9Kujj2jl6q6quqOj15/1tJy5K2j7uqjVp1eXJz8+Rt1DPqPIe+XdIXV91e0Qb5HyWN7bsk7ZZ0Ytwl7djeZPuMpPOSjlTVqMc2z6F7jY/N5/OQYLZvlXRY0rNVdWnsPa1U1fdVtUvSDkl7bI/69GueQ1+RtPOq2zskfTnSFvwAk+evhyW9XFWvj72nh6r6RtI7kvaNuWOeQz8p6R7bd9u+WdLjkt4YeRMGmrxg9ZKk5ap6Yew9LdlesH375P1bJD0s6dMxN81t6FV1RdIzkt7W6gs5r1XV2XFXtWP7FUnvSbrX9ortp8fe1NiDkp6StNf2mcnbI2OPamSbpKO2P9LqCelIVb055qC5/fYagOHm9owOYDhCBwIQOhCA0IEAhA4EmPvQbe8fe0NPHN98Wy/HN/ehS1oX/yI74vjm27o4vo0QOoApuvzAzG12LTS/17VdkrR1Rl/rP/6mWf79hIuS7pjh15s1jq+tf6jq4jV/4eumHl9qQdLvetzxOvEb/WnsCcD/8NiaH+WhOxCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCDAodNv7bH9m+5zt53qPAtDW1NBtb5L0oqRfS7pP0hO27+s9DEA7Q87oeySdq6rPq+o7Sa9KerTvLAAtDQl9u6Qvrrq9MvkYgDkx5JJM11zHSdI1F2ybXB52vyT95AZHAWhryBl9RdLOq27vkPTlf39SVR2sqsWqWpz1RQ8B/H9DQj8p6R7bd9u+WdLjkt7oOwtAS1MfulfVFdvPSHpb0iZJh6rqbPdlAJoZdNnkqnpL0ludtwDohJ+MAwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQb9uufrdYekJ3vc8Trx5OljY0/oyg88NPYENMYZHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwGmhm77kO3ztj+exSAA7Q05o/9e0r7OOwB0NDX0qnpX0sUZbAHQCc/RgQDNQre93/aS7aULre4UQBPNQq+qg1W1WFWLC63uFEATPHQHAgz59torkt6TdK/tFdtP958FoKWbpn1CVT0xiyEA+uGhOxCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCDD11z3/EP/c/QstH1vqcdfrwn1bd409oas6PfaCvvzAQ2NPmDnO6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQgwNXTbO20ftb1s+6ztA7MYBqCdIVdquSLpt1V12vZtkk7ZPlJVn3TeBqCRqWf0qvqqavUiPVX1raRlSdt7DwPQznU9R7d9l6Tdkk70GAOgj8Gh275V0mFJz1bVpTX++X7bS7aXLn59oeVGADdoUOi2N2s18per6vW1PqeqDlbVYlUt3nHnQsuNAG7QkFfdLeklSctV9UL/SQBaG3JGf1DSU5L22j4zeXuk8y4ADU399lpVHZfkGWwB0Ak/GQcEIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIMuZrqddvyl1P6+daN/Buiz409oLNjYw/o7GdjD5g5zuhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IMDV021tsf2D7Q9tnbT8/i2EA2hlypZZ/SdpbVZdtb5Z03Pafq+r9ztsANDI19KoqSZcnNzdP3qrnKABtDXqObnuT7TOSzks6UlUn1vic/baXbC9daL0SwA0ZFHpVfV9VuyTtkLTH9v1rfM7BqlqsqsWF1isB3JDretW9qr6R9I6kfV3WAOhiyKvuC7Zvn7x/i6SHJX3aexiAdoa86r5N0h9sb9LqHwyvVdWbfWcBaGnIq+4fSdo9gy0AOuEn44AAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQACvXkOx8Z3aFyT9vfkdr+1OSV/P6GuNgeObb7M+vp9W1TVXResS+izZXqqqxbF39MLxzbf1cnw8dAcCEDoQYCOEfnDsAZ1xfPNtXRzf3D9HBzDdRjijA5iC0IEAhA4EIHQgAKEDAf4NiPSlJyxmsmsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAECCAYAAADXWsr9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAIjElEQVR4nO3dQYic9R3G8efpuhKpEYnmEJJQLRWpWDDpkotQIVjYetGjQTwJOQkRevEqvXso9LJgaAuiCJEiYpEcohLQmE0arXG1BGlxUUhCkCSHKpFfDztCMIPz2n3/8+7s8/3AQmYdZp838es7M7vxdVUJwOb2k6EHAGiP0IEAhA4EIHQgAKEDAQgdCDDTodtetP2p7XO2nx16T59sH7Z93vZHQ29pwfZu28dsr9g+a/vQ0Jv6YnuL7fdtfzA6tucG3zSr30e3PSfpX5J+K2lV0klJB6rq40GH9cT2byRdlfTXqrp/6D19s71D0o6qOm17q6RTkh7bDH9+ti3pp1V11fa8pOOSDlXVe0NtmuUz+j5J56rqs6r6RtLLkh4deFNvquodSZeG3tFKVX1ZVadHv74iaUXSzmFX9aPWXB3dnB99DHpGneXQd0r6/Lrbq9ok/6KksX2XpD2STgy7pD+252yfkXRe0tGqGvTYZjl0j/ncbL4OCWb7VklHJD1TVZeH3tOXqvq2qh6QtEvSPtuDvvya5dBXJe2+7vYuSV8MtAX/h9Hr1yOSXqyqV4fe00JVfSXpLUmLQ+6Y5dBPSrrH9t22b5b0uKTXBt6EjkZvWL0gaaWqnh96T59sb7d9++jXt0h6WNInQ26a2dCr6pqkpyW9qbU3cl6pqrPDruqP7ZckvSvpXturtp8aelPPHpT0pKT9ts+MPh4ZelRPdkg6ZvtDrZ2QjlbV60MOmtlvrwHobmbP6AC6I3QgAKEDAQgdCEDoQICZD932waE3tMTxzbaNcnwzH7qkDfEb2RDHN9s2xPFthtABTNDkB2a22rW990cd77Kk26b0tb6zbYpf64Kkaf1efue/e349ta916eIFbbtzuke45R+npva1pv3n929JF6tu+AtfTUL/uV1/6P1RN44nhh7Q2Mrlzf3Tkr+8bdxffNwcFiQtjwmdp+5AAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCAAoQMBCB0IQOhAAEIHAhA6EIDQgQCEDgQgdCBAp9BtL9r+1PY528+2HgWgXxNDtz0n6U+SfifpPkkHbN/XehiA/nQ5o++TdK6qPquqbyS9LOnRtrMA9KlL6DslfX7d7dXR5wDMiJs63Gfc9WtuuGbP6PKwByXpjnWOAtCvLmf0VUm7r7u9S9IX379TVS1V1UJVLUz7oocAfliX0E9Kusf23bZvlvS4pNfazgLQp4lP3avqmu2nJb0paU7S4ao623wZgN50eY2uqnpD0huNtwBohJ+MAwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAVx1w9WV1v+g/lVJf+v9cTeKOv320BOa8t4/Dj2hsSNDD2joMVX984bLqHFGBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQICJods+bPu87Y+mMQhA/7qc0f8sabHxDgANTQy9qt6RdGkKWwA0wmt0IEBvods+aHvZ9jJPAICNpbfQq2qpqhaqakHa1tfDAugBT92BAF2+vfaSpHcl3Wt71fZT7WcB6NNNk+5QVQemMQRAOzx1BwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDASb+755xI+99aOgJTdXpoRe09vbQA5pZeOLK2M9zRgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6ECAiaHb3m37mO0V22dtH5rGMAD96XKllmuSfl9Vp21vlXTK9tGq+rjxNgA9mXhGr6ovq9Yu0lNVVyStSNrZehiA/vyo1+i275K0R9KJFmMAtNE5dNu3Sjoi6Zmqujzmnx+0vWx7WbrU50YA69QpdNvzWov8xap6ddx9qmqpqhaqakHa1udGAOvU5V13S3pB0kpVPd9+EoC+dTmjPyjpSUn7bZ8ZfTzSeBeAHk389lpVHZfkKWwB0Ag/GQcEIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwJ0uZoqwnjvQ0NPaOwXQw+YOs7oQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCEDoQABCBwIQOhCA0IEAhA4EIHQgAKEDAQgdCDAxdNtbbL9v+wPbZ20/N41hAPrT5UotX0vaX1VXbc9LOm7771X1XuNtAHoyMfSqKklXRzfnRx/VchSAfnV6jW57zvYZSeclHa2qE2Puc9D2su1l6VLfOwGsQ6fQq+rbqnpA0i5J+2zfP+Y+S1W1UFUL0ra+dwJYhx/1rntVfSXpLUmLTdYAaKLLu+7bbd8++vUtkh6W9EnrYQD60+Vd9x2S/mJ7Tmv/YXilql5vOwtAn7q86/6hpD1T2AKgEX4yDghA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBCB0IAChAwEIHQhA6EAAQgcCEDoQgNCBAIQOBPDaNRR7flD7gqT/9P7A490p6eKUvtYQOL7ZNu3j+1lVbf/+J5uEPk22l9eu97Y5cXyzbaMcH0/dgQCEDgTYDKEvDT2gMY5vtm2I45v51+gAJtsMZ3QAExA6EIDQgQCEDgQgdCDA/wBtgrRyeAjOHgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy import array, linalg, dot, zeros\n",
    "\n",
    "# matriz\n",
    "A = array([[16, -4, 12, -4],\n",
    "           [-4, 2, -1, 1],\n",
    "           [12, -1, 14, -2],\n",
    "           [-4, 1, -2, 83]])\n",
    "\n",
    "# fator de Cholesky do Scipy\n",
    "L = linalg.cholesky(A, lower=True, overwrite_a=False, check_finite=True)\n",
    "\n",
    "# fator de Cholesky implementado\n",
    "n = A.shape[0]\n",
    "G = zeros(A.shape, dtype=float)\n",
    "\n",
    "print('Matriz A = \\n', A)\n",
    "print('Matriz L = \\n', L)\n",
    "print('Matriz L^T = \\n', L.T)\n",
    "\n",
    "# prova real por produto interno\n",
    "A2 = dot(L, L.T)\n",
    "print('Matriz LL^T = \\n', A2)\n",
    "\n",
    "# prova real usando norma de Frobenius da diferenca de matrizes\n",
    "print('Norma || A - LL^T || = ', linalg.norm(A-A2))\n",
    "\n",
    "plt.matshow(A,cmap=plt.cm.flag);\n",
    "plt.matshow(L,cmap=plt.cm.flag);\n",
    "plt.matshow(L.T,cmap=plt.cm.flag);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}