Linear_Algebra-project (1).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CompMech04-Linear Algebra Project\n",
    "# Practical Linear Algebra for Finite Element Analysis\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project we will perform a linear-elastic finite element analysis (FEA) on a support structure made of 11 beams that are riveted in 7 locations to create a truss as shown in the image below. \n",
    "\n",
    "![Mesh image of truss](../images/mesh.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The triangular truss shown above can be modeled using a [direct stiffness method [1]](https://en.wikipedia.org/wiki/Direct_stiffness_method), that is detailed in the [extra-FEA_material](./extra-FEA_material.ipynb) notebook. The end result of converting this structure to a FE model. Is that each joint, labeled $n~1-7$, short for _node 1-7_ can move in the x- and y-directions, but causes a force modeled with Hooke's law. Each beam labeled $el~1-11$, short for _element 1-11_, contributes to the stiffness of the structure. We have 14 equations where the sum of the components of forces = 0, represented by the equation\n",
    "\n",
    "$\\mathbf{F-Ku}=\\mathbf{0}$\n",
    "\n",
    "Where, $\\mathbf{F}$ are externally applied forces, $\\mathbf{u}$ are x- and y- displacements of nodes, and $\\mathbf{K}$ is the stiffness matrix given in `fea_arrays.npz` as `K`, shown below\n",
    "\n",
    "_note: the array shown is 1000x(`K`). You can use units of MPa (N/mm^2), N, and mm. The array `K` is in 1/mm_\n",
    "\n",
    "$\\mathbf{K}=EA*$\n",
    "\n",
    "$  \\left[ \\begin{array}{cccccccccccccc}\n",
    " 4.2 & 1.4 & -0.8 & -1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " 1.4 & 2.5 & -1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " -0.8 & -1.4 & 5.0 & 0.0 & -0.8 & 1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " -3.3 & 0.0 & -0.8 & 1.4 & 8.3 & 0.0 & -0.8 & -1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " 0.0 & 0.0 & 1.4 & -2.5 & 0.0 & 5.0 & -1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & -1.4 & 8.3 & 0.0 & -0.8 & 1.4 & -3.3 & 0.0 & 0.0 & 0.0 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & 1.4 & 8.3 & 0.0 & -0.8 & -1.4 & -3.3 & 0.0 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.4 & -2.5 & 0.0 & 5.0 & -1.4 & -2.5 & 0.0 & 0.0 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & -1.4 & 5.0 & 0.0 & -0.8 & 1.4 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & 1.4 & 4.2 & -1.4 \\\\\n",
    " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.4 & -2.5 & -1.4 & 2.5 \\\\\n",
    "\\end{array}\\right]~\\frac{1}{m}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 373,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00416667,  0.00144338, -0.00083333, -0.00144338, -0.00333333,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.00144338,  0.0025    , -0.00144338, -0.0025    ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [-0.00083333, -0.00144338,  0.005     ,  0.        , -0.00083333,\n",
       "         0.00144338, -0.00333333,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [-0.00144338, -0.0025    ,  0.        ,  0.005     ,  0.00144338,\n",
       "        -0.0025    ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [-0.00333333,  0.        , -0.00083333,  0.00144338,  0.00833333,\n",
       "         0.        , -0.00083333, -0.00144338, -0.00333333,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.00144338, -0.0025    ,  0.        ,\n",
       "         0.005     , -0.00144338, -0.0025    ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        , -0.00333333,  0.        , -0.00083333,\n",
       "        -0.00144338,  0.00833333,  0.        , -0.00083333,  0.00144338,\n",
       "        -0.00333333,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        , -0.00144338,\n",
       "        -0.0025    ,  0.        ,  0.005     ,  0.00144338, -0.0025    ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        , -0.00333333,\n",
       "         0.        , -0.00083333,  0.00144338,  0.00833333,  0.        ,\n",
       "        -0.00083333, -0.00144338, -0.00333333,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.00144338, -0.0025    ,  0.        ,  0.005     ,\n",
       "        -0.00144338, -0.0025    ,  0.        ,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        , -0.00333333,  0.        , -0.00083333, -0.00144338,\n",
       "         0.005     ,  0.        , -0.00083333,  0.00144338],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        , -0.00144338, -0.0025    ,\n",
       "         0.        ,  0.005     ,  0.00144338, -0.0025    ],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        , -0.00333333,  0.        ,\n",
       "        -0.00083333,  0.00144338,  0.00416667, -0.00144338],\n",
       "       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "         0.00144338, -0.0025    , -0.00144338,  0.0025    ]])"
      ]
     },
     "execution_count": 373,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "fea_arrays = np.load('./fea_arrays.npz')\n",
    "K=fea_arrays['K']\n",
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project we are solving the problem, $\\mathbf{F}=\\mathbf{Ku}$, where $\\mathbf{F}$ is measured in Newtons, $\\mathbf{K}$ `=E*A*K` is the stiffness in N/mm, `E` is Young's modulus measured in MPa (N/mm^2), and `A` is the cross-sectional area of the beam measured in mm^2. \n",
    "\n",
    "There are three constraints on the motion of the joints:\n",
    "\n",
    "i. node 1 displacement in the x-direction is 0 = `u[0]`\n",
    "\n",
    "ii. node 1 displacement in the y-direction is 0 = `u[1]`\n",
    "\n",
    "iii. node 7 displacement in the y-direction is 0 = `u[13]`\n",
    "\n",
    "We can satisfy these constraints by leaving out the first, second, and last rows and columns from our linear algebra description. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Calculate the condition of `K` and the condition of `K[2:13,2:13]`. \n",
    "\n",
    "a. What error would you expect when you solve for `u` in `K*u = F`? \n",
    "\n",
    "b. Why is the condition of `K` so large?\n",
    "\n",
    "c. What error would you expect when you solve for `u[2:13]` in `K[2:13,2:13]*u=F[2:13]`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 374,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition K 1.4577532625238035e+17\n",
      "Condition K[2:13,2:13] 52.23542514351006\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\"a. There would be a linear algebra error since K is a singular matrix\\nb. Its large because the rows of K are almost dependent\\nc. There technically shouldn't be an error though the results won't be good\""
      ]
     },
     "execution_count": 374,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "condK=np.linalg.cond(K)\n",
    "condK2=np.linalg.cond(K[2:13,2:13])\n",
    "print('Condition K',condK)\n",
    "print('Condition K[2:13,2:13]',condK2)\n",
    "'''a. There would be a linear algebra error since K is a singular matrix\n",
    "b. Its large because the rows of K are almost dependent\n",
    "c. There technically shouldn't be an error though the results won't be good'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Apply a 100-N downward force to the central top node (n 4)\n",
    "\n",
    "a. Create the LU matrix for K[2:13,2:13]\n",
    "\n",
    "b. Use cross-sectional area of $0.1~mm^2$ and steel and almuminum moduli, $E=200~GPa~and~E=70~GPa,$ respectively. Solve the forward and backward substitution methods for \n",
    "\n",
    "* $\\mathbf{Ly}=\\mathbf{F}\\frac{1}{EA}$\n",
    "\n",
    "* $\\mathbf{Uu}=\\mathbf{y}$\n",
    "\n",
    "_your array `F` is zeros, except for `F[5]=-100`, to create a -100 N load at node 4._\n",
    "\n",
    "c. Plug in the values for $\\mathbf{u}$ into the full equation, $\\mathbf{Ku}=\\mathbf{F}$, to solve for the reaction forces\n",
    "\n",
    "d. Create a plot of the undeformed and deformed structure with the displacements and forces plotted as vectors (via `quiver`). Your result for aluminum should match the following result from [extra-FEA_material](./extra-FEA_material.ipynb). _note: The scale factor is applied to displacements $\\mathbf{u}$, not forces._\n",
    "\n",
    "![Deformed structure with loads applied](../images/deformed_truss.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 375,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L=\n",
      "\n",
      " [   1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [  -0.17   0.29   1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.29  -0.50   0.12   1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [  -0.67   0.00  -0.18  -0.10   1.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00  -0.19  -0.72  -0.08   1.00   0.00   0.00   0.00   0.00   0.00]\n",
      " [   0.00   0.00  -0.43   0.13  -0.24   0.33   1.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.25  -0.79   0.18   1.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00  -0.57  -0.09  -0.25  -0.22   1.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00  -0.23  -0.88  -0.33   1.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00  -0.54   0.24  -0.59   0.29   1.00] \n",
      "\n",
      "U=\n",
      "\n",
      " [   0.00   0.00  -0.00   0.00  -0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.01   0.00  -0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.01   0.00  -0.00  -0.00  -0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00  -0.00  -0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.01  -0.00  -0.00   0.00  -0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.00  -0.00  -0.00   0.00   0.00]\n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.01   0.00  -0.00  -0.00  -0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  -0.00  -0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  -0.00  -0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      " [   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import lu\n",
    "import numpy as np\n",
    "from pretty_print import array_print as AP \n",
    "import matplotlib.pyplot as plt\n",
    "'''A'''\n",
    "K_lu=K[2:13,2:13]\n",
    "P,L,U=lu(K_lu)\n",
    "print('L=\\n')\n",
    "AP([L],['[L]'])\n",
    "print('U=\\n')\n",
    "AP([U],['[U]'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 376,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u steel=\n",
      "-----------\n",
      " [   0.00] \n",
      " [   0.00] \n",
      " [   1.95] \n",
      " [  -2.13] \n",
      " [   0.43] \n",
      " [  -4.00] \n",
      " [   1.08]\n",
      " [  -5.38] \n",
      " [   1.73] \n",
      " [  -4.00] \n",
      " [   0.22] \n",
      " [  -2.13] \n",
      " [   2.17] \n",
      " [   0.00] \n",
      "\n",
      "u aluminum=\n",
      "-----------\n",
      " [   0.00] \n",
      " [   0.00] \n",
      " [   5.57] \n",
      " [  -6.07] \n",
      " [   1.24] \n",
      " [ -11.43] \n",
      " [   3.09]\n",
      " [ -15.36] \n",
      " [   4.95] \n",
      " [ -11.43] \n",
      " [   0.62] \n",
      " [  -6.07] \n",
      " [   6.19] \n",
      " [   0.00] \n",
      "\n"
     ]
    }
   ],
   "source": [
    "'''B'''\n",
    "l=300 # mm\n",
    "nodes=np.array([[1,0,0],[2,0.5,3**0.5/2],[3,1,0],[4,1.5,3**0.5/2],[5,2,0],[6,2.5,3**0.5/2],[7,3,0]])\n",
    "nodes[:,1:3]*=l\n",
    "elems=np.array([[1,1,2],[2,2,3],[3,1,3],[4,2,4],[5,3,4],[6,3,5],[7,4,5],[8,4,6],[9,5,6],[10,5,7],[11,6,7]])\n",
    "numnodes=7 #from truss\n",
    "num_nodes = 7\n",
    "A=0.1 #mm^2\n",
    "E_st=200e3 #N/mm^2\n",
    "E_al=70e3 #N/mm^2\n",
    "F = np.zeros(2*num_nodes-3)  # (N) forces on nodes\n",
    "F[5] =-100 \n",
    "#steel\n",
    "ly_st=F/(E_st*A) \n",
    "u_st_sub = solveLU(L,U,ly_st)\n",
    "u_st=np.zeros(num_nodes*2)\n",
    "u_st[2:13] = u_st_sub\n",
    "print('u steel=\\n-----------')\n",
    "AP([u_st], ['[u_st]'])\n",
    "#aluminum\n",
    "ly_al =F/(E_al*A)\n",
    "u_al_sub = solveLU(L,U,b_al)\n",
    "u_al=np.zeros(num_nodes*2)\n",
    "u_al[2:13] = u_al_sub\n",
    "print('u aluminum=\\n-----------')\n",
    "AP([u_al],['[u_al]'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 377,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reaction Force Steel=\n",
      "-------------------\n",
      " [-3.03576608e-18  2.50000000e-03  0.00000000e+00  4.33680869e-19\n",
      "  8.67361738e-19 -1.08420217e-18 -4.14613386e-19 -5.00000000e-03\n",
      " -8.67361738e-19 -1.08420217e-18  0.00000000e+00  0.00000000e+00\n",
      "  0.00000000e+00  2.50000000e-03]\n",
      "Reaction Force Aluminum=\n",
      "-------------------\n",
      " [-6.93889390e-18  7.14285714e-03  3.46944695e-18  2.16840434e-18\n",
      " -1.73472348e-18 -2.16840434e-18 -2.10100208e-18 -1.42857143e-02\n",
      "  3.46944695e-18 -1.30104261e-18 -8.67361738e-19 -5.20417043e-18\n",
      "  3.46944695e-18  7.14285714e-03]\n"
     ]
    }
   ],
   "source": [
    "'''C'''\n",
    "#F=Ku\n",
    "F_st=K@u_st\n",
    "F_al=K@u_al\n",
    "print('Reaction Force Steel=\\n-------------------\\n',F_st)\n",
    "print('Reaction Force Aluminum=\\n-------------------\\n',F_al)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''D'''\n",
    "from __future__ import print_function\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "def f(s,F,u,material):\n",
    "    plt.plot(r[ix],r[iy],'-',color=(0,0,0,1))\n",
    "    plt.plot(r[ix]+u[ix]*s,r[iy]+u[iy]*s,'-',color=(1,0,0,1))\n",
    "    #plt.quiver(r[ix],r[iy],u[ix],u[iy],color=(0,0,1,1),label='displacements')\n",
    "    plt.quiver(r[ix],r[iy],F[ix],F[iy],color=(1,0,0,1),label='applied forces')\n",
    "    plt.quiver(r[ix],r[iy],u[ix],u[iy],color=(0,0,1,1),label='displacements')\n",
    "    plt.axis(l*np.array([-0.5,3.5,-0.5,2]))\n",
    "    plt.xlabel('x (mm)')\n",
    "    plt.ylabel('y (mm)')\n",
    "    plt.title('Deformation scale = {:.1f}x'.format(s))\n",
    "    plt.legend(bbox_to_anchor=(1,0.5))\n",
    "ix=2*np.block([[np.arange(0,5)],[np.arange(1,6)],[np.arange(2,7)],[np.arange(0,5)]])\n",
    "iy=ix+1\n",
    "r=np.block([n[1:3] for n in nodes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9aff284ff71e4183a8540308d0a5ac9b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='s', max=20), Output()), _dom_classes=('widget-interact'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#steel\n",
    "interact(lambda s:f(s,F_st,u_st,'steel'),s=(0,20,1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 380,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "701ae0e92b7b47c3b6351e6960319858",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='s', max=20), Output()), _dom_classes=('widget-interact'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#aluminum\n",
    "interact(lambda s:f(s,F_al,u_al,'aluminum'), s=(0,20,1));\n",
    "#aluminum deforms more makes sense"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Determine cross-sectional area\n",
    "\n",
    "a. Using aluminum, what is the minimum cross-sectional area to keep total y-deflections $<0.2~mm$?\n",
    "\n",
    "b. Using steel, what is the minimum cross-sectional area to keep total y-deflections $<0.2~mm$?\n",
    "\n",
    "c. What are the weights of the aluminum and steel trusses with the chosed cross-sectional areas?\n",
    "\n",
    "d. The current price (2020/03) of [aluminum](https://tradingeconomics.com/commodity/aluminum) is 1545 dollars/Tonne and [steel](https://tradingeconomics.com/commodity/steel) is 476 dollars/Tonne [2]. Which material is cheaper to  create the truss?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.    0.    0.    0.    0. -100.    0.    0.    0.    0.    0.]\n",
      "displacements:\n",
      "----------------\n",
      "u_1y:0.00 mm\n",
      "u_1y:0.00 mm\n",
      "u_2y:0.17 mm\n",
      "u_2y:-0.19 mm\n",
      "u_3y:0.04 mm\n",
      "u_3y:-0.36 mm\n",
      "u_4y:0.10 mm\n",
      "u_4y:-0.48 mm\n",
      "u_5y:0.15 mm\n",
      "u_5y:-0.36 mm\n",
      "u_6y:0.02 mm\n",
      "u_6y:-0.19 mm\n",
      "u_7y:0.19 mm\n",
      "u_7y:0.00 mm\n",
      "\n",
      "\n",
      "\n",
      "3.2 mm^2 minimum of aluminum for total y-deflections<0.2mm\n"
     ]
    }
   ],
   "source": [
    "'''A'''\n",
    "#Aluminum\n",
    "E=70e3\n",
    "A=3.2\n",
    "Ff=np.zeros(2*len(nodes)-3)\n",
    "Ff[5]=-100\n",
    "print(Ff)\n",
    "# step 1 solve for uf (the joints without constraints)\n",
    "uf = np.linalg.solve(E*A*K[2:13,2:13],Ff)\n",
    "u=np.zeros(2*len(nodes))\n",
    "u[2:13]=uf\n",
    "xy={0:'y'}\n",
    "print('displacements:\\n----------------')\n",
    "for i in range(len(u)):\n",
    "    print('u_{}{}:{:.2f} mm'.format(int(i/2)+1,xy[i%1],u[i]))\n",
    "    \n",
    "print('\\n\\n\\n3.2 mm^2 minimum of aluminum for total y-deflections<0.2mm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.    0.    0.    0.    0. -100.    0.    0.    0.    0.    0.]\n",
      "displacements:\n",
      "----------------\n",
      "u_1y:0.00 mm\n",
      "u_1y:0.00 mm\n",
      "u_2y:0.17 mm\n",
      "u_2y:-0.19 mm\n",
      "u_3y:0.04 mm\n",
      "u_3y:-0.36 mm\n",
      "u_4y:0.10 mm\n",
      "u_4y:-0.48 mm\n",
      "u_5y:0.15 mm\n",
      "u_5y:-0.36 mm\n",
      "u_6y:0.02 mm\n",
      "u_6y:-0.19 mm\n",
      "u_7y:0.19 mm\n",
      "u_7y:0.00 mm\n",
      "\n",
      "\n",
      "\n",
      "1.12 mm^2 of steel minimum for total y-deflections<0.2mm\n"
     ]
    }
   ],
   "source": [
    "'''B'''\n",
    "#Steel\n",
    "E=200e3\n",
    "A=1.12\n",
    "Ff=np.zeros(2*len(nodes)-3)\n",
    "Ff[5]=-100\n",
    "print(Ff)\n",
    "# step 1 solve for uf (the joints without constraints)\n",
    "uf = np.linalg.solve(E*A*K[2:13,2:13],Ff)\n",
    "u=np.zeros(2*len(nodes))\n",
    "u[2:13]=uf\n",
    "xy={0:'y'}\n",
    "print('displacements:\\n----------------')\n",
    "for i in range(len(u)):\n",
    "    print('u_{}{}:{:.2f} mm'.format(int(i/2)+1,xy[i%1],u[i]))\n",
    "    \n",
    "print('\\n\\n\\n1.12 mm^2 of steel minimum for total y-deflections<0.2mm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.025521696000000003 N aluminum\n",
      "0.026039664000000007 N steel\n"
     ]
    }
   ],
   "source": [
    "'''C'''\n",
    "d_alum=0.00271 #g/mm^3\n",
    "d_st=0.0079 #g/mm^3\n",
    "#dv=m\n",
    "#v=a_c*L\n",
    "v_alum=3.2*l\n",
    "v_st=1.12*l\n",
    "m_alum=d_alum*v_alum\n",
    "m_st=d_st*v_st\n",
    "#w=g*m\n",
    "g=9810 #mm/s^2\n",
    "w_alum=m_alum*g\n",
    "w_st=m_st*g\n",
    "print(w_alum*1e-6,'N aluminum')\n",
    "print(w_st*1e-6, 'N steel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$ 0.04019472000000001 per aluminum beam\n",
      "$ 0.012634944000000004 per steel beam\n",
      "\n",
      "\n",
      "Steel is cheaper by \n",
      "$ 0.027559776000000008 per beam\n"
     ]
    }
   ],
   "source": [
    "'''D'''\n",
    "m_alum_tonne=m_alum*1e-5\n",
    "m_st_tonne=m_st*1e-5\n",
    "al_price=1545\n",
    "st_price=476\n",
    "price_al_beam=m_alum_tonne*al_price\n",
    "price_st_beam=m_st_tonne*st_price\n",
    "print('$',price_al_beam,'per aluminum beam')\n",
    "print('$',price_st_beam,'per steel beam\\n\\n')\n",
    "print('Steel is cheaper by \\n$',price_al_beam-price_st_beam,'per beam')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Future Predictions using past data\n",
    "\n",
    "The data from the price of aluminum and steel are in the data files `../data/al_prices.csv` and `../data/steel_price.csv`. If you're going to produce these supports for 5 years, how would you use the price history of steel and aluminum to compare how much manufacturing will cost?\n",
    "\n",
    "a. Separate the aluminum and steel data points into training and testing data (70-30%-split)\n",
    "\n",
    "b. Fit the training data to polynomial functions of order n. _Choose the highest order._\n",
    "\n",
    "c. Plot the error between your model and the training data and the error between your model and you testing data as a function of the polynomial function order, n. [Create the training-testing curves](../notebooks/03_Linear-regression-algebra.ipynb)\n",
    "\n",
    "d. Choose a polynomial function to predict the price of aluminum and steel in the year 2025. \n",
    "\n",
    "e. Based upon your price model would you change your answer in __3.b__?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "'''A'''\n",
    "al_prices_data=pd.read_csv('../data/al_price.csv')\n",
    "steel_prices_data=pd.read_csv('../data/steel_price.csv')\n",
    "al_y=al_prices_data['Year'].values\n",
    "al_y_n=(al_y-al_y.min())/(al_y.max()-al_y.min())\n",
    "st_y=steel_prices_data['Year'].values\n",
    "st_y_n=(st_y-st_y.min())/(st_y.max()-st_y.min())\n",
    "al_p=al_prices_data['dollars/MT'].values\n",
    "st_p=steel_prices_data['dollars/MT'].values\n",
    "train_per=0.7\n",
    "np.random.seed(103)\n",
    "al_rand=random.sample(range(0,len(al_y_n)),len(al_y_n)) #start,stop,size\n",
    "st_rand=random.sample(range(0,len(st_y_n)),len(st_y_n))\n",
    "#train & test\n",
    "train_al=al_rand[:int(len(al_y_n)*train_per)]\n",
    "train_st=st_rand[:int(len(st_y_n)*train_per)]\n",
    "test_al=al_rand[int(len(al_y_n)*train_per):]\n",
    "test_st=st_rand[int(len(st_y_n)*train_per):]\n",
    "\n",
    "al_train_year=al_y_n[np.sort(train_al)]\n",
    "al_test_year=al_y_n[np.sort(test_al)]\n",
    "st_train_year=st_y_n[np.sort(train_st)]\n",
    "st_test_year=st_y_n[np.sort(test_st)]\n",
    "al_train_price=al_p[np.sort(train_al)]\n",
    "al_test_price=al_p[np.sort(test_al)]\n",
    "st_train_price=st_p[np.sort(train_st)]\n",
    "st_test_price=st_p[np.sort(test_st)]\n",
    "\n",
    "Ztrain_al=np.block([[al_train_year**0]]).T\n",
    "Ztest_al=np.block([[al_test_year**0]]).T\n",
    "Ztrain_st=np.block([[st_train_year**0]]).T\n",
    "Ztest_st=np.block([[st_test_year**0]]).T\n",
    "#split 70-30 training testing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Price vs Year Data')"
      ]
     },
     "execution_count": 386,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhU5dn48e+dkJBEFiEsIpAEBVc2Bau4VNwBERWtS6PYasWm1mpftYJYlyrub1+lLj9QrGgCLlgV9wV3wUagIKBQIhJEkSWigmHJcv/+OGfIZDJnlmSWTHJ/rutcM/Oc7TlZzj3Pcp5HVBVjjDEmlLRkZ8AYY0zzZ8HCGGNMWBYsjDHGhGXBwhhjTFgWLIwxxoRlwcIYY0xYFixMiyEir4nIRcnOhzEtkQUL02yJyBoR2S4i20Rkg4j8U0TaeW2vqiNVdUYi8xiOiHQXkc0iMjwg/Z8iMisB5+8rIur+DLeJyHci8pKInBDFMX4nIu/FMZsmBViwMM3daaraDjgUOAy4IXADcTTLv2VV3QD8GXhERLIB3Bv1qcCfYnkuEWkTIh/t3J/jIcA7wBwRuSCW5zctW7P8BzMmkKp+A7wG9AcQkfdEZLKIfAxUAvu4ab/z7SMil4rIFyKyVUQ+F5FD3fS9ReQ5EdkkIl+JSNCbtogc4X4TT/dLO1NEPnPf/0JEFojIT27J5+8eeX8SWAn8zQ0YU4E/qeom9zi9ROR5v/xc7ne+YSLyiYj8ICLrRWSKiGS469q4pYY/iEgZsCKCn+N6Vf0/4FbgbhER91g3iMhq92e1XETGuOkDgAeAY9ySyWY3fYyILHa3Xysifw13bpPiVNUWW5rlAqwBTnTf9waWA7e6n98D1gIHA22ADDftd+76XwHf4JRGBOgL5ON8QVoI3AhkAvsAq4FTPPLwJXCS3+dngQnu+/nAhe77dsARIa6lF1ABvAi84JeeDiwGrnfz09e97hPc9YcBh7vXuA/wX+CP7ro2gAKvA52A7CDn7ev8mzdI38/dt5/7+Rygh/vz+TWwDejurvsd8F7A/sfjBO40YBCwGRid7L8ZW+K3WMnCNHcviMgPwEfA+8DtfuseV9XlqlqtqlUB+/0OuFtVP1VHmaqW49x8u6rq31R1l6quBh4BzvM4/yzgfAARaQ+MctMAqoC+ItJFVbep6ideF6Gq63AC1IlAkd+qI4AOqnq7m58yYLovP27+/+1e42pgGnBswOFvV9Utqrrd6/xBfOu+dnbP84w6pY5aVZ2JE7CGhried1R1mbv9EuCpIPkyLYhnHacxzcQZqvq2x7qvQ+zXG6dUECgf2NsNQD7pwIcex5kJzBORImAssMgNOgCXAH8DVojIV8AtqvpyiDwtB7ao6vqA/OQFyc97ACJyAPC/wBAgB+d/9t8Bxw31c/DS03393j3Pb3DaVvLd9HZAF6+dRWQYcAdOyS4TaEtdEDUtkAULk8pCDZn8NbCvR/pXqtovohOofi4i5cBInOqZmX7rVgHnu43rY4HZIpKrqj9HegFuflap6oEe66cCnwDnquo2EbkGGB2YzSjO53Mm8B1QJiL7AA8DJwD/VtUaEVmGU33ndfyngHuBEaq6Q0QewAkwpoWyaijTUj0KXCMiQ9zeUn1FJB8oBX4SketEJFtE0kWkv4gcFuJYM3F6Lv0Sp80CABG5QES6qmot4CsZ1ESZz/nALhG5WkSy3PwMEJEh7vr2wI/AzyJyIHBZlMevx+3K+yecXmXXqari3OQV2ORsIr8DDvDbbQPQy9ew7pev791AcQTe1XimhbBgYVokVX0WmIxzo98KvAB0VtUa4DRgMPAVTsPso0DHEIebBQwH3lHVzX7pI4DlIrINuB84T1V3RJnPapx2kF/gtBNsxilNdHA3uRq4yL2GqcDT0Rzfx/ecBfAZcAowVlWfcPPwGTAFJ5CuxwkU/lVdbwGrgA0i8p2bVgTcISJbcRrnn2lMvkzqEOeLhTHGGOPNShbGGGPCsmBhjDEmLAsWxhhjwrJgYYwxJqwW+5xFly5dtKCgINnZMMaYlLJw4cLNqto1ML3FBouCggIWLFiQ7GwYY0xKcR9CbcCqoYwxxoRlwcIYY0xYFiyMMcaE1WLbLIwxLUNVVRXr1q1jx46oRlIxYWRlZdGrVy8yMjLCb4wFC2NMM7du3Trat29PQUEB7sR+polUlYqKCtatW0efPn0i2seqoYxJMSUlUFAAaWnOa0lJsnMUXzt27CA3N9cCRQyJCLm5uVGV1qxkYUwKKSmB8eOhstL5XF7ufAYoLExevuLNAkXsRfsztZKFMSlk0qS6QOFTWemkGxNPcQsWIrK/iCz2W34SkatEpLOIvCUiq9zXTn77TBSRMhFZKSKn+KUPEZGl7ropYl8zTCu1dm106SZ2nn/+eUSEFStWALBmzRr69+8fk2MvWLCAP/3pTzE5VrzELVio6kpVHayqg3HmD64EngcmAHPdaS3nup8RkYNwZts6GGdSmYdEJN093MPAeKCfu4yIV76Nac7y8qJLb43i1aYza9Ysjj76aJ566qnYHNDP0KFDmTJlSsyPG0uJqoY6AfjSnej+dGCGmz4DOMN9fzrwlKruVNWvgDLgFyLSA+igqvPdKSCf8NvHmFZl8mTIyamflpPjpJu6Np3yclCta9NpasDYtm0bH3/8MdOnTw8aLB5//HH++Mc/7v48evRo3nvvPQDatWvHddddx5AhQzjxxBMpLS1l+PDh7LPPPsyZMweA9957j9GjnanVb775Zi6++OLd2/iCSGBJ5t577+Xmm28GYPjw4fz5z3/ml7/8JQceeCCffvopY8eOpV+/ftxwww1Nu3hXooLFeThTUwJ0V9X1AO5rNze9J87k9T7r3LSe7vvA9AZEZLyILBCRBZs2bYph9o1pHgoLYdo0yM8HEed12rSW3bgdjXi16bzwwguMGDGC/fbbj86dO7No0aKI9/35558ZPnw4CxcupH379txwww289dZbPP/889x4441B91mxYgVvvPEGpaWl3HLLLVRVVYU9T2ZmJh988AG///3vOf3003nwwQdZtmwZjz/+OBUVFRHn10vcg4WIZAJj8Jvo3mvTIGkaIr1houo0VR2qqkO7dm0waKIxLUJhIaxZA7W1zqsFijrxatOZNWsW5513HgDnnXces2bNCrNHnczMTEaMcGrOBwwYwLHHHktGRgYDBgxgzZo1Qfc59dRTadu2LV26dKFbt25s2LAh7HnGjBmz+xwHH3wwPXr0oG3btuyzzz58/fXXYfYOLxFdZ0cCi1TVd7UbRKSHqq53q5g2uunrgN5++/UCvnXTewVJN8aYevLynKqnYOmNVVFRwTvvvMOyZcsQEWpqahAR/vCHP+zepk2bNtTW1u7+7P/8QkZGxu5uqmlpabRt23b3++rq6qDn9G0DkJ6eTnV1dchz+O/jf45w54lGIqqhzqeuCgpgDnCR+/4i4EW/9PNEpK2I9MFpyC51q6q2isgRbi+ocX77GGPMbvFo05k9ezbjxo2jvLycNWvW8PXXX9OnTx/WraurHS8oKGDx4sXU1tby9ddfU1pa2vgTeujevTsbN26koqKCnTt38vLLL8f8HKHEtWQhIjnAScBlfsl3As+IyCXAWuBXAKq6XESeAT4HqoHLVbXG3acIeBzIBl5zF2OMqcdXJTdpklP1lJfnBIqmVNXNmjWLCRMm1Es766yzuP3223d/Puqoo+jTpw8DBgygf//+HHrooY0/oYeMjAxuvPFGDj/8cPr06cMBBxwQ83OEIk4Ho5Zn6NChapMfGZP6vvjiCw488MBkZ6NFCvazFZGFqjo0cFt7gtsYY0xYFiyMMcaEZcHCGGNMWBYsjDHGhGXBwhhjTFgWLIwxxoRlwcIYY0JIT09n8ODBDBo0iEMPPZR58+Y16jj33XcflYEDV7mGDx/O/vvvz+DBgxk8eDCzZ88G4MgjjwScQQRnzpzZuAuIEQsWxhgTQnZ2NosXL2bJkiXccccdTJw4sVHHCRUsAEpKSli8eDGLFy/m7LPPBtgdmCxYGGNMCvnpp5/o1Gn3fG3cc889HHbYYQwcOJCbbroJcEaZPfXUUxk0aBD9+/fn6aefZsqUKXz77bccd9xxHHfccRGfr127dgBMmDCBDz/8kMGDB/N///d/sb2oCNkc3MaYlHHVVVexePHimB5z8ODB3HfffZ7rt2/fzuDBg9mxYwfr16/nnXfeAeDNN99k1apVlJaWoqqMGTOGDz74gE2bNrH33nvzyiuvAPDjjz/SsWNH/v73v/Puu+/SpUuXoOcpLCwkOzsbgLlz55Kbm7t73Z133sm9996b8PGg/FmwMMaYEHzVUADz589n3LhxLFu2jDfffJM333yTQw45BHAmSFq1ahXHHHMM11xzDddddx2jR4/mmGOOieg8JSUlDB3aYJSNZsOChTEmZYQqASTCsGHD2Lx5M5s2bUJVmThxIpdddlmD7RYuXMirr77KxIkTOfnkkz0nOUol1mZhjDERWrFiBTU1NeTm5nLKKafw2GOPsW3bNgC++eYbNm7cyLfffktOTg4XXHAB11xzze5Z9dq3b8/WrVsbdd6m7BsrVrIwxpgQfG0WAKrKjBkzSE9P5+STT+aLL75g2LBhgNMYXVxcTFlZGddeey1paWlkZGTw8MMPAzB+/HhGjhxJjx49ePfdd6PKw8CBA2nTpg2DBg3iN7/5DX/+859je5ERsCHKjTHNmg1RHj82RLkxxpiYsmBhjDEmLAsWxphmr6VWlydTtD9TCxbGpLiSEigogLQ057WkJNk5iq2srCwqKiosYMSQqlJRUUFWVlbE+1hvKGNSWEkJjB8PviGHysudzwCFhcnLVyz16tWLdevWsWnTpmRnpUXJysqiV69eEW9vvaGMSWEFBU6ACJSfD2vWJDo3piVISm8oEdlTRGaLyAoR+UJEholIZxF5S0RWua+d/LafKCJlIrJSRE7xSx8iIkvddVNEROKZb2NSxdq10aUb01jxbrO4H3hdVQ8ABgFfABOAuaraD5jrfkZEDgLOAw4GRgAPiUi6e5yHgfFAP3cZEed8G9Ns+bdRpHn8B+flJTRLphWIW7AQkQ7AL4HpAKq6S1V/AE4HZribzQDOcN+fDjylqjtV9SugDPiFiPQAOqjqfHXqzJ7w28eYVsXXRlFeDqpQU9Nwm5wcmDw58XkzLVs8Sxb7AJuAf4rIf0TkURHZA+iuqusB3Ndu7vY9ga/99l/npvV03wemNyAi40VkgYgssMYw0xJNmlTXmO0vPR1EnLaKadNaTuO2aT7i2RuqDXAocIWq/ltE7setcvIQrB1CQ6Q3TFSdBkwDp4E7uuwa0/x5tUXU1jqLMfESz5LFOmCdqv7b/TwbJ3hscKuWcF83+m3f22//XsC3bnqvIOnGtDpebRHWRmHiLW7BQlW/A74Wkf3dpBOAz4E5wEVu2kXAi+77OcB5ItJWRPrgNGSXulVVW0XkCLcX1Di/fYxpVSZPdtok/FkbhUmEeD+UdwVQIiKZwGrgtzgB6hkRuQRYC/wKQFWXi8gzOAGlGrhcVX3Nd0XA40A28Jq7GNOqlJTUtVmkpzuN2/n5TqCwNgoTb/ZQnjEpIPBJbXBKFNaYbWLNhig3JoUF6wVVWemkG5MIFiyMSQH2pLZJNgsWxqQAr95Oqi1zpFnT/FiwMCYFBOsF5eMbadYChoknCxbGpIDCQqcxOz8/+PrKSrjyypY9r4VJLgsWxsRJrCclKix0hh33GnO5oqJuzKjycrjgAujSxYKGiQ0LFsbEQeCAf7GsKormae2KCquiMrFhwcKYOIhnV9dQ7RfBWBdbEwsWLIyJA68ureXlTa+W8m+/8I00m5vbuPwYEykLFsbEQaiqolhUS/naL2prndf77w9d2rCBBk1TWbAwJg4iqSqKZfWQr7QRrIRhAw2aWLBgYUwcBFYVeYll9VBhIWzeDMXF9auobPwoEws2kKAxCVBQ4FQ9BcrPd6qRjGkubCBBY5LI5qEwqc6ChTEJkp1d9z4316qHTGqxYGFMnPke0KuoqEvbvr3+ehumwzR31mZhTJx5tVeAU8LYuhV27apLs0mNTDJZm4UxSVBS4h0owClt+AcKsCeuTfNkwcKYOPFVPzVGsC61qsrOnTvZsWNH0zIWA1Z11vq0SXYGjGmpgo0PFalevX7ktdfm8eGHH/Lxxx9TVlbG+vXr8VUbt23bll69ejFw4EBOOukkTj75ZPbdd98Y5t5b4HzgvqfRwarOWjJrszAmTtLSnKE9IrcdeIG0tBmovoVqLW3atGHIkCEceOCB9OzZkxy3/+2WLVsoLy+ntLSUcree64ADDuCyyy7j4osvpkOHDjG/Hh97ZqRl82qzQFVb5DJkyBA1JlaKi1Xz81VFnNfi4vD75OerOuGi/iJS/3ObNps0O/sGhT0V0NzcPJ0wYYLOnTtXt23bFvIctbW1unLlSp0yZYoOGzZMAe3cubPefffdunPnzlhcegOB+fe/LpP6gAUa5J4a1xs2sAZYCiz2ZQDoDLwFrHJfO/ltPxEoA1YCp/ilD3GPUwZMwS0RhVosWJhYKS5Wzcmpf2PMyQkfMLz2KyryBZIK7dDham3bNkdFRMeOHatz587VmpqaRue1tLRUR44cqYAedthhunr16kYfy4tXEMzPj/mpTBIkM1h0CUi7G5jgvp8A3OW+PwhYArQF+gBfAunuulJgGCDAa8DIcOe2YGFipSk3x2AlkqqqKn3wwQe1c+fOmpaWphdccIHeddfyqEsuoTz33HPasWNH7dixoz733HNNO1iAxgZPkxqaU7BYCfRw3/cAVmpdqWKi33ZvuAGiB7DCL/18YGq4c1uwMLESy2qX+fPn64ABAxTQ4447Tj/77LO43XxXr16thx12mAJ6xRVX6I4dO6LaP1TVW2Oq5UxqSFaw+ApYBCwExrtpPwRss8V9fQC4wC99OnA2MBR42y/9GOBlj/ONBxYAC/Ly8uL1szStTCyqXbZv367XXnutpqWlae/evXX27NlaW1sbs+N72blzp1511VUK6JAhQ7SsrCyi/az00HolK1js7b52c6uYfhkiWDwYJFicBRwWJFi8FO7cVrIwsdLUG+eiRYv0gAMOUEDHjx+vP/74Y731iWgw/vOfX9C0tD0VOuteey0Mm3drl2i9vIJFXB/KU9Vv3deNwPPAL4ANItIDwH3d6G6+Dujtt3sv4Fs3vVeQdGPiyvfg2YUXOoMA5uZGP0fE9OnTGTZsGFu3buWNN95g6tSpDbq1es1iF6vZ7UpKYOrU06mtXQC047vvTuCSSz4N+SCd1zwbNj1rKxYsgsRiAfYA2vu9nweMAO6hfgP33e77g6nfwL2augbuT4EjqGvgHhXu/FayME3R1NJEZWWlXnzxxQroiSeeqBs3bozbucKpX0pYo9BHoYP26LEown2sZNGakOhqKGAf9+a/BFgOTHLTc4G5OF1n5wKd/faZhNMLaiV+PZ5w2i2WuesewLrOmjjxNdwGu1FGerMsKyvTwYMHK6A33HCDVldXR3Te3Ny68+Tmxi5YNKzmWqvQW6G7Z9daa7NovRIeLJK9WLAw4QT26CkqaniDjLYd4YUXXtCOHTtqp06d9OWXX44qL/G6OQcPfss1La2T9uvXz7PUYz2eWicLFsb4CXZz9mpojqRkUVtbq3/961939zr66quvospPPKt9vALRjTd+pFlZWXrEEUdE3a3WtFwWLIzxE6qqyWsJ9k2/uFg1L2+Hwq8V0GOPvVi3b98eUR78v7k3tiQTKa9SwuzZsxXQSy65ZHdXXtO6eQULG6LctErR9uoJ1gOqpAQuvbSCtWtPAmYCd1Ba+ijPPZcV9ni+kVvLy52w4MW/R1RThgUvLHQG+XvySefzhRc6x9ix4yyuv/56pk+fztSpUyM/oGl9gkWQlrBYycKEEukgf6HaDfbeu0xhP4VMhVlRVR1FUrLxP3cs2jS8jvHEE9U6cuRIzcjI0I8++sjaKlo5rBrKmDpFRcEDg2+Qv3A3yvnz5yt0Ueis8EHUVUfhqp4Czx2LNo1Qx9iyZYv27dtXO3bcS7OzN1gvqFbMgoUxLq/G7aKiyPZ//fXXNTs7W9u02VdhZaNu4F437tzc4MEqFk95hzvGkiVLFNoqjFaojShfpuXxChbWZmFanWAz2KnCq6+G3/fZZ5/l1FNPo7p6P6qrP0Zkv3rrc3Jg8uTwx5k82dnWX0YGbN1a147hm4GupCQ2T3mHOkZJCYwZMxC4C3gZeKjeNhUVwfNlWpFgEaQlLFayMF5CVQGF+sb8yCOPqEiapqUdpbClQTtHtN+4A9sG/B/KCyypxKvNwpf3up9JrcJIt4SxNGSbij3N3TJh1VDGOEI1LnvdgO+55x4FNCtrhMLPcblxhqsmikXDs/8T6t5B8zuFbgr9FSpDtq2YlscrWFg1lGl1glUB+VRWOtVUPqrK9ddfz7XXXss555zDjh0vAg13Li9verWMVzWRqtPNFZzur7W1zmskAxmGouq1pjswA1hG+/Z/ITc3+FaxGujQpIhgEaQlLFayMKEUF3uXLnzfmGtqarSoqEgBvfTSS7W6urpRpZJo8hRquJF4Hz9wad/+SgX06qtfsnGiWhGsGsqY+kJ1Ja2urtbf/va3Cui11167++nmcDfcplZHxWIgw2iv1yswPfbYdh04cKB269ZNH354k/WGaiUaHSxwJi7aw32fjTMy7J24U6M218WChQkn1ENqF154oQJ60003NRgGI5JSSVPFY0KkcGNfBWuoX7x4sWZkZOi5554bk+syzZ9XsIikzeIpnGHFAW4B+gJbcMY3MCZlFRY6Q3jk59dNavTww9W89tqFPPnkk9x6663cfPPNiEiD/fLzgx8zVvX48ZgQKdS++fnOUCCq9dtDBg0axI033sjTTz/Ns88+2/iTm9QXLIL4FuAioBwY5/f+cvf9ajd9YKhjJGuxkoWJ1q5du/RXv/qVAnrHHXeE3Dbe8z3E4/jRHNO/51Ve3i7t02eIdunSRTds2ND4DJiUQGOqoYB8YAXOPNgnAqVAnpv+qfu+Y6hjJGuxYGGisXPnTh07dqwCeu+990a0T7zHUIrH8SM5ZvB2mWUKmZqTM1affNJGp23JGhUsnP0oAjYA3wOj3bQ84LVw+yZzsWBhIrVjxw4dM2aMAnrfffclOztJ590QfqcCmpk50xq4WzCvYCHOutBEpB1Qq6qV7uc9gAxV/aHxFWDxNXToUF2wYEGys2GauR07dnD22Wfzyiuv8MADD3D55ZcnO0tJl5bm9QxGNXA08F969lzOunU9EpsxkxAislBVhwamR/RQnqpu8wUK9/PPzTlQGBOJ7du3c8YZZ/DKK68wderUJgWKpsw10dx4N4S3AR4HtvPNN5cRyRdN03LYE9ymVaqsrGTMmDG8+eabTJ8+nfHjxzf6WIETGaXaQHuBgW7UKO8n3OEAYDLwEsXFxYnKomkOgtVNtYTF2iyMl23btulxxx2nIqIzZsxo8vHiOX92vHn1kPLN6+H//IVvyc6u1v32O0pzcvbUnj3X2YN6LQxNGRtKRPYQkTT3/X4iMkZEMiLcN11E/iMiL7ufO4vIWyKyyn3t5LftRBEpE5GVInKKX/oQEVnqrpsigR3fjYnQ1q1bGTlyJO+//z7FxcWMGzeuycf0mqI12qlbkyHYcO2Vlc5w7WvWOOHhySfrP4vym9+ks3HjP6ms3Mk331yKqqZcaco0QrAIErgAC3FGT+sJfA08D5REuO//4DzA97L7+W5ggvt+AnCX+/4gYAnQFugDfAmku+tKgWGAAK8BI8Od10oWJtCPP/6oRx55pKanp+vTTz8ds+Omcski2ifF65dE7ldA4dGUumYTGk0cdVbUaeAeC/xDVc90b+6hdxLpBZwKPOqXfDrOkJa4r2f4pT+lqjtV9SugDPiFiPQAOqjqfPdCnvDbx5iI/PDDD5x88smUlpby9NNPc84558Ts2KNGRZcejXg2nJeUOMcNxquRu35J5I/AcOAqnO92qVGaMo0TcbAQkWFAIfCKm9Ymgv3uA/4C1PqldVfV9QDuazc33Vdq8VnnpvV03wemGxORLVu2cNJJJ7Fo0SJmz57NWWedFdPje82wF8nMe6HEs+Hcd+yamobrQs32Vz8YpOF8d2uDc2uosmHLW7BIg8VVwETgeVVdLiL7AO+G2kFERgMbVXVhhOcI1g6hIdKDnXO8iCwQkQWbNm2K8LSmJauoqOCEE07gs88+41//+henn356zM8RrzYLr/YE//k2YnlsgPR0Z7wsr7kyGgaD3sBU4N+0aXNbRFPKmtQU6XMW76vqGOAB9/NqVf1TmN2OAsaIyBqcwQiPF5FiYINbtYT7utHdfh3OX55PL+BbN71XkPRg+ZymqkNVdWjXrl0juTTTgm3YsIHjjz+ezz//nBdffJHRo0fH5TzxGPQP4ttw7nWM2trQkyoFnzjqHNq2HUdNzW0UFHzcop45MX6CNWQELjiNy58Da93Pg4CHItnX3X44dQ3c91C/gftu9/3B1G/gXk1dA/enwBHUNXCPCndOa+Bu3dauXav77bef5uTk6FtvvRXXc8VrUMF4Npw35djBxpf68ccftU+fPtq1a4FmZ/8YtwEWTfzRlMmPgH/jfOv/j1/askj21YbBIheYC6xyXzv7bTcJp6VsJX49noChwDJ33QM4De4WLExQZWVlmp+frx06dNCPPvooIeeM16B/gUFIxHkGIh7HbupNfd68eQrpChemZM8w42hysHBf/YPFkkj2TdZiwaJ1Wr58ufbo0UNzc3N1wYIFyc5OkxUVNezeGqtv6vEIcHCz2512VkRdcU3z4xUsIm3g/lpEjgRURDJF5Brgiwj3NSYhFi1axLHHHouq8v777zNkyJBkZ6nJXn214aB+sWrkLix0Hryrra0/4VFT5OVNwqm1/j3O9De+9KYf2yRXpMHi9ziTHvm6sQ52PxsTU16No+EaTefNm8fxxx9PTk4OH374IQcffHBiMx4nqfZ0+O23tyErqxinw+JYoJKMDNi2zRq8U16w4kZLWKwaKnX4qkOCjUPkG6coVP3622+/rTk5OdqvXyVM1bwAACAASURBVD9du3ZtUq8l1lLx6fDiYtWuXV9SEM3M/LVmZNRag3cKoYltFjOAPf0+dwIei2TfZC0WLFJD8FnZ6i/p6d43zDlz5mjbtm21f//+un79+mRfTszFe/rWeLrtttvc9ot7UirYtXZewSLSaqiB6jd/hapuAQ6JaRHHtEpeD4f5C/aUMUB5+dOMHTuWAQMGcPnl73HEEXu1uKqOwkLnITn/gfxCPTTXnFx//fXA2cB1wBv11sWiGs2e50iwYBEkcMF5/qGT3+fOwNJI9k3WYiWL1OA1kF34ksWDCml6zDHH6LRpP6bst++WrnfvrQoDFPZUWNWgdNHY31Eql7iaO5pYsvhfYJ6I3CoitwLzcEaPNaZJwvWSyclxxjCqe2pYgeuByznkkFG8/vrrTJ7cIW7DYpimueOOdmRlvYjTl+Z0YMvudU0Z6yqeQ6GY4CId7uMJ4CxgA87wHGNV9cl4Zsy0DsGGj/DNVuKrcnnoIec1L68K+A1wB3ApFRXP8/zzOSnXY6g1KSyERx/tQ/fus3EGkj4V+Hn3+sbe4O13nnghg4WIdHBfOwPf4cxLUQJ856YZ0yTB6uSffNKpWPDv+z9mzFY6dhyNM8rp34CprF3bhvHjobPHX2JamtVnNweFhfDdd8cBs3AGgxgL7Ny9vjE3+HiNx2W8hStZzHRfFwIL/BbfZ2OaLNzDYd999x3Dhw9n6dK5OFOj/BXfYMS+qohgc0bX1MR+aG/TePn5Y4FHgDeBXwO7gMbd4IOVSEMNrW6aLmSwUNXR7hSmx6rqPn5LH1XdJ0F5NCkoVE+VSB6869LFWUT+S69ew1i+fAUwB7ikwbm+/75+6SQ9vWF+rD47+Zwb/MU409z8CziD7OzKejf4cD2cfOsvvBCysyE310lPT6/7HduXgjgJ1uoduAALI9muOS3WGyp5QvVU8VoX7ME7mK+Qq9BV27Yt1dzc4L2lAvvsh+phZZKr7gHMRxRE99//GP3hhx92rwvVwynY+owM1cxM6xUVSzTxobwHgcMi2ba5LBYskifUU8de6xp2j31RIVth391dLnNzgz/Al5tb/+bgdQ4Ru4k0J88884xmZGToIYccouvXrw/7tLrX+lR7wr258woW4qwLTUQ+B/YH1uB0ZRCnUKID41PeabqhQ4fqggXWrJIMaWnOv2wgXy+n8H9yU4E/AIfizOLbbff+Tz4JV14JFRX198jJcaqiIPh6n/x8p13ENA+vv/46Z511Fp07d2bduheAhoM/ijjtWV5/V8H49jHRE5GFqjo0MD3S5yxGAvsAxwOnAaPdV2MaCNVTxWud086gwI0441aOwJm5t9vubfLynMbvdu0a7l9Z6QSJ8eO9AwVY18rmZsSIEXz88cekpaXhTK75EATMmuz7m4mmIdx6RcVeuK6zWSJyFXAtzn/vN6pa7lsSkkOTckL1VPFad8klVaSnXwLcClwMvAi0a7A/eN/wKyrCDx1iN5HmpaQEzjhjMGvXLiAj43icwazPAr4H6v/eg/3tZGRAZmb9NOsVFSfB6qZ8C/A0UAxcBrwA3B9q++a0WJtFcoWaWCdw3aOPbtWRI0cqoB073qhQq7m5TltEsP2jqbu2hs/mq2GDdY2mpf2vQobCXpqb+5g+8URNvX2Kiurat9LTnc/xmMSpNaMxDdz4jf8EtAEWhdq+OS0WLFLDhg0bdOjQoZqWlqbTpk2rt87rJuDVa8art1RTxyEy8eEV9Pfaa4EefvjhCughhxyi7733nqraeFCJ0thgsSjU5+a8WLBIvnDf+FatWqX77ruvZmdn6//8z5x624abwyLYse1mklq8ujiLqNbW1urMmTO1d+/eCuiZZ56pe++9yvOLQCArbTReY4NFDfCTu2wFqv3e/xRq32QvFiwSJ5obd1GR7xvlvzUtrYu2a5erN9/8SYNtvW4k4bpE2k0idUQysVNlZaXedtttuscee7jVU1crbGnwt+LPvjQ0TaOCRSovFiwSI9oqIScIvKaQo9BHs7JWhqw+Cra/aRmiual/++23usceFyuIQhd1hqivCvoFIhVnF2xOvIJFpF1njQnKa6hor+6rqrNwel33A+axY8d+Ibu6BrLeTC1HsEEkL7rI+ZsKHO6jR48eTJ06nayshUB/nF5TB5KZOYO//a263nFtRNo4CRZBYrEAWUApzsRJy4Fb3PTOwFvAKvfVf1KliTjjGK8ETvFLHwIsdddNAedhwlCLlSwSI5LJi+qWB9xvhr9U+CGiUoRVJbQekZQ0iotVO3euVXhBYbAC2q3bPjp9+nTdtWuXFheHnobXhEeiq6FwnvJu577PwBmb+AicSZMmuOkTgLvc9we5gaUt0Af4Ekh315UCw9xjvgaMDHd+CxaJ4VXkrz80R63CzQoonKZQGWLbupuEr33D2h9apsD2pUjG/qofUJygIXKoAtq1ax/NzHxEYWfQ4xQVJelCU0zCg0W9k0AOsAg43C019HDTewAr3fcTgYl++7zhBogewAq/9POBqeHOacEiMcINGpiXV6PwRwV0//0v0uzsKs9tLTC0HsH+biJppwr+5aRWu3Z9STMzh7pfSPIV/l+DoJHKJYtE/n8kJVgA6cBiYJtfCeKHgG22uK8PABf4pU/Hme19KPC2X/oxwMse5xuPO+dGXl5evH6WJoDXH3J1dbVedNFFCujVV1+tNTU1FhSMqjZ+UMBQ3W2dksYrCr9wg0ZvhYcUdjQIOqkk0b27kl2y2BNnoJ/+IYLFg0GCxVnAYUGCxUvhzmkli+SqqqrSwsJCBfSWW27R2traZGfJNCORtnUF3hQjG9G4VuF1hWFu0Oip8A/Ny9sel2uJ1xeguuHcgy/p6fEJGEkNFs75uQm4xqqhWr6qqio999xzFdDJkycnOzumGQp1E/QFkmA33ujmSqlVeFPhKAW0U6e99cEHH9SdO3fG7Dri9a0/0mq6eJQwktHA3RXY032fDXyIM1rtPQEN3He77w8OaOBe7dfA/SlO47ivgXtUuPNbsEiOXbt26dlnn62A3nXXXcnOjmmmwt0MQ7UvRDLumO+bN6jm5dXqxIlz9eijj1ZACwoKdMaMGVpdXd3k64jXMx3JnLsjGcFiIPAf4DNgGXCjm54LzMXpOjsX6Oy3zyScXlAr8evx5LZbLHPXPYB1nW12nMbsnQpnKqCFhX9PdpZMM1dcHLp0EWu1tbX62muvaX7+IW711IHapcvTOmNGVaOPGaoNJZhIq6yi6ZIe659V0quhEr1YsEic4mJ1ezmd7f4T3m/PRJiIJPppa+dvtUbhWYUDFFCRvfXMM2/SdevWRX28aPIfTZVVqypZJHuxYJEYxcWqaWk1Che4geLvcf+HNy1HvHv6BKuWqluq1Zm+d4SCaHp6up555pn65ptvak1NTdhje+VfJPgzHU0LLLUK/3S/kA1TuKvltFkke7FgEX/Ot7Rahd+7gWJy3KsSTMsTz95EkT7LAV/qX/7yF+3SpYsC2rdvX7333nt18+bNYc9TVNSw2ihYwAhVZeU1GKcTYH7yK7UXKPRTyNCcnC9bZm+oRC8WLOLP+UP+m/tHfF3CqhKMiURjqnJ27NihxcXFetRRTg+qtm3b6rhx43T+/Pme3b+9zuMLAuG28xrBoLhYtaysTOEghTSFe9wSxjcK2SpyQVx+bhYsTMzBI26gGOf+ETf8YzcmWRr7LIfPkiVLtKioSNu1a6eADh48WKdOnapbt26N+DzeQ5XUndtrmJMePZZo9+7dFTorvB2w/joF0c8++yzmPzcLFiYq4aoGXnrpJYV0dep7d9X7Q47Xw0LGRCOSkkUk1V4//fSTPvzwwzpw4EAFtH379nr55Zfr0qVLw54n2Fwbgf9XwYPNPIU9NT29p4p8HmT99wp76mmnnRbzn5sFCxOWf2NgqBFf58+fr9nZ2dqnz1DNzt5qJQrTLIVqs2jM32ltba3OmzdPL7zwQm3btq0CevTRR+sf/lCiviFFGlMV2zDYPKaQpdBX4StNSwt+DYcffrsC+tFHHzXmx+PJgoUJKZLGwPx81RUrVmhubq727dtXN2zYYGM9mWYtWG+oWPydbtq0Se+55x7dd999FdCsrK5u1VDd1K+hApL//01urmpmpipsV7jUrdo9XmHD7mPtsUdd/tPTncbzbdu26V577aXHHHNMTIfTsWBhQoqsMfAbzc/P127dumlZWVmys2xMo8XqS05NTY2+8cYbesYZZ6hImvoe9uvQ4Rq9+uqXdcuWLUHPXf+LWZWmp8/R9PSD3P2vV6dbb902Xj0LH3roIQX01VdfbdwFBOEVLMRZ1/IMHTpUFyxYkOxspIy0NOfP0tuPZGQcS2ZmGe+//z5DhgxJVNaMiamSEhg/vv4Mjzk5zqx9hYWNP+66deuYPXs2r7zyCu+//z5VVVWICIMGDeKoo46ioKCA7t27c801e7FxYybwLc78b3OACtq06UvnzvezceOoBsfOz4c1axqec9euXRx44IG0b9+eRYsWkZbW9MlPRWShqg5tsCJYBGkJi5UsohO6ZLFL4USFNtqt2xtW1WRSWjyfGq+r9vpZu3d/R8eOvVkPOuh4FdnDLTUELh0VChVmK+yM+iHF4mLV3NwZCmj37u/F5H8Tq4YyoXg9iep0ifU9dPdPa8Q2KS/a8ZwiFex/KCPD1x5Rq/Cjwn8VPlB4S2Gp+k/QlJ5e14aRmxu+iqzufBXu/+fkmPxvWrAwYQWrx+3UaYr7h/iXmH8LMyYZmsNIseGeAYnkpl//fAcqnBqT6/AKFk2v4DItRmGhUy9aW+u85ua+zpYtVwGnA3fU23bt2iRk0JgYmDzZaaPwl5PjpDdFNP8Tqk47hAikpzdcX1kJkyZFc74jgfmAxu1/04KFCerzzz/n3HPPJSNjAFBM4J9KXl5SsmVMkxUWOo3Zvpt1fn7TG7chuv8JX4N1ba2zBBPupl//fEcC3wP/jdv/pgWLVq6kBAoKnN5QBQXO582bNzN69Giys7O555455OS0q7dPLL6FGZNMgaXopgYKCF5iCdY5KfD/p3Pn4McLd9Ovf74j3NdStm1z/o9jLljdVEtYrM0ivGANctnZO3T//Y/Rtm3b6ieffLJ7O3vwzpjwgj9sV7+twn802uJipxE8sM0iMzOy/zOnN5QqVClk7m5bbEpDN/achQlUUADl5f4pClwMPM7MmTM5//zzk5IvY1qChv9fDv9nJry2yc2FzZujPc9goCfwSoPzRMPrOQurhmrFGtaJ3gM8Dty4O1AEq6YyxoTn1ebgn+61zfffN+Y8/XFmnw597MayYNGK1a8TfRGYAPyKvLybgLonXcvLncJxebnz2QKGMeF5tTn4p0eyTeTn6Q+sBX4CvNtCGsuCRStW10C2GCgEhpCd/Ti33+78WUyaVH9IBIisS58xJrIuurHoxjt5MmRkgBMswFe62Lo1xl/sgjVktISltTVwN7YR+r771mh6eg+Fntqz5zdaVFR3nFAPFBljwovk/zIWHUicRu6v3AdopzbpQUPsCe6Wq7GT3m/evFkPOOAA7dixoy5dujTiOYt9f9DWQ8qY5qFuaJ52Clc06YudV7CIWzWUiPQWkXdF5AsRWS4iV7rpnUXkLRFZ5b528ttnooiUichKETnFL32IiCx1100REYlXvpurUA3Njakuqqys5LTTTmP16tW8+OKL9O/fP+hxAuXkwKhR1pZhTHPitFsIgY3cMX1AL1gEicUC9AAOdd+3B/4LHATcDUxw0ycAd7nvDwKWAG2BPsCXQLq7rhQY5v40XgNGhjt/SypZhCs5RDswWlVVlY4ZM0ZFRJ999tnd6eGqnvxLFPEatdMYE726e8TvFLo26VkLkl0NhdPd5iRgJdBD6wLKSvf9RGCi3/ZvuAGiB7DCL/18YGq48zXXYNGY6ptwN+dobt61tbU6fvx4BfQf//hHVOfxideoncaYxnMG/rxPIUd79fo+5g/lJSpQFOD06eoA/BCwbov7+gBwgV/6dOBsYCjwtl/6McDLHucZDywAFuTl5TXuJxVHjWlbKC4O/W0/2uPecsstCuiECRManT8rWRjTPO3YsUNramqadIykBQugHbAQGOt+9goWDwYJFmcBhwUJFi+FO29zLFlEe5MN1+Dsv18kJZZHHnlEAR03bpznnL2R9t5oTIO6Mab5S0qwADLc6qT/8UtrtdVQ0VbfhBofP9qb85w5czQtLU1HjBihu3btavK1WG8oY1omr2ARz95Q4pYOvlDVv/utmgNc5L6/CKctw5d+noi0FZE+QD+gVFXXA1tF5Aj3mOP89kkp0T6tGepx/WiGVP7ggw8499xzOfTQQ3n22WfJcJ7gaZJ4jNppjGm+4vkE91HAhcDxIrLYXUYBdwInicgqnAbvOwFUdTnwDPA58DpwuarWuMcqAh4FynB6Sb0Wx3zHRUkJbNvWMD3U05peQSQ9HT7+OHhX2sAuttde+yqnnHIK+fn5vPLKK7Rr1y74QY0xJgQbdTYBfGMsBT7DkJsL99/v/a3ca79gcnLgootgxgz/7Z8BCsnPH8Cnn75B165dm3AVxpjWwEadTSKvh93atQtdfeOb0SvYtIuBKiudbevOMx2neecIvv76Xbp372qjxhpjGs2CRQJEMlSxl8JC72kXA9XUAFQBVwK/w6nle4Pa2o6o2pPWxpjGs2CRAE0dhjjS7dLS1gPHA1NwAsZLQP0hLW3UWGNMY1iwSICmDkMcbP9AmZkvkZU1BFgEzATuw+m53FCsJ0UxxrR8FiwSwNf2kJ8PIs5rNF1fg+1fVOS8wjfk5JzFrl1jKCjoxB13fEJ+/vmIeLd1xHRwMWNMq2DBohnz7wY7aZJTwqitdV5ffvkHysv/hsiBVFW9yu23385//vMfJkwYsPv5hxkzmj6xijHGADafRSIEGx7D9zR3NENqZGfX6tlnL9D09EsUshVQOEOzsso8n6C2J62NMdHA4wlue84iAQoKnJ5IXnJyGlZLOfvUAqtxRm5/H3gZ+Aqn0boQuBwYBDhVUmvWOPuWlDglkbVrnSqnyZPtCWtjTGS8nrNok4zMtDbhGpQrK2HixK0UFHzGkiVLWLJkCeXlnwFLgZ/drbKBE3CmADkX6Bj0HIEP8vm6y4IFDGNM41nJIgGClyw24ZQU3sDpwVQGOL+LTp06UVk5kJ07BwEDcUoP/YEs0tN9z1PU5ytZeJVi/EsexhjjxZ7gTqK6rq+1wKs4D8vtBVwMfAQMoGPHW3jppZdYu3YtFRUVTJ/+Hjk59wOX4EzpkUVOjlNKCNVo3ZQHAI0xxosFi0YKNSd2oMJCuOaaeWRmHgacCqwAJuFM8/E1OTnP8eCDf2X06NH07t0bEfHsbvvQQ6G74Tb1AUBjjAkqWKt3S1ji2Rsqmsl/amtr9bbbblMR0Z49e+rjjz+ujz++q0EPpVj1WvLqeVVU1PjrNca0HiR7Du5EL/EMFpHOePfTTz/p2LFjFdA99ihU2Bo0EMR65rmiooYTLdlMdsaYSHgFC2vgboS0NOcWHEikbtC/TZs2cdxxx/HFFytIT7+XqqorAQHqd5UtKXGGFg/VaB0ta+Q2xjSWdZ2Noby84DdjX7tAZWUlp512Gl9++SVdurzBxo0n1NvOfzC/8eODBwpofKO0NXIbY2LNGrgbYfJkyMysn5aZ6aTX1NTw61//mtLSUmbOnMmmTScEPUZ5uVOiCDWxUWMbpa2R2xgTaxYsGimwGsrX/nPllVfy4osvcv/993PmmWd63qBFvEsUPqNGNS5vTR3l1hhjAlmwaIRJk6Cqqn5aVRX86U8P8OCDD3L11VdzxRVXAMFv3CLB2zwCvfpq4/LX1FFujTEmkDVwN0LwBu7lwBDgZPLyXuDUU9N49VWnnaBzZ2eL77/3bu8Ixr/B3BhjEsGe4I6hhlVLu4ALgA7Ao6xdm8bDDztBQRUqKmD7dnjySac3kjMPRWPOY4wxyRG3YCEij4nIRhFZ5pfWWUTeEpFV7msnv3UTRaRMRFaKyCl+6UNEZKm7boqISLzyHKmGVUu3AIuBR4BuQffx7wEVycx31sZgjGlO4lmyeBwYEZA2AZirqv2Aue5nROQg4DzgYHefh0TEN8/bw8B4oJ+7BB4z4fzbBGAecCfwW+D0kPv5uq7W379uRjvfq7UxGGOam7i2WYhIAfCyqvZ3P68EhqvqehHpAbynqvuLyEQAVb3D3e4N4GZgDfCuqh7gpp/v7n9ZuHMnYtTZbdu2MXjwYNasqaGmZglONZQ3eyjOGNPcNZc2i+6quh7AffXV2fQEvvbbbp2b1tN9H5gelIiMF5EFIrJg06ZNUWcumsEBAa677jpWr17NxIkzyMkJHSisWskYk8qaSwN3sHYIDZEelKpOU9Whqjq0a9euUWXAN2mQr1HaN2lQsIChqvzjH//goYce4sorr+TWW3/ZoKtqUZF1XTXGtByJHu5jg4j08KuG2uimrwN6+23XC/jWTe8VJD3mJk2CyspqnNqvnkDR7kZp/5t8VVUVV1xxBVOnTmXMmDHcfvvtgLONBQNjTEuV6JLFHOAi9/1FwIt+6eeJSFsR6YPTkF3qVlVtFZEj3F5Q4/z2iSmn8TkdKMWZa2KLX7rj+++/Z8SIEUydOpXrrruO559/nuzs7AbHirY6yxhjmr1gQ9HGYgFmAeuBKpwSwiVALk4vqFXua2e/7ScBXwIrgZF+6UOBZe66B3Ab5cMt0Q5RXjfs+GIFUbim3rDjX3zxhfbt21czMzN1xowZnseJ9XDjxhiTSNgQ5aH52iycgf0uBkrIylrBo4/2oWvXNznnnHPIzMzkhRde4Mgjj/Q8jg0PboxJZc2lN1SzVf/Zh1sRSWfw4Ils2fIAo0aNIi8vj08//TRkoAAbHtwY0zJZsPBTWOh8+1ftyaRJV/PJJ09zxRVXcOqpp/Lxxx+TH8E4HTY8uDGmJbJg4eEvf/kLhx9+OBMnTuRf//oX7du3j2g/Gx7cGNMSWbDwMGdOe7777hPuvPN29t03PeIeTTY8uDGmJbJpVYOo39hd94AeRHbTt2cujDEtjZUsgnAe0Kuf5j9qrDHGtDYWLIKwHk3GGFOfBYsgrEeTMcbUZ8EiCOvRZIwx9VmwCMJ6NBljTH3WG8qD9Wgyxpg6VrIwxhgTlgULY4wxYVmwMMYYE5YFC2OMMWFZsDDGGBNWi538SEQ2AUGmIQqrC7A5xtlp7uyaWwe75tahqdecr6pdAxNbbLBoLBFZEGyWqJbMrrl1sGtuHeJ1zVYNZYwxJiwLFsYYY8KyYNHQtGRnIAnsmlsHu+bWIS7XbG0WxhhjwrKShTHGmLAsWBhjjAmr1QYLERkhIitFpExEJgRZLyIyxV3/mYgcmox8xlIE11zoXutnIjJPRAYlI5+xFO6a/bY7TERqROTsROYvHiK5ZhEZLiKLRWS5iLyf6DzGWgR/2x1F5CURWeJe82+Tkc9YEZHHRGSjiCzzWB/7+5eqtroFSAe+BPYBMoElwEEB24wCXgMEOAL4d7LznYBrPhLo5L4f2Rqu2W+7d4BXgbOTne8E/J73BD4H8tzP3ZKd7wRc8/XAXe77rsD3QGay896Ea/4lcCiwzGN9zO9frbVk8QugTFVXq+ou4Cng9IBtTgeeUMcnwJ4i0iPRGY2hsNesqvNUdYv78ROgV4LzGGuR/J4BrgCeAzYmMnNxEsk1/xr4l6quBVDVVL/uSK5ZgfYiIkA7nGBRndhsxo6qfoBzDV5ifv9qrcGiJ/C13+d1blq026SSaK/nEpxvJqks7DWLSE/gTOD/JTBf8RTJ73k/oJOIvCciC0VkXMJyFx+RXPMDwIHAt8BS4EpVrU1M9pIi5vev1jpTngRJC+xDHMk2qSTi6xGR43CCxdFxzVH8RXLN9wHXqWqN86Uz5UVyzW2AIcAJQDYwX0Q+UdX/xjtzcRLJNZ8CLAaOB/YF3hKRD1X1p3hnLklifv9qrcFiHdDb73MvnG8c0W6TSiK6HhEZCDwKjFTVigTlLV4iueahwFNuoOgCjBKRalV9ITFZjLlI/7Y3q+rPwM8i8gEwCEjVYBHJNf8WuFOdCv0yEfkKOAAoTUwWEy7m96/WWg31KdBPRPqISCZwHjAnYJs5wDi3V8ERwI+quj7RGY2hsNcsInnAv4ALU/hbpr+w16yqfVS1QFULgNnAH1I4UEBkf9svAseISBsRyQEOB75IcD5jKZJrXotTkkJEugP7A6sTmsvEivn9q1WWLFS1WkT+CLyB05PiMVVdLiK/d9f/P5yeMaOAMqAS55tJyorwmm8EcoGH3G/a1ZrCI3ZGeM0tSiTXrKpfiMjrwGdALfCoqgbtgpkKIvw93wo8LiJLcaporlPVlB26XERmAcOBLiKyDrgJyID43b9suA9jjDFhtdZqKGOMMVGwYGGMMSYsCxbGGGPCsmBhjDEmLAsWxhhjwrJgYUwjuX3YPxKRkX5p57jdUo1pUazrrDFNICL9gWeBQ3D6+C8GRqjql004ZhtVTdlB7kzLZMHCmCYSkbuBn4E9gK2qequIXARcjjNk9jzgj6paKyLTcIaWzgaeVtW/ucdYB0wFRuCMV9ULuBSoApaq6gUJvixj6mmVT3AbE2O3AIuAXcBQt7RxJnCk+3TxNJwhKGYCE1T1exFpA7wrIrNV9XP3OD+r6lEAIrIeyFfVXSKyZ8KvyJgAFiyMaSJV/VlEnga2qepOETkROAxY4A6bkk3dcNHni8glOP97ewMH4UxEBPC032GXA8Ui8iKQymNVmRbCgoUxsVHrLuCMPfSYqv7VfwMR6QdcCfxCVX8QkWIgy2+Tn/3enwIcizOJzQ0i0l9Va+KWe2PCsN5QxsTe28A5ItIFQERy3RF9OwBbgZ/cWctOCbaziKQDvVT1HeBanGlAcxKSc2M8WMnCmBhT1aUicgvwtoik4TRS/x5YgFPltAxneOyPPQ7RBpgppJxzKwAAAElJREFUIu1xvtDdpapb459zY7xZbyhjjDFhWTWUMcaYsCxYGGOMCcuChTHGmLAsWBhjjAnLgoUxxpiwLFgYY4wJy4KFMcaYsP4/2N0icUTShXcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "'''B'''\n",
    "max_N=30\n",
    "SSE_train_al=np.zeros(max_N)\n",
    "SSE_test_al=np.zeros(max_N)\n",
    "for i in range(1,max_N):\n",
    "    Ztrain_al=np.hstack((Ztrain_al,al_train_year.reshape(-1,1)**i))\n",
    "    Ztest_al=np.hstack((Ztest_al,al_test_year.reshape(-1,1)**i))\n",
    "    A_al = np.linalg.solve(Ztrain_al.T@Ztrain_al,Ztrain_al.T@al_train_price)\n",
    "    SSE_train_al[i]=np.sum((al_train_price-Ztrain_al@A_al)**2)/len(al_train_price)\n",
    "    SSE_test_al[i]=np.sum((al_test_price-Ztest_al@A_al)**2)/len(al_test_price)\n",
    "plt.plot(al_train_year,al_train_price,'o',color='blue',label='Aluminum')\n",
    "plt.plot(al_train_year,Ztrain_al@A_al,'-',color='black',label='Best Fit')\n",
    "plt.legend()\n",
    "plt.xlabel('Years')\n",
    "plt.ylabel('Prices $')\n",
    "plt.title('Price vs Year Data')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Price vs Year Data')"
      ]
     },
     "execution_count": 387,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhU5fnw8e+dFSiLGgLKlqAsyhohUqgbtVoEFap1QXFp1R8uuKDQKkWtC1RFRcSlimJFCYpFaxVRQcStAhp4WQREUAERhACyBgJJ7vePcxImycxkJpmTmUnuz3XNlclZnzOZnPs8u6gqxhhjDEBCtBNgjDEmdlhQMMYYU8qCgjHGmFIWFIwxxpSyoGCMMaaUBQVjjDGlLCiYuCIi74nIVdFOhzG1lQUFE1Uisk5E9ovIXhHZIiL/EpGGgbZX1f6qOqUm01gZEWkuIttEpG+55f8SkVdr4PztRETdz3CviPwsIu+IyO/COMa1IvKxh8k0ccKCgokF56lqQ6AHcBJwV/kNxBGT31dV3QLcBjwvIvUB3BvyOcAtkTyXiCQFSUdD93M8EfgIeFtELo/k+U3tF5P/ZKZuUtWfgPeALgAi8rGIjBWR/wH5wLHusmtL9hGR/xORVSKyR0RWikgPd3kLEXlDRPJE5AcR8XtzFpHe7pN1os+y80Vkmfu+l4jkishuNyczPkDaXwFWA/e7geE54BZVzXOP00pE/uOTnmE+5+sjIgtEZKeIbBaRiSKS7K5LcnMBN4rIWuCbED7Hzar6OPAAME5ExD3WXSLyvftZrRCRge7yrsBTwKluTmObu3ygiCxxt98gIndXdm5TC6iqvewVtRewDjjTfd8aWAE84P7+MbAB6AwkAcnusmvd9RcBP+HkLgRoB2TgPOwsAu4BUoBjge+BfgHS8B1wls/v/wbudN/PB65w3zcEege5llbAduC/wFs+yxOBJcDf3PS0c6/7d+76k4Bfu9d4LPAtcJO7LglQ4H3gSKC+n/O2c/6VKyzv4O7b3v39YuAY9/O5DNgLNHfXXQt8XG7/M3ACdALQHdgGnBvt74y9vH1ZTsHEgrdEZCfwOfAJ8A+fdS+p6gpVLVTVQ+X2uxYYp6pfqWOtqq7Hucmmq+r9qnpQVb8HngcGBzj/q8ClACLSCBjgLgM4BLQTkaaquldVFwS6CFXdiBOIzgRu8FnVG2isqv9w07MWmFySHjf9C91r/B6YBJxe7vD/UNVfVHV/oPP7scn9eZR7ntfVyUUUq+o0nMCUHeR6PlLVr93tlwKv+UmXqWUClk8aU4P+oKofBlj3Y5D9WuM85ZeXAbRwA02JROCzAMeZBnwhIjcAFwCL3eACcA1wP/CNiPwA3KeqM4OkaQXwi6puLpeeNn7S8zGAiBwPPAb0BBrg/F8uLHfcYJ9DIC3dnzvc8/wJp+4jw13eEGgaaGcR6QM8iJNTSwFSORwsTS1lQcHEumDD+P4IHBdg+Q+q2j6kE6iuFJH1QH+cYpVpPuvWAJe6ldwXADNEJE1V94V6AW561qjqCQHWPwcsAC5R1b0iMhI4t3wywzhfifOBn4G1InIs8E/gd8BCVS0Ska9xit0CHf814FHgbFU9ICJP4QQSU4tZ8ZGJZy8AI0Wkp9s6qZ2IZABfArtF5A4RqS8iiSLSRUROCnKsaTgthU7DqVMAQEQuF5F0VS0GSp70i8JM53zgoIiMEJF6bnq6ikhPd30jYBewT0ROAK4L8/hluE1kb8FpxXWHqirOzVyBPGcTuRY43me3LUCrkgpun3TtcANCbwIXv5laxIKCiVuq+m9gLM4NfQ/wFnCUqhYB5wFZwA84FaQvAE2CHO5VoC/wkapu81l+NrBCRPYCTwCDVfVAmOksxKmn6IVTjr8NJ3fQ2N1kBHCVew3PAdPDOX6Jkn4KwDKgH3CBqr7spmEZMBEnYG7GCQi+RVRzgDXAFhH52V12A/CgiOzBqSR/vSrpMvFFnIcIY4wxxnIKxhhjfFhQMMYYU8qCgjHGmFIWFIwxxpSK634KTZs21czMzGgnwxhj4sqiRYu2qWq6v3VxHRQyMzPJzc2NdjKMMSauuJ01/bLiI2OMMaUsKBhjjCllQcEYY0ypuK5T8OfQoUNs3LiRAwfCGomgTqlXrx6tWrUiOTm58o2NMXVKrQsKGzdupFGjRmRmZuJOOGV8qCrbt29n48aNtG3bNtrJMcbEmFpXfHTgwAHS0tIsIAQgIqSlpVlOqg7KyYHMTEhIcH7m5EQ7RSYWeRYU3CGCvxSRpe58sPe5y48SkTkissb9eaTPPqNEZK2IrBaRftU4dyQuodayzyc0tekmmpMDQ4fC+vWg6vy8/HJo2vTwddWm6zVV52XxUQFwhjtpSDLwuYi8hzNRyVxVfUhE7gTuBO4QkU4447V3BloAH4pIB3cYZGNqVMlNND/f+X39eud3gCFDopeuqho9+vC1+Nq+3bmu//0PpkypPddrqs6znII7Z+5e99dk96XAIGCKu3wK8Af3/SDgNVUtUNUfgLU448/HpbFjx9K5c2e6detGVlYWCxcuZMKECeT7+88M0Z/+9CdmzJgRwVSaQPzdRPPzneXxaMOGwOvy82HSpNp1vabqPK1TcGeYWgJsBeao6kKgecn8te7PZu7mLSk7D+1GDs8x63vMoSKSKyK5eXl51U6jF1nm+fPnM3PmTBYvXsyyZcv48MMPad26dbWDgqk5gW6iwW6usaxNm+DriwLkx+P1ek3VeRoUVLVIVbOAVkAvEekSZHN/Bd0VZgBS1Umqmq2q2enpfofuCJm/ctahQ6sfGDZv3kzTpk1JTU0FoGnTpsyYMYNNmzbx29/+lt/+9rcAzJ49mz59+tCjRw8uuugi9u51MlaLFi3i9NNPp2fPnvTr14/NmzcHPJfxRqCbaGU311g1diw0aBB4fWKi/+Xxer2m6mqk9ZGq7gQ+xpnacIuIHAPg/tzqbrYRaO2zWytgk5fp8qqI4Pe//z0//vgjHTp04MYbb+STTz7hlltuoUWLFsybN4958+axbds2xowZw4cffsjixYvJzs5m/PjxHDp0iJtvvpkZM2awaNEirr76akZbHr7G+buJNmjgLI9HQ4Y4RURpaRXXNWjgPAzVpus1Vedl66N0ETnCfV8fOBP4BngbZz5a3J//dd+/DQwWkVQRaQu0x5lP1jNeFRE0bNiQRYsWMWnSJNLT07nkkkt46aWXymyzYMECVq5cycknn0xWVhZTpkxh/fr1rF69mq+//pqzzjqLrKwsxowZw8aNG6uXIBO2kptoRgaIOD8nTYrvStchQ2DbNpg6teJ1PfNM7bteUzVetj46BpgiIok4wed1VZ0pIvOB10XkGmADcBGAqq4QkdeBlUAhMMzrlkdt2jhFRv6WV1diYiJ9+/alb9++dO3alSlTppRZr6qcddZZvPrqq2WWL1++nM6dOzN//vzqJ8JUy5AhtfOmGOi6auv1mvB42fpomaqeqKrdVLWLqt7vLt+uqr9T1fbuzx0++4xV1eNUtaOqvudV2kp4VUSwevVq1qxZU/r7kiVLyMjIoFGjRuzZsweA3r1787///Y+1a9cCkJ+fz7fffkvHjh3Jy8srDQqHDh1ixYoV1UuQMcaEqNYNcxGOkqei0aOdIqM2bZyAUN2npb1793LzzTezc+dOkpKSaNeuHZMmTeLVV1+lf//+HHPMMcybN4+XXnqJSy+9lIKCAgDGjBlDhw4dmDFjBrfccgu7du2isLCQ4cOH07lz52perTHGVE5UKzTwiRvZ2dlafpKdVatWccIJJ0QpRfHDPqe6JScn8g8/Jn6JyCJVzfa3rk7nFIypC2pb72zjrVo3IJ4xpqyqNr22sZDqJgsKxpRT226GVWl67VXHThP7LCgY46M23gyr0ju7to39ZEJnQcEYH7F8M6xqDibcptc5Of7774CNhVQXWFAwxkesDoRXnRxMOL2zS84TiI2FVPtZUPBAYmIiWVlZdO/enR49evDFF19U6TjBRlXt27cvHTt2JCsri6ysrNIhtX/zm98AsG7dOqZNm1a1C6jDYnUgvOrkYD799FP+858LGTjwFl59dTrffHMgYKujQPMugI2FVGeoaty+evbsqeWtXLmywrKa9qtf/ar0/fvvv6+nnXZalY6TkZGheXl5ftedfvrp+tVXXwXcd968eXrOOecEXB8Ln1MsmjpVtUEDVed53Hk1aOAsjyaRsmkqeYkE3mfDhg06ePBgBTQ9PV1TUxsooNBOmzf/n99rCnQeiP5nYCIHyNUA91XLKXhs9+7dHHlk6YyjPPLII5x00kl069aNv//97wDs27ePc845h+7du9OlSxemT5/OxIkTKwy1HYqGDRsCcOedd/LZZ5+RlZXF448/HtmLqsVCKWqJRuukcHIwBw4cYOzYsRx//PG89dZbXHDBvRQVraOgYBcwCyhky5ZTuPzyv5GWpmXSH+g8GRnWp6GuqNWd14YPH86SJUsiesysrCwmTJgQdJv9+/eTlZXFgQMH2Lx5Mx999BHgzJ+wZs0avvzyS1SVgQMH8umnn5KXl0eLFi149913Adi1axdNmjRh/PjxzJs3j6ZNm/o9z5AhQ6hfvz4Ac+fOJc1nXOSHHnqIRx99lJkzZ0bisuuUYAPDRasj2NixZc8LFYtzVJV33nmH2267je+//54//vGPnHLKo4wenemzX39gGTAceJAdO5oydOjtpekP5TymdrOcggfq16/PkiVL+Oabb3j//fe58sorUVVmz57N7NmzOfHEE+nRowfffPMNa9asoWvXrnz44YfccccdfPbZZzRp0iSk8+Tk5LBkyRKWLFlSJiDUZV4/xUerdVJlOZg1a9bQv39/Bg0aRL169ZgzZw4zZsxgwoRMP3UEjYAXgPOBv5Kf/0lp+uNpyPDa1p8kZgQqV4qHVzzUKaiqNmvWTLds2aK33367Pvvss3732b59u77yyit68skn63333aeqVatTKDl3XaxTqIn6gKqU7Xtt06ZNmp6erk2aNNEJEybowYMHK02v89ql0EGhmcLG6F1AFcRq3U+8wOoUouebb76hqKiItLQ0+vXrx4svvlg67eZPP/3E1q1b2bRpEw0aNODyyy9n5MiRLF68GKDMUNvhqs6+8aomnuJjrXVScXExV1xxBXv37uWLL77g1ltvJTk5OcR0NQbeBPaRmnoxBw8e9Di1kRPL/UninQUFD5TUKWRlZXHJJZcwZcoUEhMT+f3vf89ll11Gnz596Nq1KxdeeCF79uxh+fLl9OrVi6ysLMaOHctdd90FwNChQ+nfv39YFc0lunXrRlJSEt27d68zFc010ccg1qbpfOSRR5g7dy4TJ06kU6dOFdZXNjczdCYlZTIFBV8wcuRIz9IZabHan6RWCJSFiIdXrBYfxYPa+DmlpfkvJsnIiOx5pk51jini/IxWkcWCBQs0KSlJL7roIi0uLg64Xfn03nBDxfTfdtttCujUOCl/yciomb91bUWQ4iPPbthAa2AesApYAdzqLu8OzAeWA+8AjX32GQWsBVYD/So7hwWFqqttn9PUqarJyRVvEikptbOceefOndq2bVvNyMjQX375pdrHO3jwoJ566qlav359XbVqVQRS6C2rU6ieYEHBy+KjQmCEqp4A9AaGiUgnnGYPd6pqV+A/wF8A3HWDgc7A2cAz7vzOxlRq9Gg4dKji8kaNYrPlTHUNGzaMDRs2MG3aNI444ohqHy85OZnp06eTkpLCiBEjIpBCb8VTK6l44+UczZtVdbH7fg9OjqEl0BH41N1sDvBH9/0g4DVVLVDVH3ByDL2qeO7qJL3Wq42fT6Cy5B07/C+PJ+WbXv7lL7PIycnh7rvvLh3WJBKOOeYY7rrrLmbNmsXs2bMjdlyvDBkC69ZBcbHz0wJCZNRIRbOIZAInAguBr4GB7qqLcIqZwAkYP/rsttFdVv5YQ0UkV0Ry8/LyKpyrXr16bN++vVbe+CJBVdm+fTv16tWLdlIiKtZaBUVKxYHw8nnssWG0aHECo0aNivj5br75Zo477jhuv/12CgsLI358E/s879EsIg2BN4DhqrpbRK4GJorIPcDbQEk7OPGze4U7u6pOAiaBM0dz+fWtWrVi48aN+AsYxlGvXj1atWoV7WREVG3tiVux6eUYVNdRVPQJKSkpET9famoq48aN449//COTJ0/muuuui/g5TIwLVNkQiReQDHwA3B5gfQfgSz1cyTzKZ90HQJ9gx/dX0WzqrlhpFRTJtJTtfLZGIUXhSk87yxUXF+tpp52m6enpumvXLu9OVEfE0veyBFFqfSTAy8CEcsubuT8T3PVXu793BpYCqUBb4HsgMdg5LCiYWBTJljFlm16ep9BQYZPnTS+/+uorBfTee+/19kS1XKy2kopWUDgFp/hnGbDEfQ0AbgW+dV8PAeKzz2jgO5wmqf0rO4cFBROLItmG/vBNZZY77PW4GrupXHDBBdq4cWPdtm2b9yerpWK1P0VUgkJNvCwoGF+xkk2P9PhIL71UoElJHRQ6aJs2BTV2XcuXL1cR0TvvvLNmThiCWPkbhyoWx8pSDR4UbJgLUytUZ7rKSIt0S6ht256ksPBb3n33cdavT6mxppddunTh0ksvZeLEiWzZsqVmThpELP2NQxWXreICRYt4eFlOwZSIpWx6JMuRN2/erI0aNdJzzjmn9CkZVBMTD1+fl0/L3377rSYmJurw4cNVNbpP6rH0Nw6V1SlYUDBREmvZ9EjdPP/85z9rcnKyPvrotxVuLjV1k7nmmms0NTVVn3jixxq/wfl+joGGAI92UUxlYrHIK1hQEGd9fMrOztbc3NxoJ8PEgMxMpzjB8T0wHfiC5OQNtGy5m9atW5OZmUm3bt249tprIzI0hNe+/PJLfv3rX/PXv/6V6dMf9rm+ijIynF69Xli/fj3t27cnNfVq9u59tsbOXX6Wu0C8vPbaSkQWqWq2v3VWp2BqhbFjoV69lTjDZh0H/A2R7+naNZOTTz4ZEeHjjz/mL3/5C23btmXMmDExPd9EYWEh1113HS1atOCuu+6qdEhoL4eMzsjIYOjQoezdOxkn4NbMuf3NmVBebeigGHMCZSHi4WXFR0ZVtaCgQO+9915NTEzWhISjFO7Xli03+M2m/7//9/904MCBCmhaWpo+/PDDunfv3krPUdNFAJde+ojbBPUNzcgIPCx4TZWr//TTTypST+GqGjt3ZUVG4fwdYrEIJ5qwOgVTG02dqnr00YsUOiugffpcqk8/vSWkf/6FCxdqv379FNDmzZvrn/70jLZpU+x3v5quLBw//nuF+goDFYoVnGHBU1L83yBrquJywIARCgkKq2rk3JGqWI7Vyt5osqBgap2pU1VTUz9Q+JVCS4WZfm+clf3zf/bZZ3r88ae7T+V/VGfe4rL71WSrl+LiYq1X72x1ei5vKHO+tLSab33ka+vWrZqa+itt0OCSGnnijtTNPB5bLXnNgoKpdZo2fU0hWaG7wqZqFa20aVOs8IhCokIXhR/L7Bduy6bqFFW8+uqrboCaEJOtbEaPHq2ALl26tEbOF4lin3htteQlCwqmVvnXv/6lIAqnKvwSNCCE8s9/+Kb/oUIjhdYKK0r3C+dJszpPtzt27NBmzZppSkq2QmFMPtnu2LFDmzRpooMGDYp2UkIydWrgoB4Ln2e0WFAwMSvcJ8HPP/9ck5OTtV69MxXyKw0IJUUtwY5b9qa/WKG5wpHavPnnpWkM9UZfnaKK//u//9PExEQdM2ZxTJeBP/DAAwrol19+6fm5/H0/CgoKdP78+frWW2/pSy+9pGvXrg24f6C/h0jsfJ7RYEHBRFSkWnKE+1S9YcMGbd68ubZv316fe25HhX2rWhl7ww3lt/9eob0mJtbTt956K6xrrmonuk8//VQBHTFiRFjnC3fbSNi9e7empaVpv379PD1Pxe9HsaakzNBmzY5zi9gOv3r37q1PPfWU5ufnlzlGsBZMXvP9u6SlOa9Yaf1kQcFEjHfDQh9++Xuqzs/P1549e2qjRo105cqVpWkpfzOcOvVwJWyoT+v+07FVU1JO0oSEBH3uuec8uaYSW7du1czMTM3IyNA9e/aEfC7V6LWseeQRp8nsRx995Nk5yn6W8xV+o4AmJ3fWV199VRctWqQrVqzQcePGabdu3RTQTp06lX4/Kh4jtL9HJPj7u8RSzs+CgomYSP6ThfpUXVxcrEOGDFER0bfffjtix61se9irAwYMUEDvv/9+LS4urvTc4d6kCwoK9LTTTtPU1NQqFcdE66a3b98+bdu2rR599NG6bt06T87h/F0KFW5xcwRHKzyvcMjv9u+99542b95c09PTdfny5aoavaAZ6O9Sk3+jYCwomIiJ5BhDod7QHn30UQV0zJgxET2uauU5i4MHD+qVV16pgN50001aVFRU6fnL52BuuMF/8U5xcbEOHTpUAc3JyQnp2sqL5phPK1as0CZNmmiXLl10586dET9+mzb7FS50A8ItCnsqvZmuXr1aW7Rooenp6bps2TJVjU7HtWDFVr5/o2h1qrOgYCLGmwlkDr/KP8W9//77mpCQoBdeeGFIT+qBjivi3Jwr285fOoqKinTkyJEK6MUXX6z79u0rc4xg/9TBrvHpp59WwO98BaHeLKLdBn/OnDmalJSkv//97/XgwYMRO+7OnTt9+o88XulTvu/n1aLFt3rkkS21adOmumLFioilKRyh5BTS0qLXqS4qQQFoDcwDVgErgFvd5VnAApyZ2HKBXj77jALW4sy81q+yc1hQqHmRzo4Hu/mtWbNGjzjiCO3WrVvYZe033FDxaa18OgP94wZqrfTII4+oiOiJJ56o69evD+mzCHSOZs3mamJiop577rlaWFhY4TMJ9TOOhd66L7zwggJ63XXXhRy4g9mxY4dmZ2drcnKy3njjtEqDo7/PoF69NdqkSXPt2LFj2N+dSAilTiHQ0CU1EdCjFRSOAXq47xvhTL/ZCZiNO9UmzvScH7vvO1F2jubvsDmaY1JVs7zh7Ldr1y494YQTNC0tTb///vuwjxHKE3RVil5mzpypjRs31vT0dG3efF4Vz/G+wq/0hBNO0F27dlUp7b5iYVyfO++8UwF95JFHqnWcHTt2aM+ePTUlJUXfeeedkPYJ9Hk1b/6RJiQk6GWXXRaRYBWuylofRbPoLyaKj4D/AmcBHwCXuMsuBaa570cBo3y2/wDoE+yYFhTiRzhPtEVFRTpw4EBNTEzUuXPnVukYofzDVbXoZdWqVdqhQwe3aONqha0hnqNY4SmFJE1OztJNmzb5PX6szQ0RiqKiIr3oootURPSNN96o0jEWLVqkXbp00ZSUFJ05c2bI+wX7vMaMGaOAPvvss1VKk5eiWfQX9aAAZAIbgMbACe77H4GfgAx3m6eAy332mQxcGOy4FhRiT6Cn1nD+AUqaOz7xxBNllodzjFC2rU7Ry549e7Rx478qJCkc4d7sK1aEHi7GWqpwlgKakNBfJ02qmEOoynXGkvz8fO3du7fWr19fFy5cWGZdsNxMfn6+3nnnnZqYmKhHH320zp49O6zzBvu8ioqKtF+/fpqamqqLFy+u7iWGJNScWzSL/qIaFICGwCLgAvf3icAf3fcXAx+675/2ExT+6Od4Q926iNw2bdp4+LGZcAX7kof69Lt48WJNTk7W888/v0KWP5wn6FD/4apT9DJ1qmq9eisVznBzDYmakHCSDhgwQv/zn//obbf9V5OTH1c4013fUOGfev31wYsyYqGeoKq2bNmimZmZ2rx589KmqsGu59NPPy3NdV1zzTW6Y8eOsM9Z2ee1detWbdmypR533HG6e/fuSF5u2Gnxt32dan0EJLvFQLf7LNsFpTO+CbDbfW/FR3HI90sdrGlnKE+/+/bt0+OPP15btGih27Ztq3Cu6pS1e9WjdOrUkgH15mrjxqO1Y8dTNSUlxQ0CJa/jFMYq7Aj5iT8W6gmqauXKldqkSRPt3Lmz7ty5M8DfbZc2bHijApqZmalz5sypdoD2t2/JcvhUIUFPP/3qSF9uGfGSy4tWRbMALwMTyi1fBfR13/8OWOS+71yuovl7q2iObZW1sPB9kg/lCeqGG25QQOfMmRPy+UJ5gq7pJ+/9+/frF198ofCV+hvBNZbrBiJl7ty5mpSUpGeddZbCwXKfwSx1Bh0UHT58uO7du9eTv1HFY/5NAR08eFzkLrSceKkPilZQOMV9UlqG0/x0idva6BS3OGkpsBDo6bPPaLfV0WrcFkrBXhYUouPw01dor5KnpGBPgm+//bYC2rjxiEqbH4b7NBmtp7d4eWr0yuTJkxXQevX6q9MT+TOFy93c0wl69NFflG4bic+q/HejYpPPQwqXKKDjx4+v0jEr+76FWpcV7Vxg1CuavXpZUKh5oeYOwnna27x5szZq1FRFshQORPxpPlpPbzVRrxHrHnzwQU1JOUIPF6UlKdyt9esfKHOd1f0bhf69PKQlvaQnTpwY9Hj++hFU9p30l4769Yt18uR9unPnTp08eZ/Wr19cY7nWQCwomIgJJYeQmBj6Da6kdYgz/+9KT56qo/nEXp0ez7XB1Knq3gS/Vfivwrd+e5dX928UTs4VDmr9+n9QQJ955hm/aQ4WYIKlqbi4WO+770tt0uR+hX6anHyCpqb+yicoopCizjhOZyrkKOyv8dyjBQUTMZWN6RLuDe2xxx5z/1Ge8expPpZvvLW9iCnU6wt0I05LC+3vFMpYQ2VfBXreeecpoJMmTQopzcG+kwcPHtQXXnihdLRWEdHu3bvrBRdcoMOHD9eHHnrI/a4/qPBXhWsV2rrf/SMVHqsbrY+8fllQqHnB/mH8PQEGk5ubq4mJyZqY+ActmaDeq5tjrBbRxEvFZFWF24y4KkU2quHmFJztDxw4oP3791cR0RdffLHSNAf6Ts6ZM0fbt2+vgGZlZemkSZM0Ly8vhHQWKcxV6O8Gh6frRj8FL18WFGpedbLWvnbv3q3t2rXTxMRWCtsjkuuIR5ZTqN72JcKp6/L9Xu3fv1/POussFRF9+eWXg6ah/L55eXmlI+i2a9dO33nnnUqH0/CfzkMK56kzxezMutGj2auXBYXomDo18D9NqE+4V111lfy3LGoAACAASURBVCYkJCh8HPBYkQwIsZpTiOWirUgI9/qqk3MKpVWcv799fn6+nnHGGZqQkKDTpk2rtCiruLhYX375ZU1LS9OkpCQdPXp0hRnfAqXNt89M2ePvU+iq0EbhQCgfbbVYUDARV5V/vhI33pjjZpfvCXuWtKqI9RtvrAasQMIZxqHke1Lyd67s+iKRc6rKMfbu3aunn366JiQk6PTp0/1eY3Fxsf7tbx9pauppCmhqah998MHDk/kEutZA37+KgeF9BTQtbUroF1tFFhRMxIWSXfd34x0//juFRgqnuNnm0ParjtpeRFOTwmlmG62OhlU9xp49e/SUU07RxMREff3110uLgvbt26fvvfeeduhwivswc4zCPxWKtEEDpx4t2LwcgYbIrjifQrGKnKBt2/b0fFRXCwrGE1OnBh7awt+NNz8/X1NSeqgzgNz6MtuF04w1XF5X5sbbk351hBpgqxOII/F5VvUYzz+/W1NT+6gzcOERetRRrbSkKalT//WUwv4K391g/wOBXr4zrx0+zj8V0Pvv/yr8iw6DBQXjmcpaapTceIuLi91KOacyzasbtD9e5hSq8lQaz0Ek1AAbj62qDv8td6vTRPp6TUy8Qi+88AF988031bdjZSRevj39D3+HNiugycmPefq9sKBgPBNK3YKq6sSJExVwO/V4c4MOxMs6haoM0hfL9RuVqYmcQrRUluZA66uSU/D9m1c8bqbCRZ5+VhYUjGeC1S2UfPE/+eQTTUxM1IEDB+rLLxdF5abo1dN5uE/E8Xiz9OV1nUI0Vfa3DHRNweoUAr2CD/ExWKGVp7kqCwrGE/5aXJRvefHjjz9qs2bNtGPHjrpz584y+wUb+iFeilfCvcnHY7FKeeG2PoqHv6NqeIPZlf+u33BD4P+F8HNVTyigLVv+6Nm1WlAwERfKk+D+/fu1V69e2rBhQ125cmXEjhtLwk1vvOcUarNI54JC3e6GG8p/HxYqoL///QzPrtWCgom4ym5uRUVFesklzjDFb775ZsSOG4vCeSKOt6AXSLzlAkIVynWF8x2t2vEKFFK1ceMRAdP5008/6aFDh6pyiapqQcF4IFirIxHVxo1HKKDjxoU3oUltKF6pTLzfUGtLYKuqSH9H/R/vNwonB9ynS5cuOmjQoCpeQfCgkIAxVdCmTeB1qk+we/djJCXdxDHHjIzIcYOdL94MGQLr1kFxsfNzyJBopyg8o0dDfn7ZZfn5zvK6INLfUf/79UZkEQcPHqywpqCggG+++YbOnTtX7YSVsKBgqmTsWGjQwN+aHOA24A8UFk7grruk2sdt0MBZbmLDhg3hLa9t/H1HRWD9esjMhJyc6h8vJaU3qgdYunQpOTnOcRMSnJ9XXPENhYWF/OMf3ap0vkoFykJU9wW0BubhzMm8ArjVXT6dw9NzrgOW+OwzCliLMx1nv8rOYcVH0eVbDOJkeV9XSFDoq84AX4fLWsMdpiCei1dqu3is94k031ZI5Yt/qlKUVv47/8QTGxTQK6+c6Ke568sKKKyo8vmI0hzNxwA93PeNgG+BTuW2eQy4x33fCWfe5lSgLc5czYnBzmFBIXakp7+lzlSLJyvsqXDDqEtlzrVdXa9T8OVlgGzZsqU2aHCZn+P/RZ3Z2w5V+XzBgoJnxUequllVF7vv97g5hpYl60VEgIuBV91Fg4DXVLVAVX9wcwy9vEqfiZxZs2axY8dFJCT0AGYBDStsU5fKnGu7IUNg0iTIyHCKTTIynN/jrW4kErwsSuvduzf5+fP9rFkGdAaSInq+EjVSpyAimcCJwEKfxacCW1R1jft7S+BHn/Ub8QkiPscaKiK5IpKbl5fnTYJNyP7973/zhz/8ge7du/LPf75PRkbjgNvWlTLnQMqXDUe8LLgGxXtleaREutLZ9zsyd25v4AdgS7mtlgFdI3I+fzwPCiLSEHgDGK6qu31WXcrhXAKAvxpJrbBAdZKqZqtqdnp6emQTa8Ly/PPPc8kll9CrVy/mzp3L0KFHsm6d8+ToT21qQRSunBwYOtSpjFR1fg4dGt+BwUS2YUT578jOnX0ASEryfZbeBmwGulX7fAEFKleKxAtIBj4Abi+3PAkn/LXyWTYKGOXz+wdAn2DHtzqF6HnooYcU0P79++u+ffvKrLMy54qscrb2ilTDiIrfkXyFJK1X787S4zdv/pEC2qzZ7GqdjyhVNAvwMjDBz7qzgU/KLetM2Yrm77GK5phTXFysf/3rXxXQwYMHa0FBgd/trAVRWXWhU56pHv/fkWyFvqXbTJgwQQHdvHlztc4VLChUWnwkIs1E5Ffu+/oiMlpEHhKRYyrZ9WTgCuAMEVnivga46wZTtugIVV0BvA6sBN4HhqlqUWXpMzWnqKiIoUOHMm7cOK6//nqmTp1KSkqK322tzLmsutApz1SP/+9CH0S+orCwEIBly5aRnp5O8+bNPUtHKHUKrwFp7vv7gHbAL8C0YDup6ueqKqraTVWz3Ncsd92fVPVZP/uMVdXjVLWjqr4X3qUYLxUUFHDppZfywgsvMHr0aJ555hkSExOjnay4YZ3yTGUCd2Lbx4oVKwBYvnw53bp1w2m86Y2gQUFErgKOA/q67y8BcoGfgQwRuVJEugU7hol/+/btY+DAgfz73//mscceY8yYMZV+KX1bUTRt6rzKv4/3FjjhsGacpjL+viMPPdQbgPnz51NUVMTXX39N165dKzlSNQUqV3KKncgAvgFOAs4EvgTauMu/ct83CXYML19Wp+C9AwcO6JlnnqkJCQk6efLkkPYJNvGOdWozJnTFxcXarFkzveqqq3T16tUK6Isvvljt41LVOgVVXQ88AczEKe+/X1U34DQV3aaqG1R1lyfRykRdYWEhl112GR9++CGTJ0/m6quvDmk/fwOmBWKd2owJTETo3bs3CxYsYPny5QB06+Zt4UyldQqq+k+cIqRWqjrTXbwdp5+BiVOVdaRSVYYOHcqbb77J448/zp/+9KeQjx1uJ7W63qnNmGB69+7N6tWr+fjjj0lISKBTp06eni+p8k1AVfeW+32fN8kxNaGkk0zJ03xJRypwyjVVlZEjR/Kvf/2Le+65h+HDh4d1/DZtnGOGs70xxr8+fZxObFOnTqV9+/bUr1/f0/PZ0Nl1hG/O4Kqrgo+HP3bsWMaPH8/NN9/MvffeG/a5Ag+rXZG1wDEmuOzsbBISEti5c6fnRUdgQaFOKN99vihA748NG+Dpp5/m7rvv5oorrmDChAlVavpWvhVFWprzKv/eWuAYU7mGDRuWtjiqiaAQUvGR23ltv6oWi0gH4HjgPVU95GnqTESEWvF71FE53HTTTQwcOJDJkyeTkFD1Z4YhQ+xmb0yk9OnTh6VLl3rfHJXQcwqfAvVEpCUwF/gz8JJXiTKRFUpFbkrKO/zyy1X07duX6dOnk5yc7H3CjDEhOeuss0hJSSE7O9vzc4UaFERV84ELgCdV9XycSXFMHAhUkZuY6BTjNG/+CaoX06PHibz99tvUq1evZhNojAnq/PPPZ/PmzbRsWWE2gYgLOSiISB9gCPCuuyykoicTfYGGWJgyBb78Mpf8/PNo164t7733Ho0aNYpOIo0xAYkIRx11VI2cK9SgMBxnaOv/qOoKETkWZ/5lEwcCDbFw9tnbGTRoEEcddRSzZ8+madOm0U6qMSbKQgoKqvqJqg4EnnJ//15Vb/E0ZSaiyo9aetllSv/+Q9m0KY/169/ilFNa+R2HqDbNFmaMqVxIQUFE+ojISpx5lhGR7iLyjKcpM5667rqX+OqrN4GxQJbfmcBstjBj6p5Qi48mAP1whrdAVZcCp3mVKOOt7777jhdeuAXoC9xeurz8OET+mrKGOlaR5TCMiU8hVxar6o/lOjLZBDhxqLCwkMsvvxzVRJyJ8crOieDbfDVQU9bKmrhWNoyGMSZ2hZpT+FFEfgOoiKSIyEjcoiQTP3JyID39cRYsWIDIs0DrCtv4Nl+t6mxh1clhGGOiK9SgcD0wDGgJbASy3N8DEpHWIjJPRFaJyAoRudVn3c0istpdPs5n+SgRWeuu6xf+5ZhAcnLg2mt/ZOfOe4HzUB1cYZvy4xBVdbawquYwjDHRF+ooqdtw+iiEoxAYoaqLRaQRsEhE5gDNgUFAN1UtEJFmACLSCWfu5s5AC+BDEemgNk9zRIweDQcO3IYzFcYTpcsTE50WSW3aODd73+KdkvejRzs3dH/b+BNolFQbDdWY2Bdq66MpInKEz+9HisiLwfZR1c2quth9vwenuKklcAPwkKoWuOu2ursMAl5T1QJV/QFYC/QK94KMf+vXvw+8AYwG2pYuLy4+3EzV382+fFPWUOoEbD5iY+JXqMVH3VR1Z8kvqvoLcGKoJxGRTHf7hUAH4FQRWSgin4jISe5mLYEffXbb6C4rf6yhIpIrIrl5eXmhJqFOO3DgAElJN+F89CPLrPPi6d3mIzYmfoXa+ihBRI50gwEiclSo+4pIQ5xH1OGqultEkoAjgd44cz+/7vaQ9jdGs1ZYoDoJmASQnZ1dYb2p6Mknn6Sw8DtSU2dTUJBautzLp3cbJdWY+BRqTuEx4AsReUBEHgC+AMZVsg8ikowTEHJU9U138UbgTXf+6C+BYqCpu9y3OUwrYFOI6TMB/PLLL/zjH/+gf//+TJ58lj29G2OCCrWi+WURyQXOwHmiv0BVVwbbR5xODZOBVao63mfVW+5xPnbnZkgBtgFvA9NEZDxORXN74Mswr8eU8/DDD7Nr1y4efPBBune3IGCMCS5oUBCRxm6Rz1HAz8A0n3VHqeqOILufDFwBLBeRJe6yvwEvAi+KyNfAQeAqVVVghYi8DqzEabk0zFoeVc/GjRt54oknGDJkCN27d492cowxcaCynMI04FxgEWXL98X9/dhAO6rq5/ivJwC4PMA+Y3EG4zERcO+991JcXMwDDzwQ7aQYY+JE0KCgque6xUCnq6p1PYojq1ev5l//+he33HILmZmZ0U6OMSZOVFrR7Bbt/KcG0mIiaNy4caSkpNC+/SgbmM4YE7JQWx8t8OlPYGLcxo0beeWVVzj11Gv5y1+a2dDXxpiQhRoUfosTGL4TkWUislxElnmZMFN1jz/+OMXFxaxcOcIGpjPGhCXUzmv9PU2FiZgdO3bw3HPPMXjwYKZNy/S7jQ1MZ4wJpLImqfVwRkhtBywHJqtqYU0kzFTN008/zb59+7jjjjv4/HMbmM4YE57Kio+mANk4AaE/Ts9mE6Py8/OZOHEi55xzDl27drWB6YwxYaus+KiTqnYFEJHJWA/jmPbiiy+ybds27rzzTqDqQ18bY+oucVqcBlgpslhVewT6Pdqys7M1Nzc32smICYcOHaJ9+/a0atWKzz//PNrJMcbEMBFZpKrZ/tZVllPoLiK7S44D1Hd/F5wuDI0jmE5TDdOnT2f9+vU89dRT0U6KMSaOBc0pxDrLKTiKi4vp1q0bIsLSpUtJSAi1pbExpi6qTk7BxIFZs2axYsUKXnnlFQsIxphqsTtIBOTkENWhJB566CEyMjK45JJLavbExphax3IK1ZST4wwdUdJzuGQoCaiZVj6ff/45//vf/5g4cSLJycnen9AYU6tZnUI1ZWb67yCWkeFMdO+18847jwULFrB+/XoalO+UYIwxfgSrU7Dio2oKNGRETQwlsXr1ambOnMmwYcMsIBhjIsKzoCAirUVknoisEpEVInKru/xeEflJRJa4rwE++4wSkbUislpE+nmVtkgKNGRETQwlMWHCBFJTU7nhhhu8P5kxpk7wMqdQCIxQ1ROA3sAwEenkrntcVbPc1ywAd91goDNwNvCMiCR6mL6I8DeUBMDevd5WOG/bto0pU6Zw+eWX07x58wrro135bYyJT54FBVXdrKqL3fd7gFVAyyC7DAJeU9UCVf0BWAv08ip9kTJkCEyaBGlpZZdv3+7t3AXPPfcc+/fv57bbbquwrqTy2+ZRMMaEq0bqFEQkEzgRWOguusmdl+FFETnSXdYS+NFnt434CSIiMlREckUkNy8vz8NUh27IEGjYsOJyr+YuKCgo4KmnnuLss8+mc+fOFdaPHo3No2CMqRLPg4KINATeAIar6m7gn8BxQBawmcMjr4qf3Ss0jVLVSaqararZ6enpHqU6fJGqcA6l2OfVV1/l559/5vbbb/c0LcaYusfToCAiyTgBIUdV3wRQ1S2qWqSqxcDzHC4i2gi09tm9FbDJy/RVRaCbdiQqnEMp9lFVxo8fT5cuXTjzzDPDOqfNo2CMqYyXrY8EmAysUtXxPsuP8dnsfOBr9/3bwGARSRWRtkB7Ymyo7mA37UjMXRBKsc/cuXNZvnw5t99+O85HXJHNo2CMqTJV9eQFnIJT/LMMWOK+BgCv4EzaswwnEBzjs89o4DtgNdC/snP07NlTa1JGhqoTDsq+MjKc9VOnOu9FnJ9TpwY/Xvnt/R0bnPUl+vfvr82bN9cDBw6EdezK0mKMqTuAXA1wX7UezWEI8GAOOLfvcJQfHqPk+P6OU9I7euXKlXTu3Jn777+fu+++O7wTGmOMy3o0R0higF4TiYnh9wvwV1TkLyA0aAADBjjH7Nx5AiL1aNr0+iqk3hhjKmcD4gXx2Wef8fPPP3PuuedSv359ior8b1dUFP6geKG0BEpLg4svhilTID8/D3gZ1T8xcmQ6jRvbtJrGmMiznEIAxcXFXHTRRVx88cUcffTRXHPNNTRv/glQXGHbxMTw+wWE0hKoYUOYNavk2P8ECoDh1ufAGOMZCwoBLF68mC1btjBixAjOP/98Xn/9dbZs6YvTMGo08C3gFO8EykEEyw0EGh6j/P7OMQ4ATwPnAMdXemxjjKkqCwoBvPvuu4gId9xxBy+99BI///wzOTk5dO3aCXgI6MzRR89n0iSnItifYLmBkuExAu1bsr9zjEnAVmBkSMc2xpiqstZHAfz6179GRFiwYEGFdZs2baJ3794cddRR5ObmMn16UoWWRA0aODf9UMr9/bVEKtm/oCCfa645DieHMC/sYxtjTHnW+ihMW7du5auvvuKcc87xu75FixY88cQTLF26lCeffLLMU7+I8zOcm3aw/XfufBb4mebN76vSsY0xJhyWU/Dj5Zdf5qqrriI3N5eePXv63UZVOe+88/jkk09YtWoVrVq1ing69u3bR9u2benevTtz5syJ+PGNMXWT5RTC9O6773L00Udz4oknBtxGRJg4cSKFhYV+h6+OhKeffpq8vDzuu+8+T45vjDHlWVAo59ChQ3zwwQcMGDCAhITgH8+xxx7LXXfdxYwZM3j//fcjmo49e/Ywbtw4zj77bH7zm99E9NjGGBOIBYVy5s+fz65duxgwYEDlGwMjR46kY8eODBs2jP3790csHU8++STbt2/nvvvus1nUjDE1xoJCOe+++y7JycmcddZZIW2fmprKM888w/fff8+DDz4YkTTs2rWLRx99lHPPPZc1a3pVGJn1z3+Gpk0tSBhjIs+CQjmzZs3i1FNPpXHjxiHvc8YZZ3DZZZfx8MMP8+2331Y7Dffccw87d+7kvvvu8ztG0qFDznSfNtWmMSbSLCj42LBhA19//XXIRUe+HnvsMerXr8+wYcMIp0VX+aKh++//iieffJIbb7yRHj16hNRz2Ya9MMZEigUFH++++y5AwP4JwRx99NGMHTuWDz/8kNdeey2kfSpO2lPIvfdeR5MmzrEg9J7LNuyFMSYSLCj4mDVrFsceeywdO3as0v7XX389vXr14sYbb2RDCHfpikVD41H9fyQlPUGTJk2A0MZIAhv2whgTGRYUXPv372fu3LkMGDAg4DSXlUlMTCQnJ4fCwkIuv/xyCgsLA26bk+PkEA77CPgbcAHbtl1YurR8b+e0NEhJqXi8vXutXsEYU31eztHcWkTmicgqEVkhIreWWz9SRFREmvosGyUia0VktYj08ypt/nzyySfs37+/SkVHvtq1a8ezzz7LZ599xvDhw/3WL5QUGx32A3Ax0BF4iYyMskFpyBBn5rXiYti2DV580QkOvrZvtwpnY0wEBJqns7ov4Bigh/u+Ec5Y053c31sDHwDrgabusk7AUiAVaIszV3NisHNEco7mUaNGaVJSku7fvz8ixxsxYoQC+ve//73CurLzMf+k0EXhCIVvtUGD0OZTrmy+aGOMCYQgczR7NvOaqm4GNrvv94jIKqAlsBJ4HPgr8F+fXQYBr6lqAfCDiKwFegHzvUqjrx9++IGMjAzq1asXkeM98sgjpZ3Ptm3bxmOPPUZqaioA69cX4MTDVcBNwE7gLaB9yIPdBaqysApnY0x11EidgohkAicCC0VkIPCTqi4tt1lL4Eef3ze6y8ofa6iI5IpIbl5eXsTSuG7dOjIzM0t/r24vYhHh+eefZ8SIETz99NMcf/zx9O3blzZt2gD1cYqK/oAzk9tnwO/IyAh99NNAFctW4WyMqQ7Pg4KINATeAIYDhTjTlt3jb1M/yyoUyKvqJFXNVtXs9PT0iKXTNyhUbCpatfL6pKQkHn30UebMmUPr1q0pLCykb9++nH/+PaSkvAR8CnwDZNGggdPSKFT+WiWFewxjjKkgULlSJF5AMk7dwe3u711xphBb574KgQ3A0cAoYJTPvh8AfYIdP1J1Cvv371dAH3jgAVWtmfL6qVOd44k4P0OpR/DiGMaYuocgdQpetj4SYDKwSlXHuwFouao2U9VMVc3EKSLqoao/A28Dg0UkVZyJkNsDX3qVPl8lfQoy3LkxI1leH6gYyrdF0bp1VZs0JxLHMMYYX14WH50MXAGcISJL3FfA8SNUdQXwOk5F9PvAMFUt8jB9pdatWwdQWnwUqfL6SBVDGWNMTfEsKKjq56oqqtpNVbPc16xy22Sq6jaf38eq6nGq2lFV3/MqbeWVDwqRKq/3N5idjVNkjIll1qMZJygkJSXRokULIPicyaEoKTIq22P5MGs2aoyJVZ71U4gn69ato02bNiQmJpYuGzKkamX0JUVG5XMIvqzZqDEmVllOAVi/fn1pJXN1+Ssy8iUCVRiZ2xhjaoQFBSp2XKuOyoqGVGHKFKtsNsbEpjofFAoKCti0aVPEgkIoRUNW2WyMiVV1PiiU9FGIVFAIdf4Dq2w2xsSiOh8UyjdHra7yLZd86q7LsMpmY0wsqvOtjyIdFOBwq6XRo51mqSJOXUIJG6PIGBOr6nxOYf369SQmJpb2UYgE357M4ASEksncwu3zYIwxNclyCuvW0bp1a5KSIvdR+GuWquoEBDdjYowxManO5xQi2Ry1hE2AY4yJVxYUPAgKNgGOMSZe1cmgUDI2kUgBP/20id27MyN6fJsAxxgTr+pcUChbCfwjoMycmRnRHsbVHVDPGGOiRVQrzHgZN7KzszU3NzesfcqOXjoXOBOYR0ZGX6sENsbUCSKySFWz/a2rczmFspW969yfmQGHuTbGmLrEy+k4W4vIPBFZJSIrRORWd/kDIrLMnYlttoi08NlnlIisFZHVItLPi3SVrez9yf3Zwj1/2SkzjTGmrvEyp1AIjFDVE4DewDAR6QQ8UjIbGzATuAfAXTcY6AycDTwjIgEGiai6sWMPdySDPOAIIKV0vU2ZaYypy7ycjnOzqi523+8BVgEtVXW3z2a/AkoqNQYBr6lqgar+AKwFekU6XUOG+A45sRVoVmEbG8XUGFNX1UidgohkAicCC93fx4rIj8AQ3JwC0BKnOVCJje6y8scaKiK5IpKbl5dXpfQcnk8nD39BAayjmTGmbvI8KIhIQ+ANYHhJLkFVR6tqayAHuKlkUz+7V2gapaqTVDVbVbPT09OrlKbD/Qi2Av6PYR3NjDF1kadBQUSScQJCjqq+6WeTacAf3fcbgdY+61oBm7xIV0k/goQEJ6cg5cKRdTQzxtRVXrY+EmAysEpVx/ssb++z2UDgG/f928BgEUkVkbZAe+BLr9I3eHARsI277krnlVeso5kxxoC3o6SeDFwBLBeRJe6yvwHXiEhHoBhYD1wPoKorROR1YCVOy6VhqlrkVeJ27NhBcXExzZo1Y8gQCwLGGAMeBgVV/Rz/9QSzguwzFqiRgpuSSupmzfxXNBtjTF1U53o0l9i6dSsAVa2sNsaY2qjOBgXLKRhjTEV1NihYTsEYYyqqs0EhLy8PESEtLS3aSTHGmJhRZ4PC1q1bSUtLi+jczMYYE+/qdFCwoiNjjCmrzgaFvLw8q2Q2xphy6mxQsJyCMcZUVGeDguUUjDGmojoZFAoLC9m+fbsFBWOMKadOBoXt27cD1kfBGGPKq5NBoaTjmuUUjDGmrDoZFFJTU7noooto165dtJNijDExpU723OrQoQOvv/56tJNhjDExp07mFIwxxvhnQcEYY0wpCwrGGGNKeTlHc2sRmSciq0RkhYjc6i5/RES+EZFlIvIfETnCZ59RIrJWRFaLSD+v0maMMcY/L3MKhcAIVT0B6A0ME5FOwBygi6p2A74FRgG46wYDnYGzgWdEJNHD9BljjCnHs6CgqptVdbH7fg+wCmipqrNVtdDdbAHQyn0/CHhNVQtU9QdgLdDLq/QZY4ypqEbqFEQkEzgRWFhu1dXAe+77lsCPPus2usvKH2uoiOSKSG7JlJrGGGMiw/OgICINgTeA4aq622f5aJwippySRX521woLVCeparaqZtswFcYYE1medl4TkWScgJCjqm/6LL8KOBf4naqW3Pg3Aq19dm8FbAp2/EWLFm0TkfXVSGJTYFs19o83de16wa65rrBrDk9GoBVy+J4cWSIiwBRgh6oO91l+NjAeOF1V83yWdwam4dQjtADmAu1VtciTBDrnzFXVbK+OH2vq2vWCXXNdYdccOV7mFE4GrgCWi8gSd9nfgIlAKjDHiRssUNXrVXWFiLwOrMQpVhrmZUAwxhhTkWdBQVU/x389wawg+4wFxnqVJmOMMcHV9R7Nk6KdgBpW164X7JrrCrvmCPGsTsEYY0z8xq9OVQAABQRJREFUqes5BWOMMT4sKBhjjClV64OCiJztDrC3VkTu9LNeRGSiu36ZiPSIRjojKYRrHuJe6zIR+UJEukcjnZFU2TX7bHeSiBSJyIU1mT4vhHLNItJXRJa4g1J+UtNpjLQQvttNROQdEVnqXvOfo5HOSBGRF0Vkq4h8HWB95O9fqlprX0Ai8B1wLJACLAU6ldtmAM5QG4IzcN/CaKe7Bq75N8CR7vv+deGafbb7CKcF3IXRTncN/J2PwGni3cb9vVm0010D1/w34GH3fTqwA0iJdtqrcc2nAT2ArwOsj/j9q7bnFHoBa1X1e1U9CLyGM/Cer0HAy+pYABwhIsfUdEIjqNJrVtUvVPUX91ffQQnjVSh/Z4CbcXrYb63JxHkklGu+DHhTVTcAqGq8X3co16xAI7fzbEOcoFBInFLVT3GuIZCI379qe1AIZZC9kAbiiyPhXs81HB6UMF5Ves0i0hI4H3i2BtPlpVD+zh2AI0XkYxFZJCJX1ljqvBHKNT8FnIAzRM5y4FZVLa6Z5EVFxO9fno59FANCGWQvpIH44kjI1yMiv8UJCqd4miLvhXLNE4A7VLXI7Ukf70K55iSgJ/A7oD4wX0QWqOq3XifOI6Fccz9gCXAGcBzOyAmfqc9gnLVMxO9ftT0ohDLIXtgD8cW4kK5HRLoBLwD9VXV7DaXNK6FcczbwmhsQmgIDRKRQVd+qmSRGXKjf7W2qug/YJyKfAt1xJreKR6Fc85+Bh9QpcF8rIj8AxwNf1kwSa1zE71+1vfjoK6C9iLQVkRScmd3eLrfN28CVbi1+b2CXqm6u6YRGUKXXLCJtgDeBK+L4qdFXpdesqm1VNVNVM4EZwI1xHBAgtO/2f4FTRSRJRBoAv8aZ7CpehXLNG3ByRohIc6Aj8H2NprJmRfz+VatzCqpaKCI3AR/gtFx4UZ2B96531z+L0xJlAM5Mb/k4TxpxK8RrvgdIw5nyFKBQ43iEyRCvuVYJ5ZpVdZWIvA8sA4qBF1TVb9PGeBDi3/kB4CURWY5TtHKHqsbtkNoi8irQF2gqIhuBvwPJ4N39y4a5MMYYU6q2Fx8ZY4wJgwUFY4wxpSwoGGOMKWVBwRhjTCkLCsYYY0pZUDCmEm4b8M9FpL/Psovd5p7G1CrWJNWYEIhIF+DfwIk4beSXAGer6nfVOGaSqsbtYG2mdrKgYEyIRGQcsA/4FbBHVR8QkauAYThDOX8B3KSqxSIyCWfI4/rAdFW93z3GRuA54Gyc8ZhaAf8HHAKWq+rlNXxZxpRRq3s0GxNh9wGLgYNAtpt7OB/4jdvbdhLO0AvTgDtVdYeIJAHzRGSGqq50j7NPVU8GEJHNQIaqHhSRI2r8iowpx4KCMSFS1X0iMh3Yq6oFInImcBKQ6w4XUp/DwxhfKiLX4PyPtQA64Ux4AzDd57ArgKki8l8gnsdiMrWEBQVjwlPsvsAZW+dFVb3bdwMRaQ/cCvRS1Z0iMhWo57PJPp/3/YDTcSZLuUtEuqhqkWepN6YS1vrImKr7ELhYRJoCiEiaOwJtY2APsNudBaufv51FJBFopaofAX/BmT6yQY2k3JgALKdgTBWp6nIRuQ/4UEQScCqLrwdycYqKvsYZtvl/AQ6RBEwTkUY4D2gPq+oe71NuTGDW+sgYY0wpKz4yxhhTyoKCMcaYUhYUjDHGlLKgYIwxppQFBWOMMaUsKBhjjCllQcEYY0yp/w90C4x/ZjImHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_N=30\n",
    "SSE_train_st=np.zeros(max_N)\n",
    "SSE_test_st=np.zeros(max_N)\n",
    "for i in range(1,max_N):\n",
    "    Ztrain_st=np.hstack((Ztrain_st,st_train_year.reshape(-1,1)**i))\n",
    "    Ztest_st=np.hstack((Ztest_st,st_test_year.reshape(-1,1)**i))\n",
    "    A_st = np.linalg.solve(Ztrain_st.T@Ztrain_st,Ztrain_st.T@st_train_price)\n",
    "    SSE_train_st[i]=np.sum((st_train_price-Ztrain_st@A_st)**2)/len(st_train_price)\n",
    "    SSE_test_st[i]=np.sum((st_test_price-Ztest_st@A_st)**2)/len(st_test_price)\n",
    "plt.plot(st_train_year,st_train_price,'o',color='blue',label='Steel')\n",
    "plt.plot(st_train_year,Ztrain_st@A_st,'-',color='black',label='Best Fit')\n",
    "plt.legend()\n",
    "plt.xlabel('Years')\n",
    "plt.ylabel('Prices $')\n",
    "plt.title('Price vs Year Data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEWCAYAAAAgpUMxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzde5xVZd3//9d7zidmhpOiIIGKJSgijFhpamkeylJTE0vF8g4179Jf9Sg1uzXUvtLdnaWWxV14ykRulaJSEc+aykE8IKCC51GU8wAjw5w+vz/WtYfNZs/MnsOezcz+PB+P9dhrXWtd17rWbJjPXGtd67pkZjjnnHM9LSfTFXDOOZedPAA555zLCA9AzjnnMsIDkHPOuYzwAOSccy4jPAA555zLCA9AzqVA0q2SrklT2UslHZWOsp3blXkAci6OpMclbZBU2FPnNLMxZvZ4T52vNSHI1kvaEre8lOl6ub7LA5BzgaQRwOcAA76a0cpkzi/NrCxuOSjZQZLyUklrS0ePd32PByDntjsHeA64FZjc2kGSzpX0dEKaSdo3rN8q6feSHgitiH9LGiLpN6F19aqkg+Pyvi3pmLB+laRZkm6XtDncnqtKdp64c10T1o+SVC3px5JWS1ol6WRJX5L0uqT1ki7vzA9G0ohw7vMkvQs8miwtHPvVUO+NoUW5f8K1/kTSy0CtB6Hs5gHIue3OAe4My3GSdu9CWV8HrgAGAduAZ4HFYfse4Ndt5P0qMBOoBOYAN3XgvEOAImAo8F/A/wJnAROIWnf/JWnvjlxIgiOB/YHjkqVJ2g+4C7gEGAzcD/xDUkHc8WcCXwYqzayxC3VxvZwHIOcASYcDnwBmmdnzwBvAN7pQ5Gwze97M6oDZQJ2Z3W5mTcDdwMFt5H3azO4Px94BJL0N1ooG4FozayAKYoOA35rZZjNbCiwFxraR/0eh5RJbbkvYf5WZ1ZrZ1lbSzgD+ZWbzQh1+BRQDn407/gYzey+hDJeFPAA5F5kMPGRma8P2X2njNlwKPopb35pku6yNvB/GrX8MFHXgVtW6ELhi50lWl7bO/Sszq4xbEn8G7yXJE5+2J/BObMPMmsP+oe2U4bKQ3391WU9SMdEts1xJsV/+hUClpIPMLLEnWC1QEpd/SM/UFIgCUknc9hCgugfPn2z4/Pi0D4ADYxuSBOwFvN9OGS4LeQvIOTgZaAJGA+PCsj/wFNFzoUQvAWMkjZNUBFzVQ/UEeBH4hqRcSccTPX/ZlcwCvizpaEn5wA+JnoE9k9lquV2RByDnolttt5jZu2b2YWwhevj/zcTbX2b2OjAVeBhYATy9U4npczHwFWAj8E3gb91c/o8T3gNa236W7czsNaJODzcCa0Ndv2Jm9d1cT9cHyCekc845lwneAnLOOZcRHoCcc85lhAcg55xzGeEByDnnXEb4e0ApGDRokI0YMSLT1XDOuV7l+eefX2tmg1vb7wEoBSNGjGDRokWZroZzzvUqkt5pa7/fgnPOOZcRHoCcc85lhAcg55xzGeHPgJxzu7yGhgaqq6upq6vLdFVcEkVFRQwbNoz8/PwO5fMA5Jzb5VVXV9OvXz9GjBhBNMC221WYGevWraO6upqRI0d2KK/fgnPO7fLq6uoYOHCgB59dkCQGDhzYqdapByDnXK/gwWfX1dnvJm0BSNIMSaslvZKQ/j1Jr0laKumXcemXSVoZ9h0Xlz5B0pKw74YwwRWSCiXdHdLnSxoRl2eypBVhmRyXPjIcuyLkjZ+nvvu9dDcsmpHWUzjnXG+VzhbQrcDx8QmSPg+cBIw1szFE88UjaTQwCRgT8vxeUm7IdjMwBRgVlliZ5wEbzGxf4HpgWihrAHAlcCgwEbhSUv+QZxpwvZmNAjaEMtJn6WwPQM71AevWrWPcuHGMGzeOIUOGMHTo0Jbt+vrUpjr61re+xWuvvdbmMb/73e+48847u6PKvULaOiGY2ZPxrZLgQuA6M9sWjlkd0k8CZob0tyStBCZKehsoN7NnASTdTjR75QMhz1Uh/z3ATaF1dBwwz8zWhzzzgOMlzQS+AHwj5Lkt5L+5+646QXF/+OiV9o9zzu3SBg4cyIsvvgjAVVddRVlZGT/60Y92OMbMMDNycpL/XX/LLbe0e56LLrqo65XtgMbGRvLy8lrdTjVfZ/X0M6D9gM+F22BPSDokpA8F3os7rjqkDWXH+e5j6TvkMbNGoAYY2EZZA4GN4djEsnYiaYqkRZIWrVmzpsMXCkBxJWzd0Lm8zrld3sqVKznggAO44IILGD9+PKtWrWLKlClUVVUxZswYpk6d2nLs4YcfzosvvkhjYyOVlZVceumlHHTQQXzmM59h9erob/ErrriC3/zmNy3HX3rppUycOJFPfvKTPPNMNKt5bW0tp556KgcddBBnnnkmVVVVLcEx3sKFCznyyCOZMGECJ5xwAh999FFLuT/96U854ogjuOmmmzjrrLP44Q9/yOc//3kuv/xy1q5dy1e/+lXGjh3LZz/7WV555ZWWup1//vl88Ytf5Fvf+la3/Px6uht2HtAf+DRwCDBL0t5AsidY1kY6ncjTVlk77zCbDkwHqKqq6ty0scX9oX4LNDVAbsf6xzvnkvv5P5ay7INN3Vrm6D3LufIrYzqVd9myZdxyyy384Q9/AOC6665jwIABNDY28vnPf57TTjuN0aNH75CnpqaGI488kuuuu44f/OAHzJgxg0svvXSnss2MBQsWMGfOHKZOncqDDz7IjTfeyJAhQ7j33nt56aWXGD9+/E75tm3bxsUXX8ycOXMYNGgQd955Jz/72c+YPn06AJs2beLJJ58E4KyzzuKNN97gkUceIScnhwsvvJBDDz2UOXPm8NBDD3Huuee2jIX5wgsv8OSTT1JUVNSpn1Winm4BVQP3WWQB0AwMCul7xR03DPggpA9Lkk58Hkl5QAWwvo2y1gKV4djEstKjODx62roxradxzmXOPvvswyGHHNKyfddddzF+/HjGjx/P8uXLWbZs2U55iouLOeGEEwCYMGECb7/9dtKyv/a1r+10zNNPP82kSZMAOOiggxgzZufAuXz5cpYuXcoxxxzDuHHjuO6663jvve03hmL5Y04//fSWW4dPP/00Z599NgDHHnssH3zwAbW1tQCcdNJJ3RZ8oOdbQH8jeg7zuKT9gAKiwDAH+KukXwN7EnU2WGBmTZI2S/o0MB84B7gxlDUHmAw8C5wGPGpmJmku8Iu4jgfHApeFfY+FY2eGvH9P69UWVUafdRuhrNURyZ1zHdDZlkq6lJaWtqyvWLGC3/72tyxYsIDKykrOOuuspO/HFBRs74Cbm5tLY2PjTscAFBYW7nSMWfs3ZMyMsWPH8tRTT7Vb58TtxPLjtxPzdVU6u2HfRRQcPimpWtJ5wAxg79A1eyYwObSGlgKzgGXAg8BFZtYUiroQ+BOwEniDqAMCwJ+BgaHDwg+ASwFC54OrgYVhmRrrkAD8BPhByDMwlJE2f1wYTuvPgZzLCps2baJfv36Ul5ezatUq5s6d2+3nOPzww5k1axYAS5YsSdrCGj16NO+//z4LFiwAoL6+nqVLl6ZU/hFHHNHSE+/hhx9m2LBh3R54YtLZC+7MVnad1crx1wLXJklfBByQJL0OOL2VsmYQBbvE9DeJumb3iBrCl+YByLmsMH78eEaPHs0BBxzA3nvvzWGHHdbt5/je977HOeecw9ixYxk/fjwHHHAAFRUVOxxTWFjIPffcw/e//302b95MY2MjP/zhD5Perks0depUvvWtbzF27FjKyspS6r3XWUqlOZftqqqqrDMT0v3qrvv50WtnwinT4aAz0lAz57LD8uXL2X///TNdjV1CY2MjjY2NFBUVsWLFCo499lhWrFjRLd2iuyLZdyTpeTOrai2PD0aaRvllAwGwreuTdsFzzrmO2rJlC0cffTSNjY2YGX/84x8zHnw6q3fWupco7hf1g2jYsp70jvnjnMsWlZWVPP/885muRrfwwUjTqKK0iE1WQv2W9e0f7JxzWcYDUBpVFOez0Upp9ADknHM78QCURhXFBdRQSrP3gnPOuZ14AEqjypJ8NloZ8gDknHM78QCURpUl+dRQSu42H4rHud6sO6ZjAJgxYwYffvhhy3YqUzT0Zd4LLo0qiwuosTLy6ldkuirOuS5IZTqGVMyYMYPx48czZMgQILUpGrpTU1MTubm5Lds9Pf1CIm8BpVFRfg6bVUZhYw34C7/O9Um33XYbEydOZNy4cXz3u9+lubmZxsZGzj77bA488EAOOOAAbrjhBu6++25efPFFzjjjjJaWUypTNKxYsYJDDz2UiRMn8rOf/YzKysoO1aOyspIrrriCiRMnsmDBAoYNG8bVV1/NYYcdxuzZs1m8eDGHHnooY8eO5dRTT6WmpgbYedqGdPAWUBpJYlt+ObnNTdG0DIX9Ml0l53q/By6FD5d0b5lDDoQTrutwtldeeYXZs2fzzDPPkJeXx5QpU5g5cyb77LMPa9euZcmSqJ4bN26ksrKSG2+8kZtuuolx48btVFZrUzR873vf40c/+hGnn356q4GgtXp8/etfp6amhvHjx3PNNde0HF9aWsq///1vIBo3bvr06Rx++OFcfvnlXH311fzqV78Cdpy2IR28BZRmjQVhjCafksG5Pufhhx9m4cKFVFVVMW7cOJ544gneeOMN9t13X1577TUuvvhi5s6du9NYbcm0NkXD/PnzOfXUUwH4xje+kTRva/WAaOTtU045ZYfjzzgjGhps3bp11NXVcfjhhwMwefLkHQJO4rQN3c1bQGnWXFQJdUQDklbu1e7xzrl2dKKlki5mxre//W2uvvrqnfa9/PLLPPDAA9xwww3ce++9LZPBtSbVKRo6Uo/GxkaKi4uRdhwMLDa6dXtjgaZrFOwYbwGlW1GYlqjOW0DO9TXHHHMMs2bNYu3atUDUonj33XdZs2YNZsbpp5/Oz3/+cxYvXgxAv3792Lx5c4fOMXHiRGbPng3AzJkzO1SP9gwaNIji4uKW6b7vuOMOjjzyyA7Vryu8BZRmOaUDohV/F8i5PufAAw/kyiuv5JhjjqG5uZn8/Hz+8Ic/kJuby3nnnYeZIYlp06YBUbfr//iP/6C4uLhlrp723HDDDZx99tlMmzaNL33pS0lv57VWjz333LPd8u+44w4uvPBCtm7dyr777tujPfN8OoYUdHY6BoDf3vcYF798MnzltzDh3O6tmHNZIpunY6itraWkpARJ/OUvf2H27Nnce++9ma7WTnw6hl1QQZiSobF2g/+wnXMdtnDhQi655BKam5vp379/j787lE7pnJJ7hqTVYfrtxH0/kmSSBsWlXSZppaTXJB0Xlz5B0pKw7waFp2mSCiXdHdLnSxoRl2eypBVhmRyXPjIcuyLkTfssCWVl/dhmedRvXpfuUznn+qCjjjqKF198kZdffpknnniCvffeO9NV6jbp7IRwK3B8YqKkvYAvAu/GpY0GJgFjQp7fS4q9rnszMAUYFZZYmecBG8xsX+B6YFooawBwJXAo0fTbV0oKPQGYBlxvZqOADaGMtKooLWQTpTTUegByriv8ccGuq7PfTdoCkJk9CSSbh+B64MdAfI1PAmaa2TYzewtYCUyUtAdQbmbPWnSFtwMnx+W5LazfAxwdWkfHAfPMbL2ZbQDmAceHfV8IxxLyxspKm8riaEDSplrvhOBcZxUVFbFu3ToPQrsgM2PdunUUFRV1OG+PPpaQ9FXgfTN7KaFf+lDgubjt6pDWENYT02N53gMws0ZJNcDA+PSEPAOBjWbWmKSsZHWdQtTyYvjw4alfZIKK4nw2UspA7wXnXKcNGzaM6upq1qxZk+mquCSKiooYNmxYh/P1WACSVAL8FDg22e4kadZGemfytFXWzjvMpgPTIeoF19px7aksyWeFlaK6ms4W4VzWy8/PZ+TIkZmuhutmPfki6j7ASOAlSW8Dw4DFkoYQtUbihwkYBnwQ0oclSSc+j6Q8oILoll9rZa0FKsOxiWWlTWVxATWUkVfvL6I651y8HgtAZrbEzHYzsxFmNoIoUIw3sw+BOcCk0LNtJFFngwVmtgrYLOnT4RnOOcDfQ5FzgFgPt9OAR8NzornAsZL6h84HxwJzw77HwrGEvLGy0qZfUR41VkpBw6Z0n8o553qVdHbDvgt4FvikpGpJrfY4M7OlwCxgGfAgcJGZNYXdFwJ/IuqY8AbwQEj/MzBQ0krgB8Cloaz1wNXAwrBMDWkAPwF+EPIMDGWkVU6OqMsvp7CpFpoa0n0655zrNdL2DMjMzmxn/4iE7WuBa5Mctwg4IEl6HXB6K2XPAGYkSX+TqGt2j6rPr4B6oK4GSge1e7xzzmUDH4y0BzS3TMngPeGccy7GA1APaC4O78H6nEDOOdeizQAkKUfSoT1Vmb5KLQHIW0DOORfTZgAys2bgtz1Ulz4rr9QDkHPOJUrlFtw8SSelvSZ9WH4YEbvZA5BzzrVIpRfcfwIVkrYBW4lGFDAzG5DWmvUhhWVRC6h+8zo6PlqSc871TakEIO833EUVZSVssmLYst4DkHPOBe0GIDNrkvQl4IiQ9LiZPZjeavUtlcX51FgZZbXJBgd3zrns1O4zIEnXEk2f8GZYfizpmnRXrC+pLMmnhlKaP/ZnQM45F5PKLbivAAfHhsaRNANYDFyRzor1JZUl+XxopexZ5+8BOedcTKovopbHrfdLR0X6svLifDZSRu42D0DOOReTSgvol0TTJjxC1APuKOC/0lmpvqaiOJ9NVkp+vc8J5JxzMW0GoDAFwiNE0xgcShSA/svM3u+BuvUZhXm51Ob2o7BxE5iBks2N55xz2aXNAGRmJumfZjYBuK+H6tQn1eeVk9vUCPW1UFiW6eo451zGpfIMaIGk8WmvSR/XUFgZrXhHBOecA1ILQIcTBaHXJC2W9IKkxemuWF/THAtAPhyPc84BqXVCODnttcgGxZWwAQ9AzjkXtNcJIRe4z8wO6qH69Fm5JWHoPJ8TyDnngPanY2gClkka2tGCJc2QtFrSK3Fp/y3pVUkvS5otqTJu32WSVoZbfcfFpU+QtCTsuyH0zENSoaS7Q/p8SSPi8kyWtCIsk+PSR4ZjV4S8BR29rs7KLYsCkHkLyDnngNSeAQ0ClkuaK+m+2JJCvluB4xPS5gEHmNlY4HXgMgBJo4FJwJiQ5/eh9QVwMzAFGBWWWJnnARvMbF/gemBaKGsAcCVRt/GJwJWSwoQ8TAOuN7NRRDfEzkvhOrpFYZiSoXGLjwfnnHOQ2jOg6zpTsJk9Gd8qCWkPxW0+B5wW1k8CZprZNuAtSSuBiZLeBsrN7FkASbcTPZN6IOS5KuS/B7gptI6OA+aZ2fqQZx5wvKSZwBeAb4Q8t4X8N3fm+jqqrF859ZZL/ZZ15PfECZ1zbhfXagCSNMrMVpjZI5LyzKwxbt8h3XDubwN3h/WhRAEppjqkNYT1xPRYnvcAzKxRUg0wMD49Ic9AYGPcdcSXtRNJU4haXgwfPryDl7azipICaiglz1tAzjkHtH0L7u649QUJ+/7YlZNK+inQCNwZS0pymLWR3pk8bZW18w6z6WZWZWZVgwcPbu2wlMWmZPARsZ1zLtJWAFIr68m2UxY6BZwIfNPMYgGgGtgr7rBhwAchfViS9B3ySMoDKoD1bZS1FqgMxyaWlXYVJdGApN4JwTnnIm0FIGtlPdl2SiQdD/wE+KqZfRy3aw4wKfRsG0nU2WCBma0CNkv6dHi+cw7w97g8sR5upwGPhoA2FzhWUv/Q+eBYYG7Y9xjbnztNjisr7SpLCqixUnJ9JATnnAPa7oQwTNKviVo7sXXCdrvdsiXdRTRy9iBJ1UQ90y4DCoF5oTf1c2Z2gZktlTQLWEZ0a+6i2PxDwIVEPeqKiTofPBDS/wzcETosrCfqRYeZrZd0NbAwHDc11iGBKPjNDBPqvRDK6BGVxflspJTc+o966pTOObdL0/a7YAk7pDa7KJtZj/3yzrSqqipbtGhRl8owM279r0l8o+ApCn/WY3f+nHMuYyQ9b2ZVre1vtQWUTQGmJ0iiPr+cwqZaaGqE3FR6wDvnXN+V6oyorhvU51dEK3U+MZ1zznkA6kE+IrZzzm3nAagHWZHPCeScczHtBiBJ+4Zx4F4K22MlXZb+qvU9KglD0nkLyDnnUmoB/Qn4OdActpcAZ6WtRn1YbqkHIOeci0klAJWa2TOxjfBCZ0P6qtR3FZQNAqDpYx8PzjnnUglA68LoBAYg6WTgw7TWqo8q6hfNCbRt07oM18Q55zIvlZdR/pNoxIBPSXoHWEUYdcB1TEVZMZutmKYt6yjJdGWccy7DUpmS+yAz+4KkCqKRE7wLVydVFOdTQynFtf4MyDnnUpmS+5KwXuPBp2sqivOpsVIfEds550jtGdBcSZdI2kNSeWxJe836oMqSAjZaGfL3gJxzLqVnQOeHzx/GpRnQ9WlCs0xlcT7LKCV329pMV8U55zKu3QBkZnu1d4xLTXm4BZdfvzLTVXHOuYxLaUhmSZ8CRgNFsTQz+2u6KtVX5eaIrXnlFDZuBjNQpyeWdc65Xq/dACTpCqJZRT9FNNvoccDTgAegTqjPryCvoQEaPoaC0kxXxznnMiaVTghnAJ8HVpnZ2cBBpNhycjtrLAhTMnhPOOdclkslAG0N3bEbJfUjGgVh7/YySZohabWkV+LSBkiaJ2lF+Owft+8ySSslvSbpuLj0CZKWhH03KMzlLalQ0t0hfb6kEXF5JodzrJA0OS59ZDh2RchbkML1d6vm2IjYW70nnHMuu6USgF6QVAnMABYBC4DFKeS7FTg+Ie1S4BEzGwU8EraRNJpodIUxIc/vw0uwADcDU4BRYYmVeR6wwcz2Ba4HpoWyBgBXAocCE4Er4wLdNOD6cP4NoYwepWKfE8g55yCFAGRm55vZRjP7HfBl4HwzOyeFfE8CiaNungTcFtZvA06OS59pZtvM7C1gJTBR0h5AuZk9GwZBvT0hT6yse4CjQ+voOGCema03sw3APOD4sO8L4djE8/eYnJJoPDifE8g5l+1S6YTw2WRp8SNkd8DuZrYKwMxWSdotpA8Fnos7rjqkNYT1xPRYnvdCWY2SaoCB8ekJeQYCG82sMUlZO5E0hajlxfDh3ffKU37ZQACaP17vswE657JaKp0Jfha3XgRMAF4AjuzGeiTrj2xtpHcmT1tl7bzDbDowHaCqqqrV4zqqsCxqAdVvWbe9T7tzzmWhVG7BnRC3fB4YC7zfyfN9FG6rET5Xh/RqIP6F12HAByF9WJL0HfJIygMqiG75tVbWWqAyHJtYVo8p6VdBg+VSv9nnBHLOZbcO3wUys7eBAzp5vjlArFfaZODvcemTQs+2kUSdDRaE23WbJX06PMM5JyFPrKzTgEfDc6K5wLGS+ofOB8cCc8O+x8KxiefvMZUlBWyklMYtHoCcc9ktlWdA17P9VlUOcDCwNIV8dwFHAYMkVRP1TLsOmCXpPOBd4HQAM1sqaRawDGgELgpdvwEuJOpRVww8EBaI5ii6Q9JKopbPpFDWeklXAwvDcVPNLPbb/ifATEnXEN1G/HN719HdKksK2GSllH/sveCcc9lNUcOgjQOiYBHTCLxtZk+ktVa7mKqqKlu0aFG3lPX6R5vZ/LvPM3zIIAZ/98FuKdM553ZFkp43s6rW9qcyGGmPtxL6sorifKqtlBzvhu2cy3Kp3IJ7geS9xQSYmY3v9lr1YRXF+WykjPz61e0f7JxzfVgq3bAfBnKBO8L2N4HNwF/SVam+rCg/ly0qo6BhU6ar4pxzGZVKAPqsmR0Wt/2CpH+b2c/TVam+rj6vnKKmLdDcBDm57Wdwzrk+KJVu2GWSPh3bkHQoUJa+KvV9DYVhROy6msxWxDnnMiiVFtB/ALdIir24vxX4dvqq1Pc1FlRGP8WtGyA2NpxzzmWZVHrBLQQOkDQwbK9Le636uuJKqMFHxHbOZbVWb8FJ+pKk+FE4zwceknSfpE+kv2p9l0rC7BA+J5BzLou19Qzo/wHrACR9mei223eBh4A/pr9qfVd+abjt5i0g51wWaysAmZnVhvWvAX8ys/lm9gdg9/RXre+KTcnQsMXvZjrnsldbAShHUkkYBPRo4NG4fYXprVbfVlQetYC2+YCkzrks1lYnhBuJBuysAVaY2QIASQcBH/ZA3fqs8tIStlgR9Zu9BeScy16tBiAz+19Jc4luty2O27UW74bdJZXFBWykjEIfEds5l8Xa7IZtZu8STZsQn9bZyehcUFmST42VMvhjvwXnnMteHZ6QznVdRXEUgOQjYjvnslhb7wENb22f65rKkmhE7Lx6H4rHOZe92moBzQaQ9FAP1SVrlBXmsZlS8j0AOeeyWFvPgHIl/RTYX9L3E3ea2Q3pq1bfJomteeUUNW4CM5AyXSXnnOtxbbWAzgyfecDgJEunSfr/JC2V9IqkuyQVSRogaZ6kFeGzf9zxl0laKek1ScfFpU+QtCTsuyG8s4SkQkl3h/T5kkbE5ZkczrFC0uSuXEdX1OdXkGcN0LA1U1VwzrmMaqsb9nLgWkkvm9k/uuuEkoYC3wdGm9lWSbOAScBo4BEzu07SpcClwE8kjQ77xwB7Ag9L2s/MmoCbgSnAc8D9wPHAA8B5wAYz21fSJGAacIakAcCVQBXRLK/PS5pjZj3eH7qxsALqiYbjKSjp6dM751zGpdIL7glJv5T0XFimSerXxfPmAcWS8oAS4APgJOC2sP824OSwfhIw08y2mdlbwEpgoqQ9gHIze9bMDLg9IU+srHuAo0Pr6DhgnpmtD0FnHlHQ6nFWWBmteE8451yWSiUA/RloAM4JSz1wS2dPGN4j+hXR+0WrgBozewjY3cxWhWNWAbuFLEOB9+KKqA5pQ8N6YvoOecyskWg0h4FtlLUTSVMkLZK0aM2aNZ272LaU+ICkzrnslkoAGmVmPzWz18PyM2Dfzp4wPNs5CRhJdEutVNJZbWVJkmZtpHc2z46JZtPNrMrMqgYP7tIjr6RyS0ILyAOQcy5LpRKA6iR9JrYRpueu68I5jwHeMrM1ZtYA3Ad8Fvgo3FYjfK4Ox1cDe8XlH0Z0y646rCem75An3OarANa3UVaPyyuNRsRu9uF4nEgLrGAAAB84SURBVHNZKpUA9F3gT6FH2Urgf4ELunDOd4FPJ4y0vRyYA8R6pU0G/h7W5wCTQs+2kcAoYEG4TbdZ0qdDOeck5ImVdRrwaHhONBc4VlL/0BI7NqT1uNiI2HWbfEBS51x2SmVK7sXAmNCDTF2dktvM5ku6h2iA00aiEbenA2XALEnnEQWp08PxS0NPuWXh+ItCDziAC4FbgWKi3m8PhPQ/A3eEgLmeqBcdZrZe0tXAwnDcVDPLyIBspf0qabQc6resw/vAOeeykaKGgWtLVVWVLVq0qFvLfOzV1Rx41wRs/68weNLvu7Vs55zbFUh63syqWtvvg5FmSHkYkNSfATnnspUHoAypLMmnhlK01d8Dcs5lp3afAUnKIXpZc0T88T4WXNdUFufzjpUxfJsHIOdcdmo3ABH1LDNgCdCc3upkj4riqAWUvy0NL7k651wvkEoAGmFmB6a9JlkmLzeHj3P6UdC4KdNVcc65jEjlGdBcSV9Ie02y0Lb8coqbNkNzU/sHO+dcH5NKC+gp4B+SjGgcOAFmZgPSWrMs0FBQEb3ZVFezfWw455zLEqm0gK4HPgf0J5oHaBBdnA/IRZoLfTw451z2SqUFtAJ4wfyN1e5XFObc867YzrkslEoA+gB4VNL9wLZYonfD7jqVhgBU5y0g51z2SSUAVYelPM11yTp5pdFzH/t4Q9J5Ipxzri9LZTDSn/VERbJRQVkUgOq3rKMww3VxzrmelspICPNIMmmbmR2blhplkeLyaE6gus0egJxz2SeVW3BXxK0XAacS9yzIdV6/0lJqrZCGLRmZEcI55zIqlVtw8xOSnpD0RJrqk1UqS/LZSBn5tR6AnHPZJ5VbcPGdD3KACcAeaatRFqksyafGyhjo3bCdc1kolVtwS4meAYnovf23gO+ks1LZorK4gLeslEF1HoCcc9mn3ZEQzGwvMxsePkea2RfMrEu34CRVSrpH0quSlkv6jKQBkuZJWhE++8cdf5mklZJek3RcXPoESUvCvhskKaQXSro7pM+XNCIuz+RwjhWSJnflOrqqojifjZSS51MyOOeyUKsBKPxy3z1u+5uS7pX0a0mVXTzvb4EHzexTwEHAcuBS4BEzGwU8EraRNBqYBIwhmpfo95JyQzk3A1OAUWE5PqSfB2wws32JhhKaFsoaAFwJHApMBK6MD3Q9rSg/h83qR36Dj4jtnMs+bbWAphPdckPS4cCvgFlAXdjXKeGZ0hHAnwHMrN7MNgInAbeFw24DTg7rJwEzzWybmb0FrAQmStoDKDezZ8MwQbcn5ImVdQ9wdGgdHQfMM7P1ZrYBmMf2oNXjJLEtrx9FPiWDcy4LtRWA8sxsXVifBEw3s7vN7HLgk104597AGuAWSS9I+pOkUmB3M1sFED53C8cPBd6Ly18d0oaG9cT0HfKYWSNQAwxso6ydSJoiaZGkRWvWpG/SuPr8CvKtHhq2pu0czjm3K2orAOXG3eo6Gng0fl8XzpkHjAduNrODgVrC7bZWJBulxtpI72yeHRPNpptZlZlVDR6cvsG/mworohUfEds5l2XaCkCzgMck3Us0D9BTAJL2Abpyz6gaqI57v+geooD0UbitRvhcHXf8XnH5hxENkFod1hPTd8gjKQ+oANa3UVbGNLWMiO0ByDmXXVoNQGY2FbgcmAkcbmbNYVc+8P3OntDMPgTekxS7jXc0sAyYA8R6pU0G/h7W5wCTQs+2kUSdDRaE23SbJX06PN85JyFPrKzTgEfDc6K5wLGS+ofOB8eGtIxRsU/J4JzLTm2+B2RmTydJe7Ubzvs94E5JBcCbwLeIguEsSecB7wKnh/MtlTSLKEg1AheZWWwO6wuBW4Fi4IGwQNTB4Q5JK4laPpNCWeslXQ0sDMdNNbOMDkOQV+ItIOdcdkrlRdRuZ2YvAlVJdh3dyvHXAtcmSV8EHJAkvY4QwJLsmwHM6Eh90ym/XzQidkPtOvIzXBfnnOtJqUzJ7dKosHwQAHWbfTw451x28QCUYaVllTRaDg2b17V/sHPO9SGt3oKTtIHkXZQFmJkNSFutskhFSQE1lNJU68+AnHPZpa1nQIN6rBZZLBoRu5SSj/0WnHMuu7QagOJ6mgEt46gVxSVl9P2ZvqKyuIB1lFHq3bCdc1mm3WdAkr4s6XWilzjnh89H287lUlVRks9GKyXHR8R2zmWZVDohXAscBrxmZnsRDej5eDorlU36FeaxiVLy62syXRXnnOtRqQSgRjNbA+RIkpnNIxo6x3WDnBzxcW45hT4itnMuy6TyImpNGK36aeB2SauB5nbyuA7Yll9BYcMWaG6CnK6M8+qcc71HKi2gk4nmALqE6Nbb+8CJaaxT1mksqCAHgzq/Deecyx6pBKDLzKzJzBrM7M9m9mvgB+muWDZpLgwTzNZ5RwTnXPZIJQAlmzH0y91dkaxWHAKQD0jqnMsibY2EcD5wAbCfpMVxu/oBi9JdsWyS4yNiO+eyUFudEGYBjwD/jx1nLN1sZquTZ3GdkVcajWrU/PFGH5zPOZc12pqQboOZrTSz04nm2/liWNI3P3WWKug3EIC6zWszXBPnnOs5qYyEcBFRa2h4WGZJ+m66K5ZNikMAqvcRsZ1zWSSV94DOByaa2RYASb8AngF+n86KZZPyfqXUWiH1W3xAUudc9kglAAloiNtuCGmum1SW5FNDKZXvPAyzLwBrBrPwGRZsezrAZ/4TPvGZjNbbOee6oq1ecHlm1gjcATwn6d6w6xTgtq6eWFIuUW+6983sxDDa9t3ACOBt4OtmtiEcexlwHtAEfN/M5ob0CcCtRM+o7gcuNjOTVAjcDkwA1gFnmNnbIc9k4IpQjWvMrMvX0lUVxfk80jSeUxtfgbf/DRIoJ+4zLITtDW9Dw8dw9uxMV9055zqtrRbQAmC8mf1S0mPA54haPheY2cJuOPfFwHKgPGxfCjxiZtdJujRs/0TSaGASMAbYE3hY0n5huoibgSnAc0QB6HjgAaJgtcHM9pU0CZgGnBGC3JVAFdFke89LmhMLdJlSUVzAzxq/jR0xhnM+M6L9DA9fBf++AT5eDyU+L6BzrndqqxNCy202M1toZr82s//pjuAjaRjRy6x/iks+ie0tq9uIhgCKpc80s21m9hawEpgoaQ+g3MyeNTMjavGcnKSse4CjJYloJO95ZrY+BJ15JH/RtkdVFOcDUPNxQztHBqNPBmuC5f9IY62ccy692moBDZbU6pA7YUiezvoN8GOil1pjdjezVaHsVZJ2C+lDiVo4MdUhrSGsJ6bH8rwXymqUVAMMjE9PkmcHkqYQta4YPnx4By+vYwrycigtyGXj1hQD0B4HQf+RsOxvMGFyWuvmnHPp0lYAygXK6OYOB5JOBFab2fOSjkolS5I0ayO9s3l2TDSbDkwHqKqqSnpMd6osKeDBVz7knXW17R4riauGHs/QpdOhdh2UDkx39Zxzrtu1FYBWmdnUNJzzMOCrkr5ENMV3uaS/AB9J2iO0fvYAYqMtVAN7xeUfRjQdeHVYT0yPz1MtKQ+oANaH9KMS8jzefZfWeaccPJTHXlvNqpq6do99b/3HXFO2HzdbE7z6D5hwbvor6Jxz3Uxmyf+4l/SCmR2c1pNHLaAfhV5w/w2si+uEMMDMfixpDPBXYCJRJ4RHgFFm1iRpIfA9oqnC7wduNLP7w8uzB5rZBaETwtfM7OuhE8LzbJ9QbzEwwczafAGnqqrKFi3adYa/+/uL73PxzBdYOuhySnfbG875W6ar5JxzO5H0vJlVtba/rU4IR6ehPm25DviipBVEQ/5cB2BmS4lGYlgGPAhcFHrAAVxI1JFhJfAGUQ84gD8DAyWtJJo64tJQ1nrgamBhWKa2F3x2RSeO3ZORg8r4Z+NE7K0no9twzjnXy7TaAnLb7WotIID/W/Qet9w7h/sLL4ev/NZvwznndjldaQG5XdjJBw9lU8WneD93T2ypv5DqnOt9PAD1Uvm5OXz386OYve0QeOtJqPWRtJ1zvYsHoF7s1AlDmV98BLJmzF9Kdc71Mh6AerHCvFyOOepo3mweQs2iWZmujnPOdYgHoF7ujInDeSzvMMo/fM5vwznnehUPQL1cUX4u/au+Tg7NvPP0zExXxznnUuYBqA84/uijeZs92bL4/zJdFeecS5kHoD6gpDCftZ84gU/VvcQrr6/MdHWccy4lHoD6iNHHTCZXxgtzb890VZxzLiUegPqIkmFjWV/8CfZe/TBLP6jJdHWcc65dHoD6ComScafx6Zxl3PbQgkzXxjnn2uUBqA8pGncquTIKVvyL1z/anOnqOOdcmzwA9SW7jaZpwChOzFvATY96ZwTn3K7NA1BfIpF7wClM1HKee3k5b67ZkukaOedcqzwA9TVjTiGHZr6Uv5DfPfZGpmvjnHOt8gDU1+y2Pwzaj8nlL/K3F9/n3XUfZ7pGzjmXlAegvkaCMacwYssL7K4abn7CnwU553ZNPR6AJO0l6TFJyyUtlXRxSB8gaZ6kFeGzf1yeyyStlPSapOPi0idIWhL23SBJIb1Q0t0hfb6kEXF5JodzrJA0ueeuvAeNOQVhXDryde55vpo757/DA0tW8e+Va3m5eiNvr61l3ZZt1Dc2Z7qmzrks1uNTckvaA9jDzBZL6gc8D5wMnAusN7PrJF0K9Dezn0gaDdwFTAT2BB4G9jOzJkkLgIuB54D7gRvM7AFJ3wXGmtkFkiYBp5jZGZIGAIuAKsDCuSeY2Ya26rwrTsndrpsmsq1wAJ/58Iesr61v9bCi/Bz6FeVTXpRH/5ICPjdqMCcetAf7DC7rwco65/qi9qbkzuvJygCY2SpgVVjfLGk5MBQ4CTgqHHYb8Djwk5A+08y2AW9JWglMlPQ2UG5mzwJIup0okD0Q8lwVyroHuCm0jo4D5pnZ+pBnHnA8UYDrW8acQuET03jm4jGspZJNWxvZXNfA5rpGNoXPzXUNbIr7/GDjVn7zyOtc//Dr7L9HOSeO3YOvjN2T4QNLMn01zrk+qMcDULxwa+xgYD6wewhOmNkqSbuFw4YStXBiqkNaQ1hPTI/leS+U1SipBhgYn54kT2LdpgBTAIYPH96p68uoMSfDE9dRtOJfDJv4HejffhaAD2vq+NeSVfzz5Q/477mv8d9zX2PssApOHLsHXx67J0Mri9Nbb+dc1shYAJJUBtwLXGJmm8Ljm6SHJkmzNtI7m2fHRLPpwHSIbsG1Vrld1m77w+BPwdLZMPE7KWcbUlHEeYeP5LzDR1K94WP+9fIq/vnyKn5x/6v84v5XGT+8khPH7smxY3ZnSHkRebnejyUZM+PVDzfz5OtreOHVN9j7w3/xxp5f4ZBP7c0R+w1m1G5ltPFv3rmskJEAJCmfKPjcaWb3heSPJO0RWj97AKtDejWwV1z2YcAHIX1YkvT4PNWS8oAKYH1IPyohz+PddFm7njGnwOPXwTvPQFElKCduUcJ2DuQVQenAluzD+pdw/pH7cP6R+/D22trQMlrF1H8uY+o/lwHQrzCPipJ8KkvyqSwuiNaLo+2K4iitrCgvaeRPRhL5uSI3R+Tn5pCXI/Jyc3ZKy8/NITdH5OWK/Jwc8nJFXsunMvLLfd2WbTy9ci1Pvr6Wp1asYfXmbUzUcn5f9DsG2XpWfXA/33nz+1zzr70ZUl7E50YN4oj9BnP4voPoX1rQ4/V1LtMy0QlBRM941pvZJXHp/w2si+uEMMDMfixpDPBXtndCeAQYFTohLAS+R3QL737gRjO7X9JFwIFxnRC+ZmZfD50QngfGh9MuJuqEsL6tOvfKTggAa16D303sWJ4Rn4OJU+CTX4Lc5H+frFy9hWfeWMuG2gY2bq2n5uMGNm5tYOPH9Wzc2tCy3dScuYZjXghOsaDUryiPQWWFDCwtZFBZQbReVsDAskIGlRYwqF8hA0sLqCwpIDcnteDV0NTM4nc28OSKNTz5+lpe+aAGM6gsyedz+/RnSs7fOeD136H+I+CIH8Oj12C1q3l+/x9zS90XePqNddRsbUCCsUMrOGK/wRyx32DG7VVJfhpblmZGs0FTs9FsRmOzRevNRpNFn0AU4HNyyA1BPTdH5ErkpPjzcalL/E6awzZAjiAn/EGVI5Gj6A+12OeurL1OCJkIQIcDTwFLgFg/4MuJgsgsYDjwLnB6XGeBnwLfBhqJbtk9ENKrgFuBYqLOB98zM5NUBNxB9HxpPTDJzN4Meb4dzgdwrZnd0l6de20AAnj3OdiyGqw5bjGwpoS05ui4xbdDzXtQsRccch6MnwwlAzp8WjOjtr6JjR/Xs2VbY8r5msIvw4Ymo7GpmcZmo6GpeXtaczONTXFpzeG4JqOhuZmmpri0Zms5dnNdA2u31LN2yzbW1dazvrY+aYDMESnfVozVNTdHHLxXZUsAObCijtzZU+CtJ+DA0+HE66GwH3y8Hu6bAivnwQGn0XTib3hpdSNPvb6WJ1es4YV3N9Bs0S/++CAYW4v/XaMkbUpLuJuc+F/bDJrMuvyHQY621zEvJ4eM/w7swOV05MpjvxutZTtWhu2wnai17ymWbsYOQaar34dCgGrta0j2/eRILd9h/B8WuQnpOYIZ5x7CJwaWdrJuu1gA6o16dQDqqKZGeP0BmP9HePup6LbcgafDoefDkAMzXbtu09xsbNzawLot21i7pZ51tdtYuzkKTg1Nqf2fyBGMHVbJZ/cdSHlRfpT4xmNw33dg2xb40i/h4LN3/A3Q3AxP/w889gsYuC98/fboeR1Qs7WBZ99Yy8vVNS1//bYElbgqxdfOzHb4K3in3zUJvwzzcuJ/0UBuTg65OTv+Qor9tR0LsE3NsVZSc0trqTG0mBqabKfAlwmt//pNcmwHAmZi8I/9rJV4QEwb31P0CTk5agkauaElkxN+7onrsTzNIWhB9G+32aJ/G80Waz0l/w6SJVsooym0eGPfcXPLOi0t4+Zm42cnjmZIRVEqP66deADqBlkVgOJ9tAwWTIeXZkLjVvjEYdHtuU+d2OrtuazV1AiP/z946n9g8CfhtFtg99GtH//Wk3DPt6G+Fk78DRx0Rs/V1bke4gGoG2RtAIrZugFe+EsUjDa+C+XDYNyZULZ71ELKK4L8ou3ridu5BZCTCyiu80NCJwjCdk4e5PSynnU178O9/wHvPgMHnwUn/BIKUrhlsWkV3HsevPNvmHAuHD8t+rk510d4AOoGWR+AYpqb4PW5sOCP8Obj6TtPfmn0C7ywLPosKAtLWI+lo6gFUb8lfNa2vl1cCYP2i257DdoPBo2KPis/0bXW3OtzYfYF0LgtetbT0ZZMUyM8ejX8+zcwZCx8/TYYsHfn6+PcLsQDUDfwAJREfS00bI2Wxrpoaajbvh6/3bQtdHxojvsMCwnbTQ07BpFtW8J6XFCJpWFQ0C8EptK4YJW4XQK1a2HdSlj7OtSu2X4dOfkwYGRccBoF+cXhxnvT9s4aO6w3R+trX4dFM2D3A+H0W2HQvp3/eb72IMyeEt2gP+G67UGo5f9n/MOF+PUmaKyPfsaN26CpPlpi6/GfOblRd/yiimgpjq2Hz8LynVufzU1QVxO1gus2wtaN4XNDWK+B5tY6mbTy2p3Z9uuJrSd+tlVWy0OcNl7r2+H3WuLPLsn5jJ3TW/5dxnXYaU5Ia27a/u9YuTu/2pDY2s/JjX5eTQ07f0ct6w3bv08p+jeamx/lzcmP7hLkhs/49fhzkeROQ+wORE5u+EOuX8JSHpa4tKJyKKzo9F0JD0DdwAPQLij277YzXbC2boC1K2HdiiiIrF0RLevfhOaGjpVVdR4c94vuuXW24W34v3Phgxe6XlainLy4oN8aRb9wiiqjn2/dRti2qe1ycwshr3Dn9KS/V8K74FL4jDtvS1r8Z7KyWgsw7QSoHf6dtHa+hLopZ8egkpMQYFq2c0M1mpMsCYGsuSn6LvIKolvTuYVhvTAKJHmF29dzw7thzQ0haDVGn80NO67H9u30R50l2bbo+Ppa2LY5fL/txIDvzofdPtX2Ma3Y5caCc65bdKXvb3F/2OuQaInX1Ag174a/PHPDLxjFrefu+EsorzC1Zz2p6j8Cvj0X3puf0Kpo6YKVJC3UI7dg+2f8euwXWk5O9Nd7/ZYosNTVbG/BJNtGoYVUGX0W99++Hv/pz6x6N7O4YBRbanbcLt8jbaf3AORcTG5e5p+/5BXCyCPSU3ZOTmjhlKenfNf7SNEz1cIyIH2BpjW9rLuRc865vsIDkHPOuYzwAOSccy4jPAA555zLCA9AzjnnMsIDkHPOuYzwAOSccy4jPAA555zLCB+KJwWS1gDvxCUNAtZmqDrp1levza+r9+mr15ZN1/UJMxvcWgYPQJ0gaVFb4xv1Zn312vy6ep++em1+Xdv5LTjnnHMZ4QHIOedcRngA6pzpma5AGvXVa/Pr6n366rX5dQX+DMg551xGeAvIOedcRngAcs45lxEegDpI0vGSXpO0UtKlma5Pd5H0tqQlkl6U1KvnH5c0Q9JqSa/EpQ2QNE/SivDZP5N17IxWrusqSe+H7+1FSV/KZB07Q9Jekh6TtFzSUkkXh/Re/Z21cV194TsrkrRA0kvh2n4e0jv0nfkzoA6QlAu8DnwRqAYWAmea2bKMVqwbSHobqDKzXv+CnKQjgC3A7WZ2QEj7JbDezK4Lfzj0N7OfZLKeHdXKdV0FbDGzX2Wybl0haQ9gDzNbLKkf8DxwMnAuvfg7a+O6vk7v/84ElJrZFkn5wNPAxcDX6MB35i2gjpkIrDSzN82sHpgJnJThOrkEZvYksD4h+STgtrB+G9Evgl6llevq9cxslZktDuubgeXAUHr5d9bGdfV6FtkSNvPDYnTwO/MA1DFDgffitqvpI/+giP7xPCTpeUlTMl2ZNNjdzFZB9IsB2C3D9elO/ynp5XCLrlfdpkokaQRwMDCfPvSdJVwX9IHvTFKupBeB1cA8M+vwd+YBqGOUJK2v3MM8zMzGAycAF4XbPW7XdzOwDzAOWAX8T2ar03mSyoB7gUvMbFOm69NdklxXn/jOzKzJzMYBw4CJkg7oaBkegDqmGtgrbnsY8EGG6tKtzOyD8LkamE10u7Ev+Sjck4/dm1+d4fp0CzP7KPwiaAb+l176vYXnCPcCd5rZfSG5139nya6rr3xnMWa2EXgcOJ4OfmcegDpmITBK0khJBcAkYE6G69RlkkrDQ1IklQLHAq+0navXmQNMDuuTgb9nsC7dJvafPTiFXvi9hQfafwaWm9mv43b16u+stevqI9/ZYEmVYb0YOAZ4lQ5+Z94LroNCl8nfALnADDO7NsNV6jJJexO1egDygL/25uuSdBdwFNHw8B8BVwJ/A2YBw4F3gdPNrFc90G/luo4iupVjwNvA+bF78L2FpMOBp4AlQHNIvpzoeUmv/c7auK4z6f3f2ViiTga5RA2ZWWY2VdJAOvCdeQByzjmXEX4LzjnnXEZ4AHLOOZcRHoCcc85lhAcg55xzGeEByDnnXEZ4AHIuCUlNYaTiVyT9n6SSdo7f0tb+dJFUJemGdo45StI/W9l3eBjV+NWwpDwMUxhBfVBH6+xcjAcg55LbambjwqjT9cAFma5QMma2yMy+35m8koYAfwUuMLNPAYcD50v6cpJj87pW05bR5J1r4QHIufY9BewLIOkHoVX0iqRLEg+UdIekk+K275T0VUnnSrpP0oNhrpRfxh1zpqK5mF6RNC0ufYukaWGA2IclTZT0uKQ3JX01HNPSugn7n5H0Qvj8ZDvXdRFwa9yIzWuBHwOXhvJulfRrSY8B0yQNlPRQKP+PxI2NKOms0JJ6UdIfY8EmXMNUSfOBz3Tkh+76Pg9AzrUh/OV/ArBE0gTgW8ChwKeB70g6OCHLn8IxSKoAPgvcH/aNA84ADgTOUDRh2Z7ANOALYf8hkmJD2JcCj5vZBGAzcA3RXFSnAFOTVPdV4AgzOxj4L+AX7VzeGKI5auItCukx+wHHmNkPiUZeeDqUP4fobXck7R+u67AwOGUT8M24a3jFzA41s6fbqY/LMl1uVjvXRxWHoeYhagH9GbgQmG1mtQCS7gM+B7wQy2RmT0j6naTdiCbnutfMGqNhwXjEzGpC3mXAJ4CBREFmTUi/EziCaOigeuDBUPQSYJuZNUhaAoxIUucK4DZJo4iGeclv5xpF8tHc49P+z8yawvoR4Zows39J2hDSjwYmAAvDdRazfRDKJqLBOJ3biQcg55LbGv6abxEGl0zFHUQtgEnAt+PSt8WtNxH9/2urzAbbPlZWcyy/mTW38kzmauAxMztF0fwzj7dTz6VAFTsOqDsBiJ/htzYhT7KAJeA2M7ssyb66uADm3A78FpxzqXsSOFlSSRg1/BSi1lGiW4FLAMxsaTtlzgeOlDQoPDc5E3iik/WrAN4P6+emcPzvgHMljQMIA0lOA37ZyvFPEm6tSToBiE2k9ghwWmj1IWmApE905gJcdvEWkHMpMrPFkm4FFoSkP5nZC0mO+0jScqLbaO2VuUrSZcBjRC2J+82ss9MO/JLoFtwPgEdTPPdZwP+G6TgE/MbM/tFKlp8Dd0laTBQk3w3lLJN0BdGMujlAA1EHh3c6eR0uS/ho2M51s/DO0BJgfOyZj3NuZ34LzrluJCk2MdeNHnyca5u3gJxzzmWEt4Ccc85lhAcg55xzGeEByDnnXEZ4AHLOOZcRHoDc/79RMApGwSgYEAAA/QHwB+CDTxEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "'''C'''\n",
    "max_N=30\n",
    "plt.plot(np.arange(1,max_N),SSE_train_al[1:],label='Training error')\n",
    "plt.plot(np.arange(1,max_N),SSE_test_al[1:],label='Testing error')\n",
    "plt.title('Aluminum Error')\n",
    "plt.ylabel('Total Sum of Square Error')\n",
    "plt.xlabel('Polynomial Order')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3hUZfbA8e9JIZV0QKogAkpXkKYUUUHFXtFFsazY1rarq+66Kv5W11XXta0FFXsBe0dARUBRutJ7DyWUkN5mzu+PexOGkDIkmUzK+TzPPDNz55Z3mDBn3nZeUVWMMcYYXyHBLoAxxpi6x4KDMcaYQ1hwMMYYcwgLDsYYYw5hwcEYY8whLDgYY4w5hAUHY2qIiGwUkVODXQ5jaoIFB9OgichJIvKziOwXkb0i8pOInOC+dpWIzK6lcrwuIgUikuVz+602rm1MVVhwMA2WiMQBXwLPAklAa2A8kB+kIj2mqrE+t15l7SQiYf5sq8jh7m9MaRYcTEPWGUBV31NVj6rmqupUVf1dRI4FXgQGur/i0wFEJEJEnhCRzSKyU0ReFJGo4hOKyFkislhE0t0aSc/qFlJE2ouIisi1IrIZ+L6sbe6+54jIMvf6M9z3UXyejSJyt4j8DmRbgDDVYcHBNGSrAY+IvCEiZ4hIYvELqroCuAGY4/6KT3Bf+jdOUOkNHI1T27gfQESOByYC1wPJwEvA5yISUUPlHQocC4wsa5uIdAbeA24HmgFfA1+ISBOf/S8DRgEJqlpUQ+UyjZAFB9NgqWoGcBKgwMtAmoh8LiItytpfRAS4DrhDVfeqaibwCDDa3eU64CVV/dWtibyB00Q1wM8i3en+4i++vVHq9QdVNVtVc8vZdinwlapOU9VC4AkgChjks/8zqrql1DmMOWxW7TQNmltDuApARI4B3gaewvmFXVozIBpY4MQJAAQIdR8fCYwVkVt8jmkCtPKzOE+o6n0VvL6lkm2tgE3FT1TVKyJbcGo3FZ3DmMNmwcE0Gqq6UkRex2kWAqdG4Ws3kAt0U9VtZZxiC/Cwqj4cqCJWsi0V6FH8xK3ptAW2lbO/MVVmzUqmwRKRY0TkLyLSxn3eFqfG8Iu7y06gTXGbvap6cZqf/isizd1jWotIcR/Ay8ANItJfHDEiMkpEmtbSW5oMjBKRU0QkHPgLTrPWz7V0fdOIWHAwDVkm0B/4VUSycYLCUpwvVXBGAC0DdojIbnfb3cBa4BcRyQCmA10AVHU+Tr/Dc8A+d7+rDqM8fy01z2F35YccoKqrgDE4Q3N3A2cDZ6tqweGcxxh/iC32Y4wxpjSrORhjjDmEBQdjjDGHCFhwEJGJIrJLRJb6bOstIr+4M0zni0g/n9fuFZG1IrLKpwPQGGNMEASy5vA6cHqpbY8B41W1N86s08cARKQrzkSjbu4xz4tIKMYYY4IiYPMcVHWmiLQvvRmIcx/H44zbBjgXeF9V84ENIrIW6AfMqegaKSkp2r596UsYY4ypyIIFC3ararOK9qntSXC3A9+KyBM4tZbiaf+tOTD2HGArB8/6LFP79u2ZP39+jRfSGGMaMhHZVNk+td0hfSNO3pq2wB3Aq+52KWPfMsfYisg4t79iflpaWoCKaYwxjVttB4exwMfu4w9wmo7AqSm09dmvDQeanA6iqhNUta+q9m3WrMJakTHGmCqq7eCQipOCGGA4sMZ9/Dkw2s2l3wHoBMyt5bIZY4xxBazPQUTeA4YBKSKyFXgAJ/XA0+4iJHnAOABVXSYik4HlQBFws6p6AlU2Y0zNKiwsZOvWreTl5QW7KMZHZGQkbdq0ITw8/LCPrdfpM/r27avWIW1M8G3YsIGmTZuSnJyMT7pzE0Sqyp49e8jMzKRDhw4HvSYiC1S1b0XH2wxpY0y15eXlWWCoY0SE5OTkKtfmLDgYY2qEBYa6pzqfSaMMDqnpuTw5dRUbdmcHuyjGGFMnNcrgsDe7gGe+X8vqnZnBLooxpgakp6fz/PPPV+nYM888k/T09Ar3uf/++5k+fXqVzl9a+/bt6dGjB71796Z3797ceuutNXLemtYolwmNj3J67jNyC4NcEmNMTSgODjfddNMhr3k8HkJDy0/V9vXXX1d6/oceeqha5Svthx9+ICUlpdzXi4qKCAsLK/e5v8dVR6MMDnFucNhvwcGYBuGee+5h3bp19O7dm9NOO41Ro0Yxfvx4WrZsyeLFi1m+fDnnnXceW7ZsIS8vj9tuu41x48YBB9LwZGVlccYZZ3DSSSfx888/07p1az777DOioqK46qqrOOuss7joooto3749Y8eO5YsvvqCwsJAPPviAY445hrS0NC6//HL27NnDCSecwJQpU1iwYEGFQcDXsGHDGDRoED/99BPnnHMOS5YsISkpiUWLFnH88cfz97//nWuuuYb169cTHR3NhAkT6NmzJw8++CCpqals3LiRlJQU3n333Rr5N22UwaFpRBgiVnMwJhDGf7GM5akZNXrOrq3ieODsbuW+/uijj7J06VIWL14MwIwZM5g7dy5Lly4tGcY5ceJEkpKSyM3N5YQTTuDCCy8kOTn5oPOsWbOG9957j5dffplLLrmEjz76iDFjxhxyvZSUFBYuXMjzzz/PE088wSuvvML48eMZPnw49957L1OmTGHChAnllvfkk08uqc2MHTuWO+64A3BqQD/++CMAV111FatXr2b69OmEhoZyyy23cNxxx/Hpp5/y/fffc+WVV5a83wULFjB79myioqL8/SetVKMMDiEhQtOIMKs5GNOA9evX76Dx/c888wyffPIJAFu2bGHNmjWHBIcOHTrQu3dvAPr06cPGjRvLPPcFF1xQss/HHzsZgWbPnl1y/tNPP53ExMRyy1Zes9Kll1560POLL764JIjMnj2bjz76CIDhw4ezZ88e9u/fD8A555xTo4EBGmlwAIiPDrfgYEwAVPQLvzbFxMSUPJ4xYwbTp09nzpw5REdHM2zYsDLH/0dERJQ8Dg0NJTc3t8xzF+8XGhpKUVER4Ew6q8kyl35e1vmLh6qWPq4mNMrRSuB0SmfkFQW7GMaYGtC0aVMyM8sffbh//34SExOJjo5m5cqV/PLLL+XuW1UnnXQSkydPBmDq1Kns27evRs8/ZMgQ3nnnHcAJdikpKcTFxVVyVNU13ppDlNUcjGkokpOTOfHEE+nevTtnnHEGo0aNOuj1008/nRdffJGePXvSpUsXBgwYUONleOCBB7jsssuYNGkSQ4cOpWXLljRt2rTMfX37HHr27Mmbb75Z6fkffPBBrr76anr27El0dDRvvPFGjZa/tEabW+nGtxewZlcW0/88tPKdjTEVWrFiBccee2ywixFU+fn5hIaGEhYWxpw5c7jxxhtLOoyDqazPxp/cSlZzMMaYGrB582YuueQSvF4vTZo04eWXXw52karFgoMxxtSATp06sWjRomAXo8Y02g7puKhwCoq85BXashHGGFNaow0OlkLDGGPK12iDg6XQMMaY8jXa4BBvwcEYY8oVsOAgIhNFZJeILC21/RYRWSUiy0TkMZ/t94rIWve1kYEqVzELDsY0HNVJ2Q3w1FNPkZOTU/LcnzTe/ti4cSNRUVEl6bl79+7t15yGuiCQo5VeB54DSv4lRORk4Fygp6rmi0hzd3tXYDTQDWgFTBeRzqoasN7ikj6HPAsOxtR3FaXs9sdTTz3FmDFjiI6OBvxL4+2vjh07VjrfoXRa8crSjBeryRTdpQWs5qCqM4G9pTbfCDyqqvnuPrvc7ecC76tqvqpuANYC/QJVNvCpOeRYcDCmvvNN2X3XXXcB8Pjjj3PCCSfQs2dPHnjgAQCys7MZNWoUvXr1onv37kyaNIlnnnmG1NRUTj75ZE4++WTASeO9e/duNm7cyLHHHst1111Ht27dGDFiREm+pXnz5tGzZ08GDhzIXXfdRffu3Q+rzLGxsdx///3079+fOXPm0L59ex566CFOOukkPvjgAxYvXsyAAQPo2bMn559/fkk6jmHDhvG3v/2NoUOH8vTTT9fUP+EhanueQ2dgsIg8DOQBd6rqPKA14JvsZKu77RAiMg4YB9CuXbsqF6RppPPW9+dafiVjatQ398COJTV7ziN6wBmPlvty6ZTdU6dOZc2aNcydOxdV5ZxzzmHmzJmkpaXRqlUrvvrqK8DJuRQfH8+TTz5ZbqbU8tJ4X3311UyYMIFBgwZxzz33lFu24qBV7Nlnn2Xw4MFkZ2fTvXv3gxYSioyMZPbs2YCTVuPZZ59l6NCh3H///YwfP56nnnoKODi1d6DUdod0GJAIDADuAiaLk1awrFWwy8zroaoTVLWvqvZt1qxZlQsSHhpCTJNQ63MwpgGaOnUqU6dO5bjjjuP4449n5cqVrFmzhh49ejB9+nTuvvtuZs2aRXx8fKXnKiuNd3p6OpmZmQwaNAiAyy+/vNzji5uVim+DBw8GnIyuF1544UH7Fqfs3r9/P+np6Qwd6qT3GTt2LDNnzjxkv0Cq7ZrDVuBjdRI6zRURL5Dibm/rs18bIDXQhXEys1pwMKZGVfALv7aoKvfeey/XX3/9Ia8tWLCAr7/+mnvvvZcRI0Zw//33V3iustJ410ROusjIyEP6FfxNvR2IFN2l1XbN4VNgOICIdAaaALuBz4HRIhIhIh2ATsDcQBcmzlJoGNMglE7ZPXLkSCZOnEhWVhYA27ZtY9euXaSmphIdHc2YMWO48847WbhwYZnHVyYxMZGmTZuWpP5+//33a/DdQHx8PImJicyaNQuAt956q6QWUVsCVnMQkfeAYUCKiGwFHgAmAhPd4a0FwFi3FrFMRCYDy4Ei4OZAjlQqZvmVjGkYSqfsfvzxx1mxYgUDBw4EnM7ft99+m7Vr13LXXXcREhJCeHg4L7zwAgDjxo3jjDPOoGXLlvzwww9+XfPVV1/luuuuIyYmhmHDhpXbRFW6z+Gaa67h1ltvrfT8b7zxBjfccAM5OTkcddRRvPbaa36Vq6Y0zpTdGdth6Yfcs6IDizPjmHL7kJovnDGNSGNM2Z2VlUVsbCzgdIhv3749oKOHqqqqKbsb5wzprB0w9T66sNFqDsaYKvnqq6/o3bs33bt3Z9asWdx3333BLlKNapwpu6OSAEgOzbbEe8aYKrn00ktrZdRQsDTOmkO0ExySJIvsAg+FHm+QC2RM/Vefm6gbqup8Jo0zODSJhdAmxGsGYGm7jamuyMhI9uzZYwGiDlFV9uzZQ2RkZJWOb5zNSiIQlURTrzN0bX9uIcmxEZUcZIwpT5s2bdi6dStpaWnBLorxERkZSZs2bap0bOMMDgDRScR49gOWmdWY6goPD6dDhw7BLoapQY2zWQkgKonIIic4ZORZfiVjjPHVeINDdCJNCqzmYIwxZWm8wSEqibB8JwWuBQdjjDlY4w0O0UmE5O0D1EYrGWNMKY04OCQj3iKSw/Kt5mCMMaU03uDgzpJuG5lrNQdjjCmlwuAgIiEi0r+2ClOr3FnSrSPyrOZgjDGlVBgcVNUL1L00gzXBrTkcEZ5jwcEYY0rxp1lpmoicG/CS1Da35tAizIKDMcaU5s8M6T8B8SKSD+TirPesqpoU0JIFmltzSAnNYn+2BQdjjPHlT3BICXgpgiEqARCSxNJ2G2NMaZUGB1X1iMiZQPFyaTNUdUpgi1ULQkIhKoEEMsjML8LrVUJCJNilMsaYOqHSPgcReRj4K7Devf1VRP7px3ETRWSXu1506dfuFBEVkRSfbfeKyFoRWSUiIw/vbVRRVBJNvRmoQqblVzLGmBL+dEifDZyiqhNUdQIwAjjHj+NeB04vvVFE2gKnAZt9tnUFRgPd3GOeF5FQP65RPdFJxHqdNR2sU9oYYw7wdxJcnM/jpv4coKozgb1lvPRfnJqI76og5wLvq2q+qm4A1gL9/Cxb1UUlEVVkwcEYY0rzp0P6MWChiHyHM1JpGHB/VS4mIucA21T1N5GD2vdbA7/4PN/qbivrHOOAcQDt2rWrSjEOiE4ionAJABl5FhyMMaZYhcFBnG/w74AfgP44weF+Vd12uBcSkWjg7zjNUoe8XMa2MtcbdJu2JgD07du3emsSRiURnp8OWM3BGGN8VRgcVFVF5EtV7QN8XM1rdQQ6AMW1hjY4NZJ+ODWFtj77tgFSq3m9ykUnEVKUQwQFFhyMMcaHP30Oc0Xk+OpeSFWXqGpzVW2vqu1xAsLxqroD+BwYLSIRItIB6ATMre41K+XOkk4gy4KDMcb48Cc4nIQTIFaJyEIRWSQiCys7SETeA+YAXURkq4hcW96+qroMmAwsB6YAN6uqx7+3UA0+s6RtIpwxxhzgT4f0eVU5sapeVsnr7Us9fxh4uCrXqjK35tDGMrMaY8xBKuuQDgU+VtVetVSe2lWcmbVJDnssOBhjTInKUnZ7gOUiUuaw0nrPMrMaY0yZ/E28t0JE5gDZxRtV9YKAlaq2uDWHZqGWfM8YY3z5ExweDXgpgiU8EsJjSA7JIsNyKxljTIlyg4OIdFLVNar6nYiEqWqRz2sn1E7xakF0EolkWrOSMcb4qKjPYZLP49JzDl4KQFmCIyqROHWCg2r1JlwbY0xDUVFwkHIel/W8/nIzs3q8SnZB4KdWGGNMfVBRcNByHpf1vP6yzKzGGHOIijqk24jIkzi1hOLHuM8bztDW6CQii5zkexm5hbROiApygYwxJvgqCg73lvMY4G8BKEtwRCURXpBBCF6rORhjjKvc4KCqr9ZmQYImOhlBibfke8YYU8LfleAaLneWdKJYcDDGmGIWHKIOpO22WdLGGOOw4BCdCEBSSKYFB2OMcVUaHETkaBH5VkR+c5/3FJHSHdT1l1tzaBmea81Kxhjj8qfm8AowHvC6z5cAYwJWotpWnJk13DKzGmNMMX+CQ4yq/lz8RJ0cEw3nWzQiDkLCaB6WbcHBGGNc/mRl3eOu66wAInIesCOgpapNIhCVRIoNZTXGmBL+1Bz+BLwKHCMim4B7gBsqO0hEJorILhFZ6rPtcRFZKSK/i8gnIpLg89q9IrLWXat6ZBXeS9VFJ5GIpe02xphiFQYHd5nQXqo6HGjpPh6gqhv9OPfrwOmltk0DuqtqT2A17sxrEekKjAa6ucc87167dkQlEW9pu40xpoQ/y4Te7j7er6rp/p5YVWcCe0ttm+qzLsQvQBv38bnA+6qar6obgLVAP3+vVW3RScR6LTgYY0wxf5qVvhWR20WkpYjEFd9q4NrXAN+4j1sDW3xe20o5yf1EZJyIzBeR+WlpaTVQDCAqkRjPfgqKvOQVWtpuY4zxp0P6evf+Lz7bFGhX1YuKyN+BIuCd4k1l7FZmWnBVnQBMAOjbt2/NpA6PTiKycD+g7M8tJDK89lq0jDGmLqo0OKhq25q8oIiMBc4CTtEDS69tBXyv0wZIrcnrVig6mVAtJJp8MnILaREXWWuXNsaYusifmgMicgzQFSj51lTVdw/3YiJyOnA3MFRVc3xe+hx4110zohXQiUOXJg0cd5Z0kli/gzHGgB/BQUTuA0YAxwDfAiOB2UCFwUFE3gOGASkishV4AGd0UgQwTUQAflHVG1R1mYhMBpbjNDfd7HaG147o4uR7FhyMMQb8qzlcCvQGFqrqFSLSEnipsoNU9bIyNpe7RoSqPgw87Ed5al6Upe02xhhf/oxWynV/xReJSFOc2dFHBbZYtax4TQdL222MMYB/NYdF7kzmicB8IANYGNBS1bbiNR0kk/25NkvaGGP8Ga1UPJT1fyLyLRCnqg0sODhrOjQPzWav1RyMMcavDulBZW3zzdRa74WGQWQ8LQqy2WDBwRhj/GpW+ofP40igD7AIGBqQEgVLVBIpHkvbbYwx4F+z0hm+z0WkPfBIgMoTPNFJJGVnkZFnwcEYYw57DWk3I2v3mi9KkEUlEW+jlYwxBvCvz+G/HMhzFAIcBywLZKGCIjqJOF1izUrGGIN/fQ5LfR4XAZ+o6o8BKk/wRCUR48lgf74FB2OM8afPodxZzQ1KdDIR3hwKC/Ip9HgJDz3sFjdjjGkw/GlWWkTZ6bMFUFU9vsZLFQzRzlyHBDLJyC0kOTYiyAUyxpjg8adZaToQCrzlPv8DkAm8HahCBUWp/EoWHIwxjZk/wWGQqp7o83yRiPykquMDVaig8MmvZJ3SxpjGzp+G9VgRGVD8RET6A7GBK1KQHJRfyYKDMaZx86fm8EfgNREpXugnF2f954Yl2tJ2G2NMMX9GK80DuotIsvt8T8BLFQxRPmm78ywzqzGmcSu3WUlEzhSRdj6brgemisjHInJk4ItWy5pEo2FRJEimzZI2xjR6FfU5/AvYAyAio3Cakm4CpuLHSnD1kUQnkRJizUrGGFNRcFBVzXYfXwC8oqq/quqLQIvKTiwiE0Vkl4gs9dmWJCLTRGSNe5/o89q9IrJWRFaJyMiqvqFqiUqiWWg2+3MsOBhjGreKgkOIiESLiACnAN/7vObPJIDXgdNLbbsH+E5VOwHfuc8Rka7AaKCbe8zzIhLq1zuoSdGJJFvNwRhjKgwOz+Ks2/ArsEZV5wKISC+cdaQrpKozgb2lNp8LvOE+fgM4z2f7+6qar6obgLVAP3/fRI2JSiIBS9ttjDHljlZS1ZfdZUFbcPCa0bup+lDWFqq63T3/dhFp7m5vDfzis99Wd9shRGQcMA6gXbt2Ze1SddFJxKnNczDGmAonwanqZlWdp6oen23b3DUdapKUdflyyjRBVfuqat9mzZrVbCmikojxZpKRk1+z5zXGmHqmtlOP7hSRlgDu/S53+1agrc9+bYDUWi4bRCcTghdv3v5av7QxxtQlFc1zqOE2GwA+B8a6j8cCn/lsHy0iESLSAegEzA3A9SvmzpIOz9+Hx1tmxcUYYxqFimoOnwCIyNSqnFhE3gPmAF1EZKuIXAs8CpwmImuA09znqOoyYDKwHJgC3OzblFVrfGZJZ9ksaWNMI1ZR+oxQEfk7cKyI3Fr6RVV9pqITq+pl5bx0Sjn7Pww8XNE5Ay66OPmeM5w1Pjo8qMUxxphgqSg4XIYz+S0MqOGe3zoqypmTl4iNWDLGNG4VDWVdATwsIr+r6he1WKbgscysxhgD+Dda6UcReUxEfnFv/xaRpgEvWTBExKMSQoLYRDhjTOPmT3B4FSgErnRvBcBrgSxU0ISE4I1MJMmalYwxjZw/i/10UtWLfZ7/Q0QWB6pAwSbRSSRkZbLFgoMxphHzp+aQJyIDi5+4S4bmBa5IwSXRySRZn4MxppHzp+ZwE/CWiBRnYs0FrghckYJLopNIDtluwcEY06j5s0zoQqCbiCQB0mCXCS0WlUSiZNlqcMaYRs2fmgMAqlo6/XbDFJ1omVmNMY1ebSfeq/uik4mggLyczGCXxBhjgsaCQ2lufiXJaRwVJWOMKUulzUoiEoKzdGd73/0ry61Ub7mzpEPy9gW5IMYYEzz+9Dl8hrPwzhLAG9ji1AFuzSG8IB1VxVlC2xhjGhd/gkN7Ve0R8JLUFW7NIV4zyS7wEBvhd5+9McY0GP70OXwrIsMDXpK6IurgtN3GGNMY+RMcZgFfiEiWiOwVkX0i0nB7a33TdudYcDDGNE7+tJn8FxhMY+lzCGtCUXgsiUVWczDGNF7+1BzWAItUtVBVPcW36lxURO4QkWUislRE3hORSBFJEpFpIrLGvU+szjWqwxuZSKJkWtpuY0yj5U/NIRX4XkS+BvKLN1Z1KKuItAZuBbqqaq6ITAZGA12B71T1URG5B7gHuLsq16i2qCQS07PYZTUHY0wj5U/NYSswG4jDWS60+FYdYUCUiIQB0TgB6FzgDff1N4DzqnmNKpOYZBIk0/IrGWMaLX8S7/2jJi+oqttE5AlgM06G16mqOlVEWqjqdnef7SLSvCavezjCYpJJlKXW52CMabT8mSE9DWcS3EFUdURVLuj2JZwLdADSgQ9EZMxhHD8OGAfQrl27qhSh8mtEJ9k60saYRs2fPof7fB5HAhfi0/dQBacCG1Q1DUBEPgYGATtFpKVba2gJ7CrrYFWdAEwA6Nu37yFBq0ZEJxFHDlk5uQE5vTHG1HX+NCv9WmrTjyLyYzWuuRkYICLROM1KpwDzgWxgLPCoe/9ZNa5RPdHJABRlN9zpHMYYUxF/mpXifJ6GAH2AllW9oKr+KiIfAguBImARTk0gFpgsItfiBJCLyz9LgLkT4bDMrMaYRsqfZqVlOH0OgvNlvgG4rjoXVdUHgAdKbc7HqUUEX3Fm1lwLDsaYxsmfZqW2tVGQOsXNrxSanx7kghhjTHCUO89BRPqISAuf538QkY9E5EkRSaid4gWJW3OIKLTgYIxpnCqaBDcBpxkJETkJeAKYDOS5rzVcbs2hqTeTvMJqZQoxxph6qaJmpTBV3eM+Hg1MUNVJwCQR+S3wRQuiJjF4JLxkrkNkeGiwS2SMMbWqoppDqIgUfyueAnzv+1rgilQHiFAQkeik7baJcMaYRqiimsNk4AcRSQMKcNZ1QEQ6Ahm1ULag8kQmkpht+ZWMMY1TucFBVR8Ske9x5jRMUdXitRzCcbKqNmgalUiC7LOagzGmUapwKKuqzi5j28rAFafuCIlOJpEtbLXgYIxphPxJ2d0ohbppu63mYIxpjCw4lKNJ02QSyGZ/TkGwi2KMMbXOn/QZjVJITAoh4iE/e3+wi2KMMbWu3OAgIvsoYx0HnBxLqqpJAStVXeDOkvZk7Q5yQYwxpvZVVHNIqbVS1EXuLGm1zKzGmEaooqGsB+WNEJEknMV+iqUGqlB1gltzEMvMaoxphCrtkBaRUSKyGtgK/Oref1/xUQ1AcWbWvH1BLogxxtQ+f0YrPQycCKxy03ePBGYEslB1gltzaFJgmVmNMY2PP8GhyF3vOURERFWnAccHuFzBFxmPIkQU2WglY0zj489Q1v0iEgPMBt4UkV2At5Jj6r+QUPLC4ojNy6DQ4yU81KaEGGMaD3++8c7DWcPhdpzmpG3AWdW5qIgkiMiHIrJSRFaIyEARSRKRaSKyxr1PrM41akJBkwSSbJa0OVxeL3z7d0hdFOySGFNl/gSHe1XVo6qFqvqqqj4J/Lma130aJ5nfMUAvYAVwD/CdqnYCvnOfB1VRRCIJlrbbHK7URTDnOZj1ZLBLYkyV+RMcTi9j26iqXgn3CjsAACAASURBVFBE4oAhwKsAqlqgqunAucAb7m5v4NRYgkqjEkmULEvbbQ7Pyi+d+zVTIT8ruGUxpooqWkP6ehFZBHQRkYU+tzXA8mpc8yggDXhNRBaJyCtun0YLVd0O4N43L6dc40RkvojMT0tLq0Yx/BCdRIK7Gpwxflv5JcQ0h6I8WD0l2KUxpkoqqjlMBi4Gvnbvi28nqupl1bhmGM5opxdU9Tggm8NoQlLVCaraV1X7NmvWrBrFqFxoTDKJWHAwhyFtNexeDYP/Ak1bwrJPgl0iY6qk3OCgqvtUda2qXgxEAae5t+p+I28Ftqrqr+7zD3GCxU4RaQng3u+q5nWqLbxpMtGST1aWNQ0YPxU3KR17FnQ9D9ZMg7wGv3CiaYD8mSF9M04top17mywiN1X1gqq6A9giIl3cTafgNFN9Dox1t40FPqvqNWpKRFMnDhZkWvI946eVX0Gr4yC+DXQ7Hzz51rRk6iV/5jlcD/RT1SwAEXkE+Bl4vhrXvQV4R0SaAOuBq3EC1WQRuRbYjNOEFVThTZ3cg57sPUEuiakXMlJh23wY/g/neZsTIK41LP0Yel4S3LIZc5j8CQ4C+Da6F7rbqkxVFwN9y3jplOqct8a5+ZW8WRYcjB9Wfe3cH+NOAwoJcZqW5r0MuekQlRC8shlzmCoarVQcON4CfhGR+0TkPpxawxvlHdeguPmVsMysxh8rvoSkjtCsy4Ft3c4HTwGs+iZ45TKmCirqc5gLoKqPAeOAHCAXuEFVn6iFsgWfZWY1/spNh42znI5o8alYt+kL8W1t1JKpdypqVir5C1fVecC8wBenjnFrDmH5lpnVVGLNVPAWHWhSKiYC3c6DX16E3H0QFfSsMMb4paLg0ExEyk2T4abRaNjCIsgLiSKi0IKDqcTKLyG2BbQuoyut2/nw87POSKbjxtR+2YypgoqalUKBWKBpObdGIS8snmiPpe02FSjMgzXTocuZTid0aa2Oh4QjrWnJ1CsV1Ry2q+pDtVaSOio/PIHY3Aw8XiU0pFqDtExDtX4GFGY7/Q1lEXFqD3Oeg5y9BwY61EMLNu1lf24hw49pEeyimACrqOZg34RAUUQCiZJFZp6l0DDlWPklRMRB+yHl79PtfKdPYsUXtVeuGrZ2VyZXvDqXG95aSGp6brCLYwKsouBQt+YcBIkn0knb/efJv/F/Xy5n4uwNTFm6gyVb97M3uwBVDXYRTTB5Pc4w1U4jIKxJyeaf1+1m676cA/u17AWJHept01J2fhE3vL2QyPBQFOW5H9YGu0gmwMptVlJVG9wPJDU7AnbksGlPNj+v201e4cGL4EWGh9AqIYrW7q1DSgynHNuco5s3mm6Zxm3Lr5CzG445kMV+ytId3PD2AkJDhLN6tuT6IR3p2irOqT389DRk74aYlCAW+vCoKvd+vIT1aVm8fW1/vlm6g/fmbuaGIR1plxwd7OKZAPFnhnSjFpvQHLxZfHfHYFRC2JdTSGp6LtvSc9m2L5fUfTns27ebwn1ryN22jdl5Xv71TQ+OahbLyG5HMKJrC3q1SSDE+isaphVfQmgTOPpUALbuy+GvH/5Gj9bxDDgqiXd/3cxni1MZ0rkZf+5xMr31Sadpqe/VQS64/97+ZROf/5bKXSO7MOjoFDo2j2Xy/C08/d0a/nNJr2AXzwSIBYfKRCcDCqu+RvIySMrYRtL+LXTfvw32b4WMbVDgk7W1CWxpNox/hd3IyzPX88KMdbSIi2BE1yMY0a0FA45KtvWoGwpVp7/hqGEQGUehx8ut7y3Cq/Dc5cdxZHIMfzq5E2//uonXftrAeavzmR3Tmshf3yfx+KvqxQCHxVvSeejL5Qw/pjk3Du0IQIu4SK4ceCSvzt7AjcM6cnTz2CCX0gSC1Oc28759++r8+fMDe5Fln8IHYw/eFtMc4ls7SdXi2zqP49tAXBunmeG7hyAyjuyR/2Wq5zimLtvJjFVp5BZ6iIsMY/gxzRnZ7QgGd25GbITF53prxxJ48SQ4+2nocxX/nrKSF2as49nLjuPsXq0O2jWv0MPHC7dRNP0h/lDwIRfGvMZFQ47noj5tiAwPDdIbqNi+7ALOenY2AF/dehIJ0Qf6VPZk5TP4sR84+Zjm/O/y44NVRFNFIrJAVcvKb3dgHwsOlSjKd2a/RiYcCAhhERUfs3M5fHwd7FwKfa6CEQ+TFxLFrDW7mbpsB9NX7GRfjjP6qWV8JB1SYjiqWQwdUmI5qlkMR6XE0CYxul78smzUZjzq3O5czcxU4cqJc7msX1v+dUHPcg/xbF9K6Esn8kLMTfx7z0mkxDbhqkHtuW7IUUSE1Z0g4fUqV78+jznr9vDhjQPp2ebQpIFPfLuK535Yy9e3Dnb6VEy9YcEhmIry4YeH4adnIKkDnD8B2p7gvOTxMm/jPhZs2sv63dmsT8tmfVoWGXlFJYc3CQ3hyORoN3DEkhLbhNAQIUSEkBAhVIQQgRB3W2gIzmvuzR8hAv2PSiYppknlO5tDvXASNIlh1yWfcebTs0iKacJnN59EVJMKvuRV4X/90Njm/DL4TV6auY4Zq9K4cuCRPHRu99oreyWe+W4NT05bzT/P686YAUeWuc/+nEJOeux7+ndI5pWxFX7PmDrGn+BgbRqBEhYBpz3kDHH85AaYOBKG3AlD7iIsNJyBHZMZ2DG5ZHdVZW92Aet3Z7MhLZt1u7OcoLE7mx9W7aLQE5ggHhcZxp0ju/CH/kdaTeVw7NsIO5fgPe3/uGPSYrLyi3jvugEVBwYomRAnPz7GwAsLGXh1Px76YjkTf9rAaV1bMLhTDS196ymEPWuh+bGHfeisNWn8d/pqzuvdij/0b3fwi6umQF469BpNfHQ44wYfxX+mrWbxlnR6t7WU5A2J1RxqQ95++OZu+O09J5XCBS9DytF+H17k8ZJd4MHrVbyqeFRRBY/73OsFjxY/Vvz9RDNyC/nv9NX8tHYPXVvG8dC53ejbvv7O3q1Vc56Hb+/lzRM+5f5ZOfz7wh5cekK7yo8D2LUCnh8AZzwO/ceRV+hh1DOzyM738O0dQ4iPCq9e2VThoz/C0g/hvBeg9+V+H7p9fy6jnplNSmwTPr35RKKb+Px+XD8D3r7QeXzLQkg8kqz8IoY89gPdWsXx1rX9q1duU2usWamuWfYpfHm7k4tn5D+h77UHp3cOAlXl6yU7+OdXy9m+P48Ljm/NPWccQ/OmkUEtV5332pnkZOym+44HOKtnK54e3Rs5nM/yfwOcDK3XOOs8/LYlnQte+Jlze7XiyUt7V69sv06Ab+6Cpi0hOw3GfOSMqKpEQZGX0RPmsGpHJp/fchIdm/mMQtq9Bl45xRmMkb7JCThnPw3AyzPX8/DXK5g0bgD9j0ou5+ymLvEnONiYytrU7Ty4cQ4cOQi++gu8NAQ+ug6mj4d5r8Dqb2HHUmdtAH+Dttfr7L93A2xbCGu/c5LAbZwNW+c7I2p2r4X0LZCVBvmZTpODS0QY1bMl3/1lKDcN68gXv6VyyhM/8ursDRR5vBVcuBHL3o1unsP7GT1pmxTNw+d3P7zAAM6EuM1znKVFgV5tE7j55KP5eNE2pizdUfWybZkH3/4NOp8ON82B5E4w6UqntlKJf32zgoWb0/n3RT0PDgw5e+HdSyAkHMZ8CMdfCYvecf6mgDEDjqR50wj+M3W1ZQxoQIJWcxCRUGA+sE1VzxKRJGAS0B7YCFyiqhWuslPvag7FVGH+q7DkI8jY6nxBeIsO3qdJU58hsq0hLNJZD6D0LS8dtApf4hIK4dHQ/QIY8X8QGQ/A+rQsHvxiOTNXp9GlRVPGn9uNAfZr8CC68E3k81s4t/AR/nnjGHq0iT/8k6Sthv+dAKc/CgNuBKDQ4+X8539ie3oe394xhJTYSkbFlZa9x/nBERIC1890aibpW5xf/KFN4I/ToekRZR761e/bufndhVw1qD0PntPtwAtFBfDW+bB1Hoz9Atr1d875zHFOkDjLydz/5pyN3P/ZMt68ph9DOtdQv0ltK8yD8MZRY67TzUruWhF9gTg3ODwG7FXVR0XkHiBRVe+u6Bz1NjiU5vVA1i5nUt3+Lc7Euv1bD755Cpz/7L636KRDt0UmQEgoFOU5I6YKc537ouL7POc/QVEeZG53+kGatnSaCDqdBjhNTVOX7+ShL5azLT2Xc3q14u+jjqVFXOP4j1OZzc+dQ0jaMqaeOo1rBh9V9RO9cCI0iYFrp5ZsWr0zk7Oenc3Qzs2YcEUf/2skXg+8cxFs/Mk5X6veFHq85BZ6KNq6iIRJ51IQ35G1oyaRrZHkFXnJK/SQV+ghM6+If329gs5HNGXSuIE0CXMbFFTh8z/Borfhgleg58UHrvfF7bD4Hbh1McS3Jr/Iw/AnfizpqzjsmlSwLfvEGThy8RvQ5fRglybg6mxwEJE2OOtQPwz82Q0Oq4BhqrpdRFoCM1S1S0XnaTDBIZi2LYBPb4a0FdDrcjj9kZLVynILPLwwYy0vzlxPWIjQo3U8rROinFxSie59QiStEqIO7rhswJZtSOXo13syK/4sTrnjtep9Cc58HL7/J9yxzKkhuorb8J+4uBcX9WlTwQl8/PAI/PhvOPtpCntfySNfr+CNnzfidf97nxyyiFfCn+AHb2/GFf4Fb6kW5eZNI/j05hNplRB1YOPsp2D6AzD0bjj5bwdfL32zU3voew2c+TgAk+dt4a8f/c7LV/bltK71KKV3UT4819d5TxFx8MfvoFnnYJcqoOpycPgQ+BfOokF3usEhXVUTfPbZp6qHrKkoIuNw1rSmXbt2fTZt2lRbxW64ivKdL6pZT0JMMzjrv3DMmSUvb9qTzQsz1rE+LZtt6bnsyMjD4z347yYhOrwkcBwRF0lYqPOlKW7md5EDOeCLv0+Lv1gP5+s1NERKbmEhzpyPMHeuR1jJayE0CQshNiKUmIgwYt1bjHuLjQir0rDdzLxCHvvvY/xf/mNkXvopTY89+bDPcZDda+G5PjDyERh4c8lmj1e5bMIvrNiewZQ7htDa9wu7LGumO7WG3pez6+T/cPN7i5i3cR8X92lDlyOaEhkeSmR4KF02v0+P3/6P1M5XsH3QeCLCw9zXQkiJjTh4pvaKL2DSFU7fyEUTyx448fkt8NskuO03iGtJkcfLqU/+SGR4KF/fOrj+5BP75UWYcjcPh97IXWGTaBKbDNd9V9LU2hDVyeAgImcBZ6rqTSIyjMMMDr6s5lDDtv/m1CJ2LoHuF8EZj0HMof0NRR4vuzLz2ZaeW5KEMLU4EWF6Hjsz3eDh/mkplHRUFv+1aclr/v/9qYJXlSKv+t1fX56o8OLA4Xxxlg42zvOQg56npudyXdq/GBW9nPC/roXQGqgtvXiS05/0x+kHbd68J4fTn55J77YJvH1t//K/aNM3O/0Mca1ZeNpkrp+0gqy8Ih69sAfn9m596P5T73OWLC0VkA6SuhheOwOad4WrvoTwcoLTvo3wzPHQ7zo4498AfLZ4G7e9v5jnLj+Os3q2Kvu4uiQvA89TvZiX24rLC+5lQOgq3g5/hJBOp8Lo98pe2a8u2DjbSQEfX8Zn7Ie6OgnuROAcETkTiATiRORtYKeItPRpVtoVhLI1bi17wXXfw+z/OjWJDT/CmU84o6x8hIU6acpbVfaLNoC8Xme+h8fr3lTxeA5syy/0kpVfRHZBEVn5RWTlFZGd7z7OL37sISu/iPxCT0nQ8XiVIvc8uYUed5uXIo8S4i3izIjfCT/2nJoJDOD8Mv/uISflSouuJZvbJUfzj7O6cu/HS3hzzkauOrHDoccW5cPkK1Gvh4+OfoR7Jv5Gm8Qo3rq2H8ccUU46i1Mfgn2b4Nu/O3nBup5z8OsZqfDeaCfh5Oh3yw8MAIntoddlsOB1OOkOaHoEZ/Vsxf9+WMuT01ZzercjCKvjSSY9Pz1LaN5enpI7+eKWwfzj00TGbxvD+NWvw4xHYPh9wS7iofasg/cuhzZ94YqPA3cdVQ3aDRgGfOk+fhy4x318D/BYZcf36dNHTYDsWKr64hDVB+JUJ12hmrkr2CUKvrXfOf8eK76suXPu26T6z5aq45NUP71Zdffakpe8Xq+Onfirdrnva123K/PQY7+4Q/WBOH1lwtN65N1f6rWvz9P9uQWVX7MgR/XlU1T/r7nqlnkHtudnqb44WPXhVqrbl/hX/j3rVB9MVP3m3pJN3yxJ1SPv/lI/mL/Fv3MES+ZOzR/fQr+871T96vdUVVXNyivU0S/+rO/fd47zWS/7NMiFLCU/W/X5Qar/aqe6d0OVTwPM10q+X+tSWH8UOE1E1gCnuc9NsLTo5nTMnfKAs9LZ//o5vzY3zDponkSD5yly5oZk74alHzvDfzsOr7nzJ7SDP811JkQu+cDpGP3oj7BrBSLCvy/sSURYKH+e/NvB805+nwzzX+WDiAv45/qO3DmiMxOu6ENcpB+zq8Oj4LL3nWGt714Ke9c782U+ud6ZF3PRRDjCzzxPSUdBz0th/kRnxB0wstsRdG8dx9PfraagqO7Oldn62UOEePJZ3e12zuzREoCYiDBeu6Yf0zvcxULv0RR+dD3sXBbkkrpUnUm0O5fBha84NbcAshnSpnJpq2DqP2D9D86Q2sh4Z3Gbzqc799H1KOWGp9CZ7btruXPbuRyydx0Y3ltyc4cBq+fg4489Gy59OzBly9wJc56Dea9CYbZzrcF38kVac255bxF3jujMn4Z3gp3L8UwYziJPe66XB3jysr4Mrcrcgt1r4dVTnSakjqfA3JcOmnfhtz3rnKA28GYY8U8Afli1i6tfm8cdp3bmppM71rk1TNI2rSThtROZ1uQUht/13iFp0wuKvNz/9jT+vGEcEVExxN86O/h/53Nfhq/vhGF/g2EVjvKvVJ3skK5JFhxqWX6mk19n9RRYPdX5UpUQaNsfOo90gkWzY4KeEgRwfmWlb3ZmBu9a5gSBXcudwOB1az4hYc4M4rhWzq/psEjnFu7eh0VAWJRzH+7eH32qs38g5eyFX16AX1+C/P3QaQT/yTuHF9cn89kfe3LEpDPw5O7njsRn+deVp9E2qRpLdW76Gd481wn6fa+BUU9W7fP7eJwzwum23yG2GarKlRPnMmvNblrERfCH/kcyul/baqdlKSjysnJHBp1bNK3yOhger/LL4+dxfM7P7LhqDh06lJ3nrMjj5bk33+XGjbeRmnA87W/9Ggk9jLxXWWnw64uQswdGPuzMaamqLfOcQQIdT4bLJlW7o9yCgwkcrxe2L3JSfqye4ox0AqeZpONwCI9xfnV7i5wJWsX3ZW3zFjq/6L1FB+69hU6Tju9rXvdXvAglA2CLH5d8obljZnP2QUHmgfLGt3VG37ToCs27OdlKUzpDWB1OV5633/m1OOd/kLuXedKDXE8Ig2QJEzo8xTV/GFMzCwWtmuIMPjjtITicLz9fu9fAcyfAibfBaeMB50t4xqpdvP7zRmat2U14qDCqR0uuHNSe49om+D1HJDOvkBmr0pi2fCc/rNxFZn4RRzeP5bnLjyu/470C73z6BX9YPIaVR1/HMWOeqHBfr1f55NV/ceG2f/NT88sZdOPzlZd73yb4+Vl00VtODVQEjuiJXD6p3BnqFcpKc0akhTWBcTNK5iFVhwUHU3syUp1FkVZNcX6Nqtf5dRMS5qTqCAlzZm6HhB76PCTc+VIKCXNuoeHutjCf18IP/Foq+Zs9MFzWeayUbGgS6wSAFm4gqM9j1guyYf5r5M98ioi8NBZ2uYPjRj9Q92Yhf3it0z91+5JDhkCvS8virTmb+GjBVjLzi+jROp4rBx7J2b1alRngdmbkMW35TqYt38nP63ZT6FGSYppw6rHN6dE6nqe/W0tmXiH3n92Vy/u18/vfYs66PRS8cR4nhG0g+q6lEFV5mnFVZf7z13BC2sdMavsPLrr6L2XOkylMXULG9CdIXP8FXoRPvIN5oXAU7WUHL0Q8R2hsMmFjPjxoVFqlPEXw1nlO+pJrp0FLZyGpjLxCcgs8Vc5aYMHBmIakMA/dvhhp279uNN2Vtmulk4p88J/hlPvL3CUrv4hPFm3jzZ83smZXFonR4Vx6QjvGDGhHXqGHb5c5AWHxlnQAjkyOZkTXFozodgTHt0ss+VJOy8znz5MXM2vNbkb1aMkjF/SoNNX5nqx87vvv87zgeZD84eOJGHK7329Ni/LZ9swIUvYv47kO/+P2Ky5GRFixPYP1C76j3YqX6J37C9kawbueU/gx+RI6H92FQR2TWZeWxbfTv+XFkMdICC8kbPTbhBzt5wTKaffDT0+XpF7PKSjijZ838dLMdQw8KpkXxvTx+z34suBgjKldH1zlzNi+/fcKO3BVlTnr9/Dmz5uYunwHvhPue7WJ5zQ3IHRqHlturcDrVV6auZ4npq6iZXwkz152HMe1K7vJxetVrnl9Ln/edCPHxuYQfvviw0+yl5VG5rMnkpFXxB1xT9EsczlXej+hf8hK9ktT5re4BO8J19HnmI6HrK64aU82//nge25MvZdOIamkDfs3LYddV/H1ln8Gk6+EvteQN/IJ3v11M8/PWMvurAJO7tKMP5/WpWpJH7HgYIypbTuXwwsDYchdfk8g25aeyycLtxIfFc6pXVvQMr6ciXeeQlj6Efw+CY46GQbcBKFhLNi0j1vfW8TOjDzuGtmF6wYfdciM8gkz17Foyhu80ORpOOc5OP6Kqr2/1EUUvTKSQq8QRR45UUdQ1P9PxA26ptIOZ1Xli7mrSPlmHIP4jV9aX0XvsU8Q2aSMGk/aanh5ON6Uzkzq/hJPz9jMjow8BnVM5i8jOtPnyOqNnLLgYIypfZOvhHU/OLWHGug8pSAbFr7lDPPdvwViW0DWTmjZG855Flr2ZH9uIfd89DvfLN3B0M7N+M8lvUpSni/avI/RL87mx5h7aBEfg9z4c/VmuC/7xBlN1ucqJ83MYQ5q2LM/i9WvXc/A9C/5PmwwsZdOoF8nnxFw+Vnoy8MpyNjF5SGPsSA9hj5HJvKXEZ0Z1DGl6uX2YcHBGFP7diyFF0+EoffAyfdW/Tw5e53RWr++CLl7od1AJ03H0afBis/g67ucfU68DYbejYZF8M6vm3noy+XER4Xz9KW96dY6nlHPzOLswqncXfSCkxLkmFE1916rSpX1nz3MUYsfZ663C1O6PcFtZw+kaUQo2ydexhHbpnJFwT1ktBzEX0Z0YVjnZjU6AMGCgzEmON7/gzOb/vbf/RoRdJD9W53huwteh8IcZ/7MibfDkQMP3i9nr5NIcPE7kHw0nP0MtD+RFdsz+NO7C1m/O5ujm8WSunsvi+L/SpOUDnDNt3WqMz9/8QeEfnYTWz2J3BF+H6eHL+b6vFd5JXIsbc76GyO7tQjIqDQLDsaY4Nj+mzM2v8MQaHUcRKdATIqTEj462XkcnQJNfCbw7VrpjMxZMtkZltzjYqdWUNnQz3Xfwxe3OZMe+14Dpz5ITkgMD3y2jA8WbOX9rnMYsP5ZuPobZ4neumbzLxS9M5rcgkKiNYedLYfT4o8fEBrAWeUWHIwxwTPtAScHVM5uZwZ2WcJjnDkREXGwc6kzI73PWCcVR0I7/69VkA3fPwy/vgCxR8Co/8AxZ5K6I5WWr/VH2g2EP0yumfcVCHvWOet0S4iTvj3A83IsOBhjgk8V8jOc5IU5e5z77DQnaGTvce5z9kDrPtDv+jLXEPHb1gXO0qa7ljvp0CPiYOGbcONPzoTIuqw4E0BFadJrSF1dz8EY05iIOL+EI+MhuWNgr9WmD4z70WmemvmYU2PpdVndDwzgZAKoavqSALDgYIxpWMKawNC7nIWM5k90RjiZw2bBwRjTMDXrUrJ8qTl8dSvJujHGmDrBgoMxxphDWHAwxhhziFoPDiLSVkR+EJEVIrJMRG5ztyeJyDQRWePe10BSFmOMMVURjJpDEfAXVT0WGADcLCJdgXuA71S1E/Cd+9wYY0wQ1HpwUNXtqrrQfZwJrABaA+cCb7i7vQGcV9tlM8YY4whqn4OItAeOA34FWqjqdnACCNC8nGPGich8EZmflpZWW0U1xphGJWjBQURigY+A21U1w9/jVHWCqvZV1b7NmjULXAGNMaYRC8okOBEJxwkM76jqx+7mnSLSUlW3i0hLYFdl51mwYMFuEdlUanMKsLtmS1wn2Puqfxrqe7P3Vf+Ufm9HVnZArQcHcZKTvwqsUNUnfV76HBgLPOref1bZuVT1kKqDiMyvLKFUfWTvq/5pqO/N3lf9U5X3Foyaw4nAFcASEVnsbvsbTlCYLCLXApuBi4NQNmOMMQQhOKjqbKC8pY1Oqc2yGGOMKVtDnCE9IdgFCBB7X/VPQ31v9r7qn8N+b/V6sR9jjDGB0RBrDsYYY6rJgoMxxphDNJjgICKni8gqEVkrIg0qL5OIbBSRJSKyWETq7aLZIjJRRHaJyFKfbfU+4WI57+tBEdnmfmaLReTMYJaxKhpykswK3lu9/txEJFJE5orIb+77Gu9uP+zPrEH0OYhIKLAaOA3YCswDLlPV5UEtWA0RkY1AX1Wt1xN0RGQIkAW8qard3W2PAXtV9VE3qCeq6t3BLOfhKud9PQhkqeoTwSxbdbiTUVuq6kIRaQoswMl5dhX1/zMr771dQj3+3Nx5ZDGqmuVONp4N3AZcwGF+Zg2l5tAPWKuq61W1AHgfJ5GfqUNUdSawt9Tmep9wsZz3Ve815CSZFby3ek0dWe7TcPemVOEzayjBoTWwxef5VhrAB+1DgakiskBExgW7MDXMr4SL9dSfROR3t9mp3jW9+KpKksz6otR7g3r+uYlIqDvBeBcwTVWr9Jk1lOBQ1qS6+t9edsCJqno8cAbO+hdDgl0gU6kXgI5Ab2A78J/gFqfqqpoksz4o473V+89NVT2q2htoA/QTke5VOU9DCQ5bgbY+z9sAqUEqS41T1VT3fhfwCU4zWkOx023/LW4HrjThYn2gqjvd/6Re4GXqPrm7dgAABM5JREFU6WdWUZJM9/V6+5mV9d4ayucGoKrpwAzgdKrwmTWU4DAP6CQiHUSkCTAaJ5FfvSciMW6HGSISA4wAllZ8VL1SnHAR/Ey4WB8U/0d0nU89/Mz8SJIJ9fQzK++91ffPTUSaiUiC+zgKOBVYSRU+swYxWgnAHXL2FBAKTFTVh4NcpBohIkfh1BbAyYX1bn19byLyHjAMJ33wTuAB4FNgMtAON+Giqtarzt1y3tcwnKYJBTYC1xe3+dYXInISMAtYAnjdzX/DaZuv759Zee/tMurx5yYiPXE6nENxfvxPVtWHRCSZw/zMGkxwMMYYU3MaSrOSMcaYGmTBwRhjzCEsOJj/b+/+Qqsu4ziOvz9BF9rFLpSIbuyiv0Qwm2V/bEV1I4E1CEzyYgWREYQYREIEWQTbRUThhWQwEetCMuiPSGSbJoElLdg0r4SCiCiISDFd89vF8z3z535nO/NIyHY+Lxg7e37P8+x5BuP7e36H8/2amdU4OJiZWY2Dg5mZ1Tg42LwjaTIzZo5L2i1pcYv+J2e7/n+RtELSOy36PCDpsxmurcoMm8fza86pUzKT79KLXbNZg4ODzUenI6I7M6CeBTZc7gU1ExFHIuKFdsZKugb4ANgQETcDq4BnJT3SpO8l14LPzMZmUxwcbL77GrgeQNKmPE2MS9o4vaOknZIerfy8S9IaSf2S9kjal/nuByt91qnU0hiXNFBpPylpIJMhfinpTkkjkk5IWpN9pk4Fef0bSaP5/aYW+3oeGKpkDv0DeAl4OecbkvSWpGFgQNISSV/k/Nuo5BuTtD5PID9I2tYIBLmHLZIOA3dfzB/dFj4HB5u38o55NTAmqQd4ClgJ3AU8I2n5tCHbsw+SuoB7gL15rRtYC9wGrFUpBnMtMAA8mNfvkNRIdXwVMBIRPcDfwBuUeiJ9wJYmyz0O9EbEcuBV4M0W27uVUmOg6ki2N9wIPBwRL1I+lX0o5/+E8klYJN2S+7o3k7FNAk9W9jAeESsj4lCL9ViHueTjqNllsChTEkM5ObwPPAd8HBGnACTtAe4DRhuDIuKApK2SrqYUP/koIv4taXbYHxF/5dhjwDJgCSUA/J7tu4BeSsqPs8C+nHoMOBMRE5LGgOuarLkL2CHpBkpqhitb7FE0zyxcbdsdEZP5ujf3RER8LunPbH8I6AG+y30u4nzStUlK4jmzGgcHm49O513wlEykNhc7KXfOTwBPV9rPVF5PUv43ZptzIs7nnjnXGB8R52Z4D+B1YDgi+lTqB4y0WOdRYAUXJpDsAarVDU9NG9MsmAjYERGbm1z7pxJczC7gx0q2UBwEHpO0OLPX9lFOFdMNARsBIuJoizkPA/dLWprP6dcBB9pcXxfwS77un0P/rUC/pG6ATJw2AAzO0P8g+bhI0mqgUaRmP/B4npYatYSXtbMB6yw+OdiCkLWAh4Bvs2l7RIw26febpB8pj4ZazfmrpM3AMOUOfG9EtJueepDyWGkT8NUcf/d64L1M2S7g7Yj4dIYhrwEfSvqeEsB+znmOSXqFUknwCmCC8mb3T23uwzqEs7JaR8nPRIwBtzfeYzCzOj9Wso4hqVH45F0HBrPZ+eRgZmY1PjmYmVmNg4OZmdU4OJiZWY2Dg5mZ1Tg4mJlZzX/l/Qq4Tc4qQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(1,max_N),SSE_train_st[1:],label='training Error')\n",
    "plt.plot(np.arange(1,max_N),SSE_test_st[1:],label='testing Error')\n",
    "plt.title('Steel Error')\n",
    "plt.ylabel('Total Sum of Square Error')\n",
    "plt.xlabel('Polynomial Order')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2025 Aluminum Price= $ 823.5\n"
     ]
    }
   ],
   "source": [
    "'''D'''\n",
    "year_al=np.linspace(0,9,len(al_y_n))\n",
    "year_al_norm=y_n=(year_al-year_al.min())/(year_al.max()-year_al.min())\n",
    "rand_al=random.sample(range(0,len(y_n)),len(y_n))\n",
    "train_al=rand_al[:int(len(y_n)*train_per)]\n",
    "Y_al=y_n[np.sort(train_al)]\n",
    "Z_train_al=np.block([[Y_al**0]]).T\n",
    "for i in range(1,max_N):\n",
    "    Z_train_al=np.hstack((Z_train_al,Y_al.reshape(-1,1)**i))\n",
    "print('2025 Aluminum Price= $', round((Z_train_al@A_al)[1],2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2025 Steel Pirce= $ 272.1\n"
     ]
    }
   ],
   "source": [
    "year_st=np.linspace(0,9,len(st_y_n))\n",
    "year_st_norm=y_n=(year_st-year_st.min())/(year_st.max()-year_st.min())\n",
    "rand_st=random.sample(range(0,len(y_n)),len(y_n))\n",
    "train_st=rand_st[:int(len(y_n)*train_per)]\n",
    "Y_st=y_n[np.sort(train_st)]\n",
    "Z_train_st=np.block([[Y_st**0]]).T\n",
    "for i in range(1,max_N):\n",
    "    Z_train_st=np.hstack((Z_train_st,Y_st.reshape(-1,1)**i))\n",
    "print('2025 Steel Pirce= $', round((Z_train_st@A_st)[-1],1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 392,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It does not change mthe answer in problem 3d. Steel is still cheaper than aluminum\n"
     ]
    }
   ],
   "source": [
    "'''E'''\n",
    "print('It does not change mthe answer in problem 3d. Steel is still cheaper than aluminum')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "1. <https://en.wikipedia.org/wiki/Direct_stiffness_method>\n",
    "\n",
    "2. Aluminum and steel price history on <https://tradingeconomics.com>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}