Skip to content

Instantly share code, notes, and snippets.

@vladiant

vladiant/cubePCA.ipynb

Last active October 31, 2020 22:27
Show Gist options
  • Save vladiant/5f56c08c2543e75938e709340ac5766d to your computer and use it in GitHub Desktop.
Save vladiant/5f56c08c2543e75938e709340ac5766d to your computer and use it in GitHub Desktop.
Download ZIP
3D points curve fitting
Display the source blob
Display the rendered blob
Raw
cubePCA.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.decomposition import PCA"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"initial_vertices = pd.DataFrame([\n",
" [-1.0,-1.0,-1.0],\n",
" [-1.0,-1.0, 1.0],\n",
" [-1.0, 1.0, 1.0],\n",
" [1.0, 1.0,-1.0],\n",
" [-1.0,-1.0,-1.0],\n",
" [-1.0, 1.0,-1.0],\n",
" [1.0,-1.0, 1.0],\n",
" [-1.0,-1.0,-1.0],\n",
" [1.0,-1.0,-1.0],\n",
" [1.0, 1.0,-1.0],\n",
" [1.0,-1.0,-1.0],\n",
" [-1.0,-1.0,-1.0],\n",
" [-1.0,-1.0,-1.0],\n",
" [-1.0, 1.0, 1.0],\n",
" [-1.0, 1.0,-1.0],\n",
" [1.0,-1.0, 1.0],\n",
" [-1.0,-1.0, 1.0],\n",
" [-1.0,-1.0,-1.0],\n",
" [-1.0, 1.0, 1.0],\n",
" [-1.0,-1.0, 1.0],\n",
" [1.0,-1.0, 1.0],\n",
" [1.0, 1.0, 1.0],\n",
" [1.0,-1.0,-1.0],\n",
" [1.0, 1.0,-1.0],\n",
" [1.0,-1.0,-1.0],\n",
" [1.0, 1.0, 1.0],\n",
" [1.0,-1.0, 1.0],\n",
" [1.0, 1.0, 1.0],\n",
" [1.0, 1.0,-1.0],\n",
" [-1.0, 1.0,-1.0],\n",
" [1.0, 1.0, 1.0],\n",
" [-1.0, 1.0,-1.0],\n",
" [-1.0, 1.0, 1.0],\n",
" [1.0, 1.0, 1.0],\n",
" [-1.0, 1.0, 1.0],\n",
" [1.0,-1.0, 1.0]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"initial_vertices = initial_vertices.drop_duplicates()\n",
"\n",
"# [-1.0,-1.0,-1.0],\n",
"# [-1.0,-1.0, 1.0],\n",
"# [-1.0, 1.0,-1.0],\n",
"# [-1.0, 1.0, 1.0],\n",
"# [ 1.0,-1.0,-1.0],\n",
"# [ 1.0,-1.0, 1.0],\n",
"# [ 1.0, 1.0,-1.0],\n",
"# [ 1.0, 1.0,-1.0]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>0</th>\n",
" <th>1</th>\n",
" <th>2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <td>0</td>\n",
" <td>-1.0</td>\n",
" <td>-1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>1</td>\n",
" <td>-1.0</td>\n",
" <td>-1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>2</td>\n",
" <td>-1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>3</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>5</td>\n",
" <td>-1.0</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>6</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>8</td>\n",
" <td>1.0</td>\n",
" <td>-1.0</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>21</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" 0 1 2\n",
"0 -1.0 -1.0 -1.0\n",
"1 -1.0 -1.0 1.0\n",
"2 -1.0 1.0 1.0\n",
"3 1.0 1.0 -1.0\n",
"5 -1.0 1.0 -1.0\n",
"6 1.0 -1.0 1.0\n",
"8 1.0 -1.0 -1.0\n",
"21 1.0 1.0 1.0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_vertices"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"angle = np.pi * 60/180"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"rot_mat = [[np.sin(angle),-np.cos(angle),0],\n",
" [np.cos(angle),np.sin(angle),0],\n",
" [0,0,1]]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[[0.8660254037844386, -0.5000000000000001, 0],\n",
" [0.5000000000000001, 0.8660254037844386, 0],\n",
" [0, 0, 1]]"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rot_mat"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"initial_vertices = initial_vertices.dot(rot_mat)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>0</th>\n",
" <th>1</th>\n",
" <th>2</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <td>0</td>\n",
" <td>-1.366025</td>\n",
" <td>-0.366025</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>1</td>\n",
" <td>-1.366025</td>\n",
" <td>-0.366025</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>2</td>\n",
" <td>-0.366025</td>\n",
" <td>1.366025</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>3</td>\n",
" <td>1.366025</td>\n",
" <td>0.366025</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>5</td>\n",
" <td>-0.366025</td>\n",
" <td>1.366025</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>6</td>\n",
" <td>0.366025</td>\n",
" <td>-1.366025</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>8</td>\n",
" <td>0.366025</td>\n",
" <td>-1.366025</td>\n",
" <td>-1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <td>21</td>\n",
" <td>1.366025</td>\n",
" <td>0.366025</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" 0 1 2\n",
"0 -1.366025 -0.366025 -1.0\n",
"1 -1.366025 -0.366025 1.0\n",
"2 -0.366025 1.366025 1.0\n",
"3 1.366025 0.366025 -1.0\n",
"5 -0.366025 1.366025 -1.0\n",
"6 0.366025 -1.366025 1.0\n",
"8 0.366025 -1.366025 -1.0\n",
"21 1.366025 0.366025 1.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_vertices"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"initial_vertices.columns = ['x', 'y', 'z']"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"pca = PCA(n_components=3, random_state=1234)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PCA(copy=True, iterated_power='auto', n_components=3, random_state=1234,\n",
" svd_solver='auto', tol=0.0, whiten=False)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.fit(initial_vertices)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.33333333, 0.33333333, 0.33333333])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.explained_variance_ratio_"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0., 1., 0.],\n",
" [ 0., 0., 1.],\n",
" [-1., -0., -0.]])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[0].dot(pca.components_[0])"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[0].dot(pca.components_[1])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[0].dot(pca.components_[2])"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[1].dot(pca.components_[0])"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[1].dot(pca.components_[1])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[1].dot(pca.components_[2])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[2].dot(pca.components_[0])"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[2].dot(pca.components_[1])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.components_[2].dot(pca.components_[2])"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXRb5Z0+/lztki3ZWi3vexYnARKzdkonhEJapqTtdyhLF2hpf7RQvtPDnClNW8pAT3ugp+20p0M7zEzTgZkulBm+LUwXSllCgARsJ05M7CR2Ejvxql2yJGvX/f3hvi9X11ruvbpynKDnHA6xLV+9knWf9/N+Ps/n+TAsy6KKKqqooorVgeJcL6CKKqqo4t2EKulWUUUVVawiqqRbRRVVVLGKqJJuFVVUUcUqokq6VVRRRRWrCFWJn1elDVVUUUUV4sEU+kE10q2iiiqqWEVUSbeKKqqoYhVRJd0qqqiiilVElXSrqKKKKlYRpQppVVRRRRWrjlQqhZmZGcTj8XO9lKLQ6XRoaWmBWq0W/DtMCe+FqnqhiiqqWHVMTk7CaDTCarWCYQoKAc4pWJaFz+dDOBxGZ2cn/8dV9UIVVVRx/iAej69pwgUAhmFgtVpFR+NV0q2iiirWJNYy4RJIWWOVdKuooooqVhFV0q2iiiqq4GF6ehqdnZ3w+/0AgEAggM7OTkxNTZV97SrpVlFFFVXw0Nrairvvvhu7d+8GAOzevRt33XUXOjo6yr52VTJWxQUBlmXBsiwUimocUYU8uO+++9Df348f/vCHeP311/HYY4/Jct0q6VZxXoJlWWSzWaTTaaRSKaTTaajVaqhUKigUCiiVSjAMc14UY6oojof/dxRjc4uyXrOvyYR/vHFT0ceo1Wp897vfxQc+8AG88MILorS4xVANC6o4b8CyLFKpFGKxGMLhMBYXF/HWW28hk8lAoVBAoVCAZVmMjo7C6/UikUggmUwik8mgOguwCin44x//iMbGRhw9elS2a1Yj3SrWLLjRLJc8GYah0Ww6naYpBRLZZjIZeo1MJkO/ZhgGSqWyGgWfZygVkVYKhw8fxp///Ge8+eabeO9734tbb70VjY2NZV+3GulWsaaQzWaRSqWwtLSExcVFhMNhxGIxZLNZKBQKqFQqwaRJyFmhUNDHptNpJBIJxOPxHCKvRsJVcMGyLO6++2788Ic/RFtbG7785S/jH/7hH2S5dpV0qzinYFkW6XQa8Xicpgyi0SiSyWROZMolTikgJM0l4Uwmg2QyiXg8jkQigVQqhWw2WyXgKvDv//7vaGtrw3XXXQcAuOeee3Ds2DG8+uqrZV+7ml6oYtVBUgakCJbNZgFgRVRaSXAjZUKyZE1kLVyyr6Yi3l246667cNddd9GvlUolDh06JMu1q6RbRcXBsiwymQxNE5DjPIk8VaryP4axWAyxWAwmk0n07xJC5ZIwKdqR72ezWWi12lXbFKq4cFEl3SpkB5e0SPTIsixGRkawceNGaLXasokrk8kgGAwiHo9jcHAQGo0GarUakUgE2WwWer0eer0eGo1G9LXzkfCbb76Jq666CsA7EXm1IFeFFFRJtwpZQKLZVCqVkzLg5lK5/5eCpaUl+P1++Hw+JBIJ1NfXQ6VSYdu2bVAqlVRHOTw8jHg8jkOHDlEnKLvdDpPJJOm5ua+BbChkMyF559VMjVRxfqNKulVIQr7mBIJCBMQwjKgiVSaTQSAQgN/vRzAYhE6ng8Viwbp166DX6wEAg4ODUCqVOc+h0WjQ2toKk8mEZDIJn8+HM2fOYHFxEUajEXa7HTabTbYoON9rr0bBVRRClXSrEIxsNpsTzRICFRrBliJdlmVzotlUKoX6+nrYbDb09PRIavHVaDRobGxEY2MjWJZFOByGx+PB8PAwWJaF1WqFzWZDfX19WVEwWT83CgYgm/qiigsHVdKtoiBINEtIltt0QCI6KdfkIp1OIxgMwufzIRQKQa/Xw2KxYMOGDdDpdGW/Bi4YhoHJZILJZEJ3dzdSqRR8Ph9mZmYwOjqK2tpaGgVrtVpJ1+f+n6Rcqs0ZVXBRJd0qckBSBuFwGCqVKsdIRg6tLMuyiEQiNJpNp9Mwm81wOBzo7e1dVcMatVoNp9MJp9NJ1+XxeHDkyBFkMpmcKFgKSsnSUqkUdDodVCpVlYTXGFiWxdVXX42vf/3r+OAHPwgA+O///m/s2bMHzz//fFnXrpLuuxwkGiO5WRKVHT58GJdddpksJJhOpxEIBBAOhzEyMoLa2lpYLBb09fVJiihLQUpzA8MwMBqNMBqN6OrqQjqdhs/nw9zcHI4dO4ZYLIaZmRnYbDZJEXi+KPjYsWPo7u5GTU1NtSC3xsAwDB5//HF87GMfwzXXXIN0Oo2vfe1rZRMuUCXddx0KybmA3NxsOZEXP5rNZrMwm83Q6/Xo7e1FbW2tnC+pIlCpVGhoaEBDQwNYlsXrr7+OVCqFt99+G6lUikbBZrNZ0sZEon5CtNWC3NrD5s2bceONN+I73/kOotEobr/9dnR3d5d93SrpvguQT87FbU7IRxpio8VUKoVAIACfz4fFxUXU1NTAarVi06ZNNJodHR09L8mDRKGdnZ3o7OxEOp2G3++Hy+XC8ePHodfrYbPZYLfbqapCCMjfgDxHoYLcuz4K/uNuYOFtea/p3AJ88NGSD/vHf/xHbNu2DRqNBkNDQ7I8dZV0L0AUk3OVq5XlPkc4HKbRLACYzWY0NTVhw4YNea8vVjK2VqFSqeBwOOBwOKjiwuPxYHR0FMlkEhaLBTabDRaLpWgUTEx8+KjK0tYOampqcMstt6C2tla2VFiVdC8QkBtTqpyLj3yPJ9V+v9+PcDiM2tpaWK1WbNmyRZDm9UIhXS4YhkFNTQ1qamrQ0dGBTCYDv98Pj8eD8fFxaLVaGgUbDIac3+VGuqWeo1gUzPeJuOAgICKtJKQqdQqhSrrnKUg0G4vFVph0y/UhyWaziEQilGgBwGKxoKWlBUajUdINLjfpsiwLj8cDn8+HTCYDu92OTCZzzshHqVTCbrfDbrcDAI2Cjx07hng8nhMFF4p0i6GULI18T6fTVaPgNYoq6Z5HyNeccPr0aZjNZlitVlluMNLBFYvFMDg4CKPRCKvVipaWlrLHlcixPnKc9/l88Pl81HfX6XRCr9fD6/UiEAjgyJEjcDgcsNvtkhsf5IDBYEB7ezva29tph53X68XExASWlpYwPT0Nh8MBg8EgS3NGJpPB0NAQLr/8cgBVt7S1iCrprmEUknMB70SzREsr9WbKZrMIh8M0mlUoFLBYLNDpdNi2bZssDmAEUtMLhKy4DRRWqxV9fX0YGRmhFWW1Wg2j0YhIJIKWlhYkk0nMzs5idHS07PZfOaBUKmGz2WCz2QAAr7/+OpRKJU6cOIFYLAaz2UyjYCnvez5FBN8t7V1dkJOIhx56SNbrVUl3jYF0gHHlXMWaExQKBTWXEYpEIkFJNhqNwmQywWq1oq2tjd7sgUBA9lSAGNKNxWI0mk0mk6ivr4fdbhfcQKFSqWA2m6nki9/+S/KsUk1wpGBsbhH/7/AcYsksLu+ohwUM2tra0NbWhmw2S6PgU6dOQaVS0TUSHa8QkHlxQG4qgtucQVCNgs8NqqR7jsGNZpPJpCA5FxdCSDebzSIUCsHv9yMQCECpVMJisaCjo6PgDV2poleha2azWdoOHAgEoNVqYbVac8xtpD4Hv/2Xb4JjMploFJwvhSLH+zDtX8Ljr02hVquCRqXA82Nu9GgyuPovP1coFLBarbBarQCAeDyek4aoq6uD3W6H1WotGgWLUURUo+BzgyrprjLklnMVIt14PE7lXLFYDHV1dbBYLGhvbxd0dCXHUznBf13xeJxGs/F4HPX19bBarejq6spxDpMbfBOcxcVFeDwenDlzBgzD0EJYbW2tbOQz4YmCZVkYdcvvvcOoxfhCpuDjdTodWlpa0NLSQjckr9eL06dP0zRFvjUKLc5VZWnnDlXSXQXILefiQqFQIJ1O0xuTRLNqtRoWiwVdXV2SijQMw4hOWwgByR8HAgF6hO7p6YFerz8nNzbDMKirq0NdXR16enqQSCToET8SiVCXs3KhVymQ4Wxi8XQWOpWw10vy7BaLZfl3/7JZ8ddotVqRzWYlbVjvelnaKqJKuhVAPneuyclJ2Gw21NXVyRY5xGIxBINBhEIhzM/Po76+HhaLBZ2dnWVHinKlFxKJBI24A4EATCYTmpqa0NHRIWuRTi5otVo0NzejubmZpmXcbjei0SgGBwcl5VkB4JLWeuyd8GE6EIMCgErJ4L3N0l6/Tqeja2RZlkbBU1NTNKdLvIOlKiK4/ycpsIGBAfT391fd0srE2vvUn6co5TVLHlPOB5SMqCGm3hqNBhqNBmazGevXr5fldRBIKdABoMd1rhrCarWio6MDOp2OytvOBygUCpjNZhiNRgSDQWzZsiUnz2o2m2G322GxWEpucnqNEvdd242R2UUkUhn0OGoxOVr+oEOGYWA2m2E2m9Hb24uFhQXMzMxgcnIS4XCY5qutVqtk1QbZgEnUC1SHeJaDKulKhBA5FxdSSSzfiBpu3pNEkHJDTKRLOtV8Ph8ikUheNQRQmTzxaoKfZw0EAvB4PJiYmIBGo6G5YH7nGf19tRKXd5jp15MVWKNSqYTJZMK6devy5qulji7KZDI0sgXeHQW53/zmN3j44YdzvjcyMoLf//731O5RCqqkKwL8ApgYr1mlUplDzIUgZEQNF1LJvBSKESRxEfN6vTmdam1tbUWLT+drG3A+HTRfbbC0tASv14tjx44hkUjkRMGr6RFMyBFYma8uZ3QR97pcXMiytI9+9KP46Ec/Sr/+t3/7N/ziF7/Azp07y7pulXSLoFw5FxeFSLfcETVKpbIipMsvpBFnLeIiJtZ3gaBSMrRwOAyWZWmxSe7rlyILg8FANbdc/4UTJ05Ar9dTchMrfxOLYuqFckYXFSJdLi7kKHh8fBzf/OY3sX///rI30SrpclBJdy6lUkmvJ3REzaRvCc8MzyOSyODKznpcv9EOhQzNEUKxtLSEUChEPXEtFktRF7FSkPMmI5tALBbDwMAAjEYjtFotJiYmkEwmYTAYoFarZRn5I7bjj++/EI1G4fF4cPTo0Rwv3kpsQGIkY/lGF5EOPv7oIiGkm+85uP+XKkv7zsB3cNx/XNRzl8IGywZ85fKvCHpsKpXCxz/+cXz/+99HW1tb2c/9riddrpwrFothYWEBzc3NslkgkudIpVLw+/3wer2CRtQsLCbw6AunwALQKhn8amgOyXQWuy5y5jxOTtIlEx58Ph88Hg/C4TCamppyPHHLQbnpBXKEJ+Y2FosFGo0Gl19+ORiGgVqtBsMwOHToEDKZDN5++22k02lYrVY4HA6qHBGLctqsAeS4kJHNYn5+HtFoFIcPH6aKCDneY6mSsVKji2pqaqg0UWqkV0iW5nK5oFarYTabaSS8lvCNb3wDmzZtwi233CLL9d51pMuNZpPJ5Iojv9frlWU34xLY4uIilEolVCoVNm7cKOjmGp1bRDyVQWPdcqSmVDB4ZcInK+myLEvbbbmbQUNDA9RqNerq6mTRqBKIJV1+l5pOp6OeC+Q99Pl8KwhRpVKhubkZ69evRzqdhtfrxfT0NI4ePVqy+6zSIF68FosFkUgE3d3d8Hq9lNwIAUvdIDKZTNlSPIZZObpocnISbrcbBw4cgMFgoO+h1JMENwoOhUIwGo05bmnkPmUYBvdfdv85S0Xs3bsXzzzzDA4dKl9pQvCuIF2ho8MZhhFU7MqHQiNqnE4n1q1bh1AoBK/XKziaUSlzd/ssC+hUKyMAsaRLZGc+n48W6qxWKzZu3JhzAwWDwVVtAyZIJpM0mo3FYmV3qalUqpwIbnFxEW63G1NTUzlpgGK623Ij3XwgESkht87OTnrEJxuEFJOeciLRQlCpVFTt0N3dTdMlco0uSqfTUKlU9He5nxHyb/7fYDVIOBAI4DOf+Qx++ctfwmg0ynbdC5J0xcq5uD8TQzSFRtRs3rx5xU0iVL1AsLW1DrYRFxZCcaiUCmSyLD5x2coIXAjpcs1jSGW9VKGuErnifDcKKeiQ9RFVgNBOunx/r0J/Q241v7e3d4W/gcViyas4qATp5rsm/4i/uLgIr9eL4eFhAKAOZcXkXpUgXSBXMlZbW4va2lrZRhcR0iUotvlx/83PGcuNxx9/HG63G3fffXfO97/61a+WlWq4YEi3HDmXUEgZUUMglnRNOhUevKEXr074EUmkcUmLCX2NK3fbfOTIPZaTJgqr1Yre3t6CGlI+KiHv4orsiaMWUULYbDZZPHvFgK+7zac4sNvtFYt0i5Ejd4MgJj1er5fKvbgGONz3TGpOtxQKFdLkGF3EJ12hqHQU/NWvfhVf/epXZbkWF+ct6cop5yqGckfUEEiJHOv0auy6qKHoYwiREbtGrsFNOcdyuSPdpaUlBINBRCIRzM3NCd6sSkGuG0yhUNBIkmVZeoQeGRmhKalgMCg518qHWCLXaDRoampCU1MTWJZFKBSCx+PB1NQUXbvdbkc6na5YpFtqQ2QYaaOL5MhDr3YUXA7OG9IlifVkMoloNLriOCKn7o/YIMo1ogYQH+mWQjabpe220WgUY2Njoo7lpVBupEt8C8gkB51OB6VSCafTiY6OjrLWVmnwj9DBYBAnTpzA2bNnc6JMm80mmSzKVQHU19ejvr4evb29OSY9xH84lUpJNkPPh0wmI7poJnR0USqVkj06J59dUoxbSyR8XpDu7Owsfvazn+Hee+9FOp3G2NgYtm7dKtsbSDp1/H4/lpaWMDMzI/txVw7SJeskREuiWYPBgK1bt8qyTgIpLmNkfV6vl0bbNpuNRttzc3MV0xRXEkqlEnq9HhdddFFOlDk5OUmd0hwOB2pqagRfU86UBdek58iRI7BarQgEAtQMndueLPU5peh0+Sg0umhxcTFHOkeiYDneo2QyCYVCkXMfy0m8UgKT84J0GYbBvn378KUvfYmSVzlvXKERNW1tbUin0+ju7pZFVM+FlOM6ySGTdluyzg6e+fjU1JSsayXr5QrZC62P2w5Mevs7OzvzqgEuhDZgfpQZj8dpHjgWixUsxvFRqYIX6cpraWkBgBXrE2PSw4UcpMsFd3SR3+/Hpk2bVryPHo8HNputrPepEvl47rV9Pp9orjgvSNdsNiMUCpX15gkZUQMsFwbkTAMQCF07aaLw+XwIh8N0MGRra+uqFpkKEWS+IpjVasVFF10kaH2VIl2S46+EXWSxG1en06G1tRWtra20GOd2u3HixAmqZ83X+FApMuCO68m3Pq5Jj1arpWmSUgVWuUmXD71enzO6yOv1wuVyYX5+HgBy2ofFIJlM5sjRAHkjXVKMFYPzgnR1Oh2SyaSo35EyogZYJt1SEZ6cINEiSRsA8uSQywU3Mo/FYlQ7S3KFjY2NWL9+vaibQO7XQtQqY2NjCIfDUCqV0Gq1SCQSSCQSqK2tleV5hBJkoWIct/HB4XDAZDJVLNItdt18Jj0kx5pIJGiUnk9vW2nS5a+TKCKAd0YXeTweUaOLAODQoUPYuHEj9Ho9stksNBrNqr2OQjgvSJeA++HPdyOUO6IGkL/glQ9ymcdwIWfklM1mEY1GaX6WSM743hBiIUd6gWwAXq8XmUwG2WyWblAajQaxWAyHDx/GxMQExsfHc4hO6vsj5b3lF+NSqVSO5Eun04FhGMlyqUIQQ+b8HCtXb2swGGiOVafTVYx0hXhMSx1dBCzfa5XK50rFeUG6fKIlxKhQKGQdUQNUJtIl+sVkMonh4WFkMhnZJFPAO1FpOTcFt0i3tLQErVYLrVaLvr4+2W42KaRLCldEBaHRaGCz2WgH3eDgIEwmE328wWBATU0Nenp6oNVqc4iOTBQWEiHJDbVanePwNTU1BZfLhcHBwZxil5hiXD5IjaC5SgMSpXu9XupfkUwmEQ6HodVqZY3QxW46/NFFhcYrEf1yOp0+55EtH+cF6QLLOZ9YLEajA9KCKOeIGkC+SJdUZ7lOYgzDYMOGDbLb+0khXX5ag2GYnPQLiSYqJeUphnQ6TddFph/YbDbBf2PyHPwW4GAwCLfbjVOnTtF8psPhKBm9y51/JXpW0rDCLXbF43FqhiSlrZY0BJW7PhKlE5OeAwcOwOVyYXx8fIUDWTkoN9LnKjf4o4sYhkEikaATn8lrO9dY86Qbj8exb98+BINBXH311fj7v/97bNiwAS0tLVT/JyfKiXSXlpYoWaRSqRVOYocOHapIMUyoMoIr0wmFQkWLYJXqSCsE8t6RtIHYvLb2lX+EQqFA5rpvF3xuMtaGPJ/b7aaRXLE0RKU70rjFrnzH/ELFuNWCSqWCUqnE5s2bASDHgSybzdJpFFIaR+RMr3D/xr29vUgmkzhw4ACmpqYQDofR3t6O7u5uWZ6rHKxp0k2lUrjmmmtwxRVXoKmpCQ888AAuvfRSTExMVOx4qFQqkUgkBD22kHnM+vXr80azJIqWe+3FSJfrIkY2AmLCUywiqtQIdm7rZrG0gVgo3KNQMAoIPaMYDAZ0dHTQSK5YGmI124D5x3w+wZH8ZTk5aqkgz8d1IOOb9IidyVYpxQmwnNJRq9W4+OKL6aDYtQBRr/bOO+/E7373OzgcDhw9enTFz1mWxZe+9CX84Q9/gMFgwBNPPIFt27YBAJ588kl861vfAgA88MADuOOOO0o+n1qtxoEDBwAA99xzD+Lx+PKiK6gwKHXtWCxGi2DxeJy6YHV3d5c8+laqSMclXaLaIHaIpAhWaCMoBCnNEaWQzWYRiUQwOjpKZ6mJSRtUCqXSEHq9XvYNSAiR8y0WSTGORG7nMkdNkM+khzuTjWwShU4schcS+dcmJziGYaDRaM6/9MKnP/1p3Hvvvbj99tvz/vyPf/wjJiYmMDExgbfeegt333033nrrLfj9fjz88MMYGhoCwzDo7+/Hrl276FFPCMxmM4LB4PKiVaqK7Vp8YuSTmFqthtVqRU9Pj2DzGAKFQlExZQQZE0461colM7nSC9y0QSKRgEqlwrp1686pHK4Y8qUhJicn4fV6ceDAAVnUEIC03Cu/GBcMBuHxeHD69GlajDuXHX9ckx4yk83r9dLJxKTIxW2frmShK5VKnbPNqBhEreh973tf0e6nZ599FrfffjsYhsGVV16JYDCI+fl57N27F9dddx2tOF533XV4/vnncdtttwl+7vr6eoRCIQDLHz6xul2hUKlUSCQSmJubW2EeU25EJtc8M251maxRqVSivb296GBIMZBqeEM8Ifhpg76+PkSjUfj9/hy1wVoHyalqNBp0dnbKpoYQIpUqBv7mEIvF4PF4EI/HsX///qKa29UC36SHFLkmJydpGiWdTkuWSZZCKpVac3IxQOac7uzsLFpbW+nXLS0tmJ2dLfh9MbBYLHC73QCWiXFpaUmeRQM53qXcUedymccQlJNe4KshDAYD1faeOXMGNptNVqNlMZEuURt4vV5EIpGCm9TS0tJ53QZcLA2h0WjgcDgE+8hms1lZozC9Xo/W1lbMzMzgiiuuyCnG1dTU0DyxFIIrd4MAVha5SMPD7Ows0uk0IpEIbU+W630R6tO72lh7sXcB1NfX4+TJkwDkyekmk0mamyX5RavViubmZoyNjckysocPsaRLPpjEOarQXLVKGI6XKqTlUxu0traWTBucz6TLRSE1xNGjR6kaolhFXw5pFx9ENlioGEem/pKfCU3xVKIxgjQ8xGIxmEwmqNVqeDwenDx5kp6OytUtcyPdSnowiIWspNvc3Izp6Wn69czMDNXQ7d27N+f727dvF3Vti8WSk9MVS7rcCQXEspF4L3CP5KSHvxIoRbpcu0a/3w+1Wg2bzYZ169YVjZ4qNeWBe81iaQOhUqZKTF/IZDKYmppCazgCFizOHj+OZDK56uSeTw1BrCDzpSEq0QbM910AVhbj8uVZycy2QhFmJVuAiXqB2/BAUiV8EyEyuFIo+N1oawWyku6uXbvw2GOP4dZbb8Vbb72Furo6NDY2YufOnfja176GQCAAAHjhhRfwyCOPiLp2fX09JV3SaVIKhcxjilk2VtIJK18hjUhuSMRNjubnunWZkK7L5cpJG5RToJPjveX79CaTSeh0OtTU1lBvAZfLheHhYdTV1cHhcJTlewuIj5KEpCGSyaRs3hAEQoicm2cl76Xb7aYRJtcGkqCSpJtPvcA1vyFpNa4ROlljqTROKpU6Z9rmYhD1Sbztttuwd+9eeL1etLS04OGHH6Yqgi984Qu44YYb8Ic//IFW9v/jP/4DwHKU+o1vfAOXXXYZAODBBx+ku5pQ8CPdfOoFboHJ7/cjm82uCfMYAqVSiVQqldMJRqz4+BG3GMgZ6ZK0ARnBHo1GBaUNhEAq6fJzxqQK3t3djYMHD8LpdC4bVYOhR/re3l6k02m43W5MTk5CrVaLyrlyUc7RNF8awuPxwO12w+PxYHFxsazpv1yIjZ4VCkXeYhwxwCFNDyRlUQmUkoxxLSAB0PubjAMqNhQzlUrlbGzn+v4nEEW6v/rVr4r+nGEY/PjHP877szvvvBN33nmnmKfLgcVioeoFLslwR52HQqGiwyHPFchu7XK5EIlEsLi4KOsayyHdQmmDjRs3YnR0FF1dXWWvj0AM6ZKmDo/Hg0wmkzcVVOq5uEMoY7HYipxrQ0ODoM1EznwgMZkhGm+WZYumIcSgXP8NfoTp8/kwPz9P3e/m5uZETSYWArE6XTIOqL29fcVQTP5o+GohrUzU1NTQ6jfXPCadTsNisQjqshIKQg7l/JHi8TiNzpLJJOrr62E2m1FbW4uenp6y18iFEMNxLvKpDfhpAzLYc7XANWz3+Xw0ny2lQy3fuvV6PXXUIk0GJK9JcoZWqzXv56cS70M2m6UGTXKoIYD8Od1CWEpm8MIxN6KJNK7oMKPHkZvqUCqV1F7R4/FgYWEBsVgsZzKxmGJcsTVLTf3wh2JyR8On02lkMhlqOL+WcF6QbiwWwyuvvIJYLIb+/n785Cc/AQBRRRwxUCqVopPwXNkZKYLxJ/CS3JTcEBLpilUbVCIq4Ee6pI3a4/EgFArBaDTCZrOtMJaXG9wmA2Ls7Xa7Mf3TVP0AACAASURBVD4+DoPBQMmOG9FVogjIvWahNIRQNQQgPL0QTaRx+xMHMRuMI5tloVIy+P5NW3BVV/6UXzabhV6vR3d3d85k4tOnTyMajeZE6GIjbTEbRTHwrTTT6TQGBgZo84jdbqfeEecaa550jx8/jk984hO45pproFQq8dprr0Gv12NoaKhilUkyPaLU9UmhjkSMpK21UBFMruYIPgqNYS9XbSA3GIZBJpPB/Pw8naNmNpvR0NAg2ylFLLjG3kRe5Xa7cejQISgUCtjtdmQyGdlTVaUIkut1K0QNQa4phPT+cNSFmUAMapUCUDJIprN49E/jePbuK/M+nl9I4xfjyMZJinFiI/RKbPBkWsTmzZuhUCiohcBagGjSff755/GlL30JmUwGn/vc57B79+6cn99333145ZVXALyjXSQFMKVSiS1btgAA2tra8Nxzz5V8vvXr1+PgwYMAgFdffZV+UIlsrBJ5WxLp8kFSG+QITIpgQgtNlfZeEJI2OBfgFo/i8ThMJpPsjSdygCuv6u7uRiKRgNvtpgL+WCwGh8OB+vr6stctJn0ltClDaNS4GEshnWVBQgqlgsFirHB6qlirLt/flvytuYUuu92O+vr6Vd9UuZuQVqtdM581UaSbyWTwxS9+EX/+85/R0tKCyy67DLt27UJfXx99zA9+8AP673/+53+mOSBgOa92+PBhUQvkvlFENtbQ0FBR0uXOSeM7ien1ethsNklFsEqQLtkEPB4PgsGgqE2gUuA6iPn9fuh0OthsNvT29mJmZgbt7e3nZF1iodVq6WwxhmGg1WoxMzOD0dFRKkcrp+glhYSKpSFisRh0Oh2CwWDRNMRlHWZo3jiDVCYLpYJBJsvivT2F1UQk/ywE/Ajd7/djbm4OY2NjMBqNtNC1Vorc5wKiPi0DAwPo6emhFe1bb70Vzz77bA7pcvGrX/0KDz/8cPmr/AvykW4lwLIsXC4Xzpw5Q02lbTYbenp6ytqt5SDdfGmD2tpa1NfXY9OmTWVduxwQFQmRmuWLsmOx2HnbkaZSqdDQ0ICGhoa8TmQk2hRa9JNLEcEluenpaQQCAZqGKLQxXNRSh4c+tAHfe/EklpIZXLvBhq9+YH3B58hkMpLsNvmFrnA4TDvjANBZcpXuFltrnzlRpJvPQ+Gtt97K+9gzZ85gcnISO3bsoN+Lx+O49NJLoVKpsHv3bnzkIx8RtVhug4ScTmOkCEY6wVKpFHVKEuskVgxSpV2l0gaRSARnz56VbZ1CQUalEIWGxWJBc3NzUQeutXYDSAE/2oxGo3C73RgZGUE2m6UTKYrJ2yrRkcYwy+Ph29raclzI8qkhdm5qwM5NDYKuK4cTGMMwMJlMMJlMtBg3Pz+PRCKB/fv353TGyZEK46daGIY5P9MLYvDUU0/hpptuynkDz5w5g+bmZpw+fRo7duzAli1bRDm5m81m2tUmtCutEPjdatyR7PPz81AqlbISLiCOdLm541Jqg0q0AecDV5bj8/mocL1UmzJBJT702WwW09PTaIpEABaYmZhAKpWSldxLRWI1NTXo7OxEZ2cnreyfOnUK0WgUFosl7+idSkR33IIXd2NYt26dJDUE97pyq0lI95vf78fFF1+cdzS8lEYWgkr69JYLUasq5K2QD0899dSKRgny2K6uLmzfvh3Dw8OiSZc0SIhNL5AiGIkYi3WrVSp1UezDXY7aoJKky52+GgwG6ZTYiy++WLR6RI42YK6BCzGSVygUqDEstwEbjUbMzs7i0KFDsFqtNHoqJ6oUQ5D8yj4R7x87dgxGo5G2JVci0i2Wey2mhiiVn67UuHhCjIVGw4+OjiKVSuUU44T+HdaqrSMgknQvu+wyTExMYHJyEs3NzXjqqafwy1/+csXjjh8/jkAggKuuuop+LxAIwGAw0Amtb7zxBu6//35Ri+WnF0rJQIgGk1sEs1qt2LRpU1EiUyqVqyIxkUttIDfpEj+IWCyGoaEh2nZbbk5bKunyib+mpobqeYeHh9Hc3AxGsXx8dDqdmJ+fx7p16xCPx+FyuXDixAnU1tbS47XYCEhqVKpQKGgLK8lput1uTE1NIRqNYmZmBo2NjbINKhVKjvnUEIXSEEDljMYLRaP8DcLn82F2dhajo6O0GGe324tu+mvV7AYQSboqlQqPPfYYdu7ciUwmgzvvvBObNm3Cgw8+iEsvvRS7du0CsBzl3nrrrTkf1GPHjuHzn/88JYjdu3cXLMAVgsVioT68haLRRCJBiYy0WoolDK56QW5ks1mcPXuWpg2sVmvZagM5SDcWi8Hr9WLO5cGpQArW+jrUqtS49NJLZYtyxHr0+v1+eDwe6rdgt9sF/x35+ttwOAyXy4WpqSnqwyBkEjB37eWAm9Ps6enB/v37oVAoaDQnx0QKKW3ApdIQVqsV8Xh8VUmXC34Bk4wDOnjwIBiGoQTMz59zp0ZUwkazHIhOetxwww244YYbcr73zW9+M+frhx56aMXvvec978Hbb78t9ulywE0vqNVqmrsj7aN+vx9KpZLOLJOqA5VT2sUf90M+wHI2KUghXX4HnUajAWOoxwP7kwgnMmBZL5wG4Jfb0qjRrY68h1+Yk2NDAnIJr7e3l+rHuYWvhoYG1NTUFPS+lRsMw6CtrS3vYEzimyy2qCRHd1e+NMTMzAwOHjxYtkwu33rFXIc/Doh8Xk6ePImlpSWYzWZqhH7BpBfONQjpxmIxhEIh+P1+DAwMyN4+Wm5Ol9+pxk0bDA8Po7GxUdadV+jkXv4IdvK+kQ66//v0UXgiSWT+cqnpMPCzAzP4v9fIY3qTL9IlrlFer5cOMhRamCuEUu8F1/s2mUzSAs7S0hKsVisaGhpy8oeVKHpxr8k97pNUSqm25HyQO/dK1jU5OYnLL7+cRplSvSH4KDdtodVqqV83SSWSv2U2m4XRaEQ8Hl9zmuDzgnRZlsX4+Dh+85vf4I033sBNN92EH/7wh9BoNOjv75f9hpAS6fLVBoWiNHLt1TrukGo6SbcQg5X169eveN8mfTFKuACQygInvfKNRSIevaFQCB6PhzZO2O12bNmy5ZzcHBqNht64xFmL5A+JjCmTyVTMi4IPbocXUYu43W6qbSVEl8+Lt1IFL3I8z5eG4Ba7HA6HKIvKdDot62mPW4w7duwY0uk03n77bWQyGbznPe+R5XnkwHlBus888wx+8Ytf4K//+q/R0NCA3//+98hmszh06FDF+rZLRbr8tIFWq4XVai2ZNqiU/wIBt1WZRI8k3VJq9MmmxlrMheJI/YV5NQpgs1P6uBQCrhF1JBKhFoHnuj2ZD66zFrcBYm5uDsFgEMlkUvKcMSngmrh0dXUhkUhQM+94PE6JjkTl5Vo7lloLF+WoIQjS6XRZ43iKQaFQoKmpCVarFalU6vzO6QKl/ReeeOIJfPnLX6YSsXvvvRef+9znAABPPvkkvvWtbwEAHnjgAdxxxx0ln++mm27CTTfdhHQ6jSeffBKA8CO1FBSKdIulDYR+2Csxhp34LkxMTCAQCNBWZbHR49d29uC0dwlT/hhYFthkVeC2fqekNREFhMfjyTG2iUQi2Lhxo6Rr5kOlPgPcAhP3qEqMcAg5y63lLgatVouWlha0tLSsiMpNJhOi0eg5aT4Rq4YgqKSWlltIW0sbOyCBdIX4LwDALbfcgsceeyzne36/Hw8//DCGhobAMAz6+/uxa9cu2tlTCkqlclXaBrm5x3xpg3KmPMhVpCPVfa/Xi3A4jHQ6DbPZjK6uLskfsjq9Gr/+7DbMBuNQKxXwnBmHAsJvYqKAIPaRZEMiBarztRuNYRjo9Xq0tLSgu7sb8XgcbrcbY2NjsikPxIIflYdCIYyMjGBkZISmbMSoM+RCITVEvjREJZouCC6oQppY/wUu/vSnP+G6666jjkTXXXcdnn/+edx2222CnrtQZVnON5SkDRKJBAYGBqDT6QSlDYSiHNIl04G9Xi81byfNHUNDQ3SkSTlQMAxazcsRia/EaYI0KpBNSaVSrbp95GrcTPzPmE6noxMWyNF6amoKkUiERvT5xsdUCqT9V6/X4+KLL6ZjioiZNyFgKSqQcjfKfGmI6elpHD16FJlMBkajEQaDQXbyvWA60gDh/gvPPPMM9u3bh3Xr1uEHP/gBWltb8/4u0d0KBTmek04WOfJY+dIGKpUK/f39sh9NxJAun9RI2+369evzVozl3oDyRafcQhjpULPb7WhtbS35Ia+EAiCdTuP06dNoC4cBAJOjo0gkErLmzYu9r3zlAXd8TG1tLRoaGsoejCkU5F7QaDRUnZFKpaiRdzQazZGjCdkU5CzO8dMQb775JqLRKAYGBmhbsMPhkKVZ5IKKdIXgxhtvxG233QatVot//dd/xR133IGXX35ZlmubTCaEQiGapJcqOymVNhgcHKzIH6qUppZIX0j3VW1tLWw2W0lSk2sDyndNfipDrg41KeCTPhk+WFtbC4VCAafTCa/Xi0OHDsFkMslCekI3M34H2uLiYs5gzIaGBlFOZGKRTxWjVqtXtCW73W6cOHECNTU1NN9aqHurUpOAiQHN+vXroVQqZVFD5HuOtQjRn0Qh/gtEtgEAn/vc52i7b3NzM/bu3Zvzu9u3bxf1/ESra7VaaYOEkKMsV21ApErF0gakK01uUskX6fKLTkSq1NvbK/j55SbdZDKJpaUlqnk8l1OVM5kMJX0yOYGQ/sGDB+FwOOg0YKvVCqPRiI0bNyKZTMLlcuH06dPQarVoaGiAw+EQrT6QcoLgD8bkN2QkEglEIhHZx7AXWyd/UyBTMg4ePAilUkkjzdUav86NooulIcppyliLdQTRpCvEf2F+fh6NjY0AgOeee45Wq3fu3Imvfe1r1CnshRdewCOPPCLq+fn+C8WkXeWoDaTMSRMCQrpcWRfLsrBarTlFJ7GQoxWYRBtEagYsb5ROpzQFQzngb0QWiwVNTU3YsGGD6GnA69atQzQahcvlwvDwMBiGgcPhQENDg2y+B6XAbchIJBI4cOAAJiYmEIvFVki/VgP8KRmkOEjGr5PiYCXHr5N18FFIDXH69Gmo1eqSaQg+0a4lW0dAAukK8V/40Y9+hOeeew4qlQoWiwVPPPEEgGXvhG984xu47LLLAAAPPvggLaoJRX19PSXtfKQrl9pAbv8FbtstydPabLaS5jtCIUcrsFarzZGakZlXqwWiQyWFQqkbUb7opqamBl1dXejq6qIEQ/wFSBtwoc+InLlyb8yLeDoOrVaLrVu3UukXmUhBTjlShjyWA35x0OfzUVN0hmHgdrtXfU1A4UkZxdIQlYzO5YCkRFcp/4VHHnmkYAR755134s4775TytACWiZsb6SaTSeokJiRtIBRySLv4x2KTyUTzj2IsLYVAKOly2yVJK7DdbkdHR8eKD+pqyLz40bXdbi9YKJQLXIIhhSbSv08iPDnbgBeiC3h55mW8NP0SRrwj+FDHh/ABxQcA5G/IcLlcmJiYENUCLCe4JjM+nw9nz55FIBDAyZMnodPp6JrKub+kfq6EpCFqamrWrMMYcJ50pHFRX1+P+fl5DAwM0CmfJE8lZ4eTVP8FrstZIpGA2WzOORYTm0m5UYx0yVHd6/VSYxAhE3gr4dPLdYrKF12vNriFJn7USSr9ZEaaGEyHp/HSzEt4efpljPnHAAC99b24a/Nd2N64HYFTgRW/w43qSAuwy+XKachYzZQIsLxJ19TUYN26dQBAvYzJrENy1Bd7GpFDFcFPQ4RCIbjdbkz8xcj+zJkzsNvtq9rAIgTnDemePn0av/3tb7Fnzx4kk0kwDIMbb7wR6XSaaoblhNBIl9wcJKVBTFsKjfqp1ERgfnsx0fR6PB6aZuno6BB1c8gV6XIVB0tLS5ienobdbi84qv5cgRt1khOB2+2mzR7Nzc1FlRCnQ6fx8vTLeGnmJUwEJwAAfZY+3HvxvdjRsgNtxjYAyxtzSBEquhZuCzA358q1gmxoaKh4YZN/VCdrIlMy+GZB5JRQilDl1tESrTIpshIp6ujo6KrqxoVA8qsu1Qr8T//0T/jpT38KlUoFu92On/3sZ3QKrJRR7Pv27YPRaMRXvvIVDA8P44tf/CKVVlUCxSJdvmxJr9fDbrfjoosuKnmsqRTpMgyDSCRCRxCRRoWNGzdKliiVE+mS1AoZVEluhmAweE4HaAoF10AlHo+jsbERi4uLOUoIu92OyaVJvDz9Ml6eeRlTi1NgwOAi20W4b+t92NGyA401jSuuLSVdwc+5ejweTE5OIhKJwGKxyNIYkw/F8qN8syC/34/5+XmMjY3BZDLRKRn5yLXS3Wh6vR7t7e2C9OOrDUmrEdIKvHXrVgwNDcFgMOBf/uVfcP/99+PXv/41AGmj2D/96U8DAPbv349XXnllefEVnAhM1AsE/CkPUrWqlfDq9Xq9WFhYgMFgQHNzs6RROvkgNtLNpzjgD6o8ffp02es6FzCZTHA6nejp7cHB2YPYc3IPDrx5AN60FwoocIntEtzSfwuuabkGNn1xAiz3aK1SqdDY2IjGxkaqvV1YWEA0GsXIyEhRshMLoUUpIjmz2+15NcqrNY2CXHutES0XklYmpBX4mmuuof++8sor8fOf/7zMpS6jnDlpYqBSqRAOhzE9PU2Pl+X6LgDl50kLRZAqlYqK3eWCENKVS3FQLliWRRaVcW/LZDM44juCN46/gVemX4Er5oJKocLlDZfjsw2fxQb1BsT9caS9aQSZILQN2qKfETnVEER7azKZsLS0hPb2drjdbhqRk3SJ1ON1JpMR/bt8jXIsFluhFtFqtRUj3bXcjQZIJF0xo9gBYM+ePfjgBz9Ivy5nFLvFYlkxPUIucKdQuN1uZLNZtLW1lXVE50NKpMv1xE0kEnkbFaLRqOxFr0IbxLlQHORDMplEKpXCkSNHsCkSAcMwOPLmm0ilUlhaWirrb5bOpjHkHsLL0y/jxakXsTixCI1Cg6sar8I9Lffgfc3vg1FjfOcXuiFICQFUxvc2m81CpVKtaMhwuVw4cuQIWJbNGQ0vFHLIr8hRv729HalUitpARqNRAMvFOKvVKtt7kkqlcuopFwTpisHPf/5zDA0N4dVXX6XfK2cUOzfSlUvWRcTX3GkKZrMZLpcLLS0tZV2fD6FrJsTm8/kAoGhxDqiM0oBEumQzWguKAxJZezweZLNZsCyL3t5eGI8boWAU2LZtGwYHB3Hq1CkcP34cdrsdTqdT0OkkmUliwDWAl6Zfwr7ZfQglQ9Cr9OjT92HXxl24pu0aGNSFK+FClBBkZpvcRJCPyA0GQ85oeDKNgnjxNjQ0lGyzlVvzqlar0djYCJZlEY/HUVdXlzMlg7Rtl/O5uiDTC0JHsb/44ov49re/jVdffTXniFLOKHa1Wk1TClKr68lkkuZnuV6vXAlVNBqtSMGrkA9wIWLbvHmzoA+g3KSbzWaxtLSEYDCI+fl51NTUnDPFATeyVigUsNvttCI9ODiYsxFpNBrodDps2rQJKpWK+rpGo1Fa8ecSTTwdx/75/Xh55mW8NvcaoqkoatW1uLrpalzbei2udF6Jo4ePYkvLFmjVwo/ZhZQQ4+Pj9O8pJzmUalnXaDQ5XrxcfWt9fT0aGhryzmSrVKMBee3c4aGkLZlI5LhyNDG4INMLQlqBh4eH8fnPfx7PP/98Tp5RjlHsgPi8GDdyZFkWNpsNXV1dBYdXVjJfTMA3tynWqFAKCoWi7PXy88VarRZ6vR4bNmxYVWMbIsMjREscqMRE1gzD5BScSOR59uxZuINunFWdxdvxtzHoG0Q8E0edpg7XtlyLa1uvxeUNl0OtfOemLTcy5U8mnp2dxdmzZzEwMFCWJwQXYlIWSqUyZ8Iud0Mg0SYxwakU6fJzxfnakj0eD44fP45EIiGqVXotj18HJJKukFbgL3/5y4hEIvjYxz4G4B1pWLmj2EkfNbkRuP/mgoilSYurTqcTFTlWStqVSqWQSqVw9OhR2qgg1twmH6RGusUUB0SpsRqESyJ9t9sNv98PvV4Ph8OBrVu3yhINRjNRDEQH8FLwJby18BaS2STqVHXo1/Xjr+x/hfd1vg8OuyPva5UzHUAM0S0WCzZs2CCbJ4RUsyOGYXJmsvFNcEjeXG6UGtWj0+nQ2tqK1tZWGpmTVE0pAxzu1Ii1CMkrK9UK/OKLL+b9PTlGsdfW1iISicBkMtGIlKQduBaEUsbpEMhJuvxGhWw2K7pRoRTEkK5QxUEl8sRccM1MAoEAjfTL7Swk6Rt/3I+9M3vx8szLGHQNIsNm0GBowN/2/C2ubb0WW6xboGAUtPX21MlTeT1w5c7Bcq9XjicEF3IU5/JFm0NDQzh58iTGx8fLMkPnQ0xqhR+Zk9l1p06dyqvQ4EfnF0R64Vyjrq4OgUCA6j9nZ2cRCoWQTCbzakOloJxuLO7xmJiP2+12qoIYHByUXVJViiClKA4q4b3AsizdgLhDDOXy5nUvufFK8BX862v/iiO+I8iyWbTUtuAT6z+BHa07sMmyacX7zm29DYfD1A5Sp9OhoaFB9o2nEEGK9YTgohI2pDqdDhqNBlu3bgUAeL1eaoZusVjgcDgkT8iQms/mG+CQaclHjhxBNpuF3W6nwdJatHUEzkPSZVkWKpUK3/ve99Df34+enh7aG77akiUuuI0KJG9ts9nyNipUYnIrvw1YDsWBXKRL8qmkDdjv96OxsVGwTaMQ/DK5iP/VxHH8ueXTV4exA3f23YkdLTvQW98r2A7SZDLBZDKht7cXkUgELpcLkUgEw8PDNNoqt6VUSOQsVAlBCK9S49dJ1KhQKHIaMnw+HxYWFnDs2LGS3Wf5IFcRsaamJkehQWSV+/fvh9lsRktLyznlhXw4r0j3gQcewP/+7/8iHo/jxhtvxIc+9CEEAgFYrdZz8sbyC091dXWw2+3o7u4uegMQgpR7ygN31HkwGCxbcVBOeoGbK47H47BYLGhra8PS0lJFJk4Mq2vhy6SRCFyD2tQWXLJhKy6ra0RPGZMHiM+A2+1GX19fjtELIWApnzuxBFlMCUHSIalUqiIFL5ZlV6yVKAu43Wd8s/hSLmSVaAPWaDRwOp04c+YMrrjiClo0X2so61WX8l9IJBK4/fbbcfDgQVitVvz6179GR0cHgGX7xz179kCpVOJHP/oRdu7cWfL5br75Zjz44IP41re+hfXr18Nms9FJuJUAIR3uhy5fo4LYdAbJF8tRYeX2vAeDQdqhJAexiY10iQGKx+NBOp2mChFuwaSQZK5cPHjz04hnVHhl3Iun94/jqYNz+K+BWdiNGly3wYHr+xy4rL0eKqX494Rl2RUm5C6Xi+ZeSfFLaNNBOTlivhKCEN7c3BxUKhVUKlXZSggxyGcW73a76eZENgt+0axSbcAkgib3wVqEZNIV4r+wZ88emM1mnDx5Ek899RS+8pWv4Ne//jXGxsbw1FNPYXR0FHNzc3j/+9+P8fHxkn+Eiy66CEDu9AiublduEHLkTuEFlscRFWtUEHpdqcinOHA6nXTmlFwQEukSaY/H4wGAooMzgcrkiYkvBsuy+JuNDrRnatHRuxEHzoTxwpgL/+/wHH45OAOzQY33b7Dj+j4Hruy0QKMSTsBcktRqtTm5V27TASl+FSs2yZUK4BKeVqtFOp1GMpk8Z9MxgNzjPinYnjhxYkVDRqXSIWtdowuUQbpC/BeeffZZPPTQQwCAm266Cffeey9YlsWzzz6LW2+9FVqtFp2dnejp6cHAwACuuuoqQc9tsVgwMzOz/AL+YmQuJ0g+NJFI4NChQ9DpdLDb7YLlZqVAUgFiUEpxkEgksLCwUPbauChEkFwNLTn6yjUBQyiIgXUsFsPw8DCsVit0Oh1GRkawuLiIoGce7+9txo0XObGUzGDfhBcvjLnx+6Mu/PehORh1KuxYb8P1fQ68t9sKnVpa1KVWq6nTFlkTcf8iJJOvDbgSHWl6vR5NTU1lKyHkglarpQ0ZZGM8e/YsFhcXkUgk4PV6ZW3/BVbmii8o0hXiv8B9DOkL9/l8mJ2dxZVXXpnzu2JGsZvNZoyOjtLrkh7ucsBtVAiFQqitrYVGo0F3dzfq6urKvj4XYluBhSgOpBB5KXDbgIl5tc/ng1arhd1ul+RmVk6kSyr6Ho+HjmrRaDR0/JNarUZHRwcGBgagVCpx9OhRZDIZOBwOvK/TiQ9sakAilcEbp/14YcyNl4578OyRBRg0Svx1rw07Nznwvh4rarTSbgv+OHbyWSfFL9L1xbKs7Edrfo2gmBKi0GbAh5wnEu40imw2i9dffx1erxfj4+OCphILBTfSrUS7tRw4rwppBFz/hXLSC2T3JVV1/hTe8fHxiuQf+UoDgnIUB3JragnRhsNh2mrrcDjQ1tZWVgFEap6YGBDZbDb09vbS1E4+P2WFQoHm5uYcz4Fjx44hmUzC4XDg8uYG7Fi/CalMFm9NBvCnMTdePO7GH0dd0KoUuLrHiuv7HNix3g6jTtpr5RabyIbucrlw/PhxKBQKmM1mWbu9iknGxCohCCqVAiDdghs3bqSfM5fLRRsySB5YSlpkrXejAWWQrhD/BfIYcrwgo9OFejcUgpiJwHzkO6a3t7fnPXZVcsoDuS7XED0QCKC2tlaS4kAO0s1ms7RZgagfVCoVtm3bJhs5CLWLdLvdOXniDRs2CL4Judfneg6QiI+bf73Y6cRfdW/AQx/agKEzAbww5sYLxzx48bgHaiWDq7os6NaksSGahKVGWmqJX/waGxvD0tIS3nzzTdTU1NBqfzmbmVCCLKSEOHHiBIxGY05jSKVagLlr5TZk9PT05LWBFJMWWevdaEAZpCvEf2HXrl148skncdVVV+F//ud/sGPHDjAMg127duHjH/84/v7v/x5zc3OYmJjA5ZdfLvi5+cMpi5Eud5wOyUGWKvYQ8I3M5QLDMAgGg7RzTqohOv+aUkBkZm63m67Fbrejt7cXmUwGR48eXZXJhxrXpwAAIABJREFUqrFYjKYOSJQod56YG/GR6QtcM5wNTicuv2E9vv7B9TgyG8ILY278acyNfcEk/vPYa7i8w4zrNzpw3UY77EZp62IYhg5PbWhoQCQSwcLCAqampqDRaCT7MEiJSvMpIbhevGazuSLH82IaXb4NJPdvZLFYaFqk0GslUyO4r3GtQTLpCvFf+OxnP4tPfepT6OnpgcViwVNPPQUA2LRpE26++Wb09fVBpVLhxz/+sagbuxTpEt8FEj2ScTpic5ByjmHnKg4WFxdhMBjQ1dVVduecFPDTKoVkb8Q6UU5wI92lpSU6g4yMdZJiF0nWqVAoBH+O+GY4Xq8XU1NTtADW0dCA+6/vxf3X9+JXz7+BBbUTfxpz4+HfH8c3/3Ac21rrcX2fA9dvdKCpXpxvLyFIbpTX29ubM4iStL46HA5BvsDl6r75xuPRaBTT09MIhUIYGBiQVQkhtDGCu0mSHDl3HFBDQ8MK/4W1busIAEyJm2rtKYuxTKpbt27Fa6+9BpZlMTg4iP7+/pxx58R3wWw2S/4wzs/PI5VKoa2tTdLv50tl2O12hMPhsq5bCIODg7SoxAcxj/Z4PNS1yW63Fz22ZTIZHD58GP39/bKtcXh4GLW1tQgGg9BoNLSTSWwejhDt0NAQent7UVtbS5UsIyMjuPTSS6HVakVvaNlsFl6vFy6XC6FQCBaLBV6vF1dffTUAYMIdXY6Aj7kx7ooAALY0m3D9Rgd29jnQbi0tIzxx4gQsFgvsdnvBx5BjtsvlAsuylPQKyRSPHDmC7u5uUQblpRAKhTA9PY2enh643W643W6kUqmylRDkups3b5a0Lu7kX6/Xm+PUNjExgebmZtTX1yObzUKj0azKSS0PCr4xa3tLKIH5+XmwLItYLIZDhw6tGHdeLlQqFeLxuKjfEaI4iMVioq8rBXzD71JG6HzIkSfmKx8I4UtxDyNEy11Tb28vZmZmsLi4COAddypidUkiSvJfKZBR5yTnSRpP3njjDapAuOevO3DvNV2Y9L5DwN9/8SS+/+JJbHDWUgLuceQnwHxdXnxwj9n8YmA+0qtE0Ys0MMihhOCi3Fwxd/IvtyFjeHgYkUgESqUSarUaer2+ql6QA2QU+8LCAm666Sb85Cc/gVarLRjhlQMhOV0pioNKFeiA3NwoIX2p44akfmD5No1c5cPx48dF9egTos1kMnQ9SqUypymDSPuUSiV8Ph+mp6cRjUbhcDhyUiZSCNhms0Gv1+PKK6/MUSDU1dXB6XTi/3tvOz7/vk7MBGL487HlHPCPXjmNH71yGl02A67vc+ADfQ3Y4JROkPmKgSdPnkQsFoPVaoXT6ayI4U2+lIVUJQQXcqcAuA0ZxKP4+PHj0Ov1uPjii2V7Hrkg2yv3+/245ZZbMDU1hY6ODjz99NPUCYjg8OHDuPvuu7G4uAilUomvf/3ruOWWWwAsT/t99dVXqSb2iSeewCWXXJLz+6FQCH/3d3+HD3/4w2hqasIf//hHGAwGDA0NVWSnL5TTLVdxICfpsixLc6PRaJRa8K32KB1+Hp1oL/k2jULUC3yiJRIjIvgnxuYNDQ3Ytm1bTmqCOyF3bm4Ox48fp2TAnRghhYC5RadgMIiFhQWMj4/Tqv/tV7TgM+9ph2sxgRePLUfA//baFB7fN4VWsx7X9y1HwMoymiP4pOf1enHmzBkEAgGcPHkSTU1NshXASkWkYpQQXFQy75rJZNDe3o7Ozs6cjXotQbZX/uijj+Laa6/F7t278eijj+LRRx/Fd77znZzHGAwG/Od//id6e3sxNzeH/v5+7Ny5E/X19QCA7373u7jpppsKPkddXR1+97vfAQCeeeYZBINBGAwGWkyTm2S4kW6hKbxSFAflNjJwo2ufz0cLhQaDAVu2bFm1ii3X2zQYDMJkMpU0/ClEuvzUAcMwKyJaYrXY399f9KYlEarNZqNksLCwgBMnTlAtNvc4zDBMDgmXAtdekFT9FxYWcPLkSdTU1MDpdOKW/kZ84opW+KNJvHjMgz8dc+PJA2ex540zsOoVeP8GFnf8lRrddnGjaLjg+sweOHAAdrsdc3NzGBsboyN4yun4EuOPUEoJwVVlVJJ0uambtUi4gIyk++yzz2Lv3r0AgDvuuAPbt29fQbrr1q2j/25qaoLD4YDH46GkKwZmsxnBYBBNTU1QqVRIpVIVieyWlpbw9ttvr5iqUM4ftFBzRDEUiiS5o31mZ2cr7qpESIwM8iR+uFImX2Sz2RyFBCE/kiLxer3Q6/VoaGgoyymNkAG3SWF8fJyuvb6+PmcNXBIuBb7hC9ePl6z9/1zSgJsvbUYolsIrJzz47zdP4jcjbmzf6CyLdLngFtvICB7yOmtra+F0OmGz2UTlUqU6geVTQhBVBlGYVNqMZi26ixHIRroulwuNjY0AAKfTCZfLVfTxAwMDSCaTOQMpv/71r+Ob3/wmrr32Wjz66KNFNZrEyByQ1/SGqzggo3X6+vpkNR0Xml7IR3DFIklS+KpEi6nf76daXlJQWr9+vaT3JJ1O0xEwhOC4M+zInK6Ojg5ZIyI+AZOJERMTEzkETN5HoZEvQT4/3oWFBQwODtJ0yA19DnTAjca2TphNxtIXFQFu5M4dwUOcyE6dOkU3AiEtt9lsVpbuLv50jKNHj2J2dhZzc3OyekLwiVbs32+1IOoT/f73vz+vqcq3v/3tnK9Lvdj5+Xl86lOfwpNPPknJ45FHHoHT6UQymcRdd92F73znO3jwwQcLXoNEukD5QyQLKQ50Oh2GhoZkleEAxUmXpDFIjpYcE4UQnJQIuhDIOmKxGIaGhmA2myVH+dyI1mAwYGJigsrVFhcX4fV6aWdWV1fXqkh8FApFDjFxCdhkMsFqtdIGCp1Oh3Q6LSoCBpb9eHt6etDT00OjvYMHDyIWi8FoNKK+RgeoKmsSVCjqHBoaglqtpsf+fAFOqTlmUqDT6SgJG43GspUQXFSqbVluiCLdQnPPgGVT5/n5eTQ2NmJ+fj5nAjAXi4uL+Ju/+Rt8+9vfzjG9IVGyVqvFZz7zGXzve98ruhau/wJJLwiFGMVBJY4p/JwucaciVo1WqxWtra2iZ1GVmyvmN01wDWXKIVrgnYjWYrEgmUzSVIhKpUJzczOcTqckhYUcIDlao9EIt9uN2dlZ6ulQV1dHi3Nk4ofYFASQG+0NDAwgm83i8OHDYBiG5mUr/foZhqHG7N3d3VhaWoLL5aLr4DdAVHr8uhxKCC7OB1tHQMb0Amn53b17N5588kl8+MMfXvGYZDKJj370o7j99ttXFMwIYbMsi9/+9rclhdNmsxl+v3/5RQiIdLmKAzmmKpQDcnydm5vLcczKNxxSynXFgBC+2+2mGtoOztBMv98veD2EaAkxEVIiOkq/30+r2uvWrYNSqaQ+C0ePHgXLsqK6sOQAef0ulwuxWIxK7Gpra2kenRzNjUYjHA4HLBbLihREMQJmWRauxQTi6SxMumWD7ba2NvT29iIej8PlcmFkZIS+/tXywDUYDFRqxbWCJM5s8Xi8oqTLhVQlBP+6a93sBpCRdHfv3o2bb74Ze/bsQXt7O55++mkAwNDQEB5//HH89Kc/xdNPP419+/bB5/PhiSeeAPCONOwTn/gEPB4PWJbFJZdcgscff7zo89XX12NychLAck53aWlpxWPkVBzIAW4lPh6PI51O5zhmlQuhpJvPIrG7u1vSUTIf0SqVSjrKmwwQLVRs02q1dNQ2l4ABFD36loNMJgOPx5NDtPm6ubgifH5utLa2NicSK0TALMticCqIcXcECsWycsMYTeHiv7wPOp2ONkFwJ1IQ4nM6nbJ9PoqB2wBBmjECgQAikQjdCOSYAgyUloyJUUJwwTW7EdKAcq4gG+larVa89NJLK75/6aWX4qc//SkA4JOf/CQ++clP5v39l19+WdTzFfJfyDdVoRzFQb6RPWLAzRdzjVxGRkZkbwMuRrqlLBKLge9Lyk8dKBQKqFQqLC4uUnUF6Y1ft26d4PeOS8Ak8hoZGaFH8HIImGhaXS4Xna5LNhqxKgX+1GCiJLHZbCsIOBRL46Qngsa65ZbkZDqLI2fS+ECB188nPq4lpZiRQOWANGN4PB50d3cjGo3SKcCkGaOujNlzYiRjpZQQJDo2GAzvvvTCaoNbSCOpg+HhYWQyGVmO6gSk6CWUOIirmdvths/ng0ajWbVmBT7pci0SWZaF3W4XZZEI5BqZc/PFJJojRBsMBlFXVyeaaAuBG3lxCVihUAh24spHtJ2dnWVXyrkqhZ6eHtp9NzU1RdUBVqt1WWecTC2/f1kWUABqJYN0lgVbuDUfQGlLSqfTSV9HJSZRAMvvn1arhclkosZApOPv6NGj1PVL7Bj2cqJQvhLC7XZjbGwMqVQKOp0Oer1+TcvFgPOYdIPBIM6ePYubb74Z999/P7LZLPr6+mTPhZEouliuiByBSGHOYDDAbreXNPyW29leoVAgHo/j7NmzslgkEgKPRCK0j50QrdvtRigUylFXVOo4xydgl8uFI0eO5CVgQgxkdLpcRFsIXALu7u6mhtyEgOstNqgYFqGlBGq0KrjDCZjVWaiVCsF//1KWlHa7vWKky83p5su7ulwuHDt2jG64JNpfDfA9IY4dOwa/34/9+/fDarXmjA5bSzjvSHdgYAD33HMP7HY7YrEYfvazn8Fms+Ho0aMVKT4UknfxW4GNRiPsdvuKttdCkFNTS9qA5+fn6cwuqZF1NpvNeb2tra2YmJigR8J0Ok0HYW7cuHHVj3DcHGgsFoPL5cLw8DCNnojxdXt7u2w5SKHgWjX29PRQAq5fcmM8kEU0xcKiBbZvasppcRbTQZXPkvLMmTMIh8M4fvx4WZIrPoqpF/K1RRPJHenKE+OxUS7UanXO2J9AIHDh53QBYf4LwDKRbdmyBQDQ1taG5557DgAwOTmJW2+9FT6fD/39/fiv//qvFcTR19eHvXv3Qq/Xo7+/n8p5KjURmJsvJk0CxBOXGH5LKcwRMpdKuiSFQXwIHA4HmpqaoNfr4XQ6RV2LT7QkdRAIBBAIBJDJZFBXVweFQoFwOIxoNEq1pueqWpzNZhGJRBCJRJDNZun6QqEQwuEwampqoNfrz+n6YrEYlpaWUKMCtvcs+yGEQiEshRcxPz9Pp0Xwi3CAMAImbcB1dXUYHR2lA1vJv51OZ9k+DFLaorn5btK6bbfbodFoKnr0JzldpVIJi8VSsecpF7KSrhD/BWDZtu7w4cMrvv+Vr3wF9913H2699VZ84QtfwJ49e3D33XfnPIYUEsgfj0Q4lfpjMgxDBwxGo1GYzWY0NjaWbR8p1vSGb5FIJhRzLRJnZmYES8YI0ZIjLpdoSecZKUL29fXlvNZYLIaFhQUMDw9DrVbD6XSWPW5G6Jr9fj9cLhcWFxdhsVjQ2tq6okhK9KdkfUI7sORYH/F5WFxchNVqRVtb24r1ca0IueuTSsBk8+ZbUhIfBu5QzEpHf/m68rjG7JVKhQC5Ot212o0GyGxivn79euzdu5c2SGzfvh0nTpxY8bja2lpEIpHcJ/pLoWdhYQEqlQoHDhzAQw89hD/96U8Fn++SSy7Bvn37wDAMBgYGRI38KQau4Xc4HIbRaJQ9Lzg6Okr1sIVQyCLRYrHkJbi5uTlks1m0tLTkvV6+iBYAJVru2HChao9oNIqFhQV4PB4aZYvt8S+GfERLojuh6yMG3BqNRvYNgutzEAwGKcEJPeKTDcLj8UClUq2YiiuEgMPhMCYnJ3HRRRet+Bk39+r3+6klpZCmg/379+M973mPkLdBEGKxGGZnZ3HmzBkquStmzC4Ww8PDWL9+PQwGA7LZLLRa7blMMayOiblQ/4V4PI5LL70UKpUKu3fvxkc+8hH4fD7U19fTm0HIWHau0Q2pskslxUQiQYk2nU5TSRHJDRmN8vbJF4p0+cY2xDZSSK6Y5DQJ8nnRkg8hSZOQYlO+iEwIampq0N3dja6uLhrVnD59mpqsSHG54hJFKBSiLcj8iFvo+kgDALcVl+g9pRAwt2mCS2RSTj/cBgWSoz5y5AhNHRACLhYBF6sNFLKkPHHiBJX1yblJFoNer0djYyMikQg2btyYozyQQxJ3wUrG5PBfOHPmDJqbm3H69Gns2LEDW7ZsoT66YlBXV4dgMAiHw0FJTMwNxHWzApB3ykM4HEYymRS9tlLINxGYa5HocDhyjG1YlsVMMIZoIgOzQQ1HnsGI5JrkP/L+kxuK3+IrZ7GJW0Tq7u6mjQQnT56EyWSi+cVCBJyPaBsbG2Ut1nHlRmSDGBoaojKvYoUfbq7S6/XSLikp7mqFoNfr0dHRgY6ODjquh8jkSASs0WhWELBQSSM/95rPkpK8B5VK15GCLFeTTSRxExMTiMVisNlsok5b/GuvdYheoRz+C2TceldXF7Zv347h4WH87d/+LYLBIH3jhIxlJ1pdh8NBC16l3vRoNEqJlgxD3Lx586pPeVAoFLTiS6RXZApvvhto/+kAhmcWoWCALAtcs86KTY3L0TfXizYY/P/bO/OgqM50/39P0yw2INA0TS/sm/sa0aTMqKM3yR2TMPGGQQVHMLGuM/Wb3GQyE8tbuZXKZCY1y03VTGWc/Ly3pn5BI5vGEI2D3mw6jomaixHFKAgiCvTpbrpB9gZ6+f1h3jenD93Qy+kF6M9fMSL9dvc5z3ne5/0+3+c+pFIpbWMl5jnkYvalfIrAbyQg7/PWrVtUYkYOWLlb8/j4eL+pIrgeBHyZF7dEQv6uu7sbEokECoXCL6Y83HE9RI/a2NhIPRKITI5I6LjnGu42eziypPSV9MvRPcqXxHGHhJIDQVfLNWTHG8wI+lhwxX+ht7cXEokEkZGRMBgM+OKLL7B3714wDIPvf//7eP/997Ft2zan/55LfHy8nemNIwUD/wAqMjLSrcnA3jqYceFaJPb09CAmJgbp6elTOoj1DI2hoasfyrhIiBgG4xYr/tFqRHZiFMKY7/S+xBmrubmZ9s0nJye71XklNNzsirz/u3fv4tq1awAePDhTUlIEm2vnCSQAkwy4s7MTTU1NsFqtkEgkSE1NpeWwQMDXKbMsi8uXL2NsbIweZKamptLrVCQSuaWCcHT41dXVhcHBQVy+fNnjsfCOmEqxIxaLoVAooFAo6ATgrq4uaoLj6oFgMB+kCXoVueK/cPPmTezZs4fqVPft20dFzL///e+xbds2/Md//AdWrFiB559/ftLXS0hIoJ66XKcxR4bfRLvpbobibabLt2ok9cmYmBiEh4c7lNTxGbfYIMKDyrzVaoUINoybLRgbtyA6Khw2m81u0i85bCM1zKamJigUCsFuHHfhZru9vb2Ii4vD4sWLYbVaodfrcevWLchkMrsuK39jMpmg1Wqh1+shFovpAE+j0Yh79+7BYDDQGnUgpsuSpg+WZTEyMgKlUgmpVEoPWnt7e2kGHBkZSVUs7gZg4MFDKDU1FcPDw8jLy7NTH5DryNN2bHdKAKS5Jykpya781NTUNKEeHezZLZdpOYKd8Oabb0KhUKCoqAi3b9+GSCTC2NiY3egYb2Uyw8PDaGtrc2tctCOLxKSkJLv66VRKA4LVaoVp3IKayxqYrTbMjRKjd9iMxGgxHlWFobu7225CrCM1BNmC6nQ6v0m8+IdNk9V1SUDRarUYHh62M7b2JcRghmzPnWV0/HpuTEyMXauvr+BL0CZ7MBHzfZ1OR81y+G5t7gRgR4oIorTQ6/UeW1Leu3cPDMMgNTXV5X/Dh3ttGQwGSCQSyGQyaDQarFmzhu78hDZKchOnH/C0Drp/+tOf0NzcjLy8PMyfPx/R0dHIyMigQnkhGB0dRVNT05RTRbkyM2KRKJfLnW7rdTodTCYT0tPTJ/ydo1Hj/aNWnGk2oNNwH1GWYSxJZJCqlNPXcBWuxIvUKIWq3zkKtO7qQ0lNT6vV0qxdSFnR+Pg4DbRWq5UGDldvUGcHat7MIuP//r6+Pmi1WvT29iI+Ph5KpdItgxlilqPX62mHnlwutzsgnmow5/3799HV1YVFixY5fA3ug9wdS0pSMyYqJ28h5cPOzk50dXVRH16lUim4AbubzKyge+bMGRw4cACXLl1CXl4eXnnlFaSlpcFisXj1BHWExWJBQ0MDHnrooQl/x7dIJH3wrnzZRAOclZUFwHGgJaUN4ndrNpvtXJW8gQQPrVYLo9FITU3c7WDiWh4ajUZBgxD5fLVaLX3vnpidm81m6PV66HQ6O3mSt23jQr73wcFBsCxLA7lCoRCkmYG4y+l0OrcCMNmpzZ8/f8rX4O4YSJbtbNd169YtJCQkICkpyav3xae/vx93795FVlYWtFot1Gq1R4ooAfFf0HWlFfjMmTP4+c9/Tv/c1NSE6upqPPPMMy6NYr906RLmzJkDjUaDuro6vPnmmxOCmFDYbDbU19cjPz8fgGOLRDKJ1x16enpgNBqRlZXlNNCSbIxst31lbM3VbxL1g0KhcCrZ8WWgdQbJ3lzNUPmeuZMFAiHgSrBczfJJZ59er0dUVJTgjSV8yENMp9PRkhT/AU4CMGl1z83Ndes1uN+TI0vKGzdu0Ie7kHAfElarlU6mCCD+C7p79+6FVCqlrcC9vb0OW4EJPT09yMnJQWdnJyQSCcrKyvDUU09NOoqdcOnSJbzzzjvYv38/HeDInTgsFBcvXoRarUZ3dzcA0OK+J9MNSEbb19eH1tZWpKen00GJJKO1Wq00o/XHBAH++kh9lbhYKRQKREdHB6Su6QhHtdjk5GSEhYXR0gSxciQ1UH8yWZmFW94gB1P+aFHmQ8phOp0Oo6OjkMlkNAAPDg6ivb0dCQkJSE9P91gJwN2pEEvKgYEB5OTkYO7cuYK+H+Iql52dHfRBV/CTFFdGsXN5//338YMf/MCj7TLXyFzIicDAd40TpEbLMIxXFon8zrD4+HhkZWXh3r17aGpqglgsRnJysk/sKd2Be2JsNpvR0dGBhoYGjI6OUvnU6tWrA3KCT+CafQ8NDaG9vR2tra2wWq2Ij49HZmamYE5bnsAw9hMnjEYjlckRj4TFixcH9HsODw+nbmVmsxldXV1obGzEyMgI5syZA7VaTZ3MgKlrwM5eg29J2dnZiatXr9JSkafDBfhwp0YAwduNBvgg6Lo7ir26uhovv/yy3f9zdRQ7fzilt0GXWCRyGyeWLFmCq1evTqky4MMPtAzDUFkbKU8wDAOFQoHFixdjaGgILMvi6tWrSEpKglKp9MuYFj7kYIJ7Mpybm4u4uDiaAbMsC4VCgeTk5IBJ0MipPuleW758OcRiMZ2rxR01HiiJF8m6TSYT5HI5FixYQLP0hoYG2j5Mxr77m7GxMeh0Omi1WohEImRkZNBERqfTgWVZmgFzD4S5BkmuBjdiSdnZ2YnFixejv7+fNkB4OwUYeBB0uXEimIOuR+WFyVqBS0tLafYJ2Gtp+bAsi6VLl0Kj0dCtALmhySj27Oxsp6PYzWYzVq9ejb///e8wm824du0aVq5cOdn7sX9z3055IB1qxCJRJpPZbU3q6+uxcuXKKW8MZxNwx8bGaOmAO2LE2djr7u5usCwLs9lMg5uv5S/czqupjGu423uudtMfEjRyqk8ClrODv8HBQWi1WvrgEFKl4QzS/KHVajEwMGA34cHR++F345FOPV8GDFLrZlkW4+Pjk8q++JM3iCKHK1lzNwBfuHABa9asod8DKauRzkxPLSlbWloQFxdHXdYiIiICuhuD0OUFIVqBAeDIkSPYsmWLXYBzZxR7WFgYPYRytYmBZHKkK8yRRSIfkkU7yuocBdqwsDC7AzfiHrVs2bIpM0OuSfXY2Bi0Wi3dliqVSkH1tfxAm5yc7NJ0ZO72fnh4GFqtlnoYKJVKweq8fIUFOdV3ZRxQTEwMcnJykJ2dTX/H7du3PZKxTbXGvr4+sCxLXcZSUlKmlHgxDAOpVAqpVErbt4kRjbtOZa6ssaenByzL0odBbm7ulLVu0tGYnJxMtdT37t2j2alcLkdsbKzdtc8Nws7Wwv3cuUmIN5aU08XsBvBBecGVVmBCVVUVfvvb39r9P3dGsZMPloihnWXtziwSMzIyXAoO/IDuLNASgTqx6UtOTsby5cs93oJHRETY1S5JcCPmJJ6oBUiXWnd3N3Xa8qbFVSKRICsrC5mZmTS7JMFtKpMbZzjKUrOzsz0K5NwWV262TDIjT7Iq8uAmayTv1dORRSKRyC4A9/b2gmVZNDU1eRyAyTXPsix6enoQHx/v0sPAGVy/XhKAOzs7MTAwQAMw+YxdHU3v6HOQyWSQyWT0QUE60KaypJwuZjeAD9QLRqMRRUVFuHfvHm0Flkqldq3AANDe3o61a9eio6PD7kPcuHHjhFHskz2RV6xYgTNnziAsLMzOU9eZRaJUKnX75r158ybUajUkEsmE0gEJtMSv1deG2USaxLIsent7kZCQMOV0VkeB1pcdaY6yP3fWyJ024KstoqO68FRrJFm9Xq8XPKt3BL8jzRUvYe4aJRIJXaOvyipE7aLT6ajxvVwut9PIcgPwhQsX3Pbo5Uoae3p6qESRW/66fPkyFi1ahKioqGDw0gVmWnMEl/Xr1+O9995DQkICvvrqK2RnZ8NgMNhZJHqSbQHfZbRtbW2wWq1Qq9WIjo6GyWSiGa0/Au1k6+PXEEknDrdlM9jWSOqcxD9Wr9f7dcLDVGskY8ZjY2NpiUen01EzFl/XrydboyMzdyJD02q1fquxu7JG0h0WHx+P8fFxOsdv9erVbmXAXLh6aIPBgOjoaCQnJ6O9vZ2qaqxWK6KiogJdYpi5QbegoACbNm1CdnY2bDYblaJ4eiJMAi3ZIhG/Uo1GA41Gg/HxccyZMwepqalQKBRBs6WxWCzo7OxEZ2cnRkdHERkZiZSUFKhUqkDrFSkWiwUsy6KjowMjIyOIiIiAWq2GWq0OiArCERaLBTqdDh0dHRgaGqLSqrS0tKBZIzEWruDrAAAgAElEQVQK6ujowMDAAC1lpaene6Qd9wUkA7537x76+vqoFDErK4seCntSguDCbcm+c+cO9eGVyWR+H0rqgJkXdE+fPo2amhp89NFHWLt2Lfbu3QsA1B3KHfg1WhKsuSbn5KApJiaGnrZGRUVBqVT6dew0H262SG6++Ph4mrkFMvMhOOomk0qlNGOxWCxueyAIDTnVJ54P5OYlh41E3O/Lrrap4GfkpMZKfBAGBgZo4AlU0OGXv4gagWS6xDuaZMBcFYQ3AfjLL7/E0qVLodVqqQdLgPFv0D169Chef/113Lx5E1999RVWrVrl8OdOnz6NF198ERaLBbt378a+ffsAuDYV+L/+67+wdOlSVFdXY+PGjdi4cSOtvbrS7eJsAu7w8DANtBKJhHZd8QMWOUzRaDTo6emhkw6EEntPBpkqQCRok/md+rvGR3DHN8GRBM1fgyQddd85OkPgttAK6d8wFe7U8MkBF+nO8mcAHhkZAcuy0Ov1iI6OptaT/OuML5UjMi/uYaazsUSTQea5kXgWBBm/f4PuzZs3IRKJsGfPHrz11lsOg67FYkFeXh4++eQTpKSkID8/H1VVVVi4cCGKiorwL//yL3Qq8LJlyyZMBSa88cYbyMjIwJYtW9Da2kpPgR3hLNASra7RaKSB1p0eeJKBsCyLoaEh6nIk5A1JshlXAq0jHJ1mu+teNRX8bNEThzBf+xG46zPhCHd9IDxheHgYLMuiu7vbI7UK3y6TBGAh/Yq5tWSGYaBUKt3aUfF9lskZjFQqdSsAW61WXLp0CY888sjsDbqEDRs2OA26/Gm/RDq2b98+t6YCv/322wCAXbt2ob29ncrBCPxAS75A7jQJIX0EiKMVy7Kw2Wy0ucGTrI2MaSEeAySIeVtbdKTb9NSjwNHNLZTfAZFldXd3e+W8xdf7eioVc4SrnryuQA7tyLUvlC6b3+TgTQAmzQwsy8JkMiE5Odkj5zc+jvwqXA3Ao6OjaGxsxKpVq4LFSxfwp/eCq3R1ddnZMKakpODSpUtuTwVOSEjAnTt3AHw3PcJRoA0LC7PT6hLZiaf6T2eIxWLab84dreLq1p7cxHq9HsCDZpOlS5cKehGR0T6JiYk0O21paZmyQ4nAlwnJZDJkZGQIvo3lNjhw9bWuZumOvINzcnIELa1wG0XI9ImGhgaqdJiqTOKoA3HZsmWCft/8JgeDwYC2tjaMjIy4ZBjvqE6bmZkp6PfN96sgAbi1tRWxsbE0AJOJM9wAPJ18FwAvgu5krcBTzTYTkoSEBFy5coV+EaTuRXR6/f396O7utnt6CjnBdTKioqKQmZmJjIwMurVvbW2dUP8dHR2lGS0AyOVyLFmyxC9Pa1JDJa3XOp0OjY2NEIlEdttFvq5VKpV6JbZ3B+4NSbSrXV1duHnz5oRpCo6mZPhrvllUVBQyvp3mSyR7X3/9NSIjI2kAJpImbi2ZTKH2xwEdNwATw/jbt2/DZDLZZcCA4zqtK92A3sIPwMRG9Pbt24iJiYFcLqfJC3lokWuUZLrBjMdX4mStwK6gVqvR0dFB/0ym/yYmJro1FXh8fByff/45iouLkZmZidHRUVy5coVORyV6Rn8FWkdwu6JI/ffOnTvo7++HSCSisqTFixcHtBYVERFBx2KPjIxAo9Hg4sWLAB5kO0S/6o9pvc4QiUR2WbrRaERraysGBgbAMAyVoXnTCSgEEokEmZmZtFOPZVncvn2bBgaZTIa0tDS/HLw6gzsEkgSvW7duYWhoCACo7DCQQzkZxn6yNCkTtba20s+SZN5arRbHjh3D/Pnz8eSTTwZkva4QsPJCfn4+WlpacOfOHajValRXV6OyshIM495U4O9973soLi7Gnj170NvbC4vFgl/96ld45JFHMDQ0hN7eXojFYkRGRiI2NtaP73Ai/AOY1NRUiEQiOzMcT+u/QsE3KI+Pj8fcuXMxMjKCnp4euq0LpHUi8N22XKvVYnx8HCkpKQgLC4PRaKTz5wIpQSOQTjsynSMuLg4mkwlGoxEWi0WwCRHeQHYPZN4e0XYbjUZ0dHTQQ9EAj7/B2NgYenp60Nvbi9jYWDoNvLi4mGa9zz33nFO1VLDgk4O02tpavPDCC+ju7kZ8fDyWL1+O//mf/4FGo8Hu3btRV1cHAKirq8NLL70Ei8WC5557Dq+++iqAB3OUtm3bhp6eHqxYsQKHDx+e9OYpKCiASqXCunXr0NraitraWuTk5GDHjh1Yv3497t+/D5ZlMTo6Sp/s/sqC+ONSnEmNSP1Xp9P5VdoFTLRzdHZaTg7gSFuqNwdwnkC25SzL0gGWxGCdCymT6HQ6ap8pl8v99jAjdXmtVktLHHK53O7sgH9wJOThnis4qtMqlcoJdVq+VI7UgP1lO0rOHDQaDSwWCy153bx5ExUVFTh79iw2bNiAtWvX4saNGzh16hROnz4NmUzml/VNwsxrjpgMq9WKixcv4tChQzh//jyeeOIJFBcXIysri94MERER9HRY6MBGhOBTBVpH8KVdvtT/kiyMeAm4I88iBzIsy9KxLEKcYvPh15K5LbqufB6keUSn002orQoJUa2QZg93/Ib52lVXfCA8xVU9rSMcaa99oVUmDySNRoO+vj7a3j44OIiamhocO3YMarUapaWlePLJJ4OmU5DH7Aq6XEwmEz788EMcOnSIzm/70Y9+hIiICGg0GhiNRkECmy8yAq7V3fDwsCD6X77fgRDj2Pl6TW8zS/68MaGyQK4zmDdObQR+5k0ePN5+51yTG2Lw7e216Y2e1tnv9DSxcAb3gRATEwOVSgWJRIJTp06hsrISRqMR27dvR3FxseCDLX3A7A26XDQaDQ4fPoyamhqkp6ejuLgYmzZtwsDAAA1spPzgSsbGH/Qn9LhwLmazmd44NpuN3jiuBDa+3teX223S3KDT6agTl6tt0lyJlxBB0Rn8oB4fH0/9OqYKbHwHNU+aK1zFUcecq/pa/gMhOTkZSqXSJwe13JHv3GnArrwWV9cOgF4vjY2NOHz4ML744gs8/vjjKCsrw5IlS4JemcAhFHS5WK1WfP311ygvL8fZs2exceNGlJSUYN68efQCIKbh/FocOcAhA/0C0Y/vSv2XW9cEQEX7/jpYInVilmVpQ4JSqZwQ2Ii2lWTepOTjL9d/srVnWZZmlo7KF2ScEmnUIBNt/XUAxh3/42y68WR1WgBo6OzDTe0g4iXhWJ+biOgI4c/R+YfFjgIw+cw1Gg0GBwfpDu7+/fuoqqpCbW0tsrKyUFZWhscffzxoDJvcJBR0nTE6OoqTJ0/i0KFDYFkWhYWF2Lp1KyQSid1NFhUVhcHBwYAFWkfw679xcXGIiIhAf38/zGazS40O/lonN7AlJCRALBbj/v37YBiGrjPQNxfptiLddSSo9vb2BuSB4Ay+ciMhIYEeckokEqhUqgl12k9u6lH1v10IDxNh3GpDakIU/v2JXESF++69kAc/mXCdkJAAs9lMxy2pVCpERETg5MmTqKqqwuDgIEpKSrBt2zanrfzTiFDQdQW9Xo+KigpUV1dDKpUiLy8Pra2teOGFFxAeHg6r1UpH6QRykisX7g04PDxM9clKpRIKhSLo1smyLEZGRmhAIOsM9IOBQLa7Go2GToEmjSLB8AAjkJoqsRsFYOf3y88s/09NI2IixIgQP/jc2T4TXtiQieWpcQ5/v5Dr1Ol06Orqol2iNTU1CA8PR29vL65fv44nn3wSpaWlmD9//nQqH0xF8LUBByNyuRzbt2/HhQsX0NjYiKGhIXoQt2PHDixevBgGgwE3btwAAEEOJDyB30vPn3lF6r83btygATg5OTkgptbczFEul9t1XpEb8vr16wDglU+F0OtctGgRfWCRLXMwrJNfp12yZAkNsKSzkayTu4MwW2wIE30XB0QMYLb6JqfiHwCT1ubu7m5UVlbi6tWriI6OhtFohEwmQ0FBARYsWOCTtQQjoUyXx+joKJqbm2nRfnx8HKdPn8bBgwfR3t5O3c/i4uLsTlqJ9MZXT2pHfgfc9ldncA+2/KH/5Wt5ExMToVQqp1wnqe0Sn2Kh3cUcrZPvNuZIp+psnVz1B7/uL/Q6XdHTOlonqemLxWJcMETgf9kxJEgiMDJuwZzwMLz+1DzEzRHuwUGsTrmKoLCwMJw4cQKVlZWwWCzYsWMHioqK6Difrq4uMAwDlUol2DqChFB5QQiMRiOqqqpQWVmJuXPnoqSkBD/4wQ8wNjY2QVMoRL3X0YwsT0/KHfmyqlQqQUxLvFEDOIK0epKhj+TQSogHGte5jAyU9PRh6Uu1hZAPy+HhYXSxWnx0TYP2fgYqaSx2rM1GitT7a5TrjBYREUHryRcuXEBFRQUaGhpQUFCA0tJS5OTkzKTywVSEgq6Q2Gw23LhxAwcPHkRdXR0efvhh7NixAytWrKDbP7PZTLf17mxD+RmYtwHMEY62qZ7Uf4WyXnSGo8YI7mm8q3Cz6MjISCpLEtKjVwjrSFIWYlnWZ510pPOQOK+56x0NfFeO0Wg0GBsbo+UWjUaDyspKnDx5EsuWLcOuXbuwbt26gB88BohQ0PUVZrMZn376Kd59913cunULzzzzDIqLi5GYmEgzgKkyFX6m6M34cnfhDgwkN/pk9V+uybgvApgz+D6upBHB2YOCvC+ia/ZXHdaRG9tkuxOhHoCerJPbLEKsTie7Romevbe3lz4AgQdt/1VVVRCLxdixYwcKCwsD7nMSBEyvoOvquB/i4RoWFgaxWIz6+no/r9Se+/fvo6amBocPH0ZUVBS2b9+Op59+mg625F6sMTExbl30/sDZlpbbmOHr5gpXcBZQw8LCaGAeGRnxWwBzhiMLR9Kxxg1g/hz15IjJHvrj4+NUE06aXeLj43H+/HlUVFTgxo0b2LJlC3bu3ImMjIygLB8EKJ5Mr6Dryrgf4MGHVF9fHwzmFnbYbDa0tLTg4MGDOHHiBB566CGUlJRgzZo1uHDhAvXQlUgkSEtLg0KhCKotGMnW7ty5g76+PipFSktLCxrJFGFkZATt7e1UjJ+QkIDMzEy/+Py6g8ViQVdXFzo6OmAymSCRSJCeng6FQhFQhzE+5CC0vb0dfX19CAsLQ3h4OFasWIHOzk5UVFTg1KlTWL16NUpLS/Hoo48G1fodEaB4Mr0kY9NdPsIwDPLy8vDmm2/ijTfeQFVVFfbu3YvOzk6kpKTgrbfewsMPPwyj0Yiuri50d3dDpVIFdKowYN/1RLbw8+fPp/O6GhoafNpO6g6k240cti1duhTh4eHQarW4efOm4AdwnsJv387IyIBUKqUlBY1GQ6VdgTRucWQyk5eXh56eHvzyl79EfX09oqKi8Pzzz+PixYsBbwxyh2CLJ0EZdF2FYRg8/vjjYBgGe/bswb/+678GekkTMJlMqK6uxr/9279h48aN+Oyzz/Cb3/wGALB161Zs2bIFDMNAo9Hg9u3bgqoKXIEoJMi8NJlMhqysLFjEUWjo6EN75whWpMVj6dIkuq2/fv26S/VfoXEkK8vKyrLbJcTGxiInJ4faeTY3N7vtTOYtjuq0CxcutCtzpKSkICUlhUq7GhoaBDMgcge+yYxSqURubi7OnDmDX//617h9+zYKCwvx5z//Gd988w2OHDmCK1eu4NFHH/XL+vyJv+JJwMoLroz7mWywJfBA46dWq6HX6/HYY4/hz3/+M9atW+erJQuGzWZDW1sbDh06hA8//BCLFy9GSUkJ1q5dS7MNUpNUKpWC+yWQrEar1Toc660fGMW/f3gTfaYHnU5xUeH43TMLkBT73Tq4N6uz1lMhIA0UOp3OowMxfuODr+q8Qlhycq02yUw3X+x+LBYLVUlwzZNaWlpQWVmJjz/+GI8++ijKysqwevXqoC8fAEEZT6ZXTZcw1YfE5fXXX0dMTAx++ctf+mFlwmG1WvGPf/wD5eXlqK+vx+bNm1FSUoK0tDQ78x2VSuV13z9X/zqZacv/PdeOz5sNkMU82O4aBsewcV4ifrouc8Lv9FS8PxlWq5W2NruiVHAVrqOV1Wp1y/PWGb4wn+ebBXmrJya/k5jMDAwMQC6XQ6VSYXh4GEePHsWRI0eQmJiInTt34oc//GHAy0e+wM/xZHrVdF1haGgIVqsVsbGxGBoawscff4zXXnst0MtyG5FIhPXr12P9+vUYGhrCsWPH8Itf/AImkwlbt27Fs88+C7FYDJZlcefOHWoU4upBERmQqNfrnW7J+dwfHkdE2He/OyKMwf1hs8Of5c6wItvq9vZ2t+0E+VIrUuYQUnrEn9RMtvURERFumZs7stl86KGHBFNzMAyD2NhYxMbG2um2b9265XaDzPDwMDQaDbq7uxEXFwe1Wg2JRIJPP/0Ur732Gjo7O1FUVITa2loqAZuN+DOeBGWm68q4n7a2NmzZsgXAg5uguLiYjvuZ7thsNnR0dODQoUM4duwYcnNz6eghIjXiWuLxM0AyLkan09Fpv1ONAufySZMeB87dRfy3LaL3R8bxk3XpeGy+3OX34KpxNjf79vfIGgK/ycNRS7cvDOXdhaxBq9VicHDQrhWcC/nsuRalSUlJuHHjBioqKvD5559j48aNKCsrw8qVK4NK5eELAhRPpmd5QShc1emdPn0aL774IiwWC3bv3o19+/b5eaUTsVqtuHDhAg4dOoQvvviCjh7Kycmhrl3AA7MeANRGz5uts81mw7EGFieuamED8MNlCjy7XOnxzckfEZOYmAiTyYTu7m5ERUW5ZXTuS/gG5QkJCYiLi0NfX19Q6Gm5cMclEbvRqKgoGI1GDA0N0V1Gf38/ampq8P7770OtVmPnzp146qmngnXEDYDpfb9ymN1B1xWdnsViQV5eHj755BOkpKQgPz8fVVVVWLhwYQBW7JiRkRE6euj+/fvYsmULwsLCMDw8jPz8fNhsNsTFxSEtLS3gUik+4+Pj0Gq16OrqwtjYGAAgKSkJarXab6oCVzGZTNBoNHSqMLF2VCqVfhvI6CqDg4Po7OyEXq8HADrqXa1W429/+xsMBgO2bduG4uJi+mAOdmbI/Trzarru4IpO76uvvkJOTg6ysrIAANu2bcPx48eD6UvEnDlzsH37dmRlZeE///M/8ac//QkqlQrp6elYuHAhNm3aRGt4zc3NSEpKonOmAgE/G0tOTsby5csRFRXlVf3XF3AP2YhKYs2aNQ9sEb/9u6amJreHTvoCrslMeHg4VCoVcnJycPXqVXz55Ze4fv06Ojs7kZaWhp///OfYunVrQNbpKTPlfnXGrAi6rtDV1YXU1FT655SUFFy6dCmAK3KO1WrFK6+8gtWrV8Nms+Hy5csoLy/Ha6+9hn/6p39CSUkJ8vPzYTAYaKDwxHzHE/hTImQyGXJycibUHUUiEZKSkpCUlGTnqyvU4ERXIDVSlmUxNDQEuVw+QU8L2B/AjY6OQqvVUl2tv6ZJcL0nuA8wo9GIyspKfPDBB8jMzERZWRneeecdhIeHo62tDY2NjT5dV6CYTvcrnxkTdF3R6c0UHnnkEfrfDMMgPz8f+fn5GB0dxUcffYTf/OY30Gq1+NGPfoSioiLExsZCq9Xi8uXLVFObmJgo2Jbe0Tw0lUqFhQsXuvQa4eHhtFmA1H/r6+vdHhHu6lr5etq0tDSX67SRkZFIT09Heno6nZvW3t5OGwuENCnir1UqlSIrKwtisRgnT55EZWUlBgYGUFJSgo8//hiJiYl2/z4rK4tmgsHGbLpf+cyYoPvpp5969e/VajU6Ojronzs7O6FWq71dll+JjIxEYWEhCgsLodPpUFFRgaKiIiQnJ6O4uBhPPPEE1ZW2tLQgMTERKpVqQhbqKnyDHIVCgZycHK+Czpw5c5CVlYXMzEyq/21paYFUKqVr9eRh4UhPm5ub69Vao6OjkZOTg+zsbNps4omsi8/o6OgEk5ns7GzU19fjL3/5Cy5evIjNmzfjj3/8IxYsWBBU9XBXmc3366w4SCNMJo42m83Iy8vDZ599BrVajfz8fFRWVmLRokUBWKlw2Gw2XL16FeXl5fjkk0+wfv16lJSUYMmSJbRVlXiiKhSKKeuU3CnD/ioF8G0dSVfZVPVfR3VaX5dYHDmLuXIAR+rfGo0GZrOZrlWv16OyshLHjx/HvHnzsGvXLmzatMnvo5cCwTS/X2e3esEVnR4A1NXV4aWXXoLFYsFzzz3nVKfX09ODrVu3or29HRkZGThy5AgSEhIm/FxYWBiWLFkCAEhLS8OJEyd89yZdYHx8HKdOncLBgwdx9+5dPPvss9i6dSukUik9mOF65JIs0GKx0A6xsbGxgE4Znkr/66hOG6hBotxhnBaLhT4syIONL1GTyWR0bM3x48dRVVWF8fFxlJSUYOvWrYiPj/f7e3AFoe8Hoe/XADG7g67Q7N27F1KpFPv27cPvfvc79Pb24ve///2EnyOeucGIwWBAVVUVqqqqEBcXh5KSEmzevJn6pxoMBkgkEthsNphMJuoF62kpwheQ+i+ZCCESiTA8PExbkYNBT0vgNqwwDIPw8HAMDw/TZoz4+Hg64ubKlSt4+umnUVpaitzc3KB5D86YCfeDDwgFXSGZN28ezp49C6VSCZZlsWHDBjQ3N0/4uelwkdlsNnzzzTd09FBeXh4AIDU1FUVFRbBarXS76wvzHW/guo6FhYUhLCwMIyMjtFYdTNMLLBYLHZluNpsRERGBK1euoLa2FgqFAo2NjVi2bBnKysqwYcOGoPJXnoqZdD8IyOzW6QqNTqejfeoKhQI6nc7hz5lMJqxatQpisRj79u3DM888489lugTDMFi8eDFWrFiBzz77DL29vYiJicH58+chk8mwfft2JCUlQafT4dq1a4KZ73iKozrtypUraZ2W1H/b2trcqv/6Ar58Ti6XY8GCBbBYLKitrcXRo0cRGRlJuwhTU1OxadMmv6/TW2bS/eAPQkHXCZNJWrgwDON0+3f37l2o1Wq0tbVh48aNWLJkCbKzs32yXm9Zs2YN/v73v9PssLe3FzU1Ndi1axdtyigoKIDNZvPYfMdTXNXTAg/0v3K5HHK5PGD6X2L6rtfrMXfuXKhUKsybNw/nz5/Hb3/7W3zzzTfYsmUL3n33XWRmZoJhGJjNZjQ1Nfl0Xd4w2+4HXxIqL3iAq9spLmVlZXjqqadQWFjop1UKg81mQ3NzMw4ePIiTJ09i1apVKC4uxsMPP0yNwicz3/EGYu7jqT8tFxIIu7u7qabWG6tEPtxJvqRtOCkpCXfu3EFlZSXq6uqQn5+P0tJSfO973wu4z4SQzKb7wQ1CNV0heeWVV5CYmEgPDnp6evCHP/zB7md6e3shkUgQGRkJg8GARx55ZNq0KTrDYrHgzJkzKC8vx/Xr11FQUIDt27dDpVLZTRQm3W+elB/40yGENsPhKwY8HetOfhdxHuOazJhMJhw7dgw1NTWQSCT48Y9/jGeffTbofBuEYrbeD1MQCrpCYjQaUVRUhHv37iE9PR1HjhyBVCpFfX09Dhw4gL/+9a/48ssvsWfPHohEIlitVrz00kt4/vnnA710wejv78fRo0fx3nvvgWEYbNu2jRrwEP9WrlH6ZBllIPS0gOf6X+58tvj4eOpvcfbsWRw+fBitra0oLCzEzp07kZKSEvTqA28J3Q8OCQXdYGQqa7rR0VHs3LkTly9fRmJiImpqapCRkRGYxTqBO3qotrYWS5cupaOHSHmgv79/gvkOyRJJecLVJgJfwfegJWPmSf2XuKSxLEtNZmQyGZqbm+1G3JSWlmLNmjVBWT6YCdfbNCIUdIMNV6zp3nnnHVy7dg0HDhxAdXU1amtrUVNTE8BVT47VasW5c+dQXl6Or7/+Gps3b0ZxcTEyMjJok8DY2BjCw8MxOjpK67TBNi6d1H91Oh0iIiJgs9moR7FCocDg4CCOHj2Ko0ePIiEhATt37sQzzzwT1CNuZuL1FuSEJGPBhivWdMePH8frr78OACgsLMTPfvYz2Gy2oApQXEQiETZs2IANGzbQ0UMvv/wy+vv7kZGRgY6ODvzhD39AREQExsbGYDabYTY7HgMUSCwWC8xmMxiGoUH37bffxvDwMAYGBtDf34+ioiJ88MEH02bEzUy83qYrwbcHmiU4sqbr6upy+jNisRhxcXEwGo1+XaenREdHY926dQgLC4NIJKIBdv/+/dDpdFi1ahXS0tLQ3d2NixcvoqWlJaDC+dHRUdy9exeXLl1CW1sbpFIpVq9eDYZhUFlZievXr9Mas0gkwooVK6ZNwAVm/vU2nQhluiF8hlwuxzvvvEO1mGT00MGDB/Hqq6/S0UOrV6+GwWBAS0sLxsfHXTbf8RYydZhrMrNy5UqqUT569ChUKhV27tyJt956i3bj9fT0YGRkxKdrCzFzCQXdAOGKNR35mZSUFJjNZvT19U3wTA1mJBKJnfhdJBJh7dq1WLt2LUZGRlBbW4tXX30VfX19KCoqQmFhIWJiYsCyLBoaGhya73gLGRmv0WioyUxubi7Cw8Nx+vRpVFRUoLu7G9u2bcPf/vY3hyNupFKpIGvxJ7PhepsuhA7SAoQr1nR/+ctf0NjYSA82PvjgAxw5ciSAq/YNXV1deO+993D06FFkZGSgpKQEmzZtorPKiIG3N34KXD/d6OhoqFQqxMfH4+rVq6ioqMC5c+fw+OOPo6ysDEuXLp1xdczQ9eZ3QuqFYMSRNd1rr72GVatWoaCgACaTCT/+8Y9x5coVSKVSVFdXB+0kACGwWq24fPky3n33XZw7dw6bNm1CSUkJFixYQJsQTCYTLT9MZb5DTGZYloXVaqVNGwaDAdXV1Th27BgyMzNRWlqKf/7nf/a5LjjQhK43vxIKurOFqbSY5eXleOWVV+jW8mc/+xl2794diKVOyujoKE6cOIFDhw5Br9fT0UNxcXHU+5fMLpPJZLT7zWaz4f79+1QfTNqTGYZBXV0dKioq0N/fj+LiYmzfvj1ot88z5XucxSQMoGcAAAOqSURBVISC7mzAFS1meXk56uvrsX///gCu1D10Oh0OHz6MmpoaKBQKlJSU4IknnsD4+Dg0Gg0MBgNiYmIgEonQ399PTWbmzp2L+vp6VFRU0BE3paWlLs9uCxQz9XucZYR0urOB6TyWejKSk5Pxi1/8Ai+//DIaGhpQXl6OX//613j44YeRkJAAvV6P0tJSiEQinDlzBgaDAdHR0fj888+Rl5eHXbt24cCBA9NmxM1M/R5DPGB6XIUhXMLVsdTHjh3DuXPnkJeXhz/+8Y92/yaYYRgGK1asoKYpZ86cQXJyMgDg/PnziImJwaVLlzA2NgatVgu1Wk3rtdOJmf49znZCzRGzjKeffhrt7e24du0aHnvsMZSWlgZ6SW4TERGBF154Aa2trfjyyy9x4sQJDA0N4dy5c9i/fz/Onj2L5uZm/PWvf0VcXFygl+sTZsL3OFsJZbozCFe0mNyDo927d2Pv3r1+W59Q5OTkICcnh/5ZJpM5nMmVlZU1LU/fZ8v3OFsJZboziPz8fLS0tODOnTsYGxtDdXU1CgoK7H6GZVn63ydOnMCCBQv8vcwQUxD6Hmc2oUx3BiEWi7F//3488cQTVIu5aNEiOy3m22+/jRMnTkAsFkMqlaK8vDzQyw7BI/Q9zmxCkrEQU/Lcc8/h5MmTkMvluH79+oS/t9lsePHFF1FXVweJRILy8nKsXLkyACv1P6HPJoQTnErGQuWFEFNSVlaG06dPO/37U6dOoaWlBS0tLfjv//5v/PSnP/Xj6gJL6LMJ4S6hoBtiStatWzepycvx48exc+dOMAxjN7ByNhD6bEK4SyjohvAaV7xaZyuhzyYEn1DQDREiRAg/Egq6IbzGFV3pbCX02YTgEwq6IbymoKAAhw4dgs1mw8WLFxEXFzetRtn4ktBnE4JPSKcbYkq2b9+Os2fPwmAwICUlBb/61a8wPj4OAPjJT36CzZs3o66uDjk5OZBIJHj33XcDvGL/EfpsQrjLVDrdECF8AsMw/w/AUwD0NpttsYO/3wDgOIA73/6vD2w22xvT5fVChHBGKNMNESjKAewHcGiSn/mHzWZ7apq+XogQDgnVdEMEBJvNdg5Az0x9vRAhnBEKuiGCmUcYhrnKMMwphmEWTf3j0+71QsxCQuWFEMHK1wDSbTbbIMMwmwF8CCB3Br1eiFlKKNMNEZTYbLZ+m802+O1/1wEIZxhGNlNeL8TsJRR0QwQlDMMomG+nRzIMsxoPrlXjTHm9ELOX/w/y14UFOltARgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure()\n",
"ax = fig.gca(projection='3d', elev=30, azim=45)\n",
"# ax.plot(np.array([0,1]), np.array([0,2]), np.array([0,3]), label='parametric curve')\n",
"ax.plot([0,pca.components_[0][0]],\n",
" [0,pca.components_[0][1]],\n",
" [0,pca.components_[0][2]], label='X')\n",
"ax.plot([0,pca.components_[1][0]],\n",
" [0,pca.components_[1][1]],\n",
" [0,pca.components_[1][2]], label='Y')\n",
"ax.plot([0,pca.components_[2][0]],\n",
" [0,pca.components_[2][1]],\n",
" [0,pca.components_[2][2]], label='Z')\n",
"\n",
"ax.scatter(initial_vertices.x, initial_vertices.y, initial_vertices.z)\n",
"\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Raw
curve_fitting.py
#!/usr/bin/evn python
import numpy as np
import scipy.linalg
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
# some 3-dim points
mean = np.array([0.0,0.0,0.0])
cov = np.array([[1.0,-0.5,0.8], [-0.5,1.1,0.0], [0.8,0.0,1.0]])
data = np.random.multivariate_normal(mean, cov, 50)
# regular grid covering the domain of the data
X,Y = np.meshgrid(np.arange(-3.0, 3.0, 0.5), np.arange(-3.0, 3.0, 0.5))
XX = X.flatten()
YY = Y.flatten()
order = 1 # 1: linear, 2: quadratic
if order == 1:
# best-fit linear plane
A = np.c_[data[:,0], data[:,1], np.ones(data.shape[0])]
C,_,_,_ = scipy.linalg.lstsq(A, data[:,2]) # coefficients
# evaluate it on grid
Z = C[0]*X + C[1]*Y + C[2]
# or expressed using matrix/vector product
#Z = np.dot(np.c_[XX, YY, np.ones(XX.shape)], C).reshape(X.shape)
elif order == 2:
# best-fit quadratic curve
A = np.c_[np.ones(data.shape[0]), data[:,:2], np.prod(data[:,:2], axis=1), data[:,:2]**2]
C,_,_,_ = scipy.linalg.lstsq(A, data[:,2])
# evaluate it on a grid
Z = np.dot(np.c_[np.ones(XX.shape), XX, YY, XX*YY, XX**2, YY**2], C).reshape(X.shape)
# plot points and fitted surface
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2)
ax.scatter(data[:,0], data[:,1], data[:,2], c='r', s=50)
plt.xlabel('X')
plt.ylabel('Y')
ax.set_zlabel('Z')
plt.show()
Display the source blob
Display the rendered blob
Raw
pca_labe_fit.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Raw
plane_fitting.py
#!/usr/bin/env python3
import numpy as np
from sklearn.linear_model import LinearRegression
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
# some 3-dim points
mean = np.array([0.0,0.0,0.0])
cov = np.array([[1.0,-0.5,0.8], [-0.5,1.1,0.0], [0.8,0.0,1.0]])
data = np.random.multivariate_normal(mean, cov, 50)
# regular grid covering the domain of the data
X,Y = np.meshgrid(np.arange(-3.0, 3.0, 0.5), np.arange(-3.0, 3.0, 0.5))
XX = X.flatten()
YY = Y.flatten()
# best-fit linear plane
model = LinearRegression()
model.fit(data[:,:2], data[:,2])
# evaluate it on grid
Z = model.coef_[0]*X + model.coef_[1]*Y + model.intercept_
print("Model score: ", model.score(data[:,:2], data[:,2]))
# plot points and fitted surface
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.2)
ax.scatter(data[:,0], data[:,1], data[:,2], c='r', s=50)
plt.xlabel('X')
plt.ylabel('Y')
ax.set_zlabel('Z')
plt.show()
Display the source blob
Display the rendered blob
Raw
plot_pca_3d.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Raw
plot_pca_3d.py
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
=========================================================
Principal components analysis (PCA)
=========================================================
These figures aid in illustrating how a point cloud
can be very flat in one direction--which is where PCA
comes in to choose a direction that is not flat.
"""
print(__doc__)
# Authors: Gael Varoquaux
# Jaques Grobler
# Kevin Hughes
# License: BSD 3 clause
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# #############################################################################
# Create the data
e = np.exp(1)
np.random.seed(4)
def pdf(x):
return 0.5 * (stats.norm(scale=0.25 / e).pdf(x)
+ stats.norm(scale=4 / e).pdf(x))
y = np.random.normal(scale=0.5, size=(30000))
x = np.random.normal(scale=0.5, size=(30000))
z = np.random.normal(scale=0.1, size=len(x))
density = pdf(x) * pdf(y)
pdf_z = pdf(5 * z)
density *= pdf_z
a = x + y
b = 2 * y
c = a - b + z
norm = np.sqrt(a.var() + b.var())
a /= norm
b /= norm
# #############################################################################
# Plot the figures
def plot_figs(fig_num, elev, azim):
fig = plt.figure(fig_num, figsize=(4, 3))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(a[::10], b[::10], c[::10], c=density[::10], marker='+', alpha=.4)
Y = np.c_[a, b, c]
# Using SciPy's SVD, this would be:
# _, pca_score, V = scipy.linalg.svd(Y, full_matrices=False)
pca = PCA(n_components=3)
pca.fit(Y)
pca_score = pca.explained_variance_ratio_
V = pca.components_
x_pca_axis, y_pca_axis, z_pca_axis = 3 * V.T
x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]
y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]
z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]
x_pca_plane.shape = (2, 2)
y_pca_plane.shape = (2, 2)
z_pca_plane.shape = (2, 2)
ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)
ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
elev = -40
azim = -80
plot_figs(1, elev, azim)
elev = 30
azim = 20
plot_figs(2, elev, azim)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment