3D points curve fitting

cubePCA.ipynb

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()

pca_labe_fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = np.exp(1)\n",
    "np.random.seed(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pdf(x):\n",
    "    return 0.5 * (stats.norm(scale=0.25 / e).pdf(x)\n",
    "                  + stats.norm(scale=4 / e).pdf(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.4104351  -0.33153453  0.09879082 ...  0.42705574  0.04086025\n",
      "  0.50043173] [-1.04508952 -0.16008775 -0.00973006 ...  0.49366242 -0.76612609\n",
      " -0.34104571] [ 0.78999815 -0.29325959 -0.22366065 ...  0.03376915  0.99667957\n",
      "  1.11360349]\n"
     ]
    }
   ],
   "source": [
    "y = np.random.normal(scale=0.5, size=(30000))\n",
    "x = np.random.normal(scale=0.5, size=(30000))\n",
    "z = np.random.normal(scale=0.1, size=len(x))\n",
    "\n",
    "density = pdf(x) * pdf(y)\n",
    "pdf_z = pdf(5 * z)\n",
    "\n",
    "density *= pdf_z\n",
    "\n",
    "a = x + y\n",
    "b = 2 * y\n",
    "c = a - b + z\n",
    "\n",
    "norm = np.sqrt(a.var() + b.var())\n",
    "a /= norm\n",
    "b /= norm\n",
    "\n",
    "print(a, b, c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_figs(fig_num, elev, azim):\n",
    "    fig = plt.figure(fig_num, figsize=(4, 3))\n",
    "    plt.clf()\n",
    "    ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)\n",
    "\n",
    "    ax.scatter(a[::10], b[::10], c[::10], c=density[::10], marker='+', alpha=.4)\n",
    "    Y = np.c_[a, b, c]\n",
    "\n",
    "    # Using SciPy's SVD, this would be:\n",
    "    # _, pca_score, V = scipy.linalg.svd(Y, full_matrices=False)\n",
    "\n",
    "    model = LinearRegression()\n",
    "    model.fit(np.c_[a,b], c)\n",
    "    print(model.coef_)\n",
    "    print(model.intercept_)\n",
    "    \n",
    "    pca = PCA(n_components=3)\n",
    "    pca.fit(Y)\n",
    "    pca_score = pca.explained_variance_ratio_\n",
    "    print(\"pca_score\", pca_score)\n",
    "    print(\"pca.components_\", 2*pca.components_)\n",
    "    V = pca.components_\n",
    "\n",
    "    x_pca_axis, y_pca_axis, z_pca_axis = 3 * V.T\n",
    "    \n",
    "    print(\"x_pca_axis\", x_pca_axis)\n",
    "    print(\"y_pca_axis\", y_pca_axis)\n",
    "    print(\"z_pca_axis\", z_pca_axis)\n",
    "    \n",
    "    x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]\n",
    "    y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]\n",
    "    z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]\n",
    "    x_pca_plane.shape = (2, 2)\n",
    "    y_pca_plane.shape = (2, 2)\n",
    "    z_pca_plane.shape = (2, 2)\n",
    "    \n",
    "    print(\"x_pca_plane\", x_pca_plane)\n",
    "    print(\"y_pca_plane\", y_pca_plane)\n",
    "    print(\"z_pca_plane\", z_pca_plane)\n",
    "    \n",
    "    ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)\n",
    "#     ax.scatter(x_pca_plane, y_pca_plane, z_pca_plane)\n",
    "    ax.w_xaxis.set_ticklabels([])\n",
    "    ax.w_yaxis.set_ticklabels([])\n",
    "    ax.w_zaxis.set_ticklabels([])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.21707439 -1.2173145 ]\n",
      "-0.0004483572238646694\n",
      "pca_score [0.72408528 0.27427885 0.00163586]\n",
      "pca.components_ [[-0.67695451 -1.54801209  1.07022949]\n",
      " [-1.42192159 -0.3243452  -1.36855368]\n",
      " [ 1.23283072 -1.2241155  -0.99079244]]\n",
      "x_pca_axis [-1.01543176 -2.13288239  1.84924608]\n",
      "y_pca_axis [-2.32201813 -0.48651779 -1.83617326]\n",
      "z_pca_axis [ 1.60534424 -2.05283052 -1.48618865]\n",
      "x_pca_plane [[-1.01543176 -2.13288239]\n",
      " [ 2.13288239  1.01543176]]\n",
      "y_pca_plane [[-2.32201813 -0.48651779]\n",
      " [ 0.48651779  2.32201813]]\n",
      "z_pca_plane [[ 1.60534424 -2.05283052]\n",
      " [ 2.05283052 -1.60534424]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAADmCAYAAADsvYEoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29eXBb93n3+8XBvu8AAZAEF5HavWmzHVuSE1uOk968aZrkzZt4kt7kve2083aaadreJNNkfG+bppnezmTuJHN7O+lt2jR9k6aJEqexW9uSI0exLIkUKVGUuIj7BkIAwQUbsZ37B/w7Ojg4AAEQIEjg95nhkARxgAMC+OLZHwnLsqBQKJSdhqn3CVAolOaEig+FQqkLVHwoFEpdoOJDoVDqAhUfCoVSF6j4UCiUuiDb4u80D0+hULaLROxCavlQKJS6QMWHQqHUBSo+FAqlLlDxoVAodYGKD4VCqQtUfCgUSl2g4kOhUOoCFR8KhVIXqPhQKJS6QMWHQqHUha3aKyiUpuKv//qv8YUvfAEMI/65zLIswuEwFhcXsbS0xH2/fv06Tp8+jd///d/f4TPeu1DxoTQ9LMuCZVmsrKzg0qVLUKvV0Gq1WFpawtLSEnw+H3w+HzY2NgAAarUaLS0taGlpgdPpREtLCw4dOoTr16/X+ZHsLSRbzHCmjaWUPQvLskin01heXuaEhFgqPp+P+765uQkAMJlMkMvlkMlkOH36NJxOJ1wuF/ddr9dDIpFAIsnvk4xGozh37hz6+/tF/97kiP5DqOVD2VVMTU1hdHQU73//+wteh2VZxONxTkx8Ph8nKkRQ7t+/j3Q6DalUCpvNBqfTyX319PTgqaee4qwXtVrNCYbf78fly5fxkY98pKzz1mg00Ov18Pl8cLlc2/ofNAtUfCi7BpZlkUgk8NJLL0GhUORYKcRSWVlZAQDI5XLO5SHuz2OPPYaWlha4XC7Y7XbIZLKyrRCHwwGPx8MJVzmcOXMGFy5cwIsvvljWcc0KdbsoNYXEU4LBYF6Qlu/+hMNhSCQSaLVazM/P4/nnn0drayvn8hCRsVgsBV2fajE2NgaLxQKbzVbWcW+//Ta+//3v4x/+4R9qdGZ7FtEni4oPJY+33noL7e3t6OjoKHgdlmWRSqWwvLycJypEUJaXl5FIJCCRSGA2m+FwOHKCtMJ4CgBIJBK8/PLLOHDgAHp7e3foEecSCAQQDAaxf//+so5LJBI4ffo0bt68SeM+udCYD2VrWJbF2NgYzp8/j9/8zd/MERN+PCWTyXDxFL6gHDhwAGfPnuUsFaVSWfYb8ZlnnsHU1FSNHuHWmM1mTExMlH2cQqGAy+XC5OQkuru7a3BmjQUVnyaAuD7r6+tFXZ/V1VUAgFQqhd/vRyQS4YTlxIkTnJVit9shlUpr9umu0+kQiUSQyWQK1tvUEqlUCplMhng8DpVKVdaxZ86cwcWLF6n4lAB1u3Ypr732Gk6cOAGz2VzwOizLIpPJIBAIFBQVn8+HaDQKIPumFgZp+fEUs9nMxVO+/e1v49Of/jTnDu00d+/ehdPphMViqcv9T09PQ6FQwO12l3XcwMAAvv3tb+P73/8+QqEQlEoltFptjc5yz0BjPnsFlmXxjW98AyzL4r3vfW/BeEoymYREIoHFYuFEhV+bQkRFp9OVbaUsLi4iHo+jq6urRo+yOJXGXarFxsYGpqencfTo0ZzLM5kMEokENjc3ue/k5/n5eXzpS1/C8vIyOjo6YDab8dJLL+Gpp56qy2PYRVDxqTcsyyIajWJxcTHHQuFbKoFAAOQ5WVtbw9mzZ3OEhR+oVSgUNXN9kskkbty4gVOnTtXk9rcik8ng2rVrOHXq1I4Eb1OpVI6YbG5uYmpqCjabDYlEAolEAkA2IK5QKKBUKrnv5Ge5XA6GYfDZz34WX//613Ho0KGan/cegQacy+XKlStob2+Hx+MpeB0ST1ldXS2Y9fH5fFhbWwMAKJVKziIhgvL4449zv9tsNi6e8s1vfhO/8zu/A41Gs1MPmUMul0OhUCASidTFbWAYBmq1GpFIBDqdrqLbYFkWyWRS1EohP6fTaQDZOI9QUHQ6HWw2GywWC+RyeckiePr0abzxxhtUfLaAWj4FYFkW3/ve93Dr1i188pOfLBhPicViAACDwSAaTyHfjUZj2fUp9XZ9FhcXsbm5ic7Ozl11/5lMpqigJBIJznqUy+U51olQYGSywp+/i4uLSCQSRUsOxLhz5w6+9rWv4fz582U/5gaFWj7AgypavrsjtFT8fj9SqRQAcKX7brcbDocDnZ2dnKXicrmg0Whq5hY4HA709fWhs7OzLnUjdrsdAwMDOyI+pA+LLyDxeBzz8/OIRCLY3NxEMpkEkHV9hCKi1+ths9m4y6uRJbNYLBgeHi5bfA4cOICxsbGKqqSbiaYRn9/+7d/GzZs3wbIs5HI57HZ7ThzloYcewrlz5+ByueBwODgz++///u+5atudRiaTQavVYn19HUajccfvXy6XQy6XIxqNVuz6EbEvZqUIXR8iICqVCjKZDG63GzqdrizXpxqoVCokk8myRYRhGBw5cgSDg4M4duxYDc9wb9M0btfS0hI0Gk3Z/T4bGxuYnJzEww8/XMOzK0wgEEAgEMCBAwfqcv+FXJ90Ol0w6yN0fQoFaMnPxd7YU1NTUCqVZae8q8Xo6ChsNhusVmtZx333u99FOBzGn/7pn9bozPYUze12Wa1WLmNRDnq9nntjKZXKGpxZcSwWC8bHx3ek4I60TPBFJB6PY25uDuFwGIlEgnN9GIbJExWj0ZhzWTXO1263Y2Jiom7iY7VaEQwGyxafM2fO4I//+I+p+BShacQHQMUmu9vtxtLSUtm+fzVgGAYWiwXBYBB2u72i2yCuj5h1Qi7LZDIAsq4eX1DUajVUKhUXNK+kU3w7aLVaxGKxulU7V9pq0dHRgfn5eSQSCSgUihqc2d6nacSHZJoymUzZb56Wlhb09fXB6/XWJfDrcrkwPT2dJz4kQFssnkIQujxarRYWi4W7vJjrw7IsNjY2yu7yrgakKXVlZaUu909aLcq1fCUSCR577DFcu3aNFhkWoGnEB8haEZlMBizLliUiMpkMer0eoVCo5uX+fNeHn/UJhUK4efMmkskkl4ljGCYvhmI0GrnLSNHbdtnJrFeh+/f7/XURHyDr+q6srJQ9JOz06dO4cOECFZ8CNI34EMEhAlSuBePxeDA3N1ex+BQryyffietDCvz4WR+bzQa1Wo22trYdd33IaNHtZL22g8lkwujoaNkfGtXCarVidna2bPE5e/YsPve5z+GrX/0qQqEQZDIZDAZDjc5y79E04kOo1PoxGo0YGRlBMpmEXC7nLhfWpoiJCiBelq/X63MuK+b6mM1mDA8P163g0Ol0Ynl5uS7WD8Mw0Gq1CIfDdWl01ev12NjYyHvNFCt2jEQi+L3f+z0EAgE8+uijsNls+JM/+RM8//zzO37+u5WmSbWTIeGkmA0oHIAuVJYfCAS4mg9+bUqhFLJCoahq/1V/fz8OHToEtVpdldsrh2QyiYGBAZw8eXLH7xvIlkrEYrGai69Yxm9zcxMLCwvQaDTIZDJ5GT+x558871/60pfwsY99DM8++2xNz3uX07ypdr7AksDz/fv3c5oIyQutUFm+Wq2G1+vFvXv3cOzYsR13fYBs4Nvn89XF+qi362Wz2TAwMFCx+PDdXjFrpZjbq1QqYbPZIJFI0NnZWdZzT+Y6N7n4iNIU4iOEZVmMjIygo6ODK8snL7KtArRkPk49Ko6dTif6+vrQ0dFRl9iH0+mE3++vS8kBCZ7zs05iVorweyErRaFQlPXcm81m3LlzJ8flLoWnn34a3/rWt+oWr9rNNIX4CF3LVCoFnU4Hj8dTdjbI4/FgYWGhLuIjk8mg0WiwsbFRl8Clw+HAwMBAzcWnUHA+k8ngxo0bkEqleVaKsBud3zhajTe9SqXiWkHKabWwWCxIJpN1a5HZzTSF+AhJpVLcJ1i5n0gWiwVjY2NIpVJFO6JrhcvlwtLSUl3Eh7hesVis7LgTsVKKFTvyG0eFcRS9Xg+tVouFhQU88sgjdWnYNJlMWF1dLavaWSKR4Mknn8SlS5fwoQ99qIZnt/doCvGZnJyEWq3m6kRIxkoikZQtPhKJBC0tLVheXi4656dWWK1W3Lt3r24VvyTrRawfMStF+L1Q9bRCoeCsFDKMa6vnYnp6usaPsDBWqxUrKysVtVpcvHiRio+AphAfi8XCNQgC+eJTLi6XC0NDQ3URH4Zhal7xK7RS+N9jsRiCwSB8Ph8A8RICnU4Hi8VSUglBuZCCv0pbTbaD2WyuaKvGe97zHnz961+ncR8BTSE+RqMRyWSS20ZAXKZKa35UKhXkcjk2NjZy6k6uvTIAADj5gUer/hj4uFwuzM7Oli0+xEop5vrwrRSh60OslEQigUOHDtUl62W327G0tFQX8SGvmXL7tUg9VyAQqMt571aaQnxYlkVrayvm5ubQ09OTY/lsp+J5YWGBG3Vx7ZUBhHyr3M9AZSIkdqzwMoPBgEgkglQqxdUcFRKTreYPV2KltLS01C3rRYo962VFkCbfcqudn3rqKbz55pv4+Mc/XqMz23s0hfgA2TfMlStX0N3djVQqxaVrK7V+SOylVtPqMpkMrrx8HclUCoGFFaTTKfh8PvQ+2ckFZ69cucIFgcWsFHKZWCwlK2gbFQmkw+HA4OBg3br89Xp93bJHVqsV8/PzFbVa/OxnP6Piw6MpxIdlWUilUq5Bkd8iQYoOyxUfhmG423O5XDj5gUe3tHj4o0LFUslDb45gdXkdADA2Oobl8QA6Hm5DLBaHTC6DSqXitoCm02mMjo5WNClPaKUVEyCxx6RQKCrOelUDm82GQCBQF/ExGAxYX18v6/WSyWTQ1dWFq1ev4uLFi/D5fHjkkUeafsB8U4gPoa2tDcPDw1zMhsAwDNcuUQ4ejwd37tyB0+lEMpnEgae6uVJ8ocCQ2+dbKeS7VquFQqHAujcGsyYMiQQwt5jQu7+3qKixLFvRVk0+Y32TorcNFHclHQ5HTtZrJyGNnvXYCiqRSKDRaBCJRKBWq4vG0EjVNMuy+OIXv4hYLIYf/vCH6OnpqdtkzN1EU4gPyWiRYe/xeDynRqeQ9bNV9WwikUA0GsXVq1ehUqlyXJ9y5uUQnvzQibLiRaTdolwBIII21jcJe6sFId9q2XGqerpecrkcUql028Irhlhfn1BcotEobty4AZVKlVc6YDAYRKc5Xrp0Cd/4xjfQ1dWFz3zmM1U9571KU4gPn/b2dgwPD+dYPhKJBLFYDOPj40UHmiuVSmg0mpx5OcvLy9jY2EBPT09Vzk/45i8mCk6nEzdu3KhIAMjtEcum0HUK3T8JTtfL9bLb7bh//z7a2tpKun6h0gG+lUI+fEhfH3nOVSoVDAYD97wnEgmMjo7i0UfLi5edOXMG//iP/0jF512aQnz4tTx2u52rpOUTiUSgUqlw+PDhsiqXHQ4Hpqam0N3dveNFf3K5HCqVikv5l2u9lBKnKnZbwoLDncRut+POnTuw2+1FRaXYUkDyQcLfNloKZLJhuYWex44dwx/+4R/Sep93aXjxERYRkk+25eXlnHU4qVQKarW6bAEhM5YDgQAcDkdVzplQSmCYtFvc/dW9ilL9pQqV2G06HA7cvHmzquJTTrNoJBLBnTt3clzeStzdSiCtFuUMl5PL5WhtbcXExAT27dtXk/PaSzS8+Ighl8sxNzcHj8fDfQKlUqmKB317PB6Mj49XXXxKwWazYWJiAgy785s1SEyjFNdLbNWOWHAWEB9pIdaGMT4+DqPRWJf/O6m0LneyJRmtSsWnCcWHZVluD/ja2hpMJhOAbMuFWq2uKO2u0+mQSqWqHgAtxS1iGAYmkwmOXgcmrs4WvW4l8APT5Hdy+yzLwmq1cuNlxURFOG+6lOBsqdjtdiwsLNRNfCrpMzt79iy++c1v4nd/93erf1J7jIYXH7FxGjKZDO3t7ZidneXEh3S6V7rhwuPxYHFxseqT9koREpfLhfn5+arcn7AWaXV1Fevra2D9KaRSKaxEgsi8k50KKZFIIJVKEYlEuBXGarUaJpMpx+2pVXzDaDTi7t27dYmhkFEd5bZaHDlyBLdv365bY/BuouHFRwgpMDSbzRgZGeFePNvt93I6nbh27Vpd9qqTloMTz58oGOPgd59fe2UAqWQS+5/qLhqcvffODGKrcVjdFqz7N9B7ohuP/8axvOBsf38/2traqp723gqJRAKDwZBjwe4kxPVqaWkp+RipVIre3l4MDw/j6NGjNTy73U/Di4/Q8uH3dZH+rM7OzhzLp5J+L6lUCpPJVNHIhWL885//BADw4lc+kvc3fnBWp9NhbGwMGo2m4Iycqevz2d3roWzF9L13ZnDyA48W7OuKz6cRkmWD2J0HvJwVVqjg0Ov1Vu1xlwpJuddDfKxWKxYWFsoSH+BB3IeKT4MzPj4Ou90OnU4HADlDwNxuN65du4aOjg4kk0nu8kqtH4/Hg6mpKU58Km0wTafTuPLzPoxen8DyTHZo/bf+6O/x3v/+hOisaYVCgdFfTyISieDcJ58pOCNn/D/nkUISjtZsZ7XJZCr6phWLOQkzcGN9k+h6pB1yF3LEZ6c6/C0WS0VjLqqBwWCoyO07e/YsXnrpJXz+85+v4dntfhpefPR6PWZnZ7k+Gn5fl1wuh9FoRDAYzJtMWInrZDAYOKtDGAcgK4uLZXz4wdmFhUVshMNIpVJgGAkUSiU8Ho9ocPbaKwPQKwwILq5g8toc3vNf8jdMXHtlgAsa359fAQA8/9lntnxMpYhH1l3NcAH3cnrHtotMJoNcLq9LsSPDMNBoNIhGo9BqtTl/E8bOyPMcjUbx0ksvYWBgACdOnEA6ncZf/uVf4v3vf/+OnvtuoOHFx+FwYHR0lBMXoci0t7fj3r17AB4ITimrla+9MgCWZXHs+YdyRESpVOLWrVuY7ltAyLeKVCrbANrzREde+rhYcJZsqCnmdgkxGI0IhUJbXq/3eOGgeClFh+Q6ALj2jHQgBZutvq5Xe3t7ze+LlAyQ55thGIyOjnJ9XmIzkfjPu81mwx/8wR/g7/7u7/BHf/RHeOKJJ2p+zruVht/blUqlcO/ePUilUrS3t2NqagpKpRJut5u7ztWrV5FMJnPW2rIsi3A4jHA4nPNiu3XxLpLJJDYCEQCAwa7D0WcO5MRNZmdnkfJJEFmJQiqVweIyVe3Tf2LBj9BGFN1uO8yGB5+2114ZwMjVcTBmFi9+/r+KHruVsPAtFnOLqeh1hdfXWbVQuCU4ceJESfdVTeLxOIaHhyvq8AfKt0qFu7n8fj8OHDhQMHYmxj/90z9hdXUVX/ziFys65z1Gc+7tYlkWHo8HfX19aGtrQzKZ5OI/hNbWVoyOjuZcJpFIcO/ePcjlcuh0Oq5yNuhZh0wmw5ohO/rC3GLCQw89lHPs+vo6PIc9GL8yDaC6b8Cv/eMv8Pr1OwAAk06Dbo8dnW4bFBspGOUMWF8Yb/zoMgxa7Zb3u12BEMaE+vv7OddrJ0SHoFKpkE6nc1zqQpXShYasCfu5+FYpGSFSyAr2+/0wGAxlpc7PnDmDz3/+880iPqI0vPgA4LYfrKysiG6dsNlsuHPnTl7tRSqVQm9vL5RKJffCe+rDpwAUf+OSLFq134DrkRguDTwQydVwFP2jM+gfncm53v87Og6jQoGHrvehy21Hl9uGLo8dXce98NjKCzAXQyzrdeGHb8Fut9dMfArNQ8pkMujv78+Zyy2slC5lyFolGI1GrK2twWw2l3xMe3s7lpaWcvaQNRsNLz7khdje3o7JyUkwDJO3+I1lWSgUCvj9/py0KX/FjpBiby6TyYTR0dG8ve7b5T+u3kYitfXcIRbAaiKBtwbH8NbgWM7flAoZnFotXEYDrDIFnDotRpfv4zc++hTMek1ZoiMMKs/eWML8xAJkaXlOXKiYmzfWN4meY514+H2HilopYuupyXej0Qi1Wo1AIICjR4/uePGe1WpFMBgsS3wkEgmOHTuGd955B2fOnKnh2e1eGl58CEajEZubm5yJzYcsEZybm8sRHzIBsdyaH7Jex+fzlTzyoRR+fvnmtm9jM5HCbGINs6G1BxcOAl9++XWY9Rp0ubNuXLfHkbWY3A50uKxQKbYWUWJJJJNJLtuTTCYRDAbzhGXozRGMvzOFzUgCk6OTuPyLK3jhfzxTcYMoy7KYnZ3d7r+nIsxmc0WtFqdPn8bFixep+DQq/IB6e3s7RkdH89wu0tcVjUYRDodzYkKV1vy43W4MDg5WTXxC6xFcvjVeldsqeB8b4m4cw0jgtprQ5bGjo8WCdocZHqsB66koDAoZzD06DA8PQ9bCQqGT4f66H6GxFcTXNiGVyRAKhfDoc0egVCq5fq7wVBzrUzFE16PQGDXQmbTbGisqkUgqcn+qAV90y7F0z5w5g+9+97u1O7FdTkOLjzCT19LSgtu3b+eJCL/fa25uDgcPHuT+VmnFs0KhgEqlwvr6elW2i776zm2k0plt304lZDIs5u+HMH8/hLcEf1PKpbD98iqMrBQ93hY8+Z5DyMTW4MqYkFqNAcgG5YUp+Cc+dAJSmQxjfZPoPd5VlRgRSbnvtPgAD1otnE5nycc4nU5sbGwgEonk1Qk1Aw0tPkIYhoFMJoPP58tJtZPYjt1ux/j4OLeShn9cpRXPCwsLVRGff/7Zr7Z9G7VgM5nGwtoGFgDcGVrFz4ZGuL/plUq4jXo8dKQD134cQrfbjs533TqlXLZtwREGvC0WCyYmJrZ1m5VitVqxuLhYlvhIJBKcOnUKv/71r3Hu3Lkant3upKnEB8haJHNzczniQ1or+LEap9O57Ypnstd9u+t1Xvu3tzG85K/4+HqxsbmJUf8mRi8Gci5nGAk8NhM63XZ0e+zoctsRHxxDl9sOj91U0v9arIpaKpVCpVLVxZKotNXizJkzuHDhAhWfRkPodqXTaa4cn+8OkSmGQLbm58aNG7BYLJz4lFLxLIZEIuFGjfLFTniOYmX4iUQCg28MI5VK4Vczi8hUsNZ5t5LJsJjzhzDnD+Vl41QKORw6LdxGPY4/1sMrFXDArN96Q6rNZsP9+/d3XHzIjKhYLFZ0kyv/uV5ZWcHNmzfx85//HKFQCMvLy/irv/qrHLe/kWlo8RFCAoJkls+RI0e4y8naY1Jgtrq6mmP5lOt6ZTIZbG5uQq/XY3x8nPudLzCkDF84X1ihUGDq+jxkGQVUcjXGw7Ea/Dd2J/FEErMrq5hdWcU7U3M5fzPrSVGlHd1uO9KxMNxGPc4+96A73G63Y2hoaEfnSpOWC41Gg+np6ZweP/IlViqQTqeh0+mgVqvxiU98Aj09PWUvI9zLNLT4FBokZrVaMTY2xomRsJ6H1ATxP8FI4DkejxdcrSJWhq9QKJDJZJsuDQYDhi+NQS6X44n/5Tjnil17ZQBxpPHQBw5w97dg9CMTYxGKxTHiu1/Lf9OeIbQRRd/IDPpGcrNxXzj/Kjy2bDauy22HNBnFukSD3jZXyW6cEOFee7H2CyIo5LlmGIYb5k/2s2/VcnH48GEEAgGsr69XtSxjL9DQ4iNEOMtnfn6em+XDt3LMZrNop3IikcC1a9e43ebESjIajdyLTKwM3+/3IxQKweFwYFq7AAA5wiMc/A48qDa+PDSy9xvsagzfjbs0kHXjvvNa9n+pUsjR4bKi2+1Ap9uGdocJbXYT3BY9VHImT1z484+E1dEajabolEaWZfHOO+/A4/GUVeh49uxZXLhwAR/+8Ier+F/Z/TS0+JCVuvwh8cTC4c/yERunYTKZEI1Gc24vlUrBYrHg8OHDZX2a2mw2/Mc/X8DqvQhWl9cwdXsOY32TePGrv5VzPbFZyX9x6XJFj52SJZ5IYmTGh5EZX97fTFoV2p0WdLps6PbYsa/NiZ42Lzrddijl5b81SK3R+vp6WcPNnnjiCfz5n/95063UaWjxWVxcRCqVgs1mA4CcgWFyuRwmkwmBQEC0OMxgMGB6ejrnBUFEqpK97gaDAaurqxA2+PL7qcwtppwlfnP+FQyM1adqtxlYjcSxOrmIW5OLOZczjAStdjM63bZ3A97Zr26PHW5bcTeOtFqUIz46nQ4qlQr379+vyzD8etHQ4uP1ejEyMpIjPsLYztjYmOgwb5ZlodFoEAwGueP5FlK5AnTuk8/gp9/5BTKrDCRsdg4OsW749S782pX/5/yblT1wyrbIZFjMLq9gdnmFc+MIKoU8K0ou0obyrjh57DDpNLBYLBXtkX/66adx8eJFfOITn6jmQ9nVNLT46PV6ZDIZRKNRaDQapFKpnA5ivV6PdDotmkJPpVJwOByYm5vjxIek6iupeCbBa+9RD6KhGKaH5nB/fqVooV01erko1SWeSOLu9BLuTi/l/c1i0KLbbYdeARybXsPnfuNpaNWldayfOXMG58+fbyrxafjdHe3t7ZiZyWZHxNyr9vZ2bq4Ln1QqBYPBgEQigXg8zl1GxAfIz6ZtxTMfewquozbcn18BK3lg/fAhltBP/+clDE0ulHX7lPqysh7B9ZFpXLw1jf/v3y9DrSy9z+vkyZN45513yn5N7WUaWnxYloXD4eDm+IiJj81mQyqV4tKmBCI0ra2tmJuby7mMFB2W+0JxOBy4f/8+eo51ovNI4bTqtVcGcOHmWMG/U3Y/p/a3lpXxUqlUsNlsdevMrwcNLT5ANgPhcrm44LOwoz2TyUCpVGJpKdeMJtdtaWnB8vIyMplMzvGVZCWkUiksFgu6TrbB3GKCuaXweNXri4uil1P2Bo94bWUfQ0ZsNAsNLT7EMmltbcX8/DwSiYToLB+tVov5+fkcS4YIjVQqhc1mg9/vzxEfhmEqsn5Is2mxc9YdtGNhI1zW7VJ2D3qNCsf2t+aVavAhFe8bGxsIBALo7+/H3NwcvvnNb+JTn/oUnn32Wa7mqFFp2IAzXxTIXvCVlZU8y4cEoWUyWc7mS77QtLW14c6dO1Cr1RX1e/H7t5LJJCKRCFr2t4BlWQwPD3MFbuScf/TrkaK3R9ndPHW0G3qtFpOTk9DpdAWLGPlzo1mWxeHDh/H222/jy1/+MsqmrVkAACAASURBVFwuV95rtdFo7EfHw+v1wufz5QkFqf1pbW3N2d3Osizns2u1WkgkEmxubuaN2iC7mPh9PPwXm3CNilKphE6nQzgchsfj4apn+bu4/vfvCafmUPYSj3TYwLIsNjY2YLVauQFqwiWOQo4fP4433ngDcrkcFotlh89652lY8RG6Q3q9HizLcml3Aql6NplMGBkZKTjQu62tDSMjI3mV0CMjI5BKpdBoNJyFxS/LFws6plIpXL9+HVarFd8ne7nerXb+4fcuYnKR9nLtVZRyGT79X56FVqXEO++8A6fTWVbgmaxS7u3treFZ7g4aOubDhwyJF2YT+BksEhvic/UXN3D1Fzdgt9uRSCTQ9x83cfUXN3KO7+npQVdXF1pbW2G322E0GqFSqQq+6GQyGYxGIy7+668wOTiNycFpXHtlAP/8f/4YP7nQV/0HT9kxnjy6Dzq1ChKJBAaDAevr62Udf/bs2aYJOjeN+KRSKahUKm41MoGffne5XPD5fEilUnmm8fVXBzHTv4jZ8XmEfKucAPFbLsrB4/Hgwr/8CusrYayvhPHLH7yNxXs+3N7YeuMoZffy/lNHuJ+tVitWVlbKOv7QoUPcGqdGp2HFRygGyWQSCoUCLpcrJ63ODyzLZDJYLBYsLy9DJpPh6i9uIORbRci3irG+CUilUqyu5ouDVCotW3xGLk8gk0pDrpRBoZLDva8Fnvf2Yj2T2vpgyq6EYSR47uSDIfgWiwXBYLCs25BKpThw4ACGhoaqfXq7joYVHyEktkOKBolYiO1un5uby8s09B7vxv73dMHSYobSoMCpDz7G/a2SimeJRIKnP/44Wo+04NCT+/HiV38LE2x8Ow+RUmeO9XphN+m530kWi29plwIZrdroNKz4iFk+JONEtpeSy/m1P2SGD8uyOPXBx7hiQCI2z33yLGwHjDm3TQaNlSo+ZIaPJMHA0KrFp77yEVz9xQ0a79njnDt1OO8ys9lctut15swZvPlm4zcVN2y2SwhfZLxeLyYmJmC1WkWrnu12O3y+7PwXIjqknsdisWB0dDSvYLGUMavCPi4y8PyX//ZrDM8sIxgpXJRG2f3w4z0Eq9UKv99f1qiM7u5uTE1Nib42G4mGtHxYluX6sQj8J5I0jMZiMdEnmBSG8c1lflaMVCkLR3EUEx1i7ZB5PcSiet9/PY1gMIjrC/ld0pS9w/72FnS48lsqTCbTu3OcSodhGDz88MPo7++v1untShpSfCQSCZaXl7G29mAlsNgsn9nZ2ZwBY4R0Og2DwYBFXn8VGacBZKcgCnvFym02Jd3rZrMZ7ocduLG8XPHjpdQfMasHeDDfORYrbwkAqfdpZBpSfACgo6ODG6UB5IuP0+lEIBDIqWQmpFIpWK3WnH4v/iJBuVwOo9GIQCCQV/EM5MebTn7gUZx44REMXLyN/jduoeeJDiwvL2N2dhYTExO4OxfAaowGm/cyz4vEewhkuqEYmUyG+0qn09z0hSeffBLf+973anW6u4KGdSjNZnNOxbLQvWIYBi0tLXnuGfBgj5dOp0MoFILFYhHNio2MjHD7vgBwLRjz8/PchotkMolMJoPxK9NIYhOMTIpf/tuvIZPJ8Ni5ozAajRhaKK8QjbK7aLWbcbS7lfudX6PDsiyMRiOmp6e5uA/LsltayG63W3TOVCPRsOIjkUg416qnp0d0lo/H48HExERekJgIjdfrxeTkpKj4GAwGJJNJqFSqnNtcXV1FMpmEx+PJabG499oCtEodACCytAlgEx6PB6l0Gm8OjtfuH0GpOc8eP4DNzc2CgqJUKhGLxZBOp7lpCOQLQM7PBJZlodPpGnqofMOKD8uycLlcePvtt9HV1ZW3mwvIuk9SqZSzbghEaIxGIzY3NxGPx0VXHlut1rxgIokXGQyGnBdN7/EuhFcj3M9ANhA9MLeIlfVIVR87ZWdpyyjQ/5+3cPz9D3PPudCVJ6t4SoXcTiOLT0PGfMgnEMMwcDqdWFpaEg0sJ5NJaLXanNgQkD9OY25uTjQrZjAYEI1Gc8xsEpgWi/s89uxRPPbsUW6AWMi3ijeH71XnQVPqgk4hx2F3tnlUKpWCYZiyGkmLoVKpyg5U7yUaUnz4kIplMcuFxHZI2p3Ad9HIJEMx8cpkMtDpdPD7/dxl/KzYVn59KpPBgC9/nxRl73Dc24rHP/gYTrzwyJbXLbcFx2AwIBRq3F6/hhQf/pOsVCqh1WrzZjQDDywcIlDCy4FsIaDdbs/b3Q5khcZqteYcy99wwT8Pfp0PKTacScYQTdJerr3Mpz7ydEnXI0Wo5WA0GsuuEdpLNKT4CPF6vaKZA2LNOJ1O3L9/nxMooZXU3t4uKj5kBKtEIkE4HOYuK7ThYnpoDtNDD4RqJE6zXHsZjUqBpx/uKem6UqlU9AOwGFR89iBC81any2aZxNYfy+XynNgQgR/kU6vVXBqdDxEpEhfiXya8DQBgJdkvAEik0nj9+p1tPEpKvTn9cA9UitLW48hksrIbTA0GQ06hbKPRkOIjJJlMQqPRFBwkBjwILBfyyzUaDZYFVcjkeLvdjmAwiHQ6zYnPVhXPN+YWEY5tiv6Nsjc4d7JwYaEQavnk07Cpdj6pVIpbfcwXHP4GUxIbKhTgYxgGkUgkJxjNj++QOUF8y4ffbMpfkXPyA4/iH/6v4Vo+ZEqNkUkZvO/YgdKvX6Hl08ji05CWT6FBYm63O2dtjVi/18zMjGiqNJ1O5x3Pb7kgI1iF/V78wDPp54rGE7jQR12uvcxBpx1GnWbrK74LtXzyaUjxEcIfJLawsFBwkJjJZEI8Hhct6kqn02htbcXi4qLo8UqlEiqVKm8EKz/wTHp4/vPqbcQ2G3snU6Nz0tu69ZV4lCs+qVQKmUwG9+7dw+uvv46XX3653FPc9TSk21VokBjZUnH//n04HA7RlouWlpa8IfIEhULBNZTa7fa8rFhbWxsGBwdzxIfs9uKb3D+/PFiNh0mpExKJBP/bb58r6xiZTIZEIoF4PJ6zZons8iI/p1IpsCyLK1eu4Kc//SkikQgcDgc6Oztr9GjqR0OKjxC+yLS3t+Pu3btwOByiVcsmkwmTk5OiRYnk+LGxMdjt9rx5PhaLBel0GolEIqeUvu8/bmLs+gR6jnXi4DOHcOkm7eXayzzU7UGLNTvNMpVK5QmI8It0q7Msi1AoBLlczu3xUqvVMBqN3O9kZtTJkydx7tw5/M3f/A2+8Y1v1PkR14aGEx+ym4vf8MkPLOt0OjAMg42NjYK72zUaDXw+HzweT97t6/V6pNNpLm0vtHIUCgUWFxfR0dHBbbgYvHAbs3cXMHVrDm9NzCBBCwv3NIc9Jly7dg1A1p3iL35UKBTQ6XTcgkCFQgGpVIpAIIBwOFyWBdPoqfaGE590Oo3+/n488cQTnFUidK+8Xi9mZmYKtlyYzWbMzc3B7XaLrkTm1/XwYVkWcrkcCwsL8Hq9ouf3+iBdhbzX+cyH3of9XldZx1SS7TKZTA0tPg0XcJbJZLDZbNwMZiBffKxWK9bX10X3rJP9Xmq1mss0CC0kUhEtjC2R1LvRaEQwGMSpDz6GUx98DDqzFgq1AvpWI+5t0KrmvUyX21628ACVZbvUanVeYWwj0XDiw7IsvF4vZmdnC2a1yHbSZDI/48Sf5UOKEvnNosCDbnnhJxk5VmgZsSwLh9eGoFmOdJnNhZTdxTneXq5yqMTy4Y/VaEQaTnyA7CgCjUbDFQyKZbXcbjc3ZZAPf5ZPLBZDPB7PqechuFyuPPEibhx/DhAAfPLPPoJH33cEg0G6g32v054ufSYPkN10e/3VwYosn0ad40NoSPEBHsR1gHzLB3iwu12sZUIul0MikRSd5SOVSt/dYPqgCIwfQxJaP11P7sOwj4rPXsakVGKf3bL1FfFAdMgkg4HXb5ctPkD2dVauxbRXaDjxISaq0WhEMpkUzUoBDzrShf1e/Lk9LpcLfr8fiURCdJyGTqfLOZ4vPmQOUDqdBsuyePWd20g3wf7tRma/zlixNcIwTEXio9frGzbo3HDiw6ejowPT09OifyPzl5VKZY71wnexGIaBw+FAMBgsOIgsEolw4zr4x0qlUthsNm7Q2L+/favaD4+ywzze3Y6QbxXXX926SPTEC4/gxAuPcPvZ+L195dDILRYNJz784JzdbkcoFBIN2PEDy/wxqul0Oic+1NbWhkAgIGr5yGQyrqeLfxn/2Pn5edxf3cD1kelqPURKHdAq5NhvLc3l4kNEqFKMRmPDTjNsOPHhI5FI4HK5RH1mIj4mk4kLLAPIG5eqUqkglUrzhpGR40k3O8uyeXVDZO/7+V/2I5NpzIxFs/Bomxt2twXmFtO2xKRcGrnQsOHER2jl2O120awWyYDxV+wAyGuZALIvgJWVlZzLiNDIZDJYLBb4/X7RosX29nb85Jd91Xp4lDpx0tu6bSsGQEWjVKnls0dhWRZqtTqn6BDIzYC1tLTkjFEVBhVJUyB/kiH/eCJeYlmxBGQYWxTfVknZGygVMvz3zz6/7duRyWRlB531ej3W1xuzMLXhxSeZTMJoNOYUHQK54kG2l/J3s/NJp9N52035Vg6Z4xyLxfIsn3//9U00aI1Y0/DU0X3QqPLre0g6vVQqqfUxmUwNa/k0XG+X2DgNtVqNdDqN1dVVmM1m7nJhYLmvr080OJ1Op+F0OjE0NISuri4wDMOl6gnt7e0YHx9HS0tLzrEvX75ZzYdHqQPPFxmXOt4/CQAluWOlig+ZjDAzM4PJyUmMjo5yMccvf/nLJZ717qfhxEcIsXA6OjowOTnJiY/QRSLdyGImbiqVgkKhgN1uh8/ng9vtzqt6ttvtGBoaynHZppcCuHkvvwGVsneQMgyePZHbUkGsnYE3bmNuZAGTg9lsaSEBYlkWyWSSG6kRDodFZ/mQHe6kU35oaAhzc3MIh8Po7e0VnbKwl2ko8SHTAvkkk0kolcqcokONRiMan3G73QgEAnm3y+/ZunXrFtxud15aXSKRQKVSIRgMcquXX/4VHRq21zl2wAuLQZt3eSaTgW9qGbFwDMlkEjd+eQuOQ+acOT5kho9EIoFMJuP2uRsMBsjl8pw5PmR1N5+jR4/i0KFD+M53voOPfvSjO/WQd4yGEh8AXKc6sUCSySS3OocUHR46dEi034vMAIpEIjkuFT+zpVAosLa2JtrvpVQq4ff7sW/fPkgkErx8eaCWD5WyAxzrcnJuDxkMBmvWIjJ4tIjGYrB5zGg/4oZKpYLBYODEhGRTCVNTU9DpdLDb7SXffyMPkW8o8SGfMPw5ynwLx2634969e0gmk6KWTyqVgsFgwOzsLA4ePJh320DuLCCxfjGyASMYTeLO9BIoe5uOjApvf+8Gek904/HfOJbzgXP8+HHOBSsl5lPpTB8qPnsE0ohHzF2+hSORSODxeDA/Pw+WZfPqeVKpFHQ6HYLBoKhlBABmsxmjo6MAICpe+/btw8zMDN64Q3ew73W8FhPSi1EkYylEVqK48doQgFyhKafup5Jsl8FgwMbGRlnH7BUaLtXOMAykUmnOLB++iHg8npwNFHzIdT0eT86KHD6k2z0ej+e5XSzLchswfvYr6nLtdZxhFuFQFNH1GKZvZxMH4/2TZaXX+VRi+VQiWHuFhhMf4IFFQrIMfAuFTDos1nLBX7EjNu2QzPIREzCJRIK4RIWJBTo+Y69zqrMNHUdaoTFooDNnY4A2j4VrLi1XhLYjJI04UKwhxYcs6wOyAWihhdLa2sqlNvkQy4ffMlFolo9MJsubBUR4e0R89Q5l79DmsOALX3sR5hYTdGYNeo51bfs2K6lwbuRphg0X8wEeBJ6FzaAEmUwGmUyWU3QIPNjpDmQDy7dv34Zer88THyBbFzQ7OwuXK3+e77//mhYW7nWeP5Wt7RH2c5UTYBZS6WAwjUaDaDTKZW0bhYYUHyAb+yk0+IkElqenp3PEh2/laDQayGQyrK+vi8Z2GIbh0u5Go5H7ZLo5PodpH+3l2uucO5GtahaKzXYaSyt1uwwGA0KhUMOJT0O6XQC4Wh8xASItF/xJh0D+PB6v14uFhQXRWT5SqTSnG55c9jNa27Pn0SsUYKdXc8agVhpk5lNJwBlo3IFiDSs+QNY1OnXqVMHYjnDSoTA4bbFYEIlE8gSMCI3FYuFK5UnR4c9pL9eep0uhxcSNqarf7nYsHyo+ewwS+xETH5lMBrvdjtXVVW4LhdiKHavVmldnQa5HVvDMzc0hnU5jfDGEhfuN2YHcTOxX62HzZFtkyBjUagwQYximosAxtXz2KPy0O4E/SIwUHQLiWy70ej0ikUhOzxh/nIbL5YLP50MqlcKvhqdr/GgotUYOCTrUD1prqjFAbLs06jTDhhcfMeuHLzKk6DCTyYhWPWcyGRgMhpxhZPzjZTIZrFYrlpf9+NXtGVD2Nsc7W2Gy6Hd8XGoxSMC50Wh48QHAWSliG0zF1ivzSaVScDqdOcPIhONS29vb8cu+2whuNO5q22bhkNGMcCjKzekByh8athWlul4sy2JtbQ2JRAK3bt3C+fPn8YMf/KBq51FvGjbVzodhGDAMw1UrC/u2vF4vBgcHC7ZcqNVqbne72WzOc880Gg2ujIq3Y1D2DgwAYzCBjoc6YG4xcYIT8mXjLddfHdy2NSSRSBCNRpFOp3Nm+Qh/Jq/FS5cu4cKFC4jH42hvb4fX693W/e8mmkJ8AEAul3PzVITioVKpoFarRf1q/oodUhcktHzS6QyujdMO9r2OV6WDzWrgXK5SrB3SwlNMSBKJBOfWx2IxjI+PQ6lUcrN8lEol9Ho997tMJuPc/5MnT+LMmTP4wQ9+gK985Su1/hfsKE0jPqTlQkx8gGzLhdggMSI0Op0Om5ubiMViec2ql2+NIxSO1fwxUGrLfo0e+x7rxKPPHUEsFkPvk51IJBJYeW0FqWQSxm4NhoeHudIKApndwx8Mxp/ro1AouA+r27dvo6uri6ukL4VGDTg3lfiQlguxcRnkxUAmHRLEtlTIZDJu8BgAvHyZTixsBHptWoyPj0PWwuaIyeEzvTnCQqyTSlYnV1Lr06ip9qYRH+BBy4XYfi3ScjEzM5MzSIwvPk6nE5OTk7BYLDAajQCAZCqNV6/QVch7nTaNFk8/fbzmGa5KqpyNRmNDrs9pimwXgZ92F35qEfFZWVnhig6B3MwWwzBwuVw5/V5vDY5ilbpce54zh/eJCk+1M12VWD4qlQqxWOO9xppKfIDsk9/d3S1a9axQKHJ2rxP4QtXa2opwOMyJz8/okPiG4HOffq7qQiNGJeJDehQbbaxGU7ldQPaJJKtv+KJC3CuPx4OrV6/C6/XmFRwC4Pz99fV1KNUa/OfV2zt5+pQa0GoyIDA4n5NSJwy8McT9vJVLVsq4jUrcrkpiS3uBprN8gMItF2RImM1mKzgoDMgK0NLSEsKxOP7HR9+HF59/HI92u9DhskIhkxY8jrI7OeFtxXj/JDcqFciOSx3vn0R0PYboeiyn6FAMYQd8IQuq0uZSqVSaEw5oBJrO8gHEt1yk02kug0WKDh0OR9FPHSUD/MFH3wcAuH79Oh5++GHI5XIsBdfw6oW3oLU4MOdfwcL9Vcz5VzDnD8EXXEMqnSl4m5Sd55S3FZY0g3AoisDCCswtJm5yYTiUrVovNMlwp0apks52h8NR9rG7laYUHyB/ywU/q0WKDoPBoOgUQ5Zl0dHRgdnZWRw6lJ14x+90d9tMePbxh5FOp/Hhpx/Jcd9S6TSWgmuY94cw5w9h3r/C/Tw5t4yVSAyN5dnvblxWI7rtFs7lIghdJzFX6l/+4icIhyLoONLGdb8Xui5hOzN91tbWqPg0AmItF3yh6ejowNjYmKj4AIDNZsP4+DhXMyTMoJHYkXDMKiORwGHUwqSWY7/bjESiNacStv/121he34Cx1wH/WgSB9RgCGzHcX4vCFwrTYsYqc+7kIa6audRmUmLt8DdbEOHZiu1OM2wkmlZ8gOItF0ajUXQrKYE/jqOzs5O7jCCTyWA2m3Hr1i3u9knWghSq8QvZtFotFAoFnC86sbS0hMOHD4ve7+Wf9+F+OILJWT+C0Rg2mAzSWjnm/CuYvx9CiDa3lgUZl1pIdIqJUceRVkzfnofOrCm5PohOM3xAU4vPVi0XDocDS0uFe7b4mTExenp6EA6HoVQqIZVKS8paiM2M5qOUydBqMkIbzzpnwj6kg2cPYm45K0Rz77p0Q0PTWFpdRyAaQ7yCF36jYtSqcOpwZ9nHEaG5/uogAgsrZW22qNTyIW5XI9H04kOyCIUGic3MzOS0Y/D3eJFZPn6/X/T2GYaBWq0WTdkXIpPJFL0+/4XP/52gUytxsMOFgx257h65fs9TvTnClI07hXBv3gf/agTxRGNlVIrxniNdkBUR+q2oZNBYpeKj1+up5dNokMCzmPhkMhkYjcYc10p4vfb2dgwNDUGMSorCthIfgnCdy1ZjH/iXmfQaHOny5Pz97t27aG1txdXXh9F/fQwpnRyBaAxhSRr+jQjWMyksBFaRSDaO5dQp16Dv9ds4/tyRHbvP7Vg+VHwajK1aLiwWCxYXF7miQ2FfGFmxU4kfL0ap4lNtiFAa1Sq0aXWwuR/MMAay4sWyLP6Pz3wLS6F1hCVpxBUSbCCDiDSDTQWDlXgcmT1ShStlAWZyAyH7+o4KUDHXO5PJiI7jeOutt/CjH/0I8/Pz+PGPfwydToc33nhjR863ljS9+ADg1uAIBYi0XJCiQ5fLJWoheTwe3L17V/S2y61OTafTZYsPP+azncbIwddvI7a2CZvHwsUyhG6e225COhiFPpxC7/FOdBxpQ2BhBTaPBZNDc5ia9SOlk2FTKcHyegSh5CZiMmCdTSGcSe+aFgEXFJBVWDhc6v86k8kglUrlicnm5iZGRkbyZv0wDJOXjJDL5Th69CgSiQQGBwfxt3/7t1AqlZWd+C6Dig+yAuH1ekVbLtRqNVd02NLSItoRr9PpkMlksLm5mfPCqNTtKhZwLsR2u7FZlgV4j50vPHw6jrTBPxOA0W7Ao88eBQCc+1/P4vqrgzAsrKBb4kR0PQqNTg1dWysW7y3DYNWj40grdHY9ApEoJmeXMTq+gKgMWGeS2EilEYjEsBqLb+sxlMOzx3uxv7MNZquhZKuHJCbi8ThSqRT8fr/o8DD+60hMTKRSKdxud96sn0KQco2+vr6GER6Aig8HcZ341g+xckjR4erqat5iQSArGFqtFrOzs+jp6dnWedTT7Tr2/EMYvjQGoHCR3fVXB/H8554R/TuxwMb7J0WtJv7v7zmS7SIfHh7GytgGVCoVlheCCMZi2FQyULebs4Hx5RAXIK9WGYGUYfDfPv4UVLJsuUQ0Gi04gZCICXldTF6bQ2x1EzKZFNFIBI+eOwqdTpc3OKyYxbuwsACdTlfW80yzXQ0MyXyl0+k88QGyRYeTk5Nwu915n1TpdBp6vR7Ly8vo6uqqyHIhZDKZgoWNtYRYaUJRKXddsFgGqJiQzc3O4dynnoFOp8O//MVPIAfwmT/7iOhth2Ob72bospm6vmuj8IcjCEsymFkMIFpipu5AqxXB5ezGkmAwmDcojNRciQ4OC8i54L7ZaYLb7S7pPvnIZLKy3WuTyUQDzo0MeVGQTzm+i2U0GpFMJhGJRPLEgYiU05ktEGxtbQVQ22xXLRB+WguzaMD23TvhbW8Ewxh4bQgqtVr0evxB68lkEkYFC41Th9jQDN7ntqP78WMYvjSGyeUU1jObSGqkMPTYsZpIYnL6PtaSSUQkGWwgg8S7Wabfet8p9Pb2Yn19HV1dpdfoANWJr5EPOeE0zWLo9fq85ZV7HSo+PPhFh0B+Wt3r9WJmZibv0464Yq2trbhx4wY8Hk/FYxDq6XbVcnQDecMee/4hJJNJxGIxRCJRJJNJ3Lw8DAkk0Nm1SKfT+J//97+h+/F2riJcaJmMXJ5AOspCLlMgMpPAvu5uhGfi0K3HoDFo8OjJIzjxwiPv9l5F0XGkFeYWEzqf3Id5fwidbhvSm7GcRZDlsF0BriTdXum2090MFR8e/DnPwmZTIFvxfOfOnbzjiIWkVCqh1WoRCoVgsVgqOofdJD6lfMqT7Q1isZJEIoHbb44gmUphIxAGANwbv4dDZ3pgP5h1IzQZNVq8Tmg0GoSDUUilUtg9Fpw8+VjBc/WZg8Dmg4AuAOjMWujM2pxY0yf/7CN5528z6gAAq8nNisVnu2ynNKPWHxI7CRUfAWTOM8uyeZktiUQCnU6HUCjEuVZA7qhVr9eLyclJmM3mij6pKkm1VwPyoiYZHSIm3mNuJBIJTE1N5YgKefOI9arJ5XKo1WoMXxpDdDmB/Sf2YU23junbc2DlDB5++GEAwMGDBzEyMgK32w2DwZAznqLYjix+YNvcYkLIt4pwKCraY1XoNkhTcT2odJohUHk2dDdCxUcAsX7I4Cbhp4xGo0EwGMyxUEhKHngQG4pGoxVtOKj2i4uIaKFsDvl5fX0dAwMDYBgGMpksr/G1nIwOEZHN9QRavA4wDIPAwgoACWweS46w8MeDCrNjxRCKCnGtSmWviQ+Qfe1FIhEYDIYanNXOQ8VHBP6sHyGZTAZms5krOgTE1yfPzs6iu7u77Psuxe0iYlIsPcyfeieVSvPiJhqNBiaTift9aGgIR48eLSsIWg49x7ryZuZcf3UQCwsLaPlYS85lpW4I3arPrRj1FJ/tdrZT8WlgSNpdjHQ6jfb2doyNjaGlpUV0FY/T6cTExATa29tLejPzMzqxWAzBYBArKyt5okJgGCZPTNRqNYxGY47rU67VVa1YQrH6HgKX7QpEMPD6bZz96Huqct+lshctHzLTp729vQZntfNQ8SmAVCrF/v37RVsutFpt0d3tZEj97OwszGZz0RW6/IyOXC5HIpFAOp2GRqPJcXXkcnlNY0G1CGQWq+8hSCS5JQmVpLIryT7V2/KJRssvfq8lCAAABtVJREFUmGy0QkMqPgVgGAYOhyOv5YJYOaToUGx3O5Dt9xodHcXKygonIBqNJkdMxKyrwcFBtLa2QqFQ1Pwx8tnJNC5fhGKZKB559rDo32vJXrV8GqnQkIpPEYRznoEH4sMPLIuJj1wuR29vb0UB591SZFhNxKyZEy88gnv3dHWpX6m35UPHajTp6pxSYRgGUqk0781B3qSk6FBMfCql3qn2ZqGeS/jIh1q5UPFpMsR2fBEcDgdCoVDe8PlC1y8FMlqhHtRKfIQ7rYT3WQ8RqKfQbsftojGfJkLYciH8m8fjweTkZNVezM1kfQD1tUDqdb+FUu1iw8TI7xcuXMC//uu/YmFhAefPn4dSqcSVK1fqcPbVg4rPFvBbLsQggWUxt6USV6bR+neA6g0722tkMhnRlpNEIoFwOIxbt25xWU8ABYeJaTQaPPfcc3C5XHjttdfw3e9+t74PrEpQ8SkB0nIhtkpHJpNBJpPlFB1S8ikkOvW0fCqZMikUEmFhp1BMhEIil8uh1WoRCATQ29tbMOspxGq1IplM4vz58xU91t0IFZ8SINZPLBYTnbWjUCgwMzPDFR1u974oO0Mmk0E8Hi8oJuR7MTERzv8pVUwmJye59dyl0mgBZyo+JSKVSkX7rkiAmF90SC6vhEZ0u4px59IYlEolnv1vZ7d9W8IeNr41wv+diEksFsPY2FhODxtfTIig1CsBIIQWGTYppBJZOENXrOiwUvg7wZqB668OYj0QhlyeEI0JiYmJ2HcSWxOLmZCmWP7lREz6+vpw5MiRuolLuTFBpVJZMPa4F6HiUwZmsxlqtVp0zjO/6FCj0ezo8Pi9AOmu5wvHysoKotEYJJIY5maBVDoFXE/miInQzVEqldDr9VVpOyGFhvUQH3LflTzfjVKTRcWnDEjsh99yIZzlMzMzg4MHD1Z0+/Wsbi4X4aiOQiniZDJfTIhwHDqdXSetkMvxzMef3pEeNj67ocWiHPFpBMHhQ8WnTIQtF/wXkMPhwMTEBJLJJOemlUM9xYc/RGyrjA557FKpNM/NEXbXy2Syoo+JfYGFXC6HTqfbwUebpd4tFmQvXLnHJRKJhlihQ8WnTBiG4V60QvEhRYfz8/Noa2sr+7arKT5iQ8QKiQkARKNR3Lp1K8edqdaojmLU89N8N1g+5UIyXk6nswZntbNQ8akAuVyOzc1Nzlrgm84ejwdXr16taKVKMfEh97VV8JVfOUuGiAkL1oR1JxKJBNevX8djjxWem1wr6lnnU2/xqbS/KxQKUfFpVvgtF8IlgjKZDDabDX6/Hw6Ho+BtCN0cMso0EolgbGysqJjwhUOj0cBoNOaJCWVr6u12NXt/FxWfCuC3XIgFDb1eL/r7+zlREW6/JMhkshwrJJPJQKFQwOFwiC+sa0Ca2fJp9pk+VHwqpFjLhUqlgtfrRSKR4GpMShGTYDAIhmFgMpU+CJ1SOfW2fCpxu6j4UHLS7hqNJu/vDoej7HqMvZRqrxYSiaRuAlBvy6ecgkFiHTZSiwUVn21ATGehYFTqRjSj+AD1aympl/iwLMsFnMn9l/I/YBgGR44cqchd241Q8dkGEokE7e3totZNJVWo9dxWWi8aIdXO//+Rn7f6n5L+MVK6QZIY5GfhF5D9X73wwgvbPt/dAhWfbWIwGBCPx/PEppI31W5albxT7MaAs1BA+N8L/Z+ISBAB4QuKmJg0chKhVKj4bBOxlgugcstHbGRHrWmUXqFCFBITEm8Sc2OEYsL/nYpJdaDiUwXEtlxUQjqd3vGVOfWmXMunEhcHgKh709nZyVmafFEh50WpLVR8qoCw5aJS6tXVXm/Lh2VZ7ov8vhVi7s1W8RLK7oKKT5Xgt1xUSiPEfMqNl5BFioC4oFAXp3Gh4lMlim25KJXdKD7VCL4Wi5cYjUYqJk0KFZ8qQUZMkC7xSqi1+BSKl6TTaa5DX4xC7g2Nl1C2AxWfKsIPPFdCueKz3eArEQ2VSoWTJ0/SeAllR6HiU0X41k8lcRRigQirXml9CaURoeJTZfgNg+XGS+x2OzcSg9aXUBodyRZmenPtcakSyWQSGxsbogJC4yWUJkT0RU3Fp0bwq2gplCZH9E1A3a4aQUWHQilO881voFAouwIqPhQKpS5Q8aFQKHWBig+FQqkLVHwoFEpdoOJDoVDqAhUfCoVSF6j4UCiUukDFh0Kh1AUqPhQKpS5Q8aFQKHWBig+FQqkLVHwoFEpdoOJDoVDqAhUfCoVSF6j4UCiUukDFh0Kh1AUqPhQKpS5sNUaVzgKlUCg1gVo+FAqlLlDxoVAodYGKD4VCqQtUfCgUSl2g4kOhUOoCFR8KhVIX/n9R3Q+JQKq2oQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "elev = -40\n",
    "azim = -80\n",
    "plot_figs(1, elev, azim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.21707439 -1.2173145 ]\n",
      "-0.0004483572238646694\n",
      "pca_score [0.72408528 0.27427885 0.00163586]\n",
      "pca.components_ [[-0.33847725 -0.77400604  0.53511475]\n",
      " [-0.7109608  -0.1621726  -0.68427684]\n",
      " [ 0.61641536 -0.61205775 -0.49539622]]\n",
      "x_pca_axis [-1.01543176 -2.13288239  1.84924608]\n",
      "y_pca_axis [-2.32201813 -0.48651779 -1.83617326]\n",
      "z_pca_axis [ 1.60534424 -2.05283052 -1.48618865]\n",
      "x_pca_plane [[-1.01543176 -2.13288239]\n",
      " [ 2.13288239  1.01543176]]\n",
      "y_pca_plane [[-2.32201813 -0.48651779]\n",
      " [ 0.48651779  2.32201813]]\n",
      "z_pca_plane [[ 1.60534424 -2.05283052]\n",
      " [ 2.05283052 -1.60534424]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAADmCAYAAADsvYEoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXgkd30n/OlTfah1jzQjzeickUa35vCMj0AIJC+JDdlNcHZxAoFsAgaCQwKxk7xJSMyGZ3fNy7FvHoJj8waSYQkGA4kXm3WMDwJje2Y8h+7RfXS3ulutvrv6quv9Q/xK1d3Vd/VINa7P8+gZjbrr6Or6fep7fr4anuehQoUKFbca2v0+ARUqVLw5oZKPChUq9gUq+ahQoWJfoJKPChUq9gUq+ahQoWJfoJKPChUq9gX6Aq+reXgVKlRUCo3UH1XLR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxL1DJR4UKFfsClXxUqFCxLygkqaFin8BxHDiOg1arhUazq0hA/lWh4naASj4HEAzDgKIo6HQ66HQ64e8ajQYajQZa7a7BSohJJScVSoRKPgcMqVQKsVgMGo0Ger1eIBQyX43neTAMI7mtSk4qlASVfA4IeJ4XiEer1QoEQiAmDikSySSn7e1txONxdHV1pW1HyIn8rpKTiv2CSj4HADzPI5FIIJFIQKfTlUUCUuREYkbkGOTfYiwnlZxUVBsq+ewzeJ6Hx+NBKpVCQ0ODbAtco9FAPAq7VMtJ6nWfz4fW1laVnFTIApV89hEcxyEWiyEQCIDneTQ2Nsq270zyKeb9Ur8TsCyLlZUVHDp0KK/lJEVMKjmpkIJKPvsEjuMQjUbBsiz0ej2SyWTWe3ieL3vBlko+xeyP/JvPcuI4LudxVXJSIYZKPvsAlmURjUbB8zz0ej20Wq3kgq1kMVaLfIp5vVhy2tzcRH19Perr64VtCDmRWJVKTrcvVPK5xaBpGhRFQaPRCDU8Wq0WHMfJehy5yQdARfuTIqdkMgmWZdPIN5flpJLT7QeVfG4hksmkZCq9GkRRjX3KDbFbWY7lJAbZJh6Po7a2NoucVGI6eFDJ5xaA53kkk0nE43HJVLoSLJ/9Xrz5yEn8OW/cuIE777xTcnvxT2bbyn5/vjcjVPKpMnieRzweRzKZzFnDkyvmc5ACztVAJZ9PDCkyzzwO+Zf8SO1DJadbC5V8qgie50FRFGiazls8qNFoDrzlo2RkXvdclhPP81nfA8/z2NzcRHd3t0pOMkMlnyqB4zhQFAWWZQtWLSvB7aoWDsKizefScRwHt9uN7u5uSXISb6eSU2lQyacKYFkWm5ubaGxsTOtKzwUpoiA3fU1NDcxmM0wmU5Y7Ueo+DxoO+vkBKErWJJ/lRCAmJUBt+gVU8pEdJJW+urqKc+fOFbVNpuXDsixmZmZgtVqh0WgQj8eRSCTA8zwMBgMsFgvMZrPwY7FYYDQas57gSljcBx3i/rhcKKd15erVqzhz5kzadm82clLJR0aI5TBKgZgoUqkUpqen0d7ejmPHjqXd+DzPg6ZpxONxxGIxxONxBAIBxGIxpFIpaDQawVLS6XSIx+MIh8Mwm80wGAyyfla5cNAXVTHkUwhSDwWapgs2/Yrvi9uRnFTykQGFUumFQCyfeDyOmZkZ9Pb2orm5Oet9Go0GRqMRRqMR9fX1Oc8jFoshFAqBpmmsr68jHo+DYRhotVqYTCbBWhJbT8W4h3JDCZaZHOQjtc9MkTip38n/y9Fy0ul0sp+33FDJp0IUk0ovBI1GA5qmMT09jZMnT6Kurq6sc9FoNDCZTALBBAIBjI2NCa8TgiOWk9/vF/5PFoSUS2cymar2ZD3oT2y5ygHEIFXdxaIct47jONTU1FR4ptWFSj4VIF8qnTyxirlxw+EwKIrCuXPnYLFY0vZfLqRiPlqtFlarFVarVXIbhmHSyIkIkpF4UywWw8zMTBoxmc3mrHhTsThI7R+5cCssn0qh1FifSj5lolAqvVjy2d7exubmJiwWSxrxVIpybkK9Xg+bzQabzZb1Gs/zuHjxIo4ePSoQFLGcMuNNYuvJYrFAr781t1k1rJRqkE+plk85OOgWJaCST1lgWRYURWWlYcUgcZx8N5nD4YDX68X4+Dhu3Lgh6zlWo71Cq9WioaEBDQ0NWa8TNUZCTOFwGB6PJyveJCamVCol6zkqiXz2I8Z20KCST4mgaRqxWAw8z+e9gfIVDvI8j9XVVcTjcYyPj1flKXirzW+NRiOQihSk4k2RSAQzMzPQaHbF8jNjTaS+qVhCUQr5VGOfYlTjOlQDKvmUAHEqvdCTKxf5cByHhYUF6HQ6DA8P57xJrvzwBrRaLe5699myzvWg+f5S8aZkMomenh7YbDYh3kRKCDweD2KxGJLJJHieh9FozCKmzHgTz/OKIArV8tmFSj5FgKSwE4lETjcrE1Lkw7IsZmdnUVdXh66urrzEE3AHodVqcenZawCA8/edLumcDxr5FEKheBOxOAlB+Xy+rHhTTU0NEokEPB6PQFCVxpuqZU2p5KOST0GQVPr169cxOjpa9I2YST4klX7kyBEcOXKkWqcrQGnkkw/i+qZ88aZQKIRQKCQZbxK7dGLrqZBVowacqweVfPKApJdpmkYikSjpCxWTTyKRwPT0NHp6etDS0lJw2zt+ZQJXfngDGo2mZIuHQAnkI6ekBok1mc1mnDhxIu11Em8ilpPf74fD4UAikQDHcWnxJrFLZzKZFJFql4JKPgoGmSxBnpyl1O0Ae+QTjUYxNzeHgYEByarkXLjjVyaEm/7Ss9ew+MYK+s/2lU1Gbwbk+n6KqW8ixBSPxxEKhYTC0VQqBZ1OB4qiskoIDAZDWYv8Vlg+SoBKPhKQSqWXQz7hcBhutxvDw8M5b/xCuPTsNVx7YQpUKIZogCqahJTw5JM7nlLu/vR6Perq6iQryzc2NsBxHBobG/PGmwgxiQkqV7yJZdkD22t3K6GSTwYYhgFFUVmp9GLqdsSIx+Pwer04ffp00WXuSkmRHlRU4/qRTFuu+iaO44R+OmI1ud1uxGIxwcLJjDUlEomqk48S7iOVfETIl0rX6XTCjK1CcDqdiEQi6OnpKZp4cllWxMLJtHjKzYIdNBwEy6fQPvM9cMTkIgWWZQV3Lh6Pw+fzwev1Ynt7G+vr65LxJovFgpqaGkUQSCVQyQe7NxghnszJEgTFqA3yPI+1tTVQFIX29vaSbp5CAeLbMd5Tjd6uagSHKyEBnU6H2tpa1NbWpv29paUFzc3NgkQKcenE8SZidWXGmohEitLJ6U1PPiRNm0gk8nalFyIfnuexsLAAjUaDkZEROJ3OkqRRyf6lsiCZpHPh0e8gGqDQM9p521hAckApFc7igLPBYIDBYJCMN5GHIiEmqXgTUTAQExOpeTroeFOTD0mlk4xGvhtXq9WCZVnJ16SKBzWa0kThMydY0DSN2dlZMAyT1rC5+Ooa/J4gXEseeNa96DvVDeeSG4DyCEhusqjUSsm1z/1KtYuLJ3PFm8T9dKFQCC6XC7W1tRgcHJT1nKuBNy35kFR6ockSBDqdTpJMSPHg4cOH0d7eLvxdq9XmFH6Sgpiskskkpqen0d3djSNHjqQFNCmKgrWtBom5OIKOIEyHDNje8uGl7/wEPp8Pd77rjBAzeLNBiZZPJdBqtZJqCHIPI6gW3pTkQ+pvipksQSDlduUrHsxnKeXbfzwex/T0NE6cOIGmpqY0gTAA+LUPdeDCo9/ZfRI2AlszO4iFY6D8cdQ2WtFzx1GhJyoajeLq1auC1VRMGvhW4qCk2vNBaUWGB72wVIz9vwNvMViWhcPhgM/ny6qEzYdMMqEoCrOzszmLB3MNAswFjUYDiqKwtrZWUM2w/2wfogEKng0vtBpAp9Ui4qOwObmVZm5fvHgRw8PDacFMkgYmiyozy6LkTMt+ZLvKgdpesYs3FfmQyRKlxmOAdLcrFAphYWEhb/FgqbO4GIbBysoKRkdHCxYkitPvALByfR1mbxjN7Y0A9tLwmsY9q6mxsTFrPyQNTIKZ4k5yAFniYCzLgmGYA2E1SaFalk81ZFSr2V6hBOIB3kTkQ+ImWq0WBoOhJJcI2CMTr9eL9fV1jI2NCa5QvvcXA7/fj3A4jJMnTxZdCX3+vtMCCT389s8glaShNxpw4dGn0X+2t6h95EoDA+mi+LFYDJFIBKlUClevXk0rnqvEanozu11qe8WbgHykJkvkCh7ng1arhd/vRzKZxMTERMEK1WLJh8iotrS0lF31OvbzQ/Bs7AAApv59HtFAFD2jnfDMu6ELmMrKgoljTcRq2t7exvnz5wEUtprE9SnieFM1rSal6PlU4zyViNuafHJNlig1GMzzPPx+P2KxGM6ePVuUyVwM+WxtbcHtdmNiYgLr6+uSc8KLwfv/6n48/I7/CgCw1pmwOecAAASjAWgpQ1VS8IWsJlK0GY/HEYlEsiQuLBYLKIqCy+VCXV1dyaqFUlCK5VNtqG7XPiPfZAnSKlHsfhYXF8EwDFpbW4v21QuRz+bmJgKBAMbHx4XzKzdTcenZa+gdPQbPuhdJKgnHhhcMw8FyyIhogMKlZ6/d0hogcX1KvljT5OQkWJYViCmRSADItpqKFaJXSp2Pil3cluRTaLJEseTDsizm5uZQW1uLpqYmhMPhos8hF/kQ/eZEIoHR0dG0CZSV1Gf0jHYCADzrXtTodeDbrEhEoqBCMbz8zxcBHJwiRGI1GY1GHDt2DEajUXitWKspk5hMJpNq+fwMquWzT2BZFtFoNK/AezFfDk3TmJmZQWtrKzo6OhAIBMpqlxCDWFEAMDQ0VHDeEs/zRS2o8/edxqVnr6Hx8G4V7OIbK2gYaMOz3/0xEgtbaGyrR8AdvOUWUCFIfbZirSap2WJS2jsk3lRurElp8RklKSPcVuQjTqVXksokFcZdXV04dOgQgNItk8z3cxyH+fl5mM1m9PT0ZN0gUvsv5SYSd7uTJtR/+qfngUAYja31WJvehNfhO1DkUw7yxZo2NzfBMIygvRONRuH1etNE4TID4IUmslYjI6diF7cN+YjlMCp5UpHiwf7+/rR+mnIrloHdp/XMzAyamppw7NgxyffLJXsqJqE6ZxzzrQaYbjrgXvfidz773or3Lyeq8ZQ2GAxobGwsymoixEQmspJYk5iY5CYLJbpx1YLiyUcqlV4uSPHg0NBQ1pO1VMuHFDIWKxxfacxH8hwAtK7HMNtVgzsp4IdPvoinP/8DNLc34hceuEfxVlAmCrlIxWToCDFRFCWQ06uvviroRGfGm0wmU0lkcivG5qhu1y0Az/MIh8NYW1tDb29vyRdd/OTd2dnB2tpazuLBUmuDCJlMTk6muW+5UE7VdS6QCudjw4cReXkZjREesw0cIpPrMFlqoNdrce2FKQD7H4SW26Up16rI1UEeiURw9913p4nQF2s1Ed0dMVTx+D0olnzEqfRQKFTyBRfr54jrbXIV+pXqdhHzfmxsDE1NTUWdj1QXfLk3Emm96Bo5BttWAK8xKaCGRZPTj8B2CF6nH6d/aaysfcuFaoiJVWvh5ROhF88VE1tNJEMnnuaq0+mE95ZqNd1uUCT5iFPp5bRKALuWDMMwsNvtCIfDQr1NLpTiFpG4kclkKop48u2/nAVFrBnXlgtMlEdtgxXHVz2YHLWghmJgiTDQG/T7bvUA8j6lq1HnUwwKzRUTW02BQAA0TWNhYQHxeDxLrVD87+0uMq848hFPlqjEfNVqtVhZWRGUBws9gYp1u8LhMG7evInh4WHMzc0VfT5yz9k6f99pGI4Ak99bQCKaxM5WAB2LUSwMWTH4Rgg7Dh8+cPwP8Ivvfwve/1e/IdtxS4GSLJ9KILaaDAYDaJrG8PAwgOxprBRFYWdnR9CaEmtEZ2boct2zB/EaSEFR5EO+JAAVEQ/HcYhEImhubsbAwEBRX1Yx5BAIBLC8vIzR0dGcguK5IGX5eL1eOJ3OrJqVYp+KGo0Gfae78dKFn8Jg0KHFwyBwiMbyiA39VwPg/DyuPj912+hDKyEtnhlwLsVqisfj2NnZEVx68Qx7cl9YrVbFCMkphnx4nsfb3vY2PPfccxX5yQzDYHp6GiaTqWSR93zwer3Y2NjA2NhYWV9+Jrm53W5sbW1hbGwsLZ7g8/mEuhVxBkb8QzrLNRoNzrxzDK9+7w3Ut9jAMiy6FsKYPd+EjcE69MxF4Fzawtf+/FsAbn3w+aB3tR+EiuliYk3iDF0ikSja1d9vKIZ8AAgFhFIo5kYhxYOdnZ0IhUJlxYqk4HK54HK5MD4+nmWRFHsDiy0fp9MJr9eL8fFxYSCdlGBZZme52+0WGmmB3c+7srKCBz77q1i8uIr/7+FvQU/z6JwPY3WiAeYog8PrMcSjCSy+saJ466camtDVEBKTK9sltprI/aEUCVVAQeRTSNy9UAwoFothZmYGJ06cQGNjIyKRiCxflN1uh9/vlwxYlzLllLzXbrcjEAhgdHS04E1aqG7l+vXrQrFd55l2aPQa6PRaNHhTaNqKY+u4FaYoAzj9+NGFn4DjOJw8d+KWN6HKBSWQj5pq34NiyAfYq4XJNdAv15dKgsDi4sFSOtulQGZ0xWKxtAZRMUqZckrGK9M0nRYAL/dG0mg0MBgMqK+vFyRZx98yjMvPXQPA4uhCFJEmI9ZH6zBwOQC/O4jXnnsDsLHgm1OScaaD3hogdx/WQRaPvx2gqKtgtVpBUVTW3/MRic/nw8LCAkZHR9MshHLIh1gnPM9jaWkJqVQKw8PDOW+mYtPzPM/D4XCApmkMDQ3JdnNmxpH+9BsPwdZkAwDoGR5dsxFwei1WJhpAa3nEA0mcHDyJhR9uwr8YQV1dHVKpFFwuF6ampkBRFF599VVcv34dCwsL2NzcFDIzB8Hcr0bM5yC7XUqHoiyf2tpaRCKRLHH1XETicrmwtbWF8fHxNNkGsk05aoYMw2BxcREmkwknTpwoyh3MB0JkLMvCZrPJerNLZeg+9j8/iM++90vgGA51/hRa7DHsHLNgdaweumt+/Nf7v4AaqxH/8aFfgdlsTnPBXn31Vdx5551IJBJCnEmcfQGkFQwtFsstWXBy1/koTR+IfNeq21UF2Gw2RKPRrL/r9fo08uF5HpubmwgGg5iYmJC88bVaLVKpVEnH12g0mJ2dRWNjI7q6ugq+v9gppzqdDn19fVheXpZ8TyWuVyb5nL/vNH7jj9+N73/pOaQSNDqWKISbaxBtMsI+UIvOm1EkqRS+98Vn8et/dF9W/CfXrChyrqTBNxaLIRwOC0FwMude3ILAsixSqZRso3+VEPNRxeP3oCjyqa2tlSQfUq0M7N6Ay8vLoGk6ZyyGbFOK28UwDCiKQmdnZ1HEA+Qnn0yJjVQqJX9jaY7apP6zfegc7AAVisGzsYOu2TCWzjZg55gF5iiLQ4446CSD5554Ec4lN/70Gw8Vfbx8WjwMwwjERCbFTk5OgqbptBYEscVUirSqEshH7Wrfg6LIx2azIRKJZP2dEAlZ0CaTCYODg3lvxFLIJ5VKYWpqCmazuWCDqBi5mkU5jksbrwxUqas9T2HkwLnjuPbCFPR6HWxBGq2bcWx3WWAfqIWJYmAL0Ij4o1i+viabCJler0ddXZ3gNrvdbtxxxx0A0ovpiEhY5hgfKXdOvJCVQD5qzGcPtw35EIJobm7OqZmTuU0xi51MJe3r64Pb7S6JIHQ6Xdbiz6XtI3d7Rb59EiIJ7UTA0Cy8dh/al6MINRuRrNVjdaweJy/7URPn4F7bxtOf/8Huho3Va2EoVExHRh/FYjEEg0FsbW0hHo+D4zgYDAZYLBZEIhH4fD6wLCtLb5QSU+2q21Ul5Ir5kPqYvr4+tLW1FbWvYrrUSW0QmUrq9XpLctU0Gk3a+0l1dVtbW9pcd3I+t9LyAYBfeOAedJw4LMR/umfDWLijEaxxNwM2cDkAHcvj5qUlXHthCsFoADp/Dc6/64ys51kI4jE+UtW7pAI8EokgHo8jHA4jFouBZVloNJosqYti54upqfbqQvHkE4vF4HA40NjYWDTxAIXdrkgkgvn5+bTaoHKkVMniZxgGU1NTaG9vx+HDh7PeeystHyC9leLX/vBefPfzP4A1vFvx7O61IlGrx/pYPXonwwDH4eVvXkTdESsm7jx4WtCknqmmpgbd3d1pekyFqsBNJlOWK2c2m4XvrhrZLtXy2YWiyKeurg5ra2vC/0nxYGdnZ8mZq3xuVyAQwNLSUlaDaLk6zsQlzCcqVo2bphhCIyRyuLcNW8tuHF6lEDpkRNxmQKjFiK0+CzpWYohF4kglU/g///AyfvPPf132c5UDUmRRqApcXDbg9/vhcDiQSCTAcZxQ56PX62URoweqa/kc1K7+XFAU+YgtH7/fL8w2TyaT2N7eLmlfudyunZ0drK+vC31VxWyT7xjJZBKTk5Po7e1Fc3NzSedYKUqxprpHjmHH6QMXTaJrJoKF843gtRp4ui0wRxg0uRNgUixWbqwf2D6wUutyxBm2zO+GuPLRaBQ1NTWIRqPCpAyGYaDT6SQD4EajMe85qAHnPSiOfCiKgtvthtPpFIoHGYYpuVpZyu0S71cqWFlqYSLLsnA6nRgcHJRMPReDSlyxYsmHEIl71YPl6+uwRBkcWaGwdWLXWtgYsqEmxsAaZgAeeO6JFw+kDIecT34yiMBqtWbF54BdN1rszgWDQUGDB4DkHHuTyaSm2kVQHPk4HA7Mz8/jLW95i2D+ltMqIY7HAIDD4cDOzg7Gx8dzmtWluF2xWAwulwttbW1lE0+lKMXyOX/f6b2sFoC2jRiCrTWI1RvA6zRYnajHwKUAjEkOEX8Uf/eJrwvbHRTcylS7Xq+HzWaDzWaT3C6zCpzoPUejUVy7dk0gJqvVKhDUm80iUgz5cByHxx9/HMvLy7j77rvTCKKSJlGe57G+vo5oNIqxsbG8T6Vi3S4io9re3l5RfKBSlBrE/txLn8aFR7+Dp/7Hv0LDA90zYczf2QRepwFdo8PqeD363whAywE+px9Pf/4HB4585G4sLSddn68K/OLFixgdHRWIKRQKweVyCWUDxJ3LtJqKrQJXYz5VgMPhQH19Pdrb27NiMeWSD6mGZhgGIyMjBb84Iv6dDyRLNjw8jGg0KswfLwfRaBRzc3PQarUwmUxZomFmsznvOZeTQXv/X/0Gpn48D/eqBwFPCB3LUTgGdp/usXoDNobq0D0ThgbA/OuLePjtn8HnXvp02Z9RTihBTKzUKnC/3y8pqSq2mMRV4Cr5VAGdnZ34zGc+g3e+851Zr5VDPsQ01mg0OHnyZMmCX1Ig2beRkRFYLJaKur0JiY2NjcFms6WZ8T6fD3a7XZDSzEdM5cSM7v/Uu/DyP1/Exe9fxqHNOIKHdnu/ACBwxCSIkIEHFi4v4SMTj+DxG4+V9TnlxO0w4C+zCjzzfPJVgdtsNpw6deqWnm8lUAz5AIDFYhG6p8Uole1ZlhUsir6+vqK3z+d2BYPBrPR8uYWDhHhGRkZgtVoLZmXE1b8kXRyPx4XZ5ZFIJIuc8i2q8/edxvn7TuMDx/8APqcfXXO77hen392GiJA17KTAsTxcKx48/PbP4P5PvWvf3bCD3l5RCYqpAlcSFEU+ctwIpMq4tbUVyWSyJNM6V7aLpP0z9ZvLIR8yNbVYEfp81b92ux2pVApNTU1CfEHcliAedCclf/GL738Lvv3YM6iJc+hYjMI+VEcOKoiQmSkWLMNi6drqgU3Blwu5yaeaYmzEnVPdriqjXF+cpmlMTU3h6NGjaGtrw/b2dklFX7kmTGxsbEhqBpVKPhzHlUQ8xZyvTqeTnF2eOeguEonA4/EIbQkGgwE1h/Ro7W6GZ3UHLc4Egq01iLTskisRITt52Q89zYNJMvjB4y+g/2wfgIOVBSsXcpNPtS2pg640mQlFkU+hGEY+Ukomk5iamkor9iOxomIzGplul8fjgcPhyFkXVAr5hEIhJBIJnDlzJo14qlXnU2hkC03TGBgYgG8tiIvffQNBVwhdcxHM32UAa9hdQCmLDqtj9ThxLQgND1DBGJ7+/A8w9vODtwX5yJ09UwsM06Eo8gH2tHsyF3s+EXnSINrf35+20EqtWBa7XS6XC263W5a6IBIvstlsWdZTJajEBCf9Uh/93O8g5Iri6vOTQDSBowtRbIzsBUOjTUbYB+vRORcCsJsB87p2YOuw4Ny9p9Iqf5UGpVk+gJrtqiqsViui0WiWG5FLRJ4EbwcHB7MKwkqtWCZk4nA44PP5MDY2lvdJVsqgwbGxMSwsLMhqOufSEyoVREzs3db3ocmVQPCQEaG2vebNnY4amKlaHHKlgFQKvs0Anv2fL2L8HUMIBAJCqljcYW6xWMAwDJLJZMGWhP2C3DKqquWTDsWRD9H0yUU+YgSDQSwuLgqp70yUmqLXarWgKAparTavSqL4/fkWf2agWi6yqBaODh+BfcaF7rkIblr1SNbu3T7242aYYhxs3hR4jodreRuPnPss/vjrH8P5+3YFwzI7zEkMLpVKZRFTKdIX1UI1LB+VfPagSPIpRsfZ5/NhdXVVskGUoBTyIY2GDMPknVghRj7yIcQjDlTLrekj9/7+y5f/E1744qt49ftX0DcZwsK5RiH+A60Gq0MWnLycQE1895jxaAJf+L3H8cmvfgTn7zud1WG+tbWVU8mQBL9J+lhqMmspEqvlQG7yqXZHu9KgOPIpRseZBIInJibyBpNLGW2zsrICmqZhMpkqyo4Bu8S4traWlSGrhqaP3Pjjr38Uf+4OYf61RfRMhbB8ugH4GQFkipABABWKFZWCz1fDktkr5fV6hV4pAEKRZSqVws7OTlHV38WgGuRTbS2fg+i+5oLiyKeQjrPT6cT29nbeQHDmNvnA8zwWFxeFSug33nij6HOVIh8i2TE2NlZUar7SbFc13Lixnx/E6uQG6vxJHFuIwn5yL5aWqNVjfaQOvVO7LRjgeTiX3BUdr9DEDEJMbre76OrvYkhFiQFnJUGR5JPL8iEay4UCweJt8pEPz/O4efMmjMvzPnIAACAASURBVEYjent7S36qZFoyYq0gKYtMbsunGn1JwG7/l3PJjZ9+9xJa7HHErTrsHNsjhlBrDbb6behYpgCOw2vPFE/Y5ZwTqf42Go0YGBgQXsus/g4EAnA6nUKRJRGlz1X9rabaq4vbgnx4nkcgEADHcTh16lRJbhFx1TLBcRzm5uZgtVrR3d1d1kIWb+P1erG5uZmTeMj5yGmpVMPyIeT4p994CA+8PIuIP4pjC1EkrHqh/wsAPJ0mmGMcmpwxsAyHBzo+gns//A68/69+Q9bzyYd81d+ZM8ZIdznpxzMajaAoChsbG7BarbIMP6y25aMklwtQKPmI3S5inQBAa2trSV+uTqeT7IfhOA4zMzNoaGhAZ2dnxee8vb0Nu91e0BXM1BiqFNW+Ge/98DvwnceeAcfx6J0K4ea5JqQse4tz44QZxjiL2p04IsEYXvpfPwWAqhBQqdetUHd5KpXC66+/DpPJJFn9LWUxFXLzVcsnHYojn7q6OrhcLgB7869qa2vR1NQk2XSaD1JuF8uymJ6exqFDh9DR0SG5XSntHTRNC1XQhW5OuS2VartxhES++/kfADSLvhtBLJxrFBpQeZ0Ga0MWDFyhYUzx8Np9+NGFn1RFBVFu+Quj0Qi9Xi+pYihuS6EoSgiAE3KRIiaDwSBMba0WVMunyiBuFyGJlpYWHD16VJjXVAoyyYdMmDhy5AiOHDkiuQ1Z0MV80W63GzRN49y5c0XddNWwfOTOnmXur/9sH+781bP46fcuwUyx6JkOY2WiXsiA0TU6rI7a0H8tBK3BAL8rgK/84derQj63KphLqr/r6+uzXsvU4/H5fEJNE03TMJvNSCQSWbrPleKgZ0mloEjyicVimJycTBtDU46mj3gbUvB27NgxtLa25tyGxGUK3eik/cJsNhf9tJOyfMLhMJLJZFk36q1I3RP5jQc6duM/9TspdCxRcPbvTYuI1RuwOWxD1ywFcDx2nH789/f9Ld76Mfnmfx2UyQ359HgWFhaEGBQJgOeq/i5WkF6Mg/D5S4HiyCeRSOD555/HBz7wgbT5V+XqOItH23R3d6OlpSXvNuQ4+Qhla2sLHo8HY2NjuHbtWsnnQ0B6vtrb2xEIBEBRFGiazjLtSUA0M5B9K+uG/vCJD+OfP/s9LF9fR+tGDPFaHfztew2y/tYamOM02hZZgOfx0+9dAm9icPfdd8ty/INCPoVgtVol7zGO49IsJhL8JtXfUkWWSpPQyISiyGdrawuf+tSncPToUQwPD6e9Vq7lk0qlMDk5ib6+PslpmJkolJHa2trC9va2kO4n1kwxLoGYLEhryPj4eFZPGsuyQryBCIhRFJU20sVqtUKr1SKVSoFhGFliDYXIbODccbjWtkEFY+icjyBp0YNq2CNEZ6cVmnobWq96AY7DjefmcOHR78gSgFYC+eQLOGu12pzzxYqp/jabzRgcHKzq+csNRZHPoUOH8Pjjj+Ov//qvs14rh3xomkY4HMb4+LikrIQU8pGP0+nEzs4ORkdHhZusFPIh+w6FQlhcXMwSJyPQ6XQ5JyeIYw6BQADRaBRXr14VCEhsKZF/5cjAkBjO6tQmlq6ugkkxP8uANYI2/Wz/Gg22bByMI81omPKCCsRx8V/ekIV85G4CrQbKTbUXU/1Neg6VBEWRj8FgwJkzZ/JWOBeLWCyGubk5mEymookHyC3DQTrdR0ZG0hZzKUFkjUaDeDyOra0tjI2NpY39LRbimIPNZkMymcTExASAXWIi1hIZgpeZPhbXtJRKTOfvO43FN1YQ9AThWt2GgQb6JsNYPNsATrdLDJxOA3sTj1qjFvoUB+fiFi49e03YvlxUQzxeblQj1U6qv6vd51YNKIp8AKCmpkZyNHIp5EOmQgwNDQk1QsVCSobDbrcjEAhIdrqXUjiYTCbh8Xhw5syZNOIpd2Flukl6vT5nloamaYGYwuEw3G53WsEdEcTf2dlBY2NjTh3o/rN9iAYo/Ogb/44klYIlzKFrPoq1kT0rjTZosHa6GccveQGOx3/7zf8X7/nkfQDKJyAlTK5Q2yvSoTjyyfXlFXujiEfbSJmxxRxfTCabm5sIhUIYGRmRPLdiyScSiWBrawttbW2ySKgCpQWcDQYDGhoasqxAIrdKURRCoRAikQh8Pl+aDrTYjRt52wB4nsfkK3Nwr22DSTFo9CQRt2jh7t273hGbBo6JQzg26QeTYvCjCz8RXiuHgG7lwMByoRYZpkNx5ENQzs1Wqji7FMRu18bGBiKRSF6JjWLIhxBiV1eXrBMI5Mh2ieVWzWYzuru7hQZPcYsCRVEIBAJwOBzgmpJ4/5f/A5757I8w+8oSwHE4skIhUV+DYPPeLedt0cDUUYND9hh8Tj/+7WuvwLnkLpt8DnoT6K3oalcSFEc+5V5gsWJgObEUAuJ2kSmnQ0NDBaec5iMf8ZicRCJR0ZDBTFQ71Z6vRYHnefAf0qG2zobLz00BPNB9w4+FOxoRr9vLgNn7rTBFadgCNILeENyr27j07LWSCagabpfc5KO6XelQ5JUwGo05LQSpxebz+QThrkqIB9glE4/HA4qiChIPeX8u8olGowLxWCyWnO0V5S6qarRXlBI8/7lfuxPv/OA7YG2shVarhZYD+m6EoE+KYnNaDVbH6pE0acFzwOr0Bn78vy/iu3//DOx2O3w+HxKJRMHjKsHtupVV2EqAIq9EbW2tZMZLaqFvb28LMhZS1cGl9FPxPA+/349UKlUU8eQ6J2Av6D08PCy4MbkyY+USyEEQJzv3y+P45N9/CE2Hd4PcxiSHvhshaNi98yIiZKxOA47h4Lm5g7q6OvA8j+3tbczOzuK1117Dq6++iuvXr2NhYQEOhwN+v18gpmqQj9LcGKWdr+LcLmCvv+vQoUNpfydSqsSvdrvd2NraKmrCRCEi4Xkeq6urYFkWhw4dKmnKaSb5UBQlEI846H3QG0vLBc+yeO8j78Z3v/AsXKseWMMMOucjaVMwErafiZBNhrByfQP/9uS/45ceeFvafsQ1LUQ8jBTbcRwnCL+JSwbKFadXXaTqQ5HkQyZYZEKcbieVxuPj43mDfMW0SxAZVYZh0NHRUVL3fCb5UBSF2dlZyWxbNfR8qt1YWgqsDXtB/mZXAgmrDp6evWsQaq3BVp8VHSsU5l9bxH9/398KkzOA/IqGOzs72NraQkNDg9CeQFEUUqkUtFotzGZzVnFlvj45lXyqD0WSTz41Q4ZhYLfb4ff70yqNc6EYNcPl5WVwHLc7RK/E7nkxocRiMczOzmJoaEgyzb+fMZpi91cOzt93GpeevYYH/u9fx9f+4ik4bjoBAO3LFBK1eoQO7VVxe3qtMFMMmtxJ/PR7l3Dh0cNFVUCTrJxUU7C4b4qiKASDQaFPjhBaJjEpiXwOgnVbDhRLPrmqnJ1OJxiGKWq0DZB/cCDP81haWgIA9Pf3Q6PRlGydkPeTwYVDQ0OS/Tvi98qFg+J2AXu1Oy//80UkYwn4tgLgGA7d02EsnGtEQjSGZ2OoDiYqAEuEwXNPvFgU+eSL+eTrmyJ9coSYSJ9cMpkEy7JgGCat6ttqtZbVJ1ft70Fp4vGAQsmnrq5OUko1HA5Dp9NhfHy86KdWrsGBJH6g1Wpx/Phx4YstZ9BgIpGA3W7H4OBgTuIh7828SWmaBs/zQpMqkH6j5bvhDuLN+KffeAgXHv0Ofvr9K9je9AJxGn2TYdy8owGsUSRCNlaHk6/5EfFHhfYLIHcBYrmZpFx9cjs7O/B6vejo6BCIaWdnB7FYLK1PLrMlJRcxKcmSulVQJPlkjs8hrhHP82hvby9ZSjXT8uF5HgsLC9DpdGnEA5Q+YplhGDidTsnu9ExkBpxpmhb+DiAn6YnJSPzZqyEgL8cT/P1/9RvonejCT597HT/9X2+gJsagdyqMpdP1gHb3nJMWPZz9tei8GcVn//OXcPoXR3Hnu8/krAGqRrYrnzZPvj45vV6f5cbp9Xq1ujkDiiSfuro67OzsAEgnira2tpLdlkzyIZrQBoMBfX19WTd0Ka5RPB6H0+lEa2trQeLJ3DchHqPRKBBK5rmIiYBsJz43hmHAcRxYls0yy/fbKjr7y+Oo6dBiZymI+dcWYQvS6LwZwebQ3kLfOWZB/U4K9TspXH9xBkFvCGf+r/FbRj75HmKl9slFo1HEYjFcuXJF0mKq1Cra7++zHCiSfGw2G9bX18FxHObn52E2m9HT0wOXy1WRmiHP85ifn0dNTU3OUTnFul3xeBzT09Po6OgomqyIZSFFPLneL/U7sEs8k5OT6O3tzVm4KOXGSe2rmrj/U+/C//PBv0M8mkCLM4F4rR7ezr1s1sZQHYZe8wE0i5UbGwh4wpIxoINUZCjVJ0dRFBYXFzE8PCwQUygUwtbWVs4+uVJmjCkRiiUfkrK22Wzo7u4GsPs0isViJe2LWBuEeEwmE3p6evIGLwsRXCKRwPT0NAYGBgTNoGLPhey7EPHkA8dxmJqaQkdHB9ra2rJeJxZTMW5ctQKZYrJo7miCY3EL4IGji7tjeCLNu2lwpkaLjUEbeqfCAMcjFopJul5yFwVWa0476ZOTakeR6pMjww9ramqyAt/iqayq5XOLYDAYcPHiRbznPe8RiAcoX82QpmnMzc3BYrGgp6cn7/sLuV2JRAJTU1MYGBhAfX09/H5/0ZYPIcGpqSlBPEpqmF2hfZDpG1KTFwDkvWFzERPP84LEBrnGxQa+84FoAP3b115BwBOChgd6pkJYON+E5M/G8ITaTPC3p9C8lUAilsQXfu9xvOsjv5RmAR0ky0cKhZpKC/XJkeGHJCNHiAnYfVCNjo7KIkR/K6E48olGo/j0pz+NmpoajI+Pp71Wro6zy+VCW1tbGpHlQj63K5lMYmpqCv39/UIsoNgYEcuy4DgO58+fFxY6RVFpM6N4nhe0dcTERISkeJ7H3Nwc6urqcOzYsZKuA4EUmZDMX0NDA6xWa97+M0JIpbhxZALq9RenQQVj0DM8+q7vjuFhDbsEYB+oRa0/hZoEByoYyxrBXI3GUjkDxCzLlk1mhYYfxuNxScXLgw7FkQ9FUfjt3/5tPPvss1mvlUo+HMfB4XDAaDQWRTxA7oyPmHjEvn4x5EOIx2AwCDe8VJYl0zQnc8lJfxNRJGxqakIwGBRE5StdlCS+NjAwIOmGSQW+xZAipkw36U+/8RAefvtncPPSEniOhynGomcqjOXTu2N4OL0WGyN1OPFGEBoAr/7rFXxk4hE8fuMx4RwOcm8XcbvkBiEmJcaFFEc+bW1tuP/++/HUU09lvVYK+ZCppLW1tRV/cYR4jh8/niXGVYh8pIgnF/KZ5ktLS0gkEjh8+DDi8Xja9APxtItSR/86nU6EQiGMjY3lXIz5rJxcxBQMBoWKdI1Gg8vPXccv/OY92Jx3gArFAB6oC9A4ukjBMbBbGxVtNGK7y4K2jRg4hkM0SAk1QEdGWw60no8qJJYNxZEPIF1kCBRPPoR4GhsbUV9fj62trbLPhYzdOX78uOTY3XzkQ4in0hqQ9fV1JJNJjIyMSBKEuIo3FoulTdgkbpyYnMiTdHt7Gy6XC6dOnarIZcj8PRKJYG1tDadOnQKwF+viOQ4nzx2HY8kN35YfTJLBoc3dMTy+jt2+sK3jVth8KViiDFKxFBbfWEH/2b4DH/NR57RnQ5HkYzQawTBM1t+LIR+WZTEzM4Pm5mYcPXpUmH5aDsRjd6SIB8hNPmLiqWSsjcPhQDAYzGuZ5Jt2Qdy4zNQvwzCgaRpHjhyBy+USiKncLnGCeDyOmZkZjI+PC3EKjUaD8/edxuXnruOXf+/teOVbr6K23gL7whZS8RSOzUeQtOgQbTSC12qwPlKHk5f9iEXi+NGFf0f/mV5MvjiHM+8cK/u8MqE0y0cln1uMzKddIfIhxNPS0iLMYS+1XQLY/aKJxVNo3pcU+chFPC6XS+jcL3ehkNSv2F2MRCKYmZnB6OioMIrH6/VifX09bYhdZjYuc2hhJshU2KGhIcnO9HP3nhJ+D7iDuPSDa7j2oyloOaBncjcDljLrkLDpsdVnxdElCj5nAE888g3UtVlRU1ODtv+8W1pQTuBbjGql2lXsQZHkU0zsIRNktntra2taCrqcDBlJh/f09BQcNJhJPhzHCTdiJcTj9XrhcDhw6tQpWW9qYpmMjY0JnffNzc1p7yFD7EixHJmmSp7uUiN4AAjXTKoqWIxz954SrKBf/r2347P/6Usw0Dz6bgSxcEcjOL0W21271c+2AI3tjR1Ay8O96E07x0yUUlgpt+ogcXFV7EGR5ENQrJ/PsiympqZw+PBhHDlyJO21UsmHpmnE43EMDg5mLUopiMmH4zihKbGQlZAPfr9fiJnIMYmUgLiRuSQ/CPINsSOWkrhMgKIoxONxGI1G+P3+tNnzxcybamirR9ATgjnKonsmjNXx3QzYxnAdBl/3Q8fwiIcSmPnJTVz54Q2cu/dUzsB3sYWVpCVFLqhuVzYUSz5msxnxeLzg+BuGYTA1NYX29va02e4EpfRq0TSNyclJmEymnDGeTJDUvFzEQ6aZnjp1qqL9ZIJhGNy4cQMnTpwoaJnkg1Qz5tLSEnieR1dXl2SZAABhFjmxmCZ+cRgGgwFXfngDvRNduP7CFHgOaPCm0L5MYetELVJmHewDteiejSC8E4WhZifneZVSWMkwjKCQQNO0LI27ald7NhRLPkTHWYp8iEVUiHiA4m8gEq/o6uoSesiKWfyEfEi3cyWEQSZdiIO1coC0Y3R2dhZlzZUCUok7OjoqlApkEjeRRyWBb1ImQNM0tC1aDP1iH268OAP+Z+TQth5DvFaPwBET/O1m1HtTaNxOIuAJ4f989SUA6fGjQhATEynU7O7uThs2INW4W4obp6bas6FY8smlZiju1crX31QKCIl1dnbi0KFD8Hg8JbVMMAyDra0t2Gy2spsFKYoSYjFyDRUEdomaBOFzEXS58Hq9cLvdOHUq2w0SI588KsuyiNsZtHa1ILQdRjySgAZA11wYSYsOsXoDNgdtqA3SMKQ4LF1bhelbe6RRCgkBwMrKCqxWa9a1yOXGEfH6TGSSEqlwrsYkVKnzUwIUTT5SaoZ6vR7JZBLz8/M4duyYpKxmKSDEc+zYMUGwvlhXjRDP6OhoViobAEwmU5aEp1RFMmlUHRkZKWvKai4QORKLxYLOzk7Z9gvsuoerq6s4ffp0RU98nU6Hn3/P3TCbzbj2o2lc+eF1BL1hMAkGfTdCuHm+EbRJh43hOvRdDyLip3D5h9ewOruOu+8/A9/ODs6/6wwsFgtqamryLlKPx4NIJCLMts+HYgorCTExDINQKAS9Xi8UVR4ERYH9hmLJJ9f4HI1Gg5mZGfT09GRNtygVhHiOHj2aRmLFBKkJ8eh0OtTV1UmOISauBkVRQmBWXJFM6mo2NjZw8uTJojSBSsHa2ho4jkNfX5+s+43H45ibm8PExIRscalz957C4hsraD7WAB48IjsUEE+hdzKMxbMNCLcYsXPUjEPOBFiaQypKY+nVdazW2NF7vlO4trnKBBKJBNbW1nDmzJmKCUDKjevt7YXJZCpbUeB2JCXFko+U20XTNEKhELq7u0smnkxzmBBPR0dHlvVUyPIRE49er5d0scgiMJvNWXEWUpEcDoexsrICi8WC5eVlsCyLmpqaLGup0BNdCg6HA5FIJG9xYjkgGbPh4WFZ3UMAeOsHz8HUasDqq3bYb25hc94Ja4RB11wE66N1cPTXwhagYaIY+LYCCO9EMPbzQwitxHDu3l0JDqkygUgkgmg0CqvVisXFxbLUBHJhfX0dFotFcP1zkUk5/XFKx21DPqTor66urmQLIXN8DqkJam9vl4wX5SOfYoinmPMxm824efMmTp48KZAfaSwlC2dnZwebm5tIJpOST/RcYufb29vweDyYmJiQPZ1MCi+lpEcrgdfrhdfrxa9/+N3QfkSLx377ywj7ogh5w2jyJJGwUnD3WrE+bMPAlQA0PMCyHILedC2lzDIBjuNw/fp1jI+PC2N3MtUEOI7LIv1iygR8Ph/8fr/QRpIP5fTHkWP85V/+Jf7xH/+x4DEOGhRLPmIdZ0I8PT09CAaDZWn6kC+VLKAjR47kDMDmEhQjxKPVassmHvE5ZLp74sbSzOLGzCe63+8XxM4NBoOwcDiOg9vtxpkzZ2TNvvA8j9nZWRw+fLhidzcTkUgEKysrOHPmjHBNH/mn38djv/1l3Ly0hASVxJEVCgmrDsE2E9w9FhxZ3W0+dSxu4Xtf3FVAkAo+LywsoLm5GS0tLQBKVxMA9soExOTEMAwWFxdx+vRpWSVSxb9zHIff//3fxwc+8AFFWkKKJZ+6ujo4nc60/qqmpiZEIpGyNH1YlhUsnsOHD+fN/Ei1ZBCRLa1WC4PBUJEKIREDyyyILPQZchX+EU1hv98Pu90Om82Ga9d2u8FNJlNJw/RyYWlpCRaLBUePHi1523xIJpOYnZ3F6OhoVvzobe+9G29779145Vuv4uL3r6BrZjcD5uqxom4nBWuYgVarAxWM4dqPpgGkE5DD4QDLsujq6sp7DvnUBKTKBCiKQigUgtlsFq6L+PrKRfpf/vKX0d3djfvvv1+W/d1qKJZ8amtrkUgkMDk5mdZRXoma4cLCAlpbWwsueqmWiUpdLQBCcLK+vr5sMTApGAwGGI1GeDwenD17ViAoEvQm1lJafc3PyCxz4Uh9ts3NTaRSKQwPD8t2zsCeBdjf359Fqpefu46AOwgAaD9+GGd+aRRXfzSN3hshLJxr2m0+fd0PlmFga8l2wwOBALa2tioOMEuVCczNzaGtrQ3t7e1pagJk9A6RUBHHlUotwbhy5Qr+5V/+Ba+88ooirR5AweSTTCbx1FNP4YEHHkh7GpWrZri4uIj29vac0qOZ7yfHkMvVIlMziIa0nCBuqdRseBL0zgRpk8gcDUPUFAkZ0TQNn8+H06dPy64kODs7i/b29oL9cwSHuw7Bvb6N3skgls42wtlfi65FCqHtEMKH6rD4xirO3XsK8XgcN2/elL0vDtjVP+I4DkePHoVGoylZTQBAll5zZplAIBDAH/zBH+Dpp59WpIIhgSLJx+l04m/+5m9w8uTJrPqUUkXkOY6D3+/Pq3mcCWIpiYmnElcLAJaXl6HRaGRPe4vbJkoJAueaWSXWE97Z2YHL5UJtbS2uXLmSFfQWz6wqFSsrKzCZTIL6QCZI8ykAvPO//AK+8ZnvwrPhhbnWDITi6JyLYGOkDvXeFDQrHiSpFELbYVx49DsY+JUuDA4OplUwy4FQKASn01mUNSWlJgDsXV9ijXq9XmGCKs/z+NznPge/34+3vvWtCAaDSCaTiiUgRZJPbW0tHn30UTz99NNZr5WjZmixWIru1QLSY0RyEM/a2hqSySSGh4dlF8SanJxEV1eXbG0TRLaT4zj4fD6cP39eWMSZM9Ezg95SY2GkPq/L5UI0Gs3S6M6EOH7zvk+/BwDw70+/jngkjmbX7hiejWEbrK/5EfSGkIinEA6H0dLbiIZzDbl2WxZSqRTm5+cxNjZWkTUl1mvO/M4YhsGZM2cwOzuLzs5OPPnkk/iLv/gL2QtEbxUUST719fW488478fWvfz3rtVLVDJuamgQiKRakb0yj0VRMPHa7HeFwWOh9kgukbaK1tbXi9pJMEDduZGQkzXrINxM9cyyM0+kU3AxxUynLstja2sLZs2dLvh7v+/R70H+2F//4l0/BsehCx1IUiVodNgdtOD5PIRGNg6Ut8K2FKrsAGSDXuq+vT7JFRC5MTk7ipZdewiuvvCK71bYfUCT5ALl7u4qtPp6dnUVjYyOOHj0Ku91ekrVkMBiws7MjdNWLx9yUovTncrng9XoxMTEha8cziR9ZrVZZA9fAbhB4cnISJ06cKKmeKp+bQUoEgsEgnE4nLBYLLl++nHMeer5rRayhxTdWMfOTeWgm1zF/qg7BDjMat+IIuqNYvbGOb3zmu4K1VClWVlZQX18ve4mBGMFgEA899BCeeuqp24J4AIWTD0VRWX8vRD6EeMQZpVLUDBmGgcViwd133502FpcUpCWTybT2iFyLZnt7Gw6HQ5Y6kEysrq6C53n09vbKul/yhO/o6JDVjSMqiKQXjMSZMrWB3G63MN1Tqi/OaDTiyg9vIOAO4tDRJoy8ZRAWmxmGuXXc6K6B1ZeCLknDueRGLBKX5fy3t7eL7gcrFxzH4eMf/zgeeeQRDAwMVO04txqKJR+9Xi9JGHq9Pif5cBwnzLUS+8kkgFwIRNuFTBPV6XSSvjlpj8hcNGTypFarRSQSwcDAgOxSC3a7HdFoVPa2CdKEWltbW3RgvliQ2qaenp60AHehoLc4KEtE9FcWNpGK0jAajVidscO9sQ1dKIHeZRpr/RacnImCTtEYectgxedNURRWV1dl6QfLh69+9atobm7GAw88ULVj7AcUSz4Exeo4k7nutbW1WUVlxbhqxc5PJ/uTSrHyPI/t7W0sLy+jvb0dPp9PqJHR6XRplhIZh1uKVeTxeLC9vS172wQAbGxsgGVZ2Z+8hNSampqKViDIF5Q9d+4cLv7rZSQSCdQdsuEHf+uEwWxAU5wF5U3C02lGQxK4eXkRP/2XS7jnP5wr61oRPfChoSFZRd0ycf36dXzzm9/Ej3/8Y8XW8+SCYskn1xch1fpA5rBbrVZ0SwwHLNQoWgrx5EM0GhU6pzP9drGLEQ6H4Xa7hZKBzCpkIr0hht/vx8bGRsUSFlLweDzw+/1VITW73Q6O4wpWGRcLnU6Ht/76XQCAxcVFbE47EPOlkIgmgNcWsdalhS/Jg5p04Kn/8X2sra3hxF3dkn1xuUiFFIMeO3ZM9h42MUKhED7+8Y/jm9/8puxNugcBiiUfYJeAMqcCSDXlzc/Pw2KxSBIPkN/ykYt4xGJgUgHDfC5GZhUyRVGgaRp6vR5WqxU6nQ5erxcjIyOyx48CcmrWCgAAFcRJREFUgYBAanLvmzSLFhIbKwculwvxeBy//4Xfg0ajwTc+8104l9044Y5joU2HYz4GPnsYlD0FTb8RJ37hBC7+62Ukk0n0njsmVHqToLeY/Hd2dqDX62V3P8XgOA4PPfQQPvnJT2JwsHIX8SBC0eRjtVpBUVTOpw8hnkJVw7nIRy7iicfjZYuB5atCpmkafr8fi4uLaGlpwcbGhmT6miycUt0DiqKwsLCAiYkJWYXqgb1m0WqQWigUwubmZlospv9sL66/MA3vVgiH4xpsH9LjsJfB7E9uorbBivP3nRa0q4eGhoR9MQwjEH84HMbm5ibC4TDMZjOuX7+eZS1VOteM4B/+4R9QV1eH973vfRXv66BC0eRDBMWkyIekm2tqagq2K0i5amQoYaV1PMlkEpOTkxgcHJRdDIzjOKytrWF8fDztGojT17FYDE6nExRFZXW4kx8paYhkMikQptypXXGzqNzjZJLJpCBklkmYHHjo9FpoOB62IIuoWQNDjQH9Z3vTesUuP3ddSNnr9XrU19ejvr4eyWQSXq8Xd999N2pqatKC3qT9RJztLLehdGpqChcuXLgt4zxiKJp8ctX6kCCm0WhEb29vwS8wM9XOsix4ni9qfno+kGkX/f39FU2EkALDMMK+M8mXpK+lCt7E5QFSxX6EjOx2O44fPy5ZMFgJ8jWLVgoihD8wMCBpKZ4824fwThgpmoeOSiJp0mD4HcNprRr59j09PY2BgQGBjHMFvcXjqSmKEhpKyeyuzEpvMflHIhF87GMfw4ULF6pasHgQoGjyEWv6EBDtFa1WWxTxAOluF5kmWinxkJ6qYgYLlgpS6NfV1VXyvg0GAxoaGrKK/cR6QMvLyzAYDFhfX8fy8rIwz13swhUzbysT5TSLlrLv+fl5tLW1Se6bWDIhXwR+VwBavQ4t7Y1geE54nRCQlO7P0tISWlpaijrvfNlOQv5EF4iIwfE8j8ceewzhcBhnz55FLBZDPB6/LQPNBIomn0wReZ7nsbi4CK1Wi66urqIXByEfuYiHkMPRo0dlr3qtVtsEkYbY2NjAkSNH0lxVcWsE0QRKJBKSjaQkAC6FQs2ilcButwNA3oruc/eeEkjm+gu7+j4T5/sL7tvtdiORSKC/v/B784HUiBmNxqxeQpZlcc899+D1119Hf38/nnzySXziE5+QXabkIOG2IR+e57G0tAStVgubzVZyrxYZ7KfX6ysiHmKet7W1lSQGVgxIHKu2tlb2tglgV28YQFZWMFdrhNhaIsRExiYT94IQE0VRRTWLlgO/3w+Px1N0sZ/YshFbPFIxn2g0ivX19bJ6zUrB/Pw8nnnmGfz4xz+W3R09qFA8+USjUfA8j+XlZQDA8ePHcfPmzZLIh1g8ZGGLg7GlEBFxKxoaGmRX9AN22yYAyN42AeympoPBIMbHx4teZLnUEzPdC6fTCZ/PB5PJhMuXL6cFYitV94vFYlhYWCg5a1bMPK/XnrmC5ZUV/NqH3iV7tk+MaDSKj370o/j617/+piEe4DYhn5WVFXAch/7+fmg0mpJkNQjx3HHHHWnBWPFTnIg7kYVSW1ubNV+LWCVmszlnPVElsNvtoChK9u53AIIrJVfaW+xe1NTUwG6346677hKkOIgLl6nulxmMLZS6ZhgG09PTGBoaqkjThhCROObD8zw2NzfR1tYme9BdDJ7n8Ud/9Ef46Ec/itHR0aod5yBC8eRz9epV3HXXXWlaOKXIapAiRSI1KvUUJxMjKIrC9va20EckbokIhULQ6XQVxwWkQNomqlGMF41Ghdnvcj/daZrG9PQ0hoeHhQxRLtmNzMkcpF+LJA8yrSWTySTo2sidSQSAH3z9eSQjNHhLuhsmNy5cuACNRoPf/d3frcr+DzIUSz48z+OZZ56B0+nE4OBgUf1dYhAVwkLz0/NNjCAtEevr64jFYjCbzbhy5QqAvbR1uS4cgbhtQu5ivEQigZmZGYyOjsquhperWTQX8l3nzNT19vY2AoGA8PAIh8Nly5qIQQjG5/MhEolUPO22EObm5vDEE0/c9vU8uaBY8nnjjTcQi8UwMTGRtaiLkdUohngKQa/XIxQKged5nD9/XiAHMtFAHIglin7i+U+FXItwOFw1q4QMRRwYGJA9zlBOs2g+ZKauPR4PGIbB+Ph4WiaOyJqIraV8siZSiMfjWFxcxHse/FVMvjgHoPR578WAoig8+OCD+NrXviZ78alSoJEaci9C3hf3G6+//jr+7u/+Dl/+8pfT/u5yucCyrGTQl4y4Ia5WJdja2oLb7S5aDCzTtSDkRBaLeKFotVosLS1hfHxc9mIzjuNw48YNtLe35x0RVC42NzcRiUQwNDQk+xM9EolgdnYWZ86cyfv9iWVNyPUmAviZMTzyAGBZFteuXcOJEyeyMntygud5fPSjH8U999yDBx98sGrHOUCQvAkUa/kAu7O7cqkZJpPJrL/LZfEAuyJSW1tbOHXqVNHuUCHXQmwpuVwuGI1GTE1NZblw5YqyA3vFeE1NTVUhnmo2i6ZSqZwzvDKRr9CPtEUQa4k8AFKpFMxmM4LBIFKpVFmyJsXgm9/8Jmiaxoc+9CFZ96s0KJ58xEWGBFJul5zE4/P5sL6+LuvoFZ1Oh7q6OpjNZtjtdoyPj6OpqSmtTyvThRPX0hSTHQJ20/U6nU42CQsxqtksSmJIfX19FbmJubSASDlAV1cXYrFYlghcMbImxWB+fh5f+cpX8PLLL8t+jZQGxZNPMTrO4qF+lRJPMBjE0tISTp8+LbuIFOl7ErdNiPu0xNXSmbU0UtkhMSmZzWZhKoTcKodAdZtFgV1tnqampqroJIfDYWHkjU6ny8qe5ZI1IfdUJimZTCZJYonFYnjwwQfx1a9+tSoZOqVB0eRjtVolZ3SJyUcc46k0aBuJRHDz5k1MTEzIvsBI20RbW1tRbROFSvUzZVwjkQiSySSampqwsrKStlgqvS7VbBYFdq0Smqarol+cSqUwNzeXd+RNoeGKhJQyh/8RWROapuF2u/Htb38bv/M7v4PTp0/L/jmUCEWTTy6zlZCPXGOMgT0xsPHxcdklJkgcxmazyVIZnRnvIEHae+65Jy22FAwGhSd4pguXOSUz37lXq1kUgDDRoho6yXKMvBFLbmTum7jLc3Nz+NKXvoS5uTlMT0/j6aefxvPPP18VC1FJUDT5EGTqOOv1ejAMI8sYY2A3/To1NYXR0dGqyBysrKxAq9XKPiYZ2KvlGR8fF2p5Ml04AGmFlF6vV9CmyeXCketZzWbRRCKB+fn5qow1Bqo78kbsLh8+fBg+nw9zc3NoaGhAKBR60xMPoHDyyfUk1Gg0SCaTSCaTsNlssoiBDQ0NVaXMfnNzE7FYrCptE0RPaHBwsCBpluLCkSAseb27uxuhUEgWF0583KmpKZw8ebIqc6q8Xm/VR94Auw+uD3/4w3jyySeF9L0a79mFoskH2HUxiEIfsBvj4Xkera2tabKimU/vYpoZaZrGjRs3qiIGBuxKNezs7FRFmJ0Ia/X09FRUs5IrZR0IBLCwsIDjx48jkUgIaoksy2aNRrZarUW5cAREoL29vb2kMdbFIhaLYWVlpeojb3iex5/8yZ/g/e9/P86ePVu14ygViicfIqXa1NQkxHiIkJi44licrhY3M9bU1GR1shO37caNG+jt7a1KLIMISVUjLU0W76FDh6rSIhCLxXDz5k2cOnVK0ioRVx2LBbPELpz438zPv7GxAYPBUBVlAJZlhWbUao68AYCnn34agUAAH//4x6t6HKVC8eRDNH0aGhoE4snUXc4n/SBugyBPb5qmkUwmUV9fj0QiAb/fL6s4eCgUEtL11ZBqWFlZgdFoTBuMKBekmkUzkUv/R+zCRaNRQfdYXEfDcRyCwWBV3KFbNfIGAJaXl/HFL35RrefJg9uCfKLRKFiWlSSefBCnUFtaWgDsuSsdHR2oq6tLq6FJJpPCuJpCAuy5QLIf1UjXA4DD4RBiSHKj1GbRTOSrOk4kEvD5fFhdXUVjYyOmpqbAsqzk9S7FhRPDbrdXfeQNsBso/9CHPoQnnniiKm7j7QLFk09tbS3C4TA0Gk3FkyZI2rixsVGoAM68eUhdB0VRCAQCcDgcgqSoVFxJfD6JREKYCFENbV6v1wu3212V1gbSLNrY2Ci7K6fRaKDX6+FwOHDq1Kk0csqleazRaLLiSvkaR4PBoKB2WE3wPI8/+7M/w3vf+16cO3euqsdSOhRPPi6XC3//93+P9fV1DA4OCl3a5Yibk+GC+VoPctV1iEWyotEoPB6PkBUym80wmUzwer3o7u6uSro+FAoJQdRqpKXJZNFqCKWRepuenp4sqyiX4D3LsmlxPCkXjhCTXq8XUvbVdoG+//3vw+Px4Ctf+UpVj3M7QNFd7cDuLOuFhQVsbGxgbm4Oi4uLoCgKhw8fxsDAAAYGBnDy5EmcPHkSjY2NkqRE9J95nhfUEOUCz/OIRqOYmZmBzWaDRqORVEgkQe9ygqDxeBw3btzAxMRE1Syqzc3Nqi3epaUlaDQaHD9+vOJ9ERcu80FgMBgEUqp0CkcurK6u4rd+67fw8ssvVyVJoWBIXmDFk48UOI6Dy+XC7Ows5ubmMD8/j/n5eYRCITQ2NmaR0re+9S3cfffduOOOO6qW8m5paUnL3mR2V5MfmqaFVHUxcQ6apnH16lUMDg5WpRyAVEefPn26KjEql8sFj8dTknZ0KVhYWEBNTQ26u7vTXLhYLIZoNCq4cIVc5kJIJpO499578YUvfAF33XWX7J9D4XjzkE8u8DyPnZ0dzM7OYn5+HnNzc3jhhRdA0zTa2trQ39+fRkrHjh2rOIY0NzcHs9lckui7eJGQn0QikSbbSiqNb968ia6urqpU6SaTSVy7dg1jY2NV6dkKh8OYn5/HmTNnqpL1c7vd8Hg8BRtpxS6zmJwytX9ydbPzPI9HHnkE3d3dePjhh2X/HLcBVPLJhMPhwEMPPYRvfetbQin/3NycYC05HA4YjUYcP34cJ0+eFIipp6cHer2+4JN6eXkZDMNgYGBAlqe6uC8rGo3C5XIBQNZQP/JTCXESYa2+vr6quBCE2KohlgZAcHXPnj1bkfaRlHVKpFksFgtefPFFhEIhXLp0Cc8//3xVp1woGCr5SCGzLyzztXg8joWFBcFamp+fx9raGjQaDXp6etIspePHjwvu0Y0bN6DVaqvSNgHsSkxoNBqcOHEiq4iS/IiDr5lFlIWuyfT0NJqbm6vSs8VxHK5du4aenp6sUcNygLiiIyMjVZs8QdM0YrEYnnjiCXz7299GZ2cn3G43PvGJT+CDH/xgVY6pYKjkIxeIls7y8nIaKS0vLwsKeKlUCg8++KBgMZWTgcuFzc1NhEIhjIyM5N1nZhGl+Mmd2cUuHge0vLwsjCKSGySrWFtbW5UiSJ7nMTU1hba2tqooNYqRTCZx33334bHHHsPP/dzPVfVYCodKPrcCL730Eh5++GE8/PDDWRm4I0eOoL+/P82Fy5WBy4Xt7W3Y7faKMk+EPKPR/7+9+3dpI4zDAP4EkaKN5RQtaJOlGJPgILp06FBxyVJBSktp/4mOFSJWDILg0kFxcHBzqVuxBkVHEYqWFnpGU4jgj0oEozEeGBOvQ7nrXUzU2Ht7NXk+4KjRwYe79/vjTZlCKZ1O6zdCuN1uvQJnZUVoa2sLyWRSyH5nAIjFYkin00J2/xhp/TxNTU3o7e0V+lklgOHzLxwcHCCbzV44AD4/P8fu7i5kWcb3798RiUQgyzKOjo5QV1d3oQJ3//79C+FyeHio384pYi4pkUhgY2MDHo/H9MSkNVHmvr4Vu9/44OBA70USUbI3rrcV3c8zMzODyclJfPz4UUhfVYlh+PyPjBU4Y1tAPB7H3bt39VCSJAmLi4sYHx8XckCrKAq+fv1acFg09zaIk5MTUxPlVXeUab1IHR0dlt8R9i9+vtHW1hZevHiBhYUFIVXGEsTwuU1UVdVL0cvLyxgeHkZ7ezt+/vyJyspKeDwe/UnJ6/Xi4cOHN660nJ2dYXV1FX6/v+iZrdwF99pXNpvVD7u1pfher1fIAfP5+TlWVlaEX3kD/J7Yf/r0KYaGhvDkyROhn1VCGD631adPn1BRUYFAIFCwAre5uQkAegVOGzUxVuDy0e7wcrlcls5saWXqVCqFaDSqvyYW20R5HbIso6amBm6327LfPx9VVdHX14f6+noEg0Ghn1ViGD6lLLcCJ8syIpEIotEoMpkM3G636VzJ6/WiqqoKHz58wKNHj4SscAV+jxzkLn/XmiiNB96np6cXmii1p6bLQmlnZweJRAKtra3CrxyenZ3FxMQEZmZmeM5THIZPucpms4jFYlhbW9OfljY2NrC3twdJkvD48WP4/X79sFuSJEv+kePxOLa3t699o6uxidJ4rgT82USpVeCqq6uRSqUQiUSEDdMabW9v4/nz55ifn7/W7SJkwvChPz5//oze3l5MTEzoT0taBS6ZTKK2thY+nw8tLS36K1y+ClwhWofxVdcaX0e+8YdUKgVFUSBJEu7du1dUE2Wxzs7O0N3djYGBAXR1dVn6s8sEw4f+0P6h83UAX1aBczqd8Hq9eij5fD64XC5TKKXTaayurgrrMFZVFV++fIHL5UJNTU3eJsp8M1k3GYxVVRUDAwNwOp149+6d5X9LmWD40N9RVRVHR0emGbhIJKLPwHk8HjQ3NyMcDqOvrw+dnZ1CZp1+/PgBAAVXcKiqaroKyNhEWewmyrm5OYyNjSEcDvOc5+YYPiSGqqpQFAXr6+sIBoM4PT2FJEmmGTjt1U2bgbvpPuz9/X29w/sm36/NZBlf34ybEbWvZDKJyspKvHz5EnNzc8JHNUocw4fEymQyGB0dxZs3b+BwOPQKXDQaNW0L0Kb93W636VxJu265UKgoioJv375Zco6UK7eJMhQKYWlpCdXV1Whra8OrV6/w7NkzSz+zjDB86P+RrwK3vr4ORVHQ2Nh4YdykoqICU1NTeP36tfCbJ1RVxeDgIO7cuYP+/n79Cc6KTYtliuFD/z9tBs542C3LMmKxGDwej+lJyefzoaGhwfI5roWFBbx//577eazD8KHbaXp6GvPz8wiFQvrrm3bYvb+/b5qB8/l88Pv9ePDgwY1CaW9vDz09PQiHw8Kv2CkjDB+6nTKZjL5wP9d1KnD5tlAW+pyenh68ffsWgUBA9J9VThg+VD6MFTjjYffm5mbBCtzIyAgcDgeGhobs/vVLDcOHyFiBM14ksLa2huPjY0SjUZ7zWI/hQ3SZbDbLRkIxGD5EZIu84SN21yQRUQEMHyKyBcOHiGzB8CEiWzB8iMgWDB8isgXDh4hswfAhIlswfIjIFgwfIrIFw4eIbMHwISJbMHyIyBYMHyKyxVVbk/7+wm4iojz45ENEtmD4EJEtGD5EZAuGDxHZguFDRLZg+BCRLX4Be5N1HKf++fcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "elev = 30\n",
    "azim = 20\n",
    "plot_figs(2, elev, azim)"
   ]
  },
  {
   "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
}

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()

plot_pca_3d.ipynb

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()