{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Polynomial Chaos Expansions (I)\n",
"\n",
"The main idea behind [PCE](https://econpapers.repec.org/article/eeereensy/v_3a93_3ay_3a2008_3ai_3a7_3ap_3a964-979.htm) is to search an interpolator that lives in the subspace of low-degree polynomials (or rather, in the tensor product of such subspaces).\n",
"\n",
"In this notebook we will tackle a regression problem: we have $N$ features that are used to train an $N$-dimensional TT model. Each feature $x_n$ is mapped to the $n$-th entry of the tensor; in other words, we learn a function as a discretized tensor:\n",
"\n",
"$$f(x) = f(x_0, \\dots, x_{N-1}) \\approx \\mathcal{T}[i_0, \\dots, i_{N-1}]$$\n",
"\n",
"Here we will use a noisy 5D synthetic function: $f(x) := \\sum w_n x_n^2 + \\epsilon$ where the $w_n$ are random weights and $\\epsilon$ is Gaussian noise added to every observation."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import torch\n",
"torch.set_default_dtype(torch.float64)\n",
"import tntorch as tn\n",
"\n",
"P = 200\n",
"ntrain = int(P*0.75)\n",
"N = 5\n",
"ticks = 32 # We will use a 32^5 tensor\n",
"\n",
"X = torch.randint(0, ticks, (P, N)).double() # Make features be between 0 and ticks-1 (will be used directly as tensor indices)\n",
"ws = torch.rand(N)\n",
"y = torch.matmul(X**2, ws)\n",
"y += torch.randn(y.shape)*torch.std(y)/10 # Gaussian noise: 1/10th of the clean signal's sigma\n",
"\n",
"# Split into train/test\n",
"X_train = X[:ntrain]\n",
"y_train = y[:ntrain]\n",
"X_test = X[ntrain:]\n",
"y_test = y[ntrain:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our first attempt will be to learn this problem using only the low rank assumption, i.e. plain [tensor completion](completion.ipynb):"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter: 0 | loss: 0.999345 | total time: 0.0020\n",
"iter: 500 | loss: 0.884514 | total time: 0.8227\n",
"iter: 1000 | loss: 0.130079 | total time: 1.7812\n",
"iter: 1500 | loss: 0.019054 | total time: 2.7306\n",
"iter: 2000 | loss: 0.009741 | total time: 3.7072\n",
"iter: 2500 | loss: 0.005941 | total time: 4.7091\n",
"iter: 3000 | loss: 0.003692 | total time: 5.7308\n",
"iter: 3500 | loss: 0.002271 | total time: 6.8586\n",
"iter: 4000 | loss: 0.001340 | total time: 7.8531\n",
"iter: 4500 | loss: 0.000724 | total time: 8.8624\n",
"iter: 5000 | loss: 0.000347 | total time: 9.8511\n",
"iter: 5500 | loss: 0.000153 | total time: 10.8660\n",
"iter: 5743 | loss: 0.000100 | total time: 11.3829 <- converged (tol=0.0001)\n"
]
}
],
"source": [
"t = tn.rand(shape=[ticks]*N, ranks_tt=2, requires_grad=True)\n",
"\n",
"def loss(t):\n",
" return tn.relative_error(y_train, t[X_train])**2\n",
"tn.optimize(t, loss)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This TT did very well for the training data, but clearly overfitted:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Test relative error: tensor(0.4168, grad_fn=<DivBackward1>)\n",
"The model overfitted: it has 512 degrees of freedom, and there are 150 training instances\n"
]
}
],
"source": [
"print('Test relative error:', tn.relative_error(y_test, t[X_test]))\n",
"print('The model overfitted: it has {} degrees of freedom, and there are {} training instances'.format(tn.dof(t), len(X_train)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can look at the groundtruth vs. prediction for the training and test splits:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X+8VGW59/HPxXZjW9Q2JHpwK4I+ptmDv85O6MHHzFLUjoo+x8yw0Cx7PNrRc5TEpJTKR5STv069LDlpdjL80dEJUyNKOxYpAm5wS8oBFYGtKR5FTXa5gev5Y9bAsJlZM7Nn1loza77v12u/9sw9a2bugTX7Wvev6zZ3R0REpBqDkq6AiIg0PgUTERGpmoKJiIhUTcFERESqpmAiIiJVUzAREZGqKZhIQWbWYmZ/NrORSddFROqfgklKBH/4cz+bzaw37/6kSl/P3Te5+87uvjqK+opUqtbneN7rPmFmZ9Wyrs1oh6QrILXh7jvnbpvZKuCL7v7rYseb2Q7uvjGOuonUQqXnuMRLLZMmYWbfNrO7zWy2mb0DnGVmHw2uytab2StmdrOZtQbH72Bmbmajgvs/CR5/2MzeMbPHzWx0gh9JZBtB1+zXzewFM3vdzO40s/bgsSFmdpeZvRGc7wvMbKiZfQf4CPBvQQvnO8l+isalYNJcTgV+CrwfuBvYCFwE7AaMB44Hvhzy/M8CXweGAauBb0VZWZEKTQGOA44E9gL6gBuCx75Itiemg+z5fiHwnrtfAiwk28rZObgvA6Bg0lx+7+4PuPtmd+9194XuvsDdN7r7C8CtwMdCnv8zd1/k7n3AncChsdRapDxfBqa6+8vu/hdgOnCGmRnZwDIc2C843xe6+7tJVjZtNGbSXNbk3zGzA4HvAH8L7ET2fFgQ8vw/5d3eAOxc7ECROAUBY2/gITPLz147CPgA8EPgb4CfmdnOwI+Br7v7ptgrm1JqmTSX/imifwA8A/wPd98V+AZgsddKpEqeTX/eAxzj7u15P+9z99fd/a/u/g13PxA4Cjgd+Ezu6UnVO00UTJrbLsBbwLtm9iHCx0tE6t33gRlmtjeAme1uZicFtz9pZgeZ2SDgbbLjhblWyavAvklUOE0UTJrbJcBk4B2yrZS7k62OSFWuA34NPBLMWPwDcHjwWAfwc7Ln+jPAQ8A9wWM3AJ83szfN7Lp4q5weps2xRESkWmqZiIhI1RRMRESkagomIiJSNQUTERGpWioXLe62224+atSopKshKbZ48eLX3X143O+rc1uiVM15ncpgMmrUKBYtWpR0NSTFzOylJN5X57ZEqZrzWt1cIiJSNQUTERGpmoKJiIhUTcFERESqpmAiIiJVS+VsLpFqZLp6mDl3OS+v72XP9jamTDiAiYd1JF0tSYDOhfIpmIjkyXT1cPl93fT2ZbOT96zv5fL7ugH0R6TJ6FyojLq5RPLMnLt8yx+PnN6+TcycuzyhGklSdC5URsFEJM/L63srKpf00rlQGQUTkTx7trdVVC7ppXOhMgomInmmTDiAttaWbcraWluYMuGAhGokSdG5UBkNwIvkyQ2sagaP6FyojIKJSD8TD+vQHwwBdC5UQt1cIiJSNQUTERGpWmTBxMz2NrNHzexZM1tmZhcF5VeZWY+ZLQl+Tsx7zuVmttLMlpvZhLzy44OylWY2Nao6i4jIwEQ5ZrIRuMTdnzKzXYDFZjYveOwGd/+X/IPN7CDgM8CHgT2BX5vZB4OHvwccC6wFFprZHHf/Y4R1FxGRCkQWTNz9FeCV4PY7ZvYsEDaSdQpwl7v/FXjRzFYCRwSPrXT3FwDM7K7gWAUTEZE6EcuYiZmNAg4DFgRFF5rZ02Z2m5kNDco6gDV5T1sblBUr7/8e55nZIjNbtG7duhp/AhERCRN5MDGznYH/AC5297eBW4D9gEPJtly+kzu0wNM9pHzbAvdb3b3T3TuHDx9ek7qLiEh5Il1nYmatZAPJne5+H4C7v5r3+CzgF8HdtcDeeU/fC3g5uF2sXERE6kCUs7kM+CHwrLtfn1c+Iu+wU4FngttzgM+Y2Y5mNhrYH3gSWAjsb2ajzWww2UH6OVHVW0REKhdly2Q88Dmg28yWBGVfA840s0PJdlWtAr4M4O7LzOwesgPrG4EL3H0TgJldCMwFWoDb3H1ZhPUWEZEKRTmb6/cUHu94KOQ5VwNXFyh/KOx5IiKSLK2AFxGRqimYiIhI1RRMJLUyXT2Mn/EIo6c+yPgZj5Dp6qnJ665Zswbgg0oVJLKVUtBLKmW6erj8vu4te3j3rO/l8vu6AapOKb7DDjsArHX3g5QqSCRLLRNJpZlzl28JJDm9fZuYOXd51a89YsQIgA2QTRUElJ0qyN1fBHKpgo4gSBXk7u8BuVRBIg1HwURS6eX1vRWVD1QcqYJkYKLq5pTCFEwklfZsb6uofCDiShWkvHOVy3Vz9qzvxdnazamAEh0FE0mlKRMOoK21ZZuyttYWpkw4oFZvYRRIFeTum9x9MzCLrVmvi6UKCkshtIXyzlUuym5OKUwD8JJKuUH2mXOX8/L6XvZsb2PKhANqsp+3uwPsA/yuf6qgYOsF2D5V0E/N7HqyA/C5VEFGkCoI6CE7SP/ZqisosXVzxiXT1RPJuVxLCiaSWhMP64jkCzd//nyADwDHKFVQfdqzvY2eAoGjlt2ccYlyZmItKZiIVOjII48EWOzunf0eUqqgOjFlwgHb/AGGmndzxiasy07BREQkQlF2c8atUbrsFExEJJWi6uaMW6N02Wk2lzQsrSOQZhDDzMSaUMtEGlKjDEqKVKtRuuwUTKQhNcqgpEgtNEKXnbq5pCE1yqCkSLNQMJGGFEe6FBEpn4KJNKRGGZQUaRYaM5GG1CiDkiLNQsFEGlYjDEqKxCnJHF4KJiIiKZD0dHkFE6l7jZAxVSRpSU+XVzCRujYt082dT6zesmOUFieKFJb0dHnN5pK6lenq2SaQ5GiTI5HtJT1dXsFE6tbMucu338M2oMWJIttKerq8urmkboUFDC1OFNlW0tPlFUykLmW6ehhkxibfvm1ioMWJIgUkOV1ewUTqSnZ649P09m0u+LgBk8aN1OC7SJ2JbMzEzPY2s0fN7FkzW2ZmFwXlw8xsnpmtCH4PDcrNzG42s5Vm9rSZHZ73WpOD41eY2eSo6izJynT1MOXepUUDSYsZN5xxKN+eOCbmmolIKVEOwG8ELnH3DwHjgAvM7CBgKvAbd98f+E1wH+AEYP/g5zzgFsgGH+BKYCxwBHBlLgBJusycu5y+zcWG3GGzu1okInUqsmDi7q+4+1PB7XeAZ4EO4BTgjuCwO4CJwe1TgB971hNAu5mNACYA89z9DXd/E5gHHB9VvSU5hbYmzadBd5H6FcvUYDMbBRwGLAD2cPdXIBtwgN2DwzqANXlPWxuUFSvv/x7nmdkiM1u0bt26Wn8EiVimqwcLeVyD7iL1LfJgYmY7A/8BXOzub4cdWqDMQ8q3LXC/1d073b1z+PDhA6usJCLT1cMl9ywtuqYENOguUu8iDSZm1ko2kNzp7vcFxa8G3VcEv18LytcCe+c9fS/g5ZBySYFccrpCU4BzbtSgu0jdi3I2lwE/BJ519+vzHpoD5GZkTQZ+nlf++WBW1zjgraAbbC5wnJkNDQbejwvKJAUKJafL19HephaJSAOIcp3JeOBzQLeZLQnKvgbMAO4xs3OB1cDpwWMPAScCK4ENwDkA7v6GmX0LWBgc9013fyPCekuMwla5a+dEkcYRWTBx999TeLwD4BMFjnfggiKvdRtwW+1qJ/Viz/a2grO4Wsy45rQxapWINAitgJfY5e9P0r5TK62DbJv1JW2tLQokIkXU6/4+CiYSq0xXD1N+tpS+Tdng8eaGPloGGe1trbzV21dXXw6RepP0bophlIJeYjX9gWVbAknOps2OGbw441PMn3pM4l+KUtasWQPwQaUKkriF7aaYNAUTidWbG/oqKq9HO+ywA8BapQqSuCW9m2IYBROJTaarJ+kq1MSIESMgO+NQqYIkVknvphhGwURikevrLaa9rTXG2tSOUgVJnJLeTTGMgonEImxxYusg46qTPxxzjaqnVEESt4mHdXDNaWPoaG/DyC7qrZeZj5rNJbEI69OdefohdfFlqJBRJFWQu79SQaqgo/uV/zbKSkvjS3I3xTBqmUikMl09jJ/xSNEkjo2YLiW7vpZ9UKogkS3UMpHITJr1OPOfL575pl76eis1f/58gA8AxyhVkEiWgolEYlqmOzSQdDTw4sQjjzwSYLG7dxZ4WKmCpCkpmEgk7nxidejj86ceE1NNRCQOGjORmpuW6Q7d6EpE0kfBRGoq09VTslUydKfGXFMiIsUpmEhNTX9gWclWyZUnNd6aEhEJp2AiNZPp6imZY+ss7eUukkoKJlIz0x9YFvr4WeNGai93kZRSMJGamDTr8dBWiQKJSLpparBU7djrf8uK194t+nh7W6sCiUjKqWUiVZmW6Q4NJEBDJnEUkcqoZSJVmb1gTejj7W2tGnAXCdTr/u21oGAiVdnk4ROB1SoRyarn/dtrQcFEqtJiVjSgDG6xVHxJpDppvhqvRNj+7Wn491AwkQGZlulm9oI1oS2T6/7+kBhrJPUo7Vfjlajn/dtrQQPwUrFpmW5+8sTqooFkp9ZB3HjGoU33x0K2F3Y13mzqef/2WghtmZjZP4c93m9jIGkCma4eflIk91aLGc9fc2LMNRqY668PP3X/+Z9DT30pU9qvxisxZcIB27TSoHH39CmkVDfXLsHvA4CPkN0xDuAk4LGoKiX1aVqmOzSJY6nB+HryzjvvALB8+XIWLlzIySefDMADDzzAUUcdlWTVUmXP9jZ6CgSOtFyNVyLXUk/r+FFoMHH36QBm9ivgcHd/J7h/FXBv5LWTupHLBhwWLlrMYqtPta688koAjjvuOJ566il22SV73XTVVVdx+umnhz1VKpD2q/FK1ev+7bVQ7gD8SOC9vPvvAaNqXhupW1+77+mS2YDPHLt3LHWppdWrVzN48OAt9wcPHsyqVauSq1DKpP1qXLYqN5j8O/Ckmd0POHAq8OOwJ5jZbcDfAa+5+/8Myq4CvgSsCw77mrs/FDx2OXAusAn4R3efG5QfD9wEtAD/5u4zyv50UhPTMt1s6Nscekyj5t763Oc+xxFHHMGpp56KmXH//ffz+c9/PulqpUqar8Zlq7KCibtfbWYPA/87KDrH3btKPO1HwHfZPujc4O7/kl9gZgcBnwE+DOwJ/NrMPhg8/D3gWGAtsNDM5rj7H8upt9RGsQF3AANuaOCZW1dccQUnnHACv/vd7wC4/fbbOeywwxKulUjjqWSdyU7A2+5+u5kNN7PR7v5isYPd/TEzG1Xma58C3OXufwVeNLOVwBHBYyvd/QUAM7srOFbBJCaTZj0e/ngK9ifZsGEDu+66K+eccw7r1q3jxRdfZPTo0UlXS2pEiybjUdY6EzO7ErgMuDwoagV+MsD3vNDMnjaz28xsaFDWAeQneVoblBUrL1TH88xskZktWrduXaFDZADmP/9G0cfMaMiurXzTp0/n2muv5ZprrgGgr6+Ps846K+FaSa3kFk32rO/F2bpoMtPVk3TVUqfcRYunAicD7wK4+8tsnTZciVuA/YBDgVeA7wTlhaYBeUj59oXut7p7p7t3Dh8+fABVk/7GXj0v9PFJY0fGVJPo3H///cyZM4chQ4YAsOeee26ZNiyNT4sm41NuN9d77u5m5gBmNmQgb+bur+Zum9ks4BfB3bVA/lSgvYCXg9vFyiVC0zLdvPrOe6HHNHqrBLKzt8wMC6Y1v/tueDp9aSxaNBmfclsm95jZD4B2M/sS8Gvg3yp9MzMbkXf3VOCZ4PYc4DNmtqOZjQb2B54EFgL7m9loMxtMdpB+DhK5Uqnlx+83LKaaROvTn/40X/7yl1m/fj2zZs3ik5/8JF/84heTrpbUSNpTmNSTcmdz/YuZHQu8TXY1/DfcPbQPxMxmA0cDu5nZWuBK4GgzO5RsV9Uq4MvB6y8zs3vIDqxvBC5w903B61wIzCU7Nfg2dw/faFyqlunqKbma/c4vfTSm2kTr0ksvZd68eey6664sX76cb37zmxx77LFJV0tqRIsm41NWMDGza939MmBegbKC3P3MAsU/DDn+auDqAuUPAQ+VU0+pXqkteCE9rRKAyy67jGuvvXabAJIrk8anRZPxMS8jn5KZPeXuh/cre9rdD46sZlXo7Oz0RYsWJV2NhjP26nklx0n22GUwC65Iz5X74YcfzlNPPbVN2cEHH8zTTz8d+jwzW+zunVHWrRCd2xKlas7r0DETMzvfzLqBA4PpvLmfF4Hugbyh1KdyBtzPGjcyNYHklltuYcyYMTz33HMcfPDBW35Gjx7NmDHhEwu+8IUvABxiZrkxP8zsKjPrMbMlwc+JeY9dbmYrzWy5mU3IKz8+KFtpZlMj+JiSgExXD+NnPMLoqQ8yfsYjTTMNuVQ310+Bh4FrgPyT/R13L74AQRpOWDbgnDTM3sr57Gc/ywknnMDll1/OjBlbM/TssssuDBsW3o139tlnc/vtt68o8JCyOzS5Zt4MLLRl4u5vufsqsrmx3nD3l9z9JaDPzMbGUUGJXqarp2QSx6E7tcZSl7i8//3vZ9SoUVx00UUMGzaMffbZh3322YfW1lYWLFgQ+twgRf3GMt9qS3aHIGNELrvDEQTZHdz9PSCX3UEaWDOvayl3avAtwJ/z7r8blEkKXDUnfIKcAVee9OF4KhOz888/n5133nnL/SFDhnD++ecP9OWU3aHJNfO6lnKDiXneSL27b0b7x6fCtEw363v7Qo9p5ESOpbj7lgWLAIMGDWLjxnIbHdtQdgdp6nUt5QaTF8zsH82sNfi5CHghyopJ9HJ7uYdpb2tNbSAB2Hfffbn55pvp6+ujr6+Pm266iX333bfi13H3V919U3ChNYutiUqLZXcIy/ogBTTCwPaUCQfQ1tqyTVmzrGspN5j8X+B/AT1kvwRjgfOiqpTEo1QgAbjq5HR2b+V8//vf5w9/+AMdHR3stddeLFiwgFtvvbXi11F2h2g1SsLGiYd1cM1pY+hob8OAjvY2rjltTKovyHLKXQH/GtmTXVKiVGp5yE4FTvuXYPfdd+euu+6q6DlnnnkmwIGAKbtDPMIGtmt1jtYqVX2zbgYWGkzM7Kvufp2Z/SsF+nPd/R8jq5lEJtPVE5paHuDGFI+TAFx33XV89atf5Stf+co2YyY5N998c9Hnzp49m7vuuuvpfou7lN0hQlEPbDfzlN5aKdUyeTb4rSW3KTL9gdIXwGn/An3oQx8CoLMz9kXsMgB7trfRUyBw1GpgO46WT9qFBhN3fyD4fUc81ZGoTct08+aG8NlbZ41r/H1KSjnppJMAmDx5csI1kXJEnbCxmaf01kqpbq4HKDJdEcDdT655jSQyk2Y9XrJ7a//dh6RqpXsxJ510UsHurZw5czQWXk+iTtgYdcunGZTq5sqlhjgN+Bu2btV7JtlBRmkQ5QSS8fsNS01q+VIuvfRSAO677z7+9Kc/bdmqd/bs2YwaNSrBmkmxgfAoB7anTDiAKfcupW/z1mvn1kHWFFN6a6VUN9d/ApjZt9z9qLyHHjCzxyKtmdRMOQPu7W2tTRNIAD72sY8B8PWvf53HHtt6Kp900km5dCmSgEQHwvs3VIs3XKWActeZDDezLSu5gvnyWorbIC6/LzydOqR/PUkx69at44UXtq6/ffHFF1HKkuQkldtq5tzl9G3atke/b5M3RU6tWik3Jco/Ab81s9y3bhTBPHqpb9My3fT2bQ49ZsjglqadsXLDDTdw9NFHb1n1vmrVKn7wgx8kXKvmldRAuAbgq1fuosVfmtn+ZBdqATzn7n+NrlpSC5munrJWuV99avoH3Is5/vjjWbFiBc899xwABx54IDvuuGPCtWpeSQ2EawC+emV1c5nZTsAU4EJ3XwqMNLO/i7RmUrWvldG9lfbFiaVs2LCBmTNn8t3vfpdDDjmE1atX84tf/CLpajWtpHJbNXNOrVopd8zkduA9IDdCuxb4diQ1kpqYlulmQ4nurWYPJADnnHMOgwcP5vHHs+ll9tprL6ZNm5ZwrZpXUrmtmjmnVq2UO2ayn7ufYWZnArh7r4VN0pdEldO91dY6SF8U4Pnnn+fuu+9m9uzZALS1tZG324IkIKncVs2aU6tWyg0m75lZG8ECRjPbD9CYSZ366s+WljzmmtMOjqEm9W/w4MH09vZuWcD4/PPPa8ykQYQlZqxV0kYpX7nB5Ergl8DeZnYnMB44O6pKycBNmvU4720Kv7JuhmzA5Zo+fTrHH388a9asYdKkScyfP58f/ehHSVdLSghbjwIoaWMCSgaToDvrObKr4MeRXcpzkbu/HnHdpELTMt0lFyeeNW5kU6RLKYe7c+CBB3LffffxxBNP4O7cdNNN7LbbbklXLTWiaiGUWo+ipI3xKxlM3N3NLOPufws8GEOdZIDuLGMasALJVmbGxIkTWbx4MZ/61KeSrk7qRLmafSDrQrRmJFrlzuZ6wsw+EmlNpCqZrp7iGTkDzZANuFLjxo1j4cKFSVcjlaJczR6213oz78OepHKDycfJBpTnzexpM+s2s9KLGCQ2pfYo2XGHQWqVFPDoo48ybtw49ttvPw4++GDGjBnDwQdrckItRLmqPGxdSKHHjGzLqF73jk+DcgfgT4i0FlKVTFdPyT1Krv0/+gNZyMMPP5x0FVIrylXl+Snpe9b30mK2pdUzZcIBXHPamC2PGVv30dBgfHRCWyZm9j4zu5js6vfjgR53fyn3E0sNJVSmq4dL7gmfCqzZW9v7y1/+wo033sjMmTP55S9/SUdHB/vss8+WH6letavKM109jJ/xCKOnPliwRTHxsI4t77EpWBuUHyzmTz2Gjva27bp/40gc2YxKdXPdAXQC3WRbJ98p94XN7DYze83MnskrG2Zm88xsRfB7aFBuZnazma0MutEOz3vO5OD4FWambfHyTMt0c/HdS7Z8kQrR7K3CJk+ezKJFixgzZgwPP/wwl1xySdJVSp1qVpXnBu971vfibA0S/QNKqXEZJXCMT6luroPcfQyAmf0QeLKC1/4R8F3gx3llU4HfuPsMM5sa3L+MbKDaP/gZC9wCjDWzYWTXuHSSbakuNrM57v5mBfVIpXJWube3tSqQFPHHP/6R7u7sFey5557LEUcckXCN0mmgq8rL3ZO9VLBQAsf4lGqZbOmId/eNlbywuz8G9F/0cArZ1g7B74l55T/2rCeAdjMbAUwA5rn7G0EAmUe2u63pXXF/d+jjba0tTbtHSTlaW1u33N5hh3KHDiUu5bYoSs3cUgLH+JT6Fh1iZm8Htw1oC+4b2SUou1b4fnu4+ytkn/yKme0elHcAa/KOWxuUFStvWpmuHqbcu4QSORyVpK6EpUuXsuuu2dPX3ent7WXXXXfF3TEz3n777RKvIFEqt0UxZcIB26xlgW2DRdR7x8tWpbbtbQl7vIYKJY30kPLtX8DsPOA8gJEj07meItPVw8V3Lyl53NCdWvVlKWHTpk2lD5LElAoSOeUECyVwjEfc7ftXzWxE0CoZAbwWlK8F9s47bi/g5aD86H7lvy30wu5+K3ArQGdnZyrTvpY7A+XKk9S9JY2tkhaFgkV9iDuYzAEmAzOC3z/PK7/QzO4iOwD/VhBw5gL/LzfrCzgOuDzmOteNQs3+/sbvN0xfLEmFsCChrMD1J7JgYmazybYqdjOztWRnZc0A7jGzc4HVwOnB4Q8BJwIrgQ3AOQDu/oaZfQvI5bv4pruHZzJMqXJW7WoasDSDKHN+ycBFFkzc/cwiD32iwLEOXFDkdW4Dbqth1RpSOV1cCiRSK9Ve+UfZcih32rDEq9zcXJKgaZnukl1cg7TvZWy+8IUvQHamYyoX5Ja7YLCS5//T3UsYVWQle6W0ELE+KZjUuWmZ7pKLEwE+OzadM9jq0dlnnw2wol9xbkHu/sBvgvuw7YLc88guyCVvQe5Y4AjgyryxwURVk+03l96n//P758aqJqAUW1vSvlNrwXKJh4JJHSs3kIzfb5i6uGJ01FFHAfRfxJuaBbkDvfLPtUjC0vtA9bmxpkw4gNaW7Zvif/7LRmUETpCCSZ2aNOvxstKl3HjGodz5pY/GVCsJsc2CXKBmC3LN7DwzW2Rmi9atW1fzivc30P1ACrVoiqmmS2riYR0MGbz9cG/fZlcCxwQpmNShTFdPye13W8xYcuVxGnCsf1UvyHX3W9290907hw8fXtPKFTLQFCSVBIhqc2O91Vt4ywWNmyRHSYnqUDlXV2eO3bvkMRKryBbkxm2gKUiKpUDpr5zAVGo2mBI41h8FkzpUztWVxkjqTqoW5BZaMFjqD3yhFCj9tbe1ctXJHw4NTOWsIyk33YrER8GkzmS6ehhkFjqIOX6/YTHWSPo788wzAQ4kO/O3KRbklvMHvv/uh4UM2XGHki2cctaRKIFj/TEvMfOiEXV2dvqiRYuSrkbFpmW6ufOJ1YU7zvOsmvGpWOojxZnZYnfvjPt9kzq3x894pGCA6GhvY/7UY7YrHz31wYLnsQEvljh/B/JcpVepjWrOaw3A14ncZldhgWQQcOMZh8ZVJZEtKp0uPNAZYQN5brWLLKU2FEzqxJR7w1PLd7S3cf0Zh+pqSxJR6R/4SmeE5e/3/u5fN263jiTsucW6xS65Z2nR/eOl9hRM6sC0THfoZle5rgQFEklKpcGhkv3f+7cs1vf2gWf35Sln7/hiraNN7mqpxEgD8AkrZy93zVCRpA1kwLvcfUYKtSz6Njs7Dd6Brm8cV/L55UxJHkgiSI3DVEbBJGFXzVlW8hidwFIPotqEqtrEjeVMSa7k9UBp7gdC3VwJW19kJW9OW6v+iyQd8sdF8scxqhmsh+271FqscArtShY0VpPsslmpZZKgcvpwrznt4BhqIhKtsCv9WixAzG819X+vgbye0txXTsEkIbk1JcUYcINmb0kDKGdsIexKP7dOpVbjE7VY0Kh0LZVTMElApqsndHFia4sx8+8PUSCRulfu2EKpK/1aj8dU+3pK11I5dcjHLNPVwz/dvSR0caICiVSi2FhEHModW6h2XCRulUwYt2vqAAAOiUlEQVRtliy1TGKU6ephyr1LQwNJi5lOWClb0rOOyh1baMQr/ahmr6WVWiYxmjl3OX2bwzNvldqlTiRf0rOOym1x6Eo//dQyiVE5M0E66rTZL/Up6VlHlbQ4dKWfbgomMSkntXy9N/ul/sQ166jYjC2lgpccBZMY5Pq1wwJJOZsGifQXx1hEqXEZtTgEFExiMf2BZaGpHobu1FpWDiKR/uJoGZSzWZWIgkmEMl09TH9gGW9uCE+Zsr7E4yJhom4ZJD0uI41BwSQima4eLrl3KZtKzN6C+p1rLwKlx2WUXVdAwSQyV9zfXVYgMZRiXupb2LhMsfGURS+9waPPrYskwCh41ScFk4i8+154OmzIBpJJ40bqiyB1LWxcZvyMRwqOp+SnC6rlQsqkF2lKcQomCTGDGz6tRI7SGIqNyxQbN+nfJq/VgL0mA9SvRFbAm9kqM+s2syVmtigoG2Zm88xsRfB7aFBuZnazma00s6fN7PAk6lyJUrmRWltMgURSoZLxvloM2GsyQP1KMp3Kx939UHfvDO5PBX7j7vsDvwnuA5wA7B/8nAfcEntNyzQt082oqQ9y8d1Lih4zZHCLEjlKahTaG76YWkw0abSEkc2knnJznQLcEdy+A5iYV/5jz3oCaDezEUlUMMy0THfJvdzPGjeSZd88XoFEUqN/zq32tlYGFdjosLXFajLRpFDwUuaI+pDUmIkDvzIzB37g7rcCe7j7KwDu/oqZ7R4c2wGsyXvu2qDslfwXNLPzyLZcGDlyZMTV317YRleQHWz/9sQx8VRGZIAGMlMqfzxl/IxHCm5FPWTwDjW5iFL6lvqVVDAZ7+4vBwFjnpk9F3JsoQ2dt5tzGwSkWwE6OztjTb07adbjoWnlQc1wqX/VzJTKBaFC61EA3ioQYAZK6VvqUyLBxN1fDn6/Zmb3A0cAr5rZiKBVMgJ4LTh8LbB33tP3Al6OtcIhJs16nPnPv1HyODXDpZ4UaoEMdKZUoT3X+9PFVPrFPmZiZkPMbJfcbeA44BlgDjA5OGwy8PPg9hzg88GsrnHAW7nusKRlunrKCiRtrYN0JSV1I/fHv2d9L87WFkixVkWpmVKFglA+jWk0hyRaJnsA95tZ7v1/6u6/NLOFwD1mdi6wGjg9OP4h4ERgJbABOCf+Km8v94UsZZDBNacdHEONpB6Y2SrgHWATsNHdO81sGHA3MApYBXza3d+07JfgJrLn9wbgbHd/Kuo6FmuBtBTZIqFUqyIs2HRoTKNpxB5M3P0F4JAC5f8NfKJAuQMXxFC1ipS6GoPswsTrtZ6kGX3c3V/Pu5+b9j7DzKYG9y9j22nvY8lOex8bdeWK/fHf5E5ba0vF6eyL5e7qaG9j/tRjqqusNIx6mhrcMDJdPUW7BPJpYaIE6mrae7GWRm4r3Uq31q1kum6mq4fxMx5h9NQHGT/jkZILfKVxKJ1Khcrt3rrxDAWSJlXzae+1Fpa4cSAzpcqdrltottjFdy9h+gPLuPIkbQzX6BRMKlRO91ZHe5u+GM2r5tPea72GKoq1GuUEoWLfnTc39ClZYwyizrasYFKhcnIAaeZK84pi2nsUa6iSWKsR9t1RssZoxZFtWWMmFSo1s0XTgJtXmqa9R6GaWWFSnbA1RLWiYFKhsMR2rYNM04Cb2x7A781sKfAk8KC7/xKYARxrZiuAY4P7kJ32/gLZae+zgH+Iv8q1FTbAXioppBY2RieObMvq5qpQfn9zz/reLXPzNZ9e0jLtfaBKdaXkvhtXzVm2Xf4uLWyMVqmtl2tBwWQAlBtIZHvlpGPJfXeaYevdevqMYTP4akXBRERqopKulLRfkNXb9sJxZFtWMMlTT1cSIo0mjq6URlGP2wtHHcA1AB8olvxOK3RFyqONq7Zqxu2FFUwCcUydE0mz/rsulpuOJY2acXthdXMFmvFKQqTW0j4WUq44BrzrjVomgWJXDIPM1NUlIhVpxlaaWiaBQlcSkE3LrbxBIlKpZmulqWUSyF1JtNj2ufc0diIiEk7BJM/EwzrYXGCnOdDYiYhIGHVzBXJrTIqlZE3zLAwRkWopmLD9atX+0j4LQ0SkWgomhG94pQSOIiKlKZhA6H7u86ceE2NNREQaU9MPwGsNiYhI9Zo6mOTGSkREpDpN3c0VNlYCFFxzIiL1QVm+60tTB5NSa0fOHLt3TDURkUrU234h0sTdXJmuHgaFtDzOGjeSb08cE2ONRKRcyvJdf5qmZZLfJH5/WyvvvreRTQVWu7e1tqQ+IZs0j7R2BSnLd/1pimDSv0m8vrev4HEtZgokkhpp7grSro71pym6uUoNtOdsdm/4L5lITpq7grSrY/1pmGBiZseb2XIzW2lmUyt5brlNX13VSJqkuSuoGfcLqXcN0c1lZi3A94BjgbXAQjOb4+5/LOf5xZrE+XRVI2mT9q6gZtsvpN41SsvkCGClu7/g7u8BdwGnlPvkQk3i1kHG0J1adVUjqaWuIIlTQ7RMgA5gTd79tcDY/APM7DzgPICRI0du8+RckEjjrBaRYnTeS5waJZgUWhCyzbxed78VuBWgs7Nzuzm/ahJLM9J5L3FplG6utUD+cvS9gJcTqouIiPTTKMFkIbC/mY02s8HAZ4A5CddJREQCDRFM3H0jcCEwF3gWuMfdlyVbK5HqVTPlXaSeNMqYCe7+EPBQ0vUQqZVqp7yL1JOGaJmIpFRVU95F6knDtEwqsXjx4tfN7KUiD+8GvB5nfWKgzxS/fWrwGiWnvMO2096BP5tZFPlQkv73TvL99dm3GvB5ncpg4u7Diz1mZovcvTPO+kRNn6lhlZzyDttOe4+sIgn/eyf5/vrstXlvdXOJJEdT3iU1FExEkqMp75IaqezmKiHS7oKE6DM1IHffaGa5Ke8twG0JTnlP+t87yffXZ68B8wK7DYqIiFRC3VwiIlI1BRMREalaUwWTRkpdYWa3mdlrZvZMXtkwM5tnZiuC30ODcjOzm4PP9bSZHZ73nMnB8SvMbHISnyWox95m9qiZPWtmy8zsokb/TI3GzFaZWbeZLTGzRUFZxf/+Zb5Xoudvkfe/ysx6gs+/xMxOzHvs8uD9l5vZhLzyiv9mJHmuh7x39J/d3Zvih+wA5/PAvsBgYClwUNL1CqnvUcDhwDN5ZdcBU4PbU4Frg9snAg+TXbcwDlgQlA8DXgh+Dw1uD03o84wADg9u7wL8F3BQI3+mRvsBVgG79Sur6N+/Uc7fIu9/FXBpgWMPCv4e7AiMDv5OtAz0b0aS53rIe0f+2ZupZdJQqSvc/THgjX7FpwB3BLfvACbmlf/Ys54A2s1sBDABmOfub7j7m8A84Pjoa789d3/F3Z8Kbr9DNmFnBw38mVKi0n//siR9/hZ5/2JOAe5y97+6+4vASrJ/Lwb0NyPJcz3kvSP/7M0UTAqlrmi0XYP2cPdXIHvSALsH5cU+W11+ZjMbBRwGLCAln6lBOPArM1ts2RQtUPm/fzXq4f/6wqAr6bZcN1OU75/kud7vvSHiz95MwaSs1BUNqthnq7vPbGY7A/8BXOzub4cdWqCsLj9TAxnv7ocDJwAXmNlRIcfG+e8c1//1LcB+wKHAK8B3onz/JM/1Au8d+WdvpmCShtQVr+a6GoLfrwXlxT5bXX1mM2sle4Lf6e73BcUN/Zkaibu/HPx+DbifbFdGpf/+1Uj0/9rdX3X3Te6+GZhF9vNH8v5JnuuF3juOz95MwSQNqSvmALkZHZOBn+eVfz6YFTIOeCtoRs8FjjOzoUGz9rigLHZmZsAPgWfd/fq8hxr2MzUSMxtiZrvkbpP9d3uGyv/9q5Ho/3W/MZ9TyX7+3Pt/xsx2NLPRwP7Akwzwb0aS53qx947ls5eamZCmH7KzJv6L7CyFK5KuT4m6zibbHO0je5VwLvAB4DfAiuD3sOBYI7vJ0vNAN9CZ9zpfIDuothI4J8HPcyTZZvLTwJLg58RG/kyN9EN2Vs7S4GdZ7vwfyL9/I5y/Rd7/34PXfzr4wzgi7/grgvdfDpyQV17x34wkz/WQ9478syudioiIVK2ZurlERCQiCiYiIlI1BRMREamagomIiFRNwURERKqmYFJnzGwPM/upmb0QpL143MxOjbkOq8xsNzNrN7N/GOBrnG1me/Z/zdrVUhqNzu10UzCpI8GCowzwmLvv6+5/S3ax0F79jotru+V2oOAXzsxaSjz3bGDPEsdIk9C5nX4KJvXlGOA9d/9+rsDdX3L3fw2uhu41swfIJuszM5tpZs9Ydo+KMwDM7Ggz+0Xu+Wb2XTM7O7i9ysymm9lTwXMODMo/YGa/MrMuM/sBW/PyzAD2s+z+BzOD137UzH4KdJvZKNt2v4hLLbtvwt8DncCdwXPbgkO+0v+9pWno3E45BZP68mHgqZDHPwpMdvdjgNPIJm07BPgkMNPKSxH+umeT/d0CXBqUXQn83t0PI7s6dmRQPhV43t0PdfcpQdkRZFfDHlTsDdz9Z8AiYFLw3N6Q95bmoHM75RRM6piZfc/MlprZwqBonrvn9mg4Epjt2eRtrwL/CXykjJfNJZ1bDIwKbh8F/ATA3R8E3gx5/pOe3fdgIAq9tzQhndvpo2BSX5aR3R0OAHe/APgEMDwoejfv2EIpogE2su3/6/v6Pf7X4PcmIL9/uty8Ovl1KPVe/RV7b0k/ndspp2BSXx4B3mdm5+eV7VTk2MeAM8ysxcyGk70CexJ4CTjIsllA30/2C1vKY8AkADM7gewWoQDvkN36s5hXgd2Dfukdgb/Le6zUc6W56NxOuaaMoPXK3d3MJgI3mNlXgXVkr5YuA9r6HX4/2X7mpWSvvL7q7n8CMLN7yGYHXQF0lfHW04HZZvYU2S6F1UF9/tvM5gcDkQ8DD/arb5+ZfZPsTm4vAs/lPfwj4Ptm1hvUU5qYzu30U9ZgERGpmrq5RESkagomIiJSNQUTERGpmoKJiIhUTcFERESqpmAiIiJVUzAREZGq/X/ig4ehyH60bgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"def show():\n",
" fig = plt.figure()\n",
" fig.add_subplot(121)\n",
" plt.scatter(y_train, t[X_train].torch().detach().numpy())\n",
" plt.xlabel('Groundtruth')\n",
" plt.ylabel('Predicted')\n",
" plt.title('Train')\n",
" fig.add_subplot(122)\n",
" plt.scatter(y_test, t[X_test].torch().detach().numpy())\n",
" plt.xlabel('Groundtruth')\n",
" plt.ylabel('Predicted')\n",
" plt.title('Test')\n",
" plt.show()\n",
" \n",
"show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Seeking a Low-degree Polynomial Solution\n",
"\n",
"Next we will try a PCE-like solution. The good news is that PCE is essentially a [Tucker decomposition](https://epubs.siam.org/doi/pdf/10.1137/07070111X) with certain custom factors, namely polynomial, and we can emulate this in *tntorch* via the TT-Tucker model. Essentially we will fix polynomial bases and impose low TT-rank structure on the learnable coefficients."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5D TT-Tucker tensor:\n",
"\n",
" 32 32 32 32 32\n",
" | | | | |\n",
" 3 3 3 3 3\n",
" (0) (1) (2) (3) (4)\n",
" / \\ / \\ / \\ / \\ / \\\n",
"1 2 2 2 2 1"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = tn.rand(shape=[ticks]*N, ranks_tt=2, ranks_tucker=3, requires_grad=True) # There are both TT-ranks *and* Tucker-ranks\n",
"t"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the shape of that tensor network: cores are no longer $2 \\times 32 \\times 2$, but $2 \\times 3 \\times 2$. Their middle dimension is now *compressed* using a so-called *factor matrix* of size $32 \\times 3$. In other words, we are expressing each slice of a TT model as a linear combination of three \"master\" slices only.\n",
"\n",
"Furthermore, we want those factors to be fixed bases and contain nicely chosen 1D polynomials along their columns:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"t.set_factors('legendre', requires_grad=False) # We set the factors to Legendre polynomials, and fix them (won't be changed during optimization) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected, the factors' columns are Legendre polynomials of increasing order:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XdcleX/x/HXxVCcKK5QxJGm4ogARc3MmSNHZUNtaK6cWWo56qtZuVLTypGampVmpmWu3LknqLlXThAFF4qDef3+uI/90NgcuM/hfJ6PBw8497nH+3D0fLjv+xpKa40QQgjxgJPZAYQQQtgWKQxCCCEeIoVBCCHEQ6QwCCGEeIgUBiGEEA+RwiCEEOIhUhiEEEI8RAqDEEKIh0hhEEII8RAXswNkRNGiRXXZsmXNjiGEEHYlODj4qta6WGrr2WVhKFu2LEFBQWbHEEIIu6KUOp+W9eRSkhBCiIdIYRBCCPEQKQxCCCEeIoVBCCHEQ6QwCCGEeIhVCoNSao5SKlwpdTiZ55VS6mul1Gml1EGllF+i5zoppU5ZvjpZI48QQoiMs9YZw/dA8xSebwFUtHz1AKYDKKU8gBFAIFALGKGUKmylTEIIITLAKv0YtNZblFJlU1ilLfCDNuYR3aWUKqSU8gQaAOu01tcBlFLrMArMz9bI9aiRy49w9NKtrNi1EEJkOZ+SBRnRumqWHye77jGUAi4mehxiWZbc8v9QSvVQSgUppYIiIiKyLKgQQji67Or5rJJYplNY/t+FWs8EZgIEBAQkuU5qsqPSCiGEvcuuM4YQoHSix17ApRSWCyGESCz6NuyeATpDfxenS3YVhmXAW5bWSbWBSK11GLAGeE4pVdhy0/k5yzIhhBAPnFgNU2vDn4MhNDjLD2eVS0lKqZ8xbiQXVUqFYLQ0cgXQWn8LrAJaAqeBu8DblueuK6U+A/ZadvXpgxvRQgjh8KLCjWJw5DcoVgW6rgOvgCw/rLVaJXVI5XkN9EnmuTnAHGvkEEKIHEFrODAf1nwEsXeh4Ufw9HvgkitbDm+Xw24LIUSOdf0MLH8Pzm4G7zrQ+mso9kS2RpDCIIQQtiA+DnZOgU1jwDkXPP8l+L8NTtk/cpFDjZUUcTeCgxEHzY4hhBAPu3QAZjWE9SPg8cbQZzfU7PpQUYiJj2FryFZ0DmqVZBMGbh7Ih1s+JDo+2uwoQggBMXdh7cdGUYi6Aq/+AO3nQ8GS/1n15+M/03tDbw5dPZTlsRyqMPR6shehUaH8fCxLRtwQQoi0++cvmF4HdnwDT70BffaAT1tQ/+33e/P+TWYcnEG9UvWoUaxGlkdzqMJQp2Qd6pWqx8yDM7l5/6bZcYQQjujudfi9F/z4Aihn6LwS2nwDeQolu8mMgzO4E3uHgf4DsyWiQxUGgIH+A7kTd4dvD35rdhQhhCPRGg4thik14dAieGYg9NoBZeuluNn5W+dZeHwhL1V8iQqFK2RLVIcrDBUKV+Clii/xy/FfOH/rvNlxhBCO4OZFWPAqLOkKhbyhx2ZoPBxc3VLddHLwZFydXenjm2RXsCzhcIUBoI9vH3I552JS8CSzowghcrKEeGN8o6mBcG47NBsD3dbDY9XStPm+K/tYf2E9Xap1oWieolkc9v85ZGEomqcoXap1YcOFDQRfyfpxR4QQDujKUZjTDP78EMrUgd47oU5vcHJO0+ZaayYETaB4nuJ0qpq9k1s6ZGEAeKvqWxTPU5wJeyeQoBPMjiOEyCli78PGUTDjGaMX80uz4PXFULhMunaz+txqDl09RD+/fuRxyZNFYZPmsIUhj0se+vn14/C1w6w+u9rsOEKInOD8Dvi2Hmz5Aqq9DH32Qo1Xk2yCmpLo+GgmB0+mUuFKtC7fOovCJs9hCwNA6/KtqexRma/2fSWd3oQQGXc/0hjfaG4LiI+GN36Dl2ZAviIZ2t2CYwu4dOcSg2oOwjmNl56syaELg7OTMwMDBnLpziUWHFtgdhwhhD06thym1IJ986BOX+i9Cyo0zvDubty/wayDs3im1DPU9qxtxaBp59CFAaC2Z23qe9Vn1sFZ3Lh/w+w4Qgh7cSsMfnnD+MpXDLptgGajIFe+TO3227+/5U7cHQYGZE9ntqQ4fGEAGOA/wOj09rd0ehNCpCIhAYLmGk1QT62DJp9Aj7+glF+md30u8hyLTiyiXcV2PF7o8UzvL6OsUhiUUs2VUieUUqeVUkOSeH6SUuqA5eukUupmoufiEz23zBp50uvxQo/TrmI7Fp1YxLnIc2ZEEELYg6unYF4rWPEeeNYwei7Xex+cXa2y+8n7JpPLORe9fXtbZX8ZlenCoJRyBqYCLQAfoINSyifxOlrr97XWvlprX+Ab4LdET9978JzWuk1m82RUb9/e0ulNCJG0uBjYMh6mPw1XDkObKdBpORSx3l/1QZeD2HBhA12rd83WzmxJscYZQy3gtNb6jNY6BlgItE1h/Q6AzQ1vWjRPUbpW78rGixsJuhxkdhwhhK0ICYKZz8LGz6FSC6MJqt+b6W6CmpIEncDEoIkUz1ucN33etNp+M8oahaEUcDHR4xDLsv9QSpUBygEbEy12U0oFKaV2KaVeSO4gSqkelvWCIiIirBD7v970eZMSeUswIUg6vQnh8KKj4M8h8F0TuHcT2v8Mr86DAiWsfqg/z/7J4WuH6e/XP9s7syXFGoUhqbKZ3BRD7YHFWuv4RMu8tdYBQEdgslIqyXMzrfVMrXWA1jqgWLFimUucjDwueXjX712OXDvCqrOrsuQYQgg7cGodTKsNu781ZlLrsxsqt8ySQ0XHR/PVvq+o7FGZVuVbZckx0ssahSEEKJ3osRdwKZl12/PIZSSt9SXL9zPAJuApK2TKsFblW1HFowqTgidxN/aumVGEENntzlVY0g3mvwyueaHLGnh+IrgVzLJDzjk8h7A7YQwKGISTso2GotZIsReoqJQqp5TKhfHh/5/WRUqpSkBhYGeiZYWVUrktPxcFngaOWiFThjkpJ4YFDiP8bjgzDs4wM4oQIrtoDQd+NuZKOLIUGgyFnlvBOzBLDxsaFcrsQ7N5rsxzBHpm7bHSI9OFQWsdB/QF1gDHgEVa6yNKqU+VUolbGXUAFuqHZ7KuAgQppf4G/gLGaq1NLQwAvsV9aft4W344+gNnI8+aHUcIkZVunIMfX4SlPaFoRei5DRoMAZfcWX7oL/Z8gZNy4oOaH2T5sdJDPfw5bR8CAgJ0UFDWthy6eu8qbX5vQ7Wi1ZjRdAbKii0QhBA2ID4Odk+Hv0YbU2w2GQEBXcEpey7nbAvdRq/1vejv159u1btlyzGVUsGWe7opso0LWjaoaJ6i9HmqDzvDdrLhwgaz4wghrCnsb/iuMaz9GMo9a9xcrtU924pCTHwMY/eMpUzBMrzl81a2HDM9pDCk4LVKr/FE4Sf4Yu8X3Iu7Z3YcIURmxdyFdcNhZkO4dQle+R46/AzuSbawzzI/HP2B87fOM7TWUHI558rWY6eFFIYUuDi5MCxwGGF3wvju0HdmxxFCZMaZzTC9Lmz/Cnw7Qt89UPVFq3ZUS4vLdy4z8+BMGpVuxNOlns7WY6eVFIZU+Jfw5/nyzzP38Fwu3LpgdhwhRHrdvQ5L+8APbYwi0Gk5tJ0CeQqbEmf83vEk6AQ+rPWhKcdPCykMaTDQfyCuTq6M2zvO7ChCiLTSGg7/BlNrwd8/w9PvGYPelatvWqRdYbtYe34tXat3pVT+7L18lR5SGNKgWN5i9PbtzZaQLWy6uMnsOEKI1ESGwM8dYPHbULAU9NgETUeCq3nDTcTGxzJm9xi88nvRpVoX03KkhRSGNOpYpSOPuz/O2D1jZRpQIWxVQgLsmWXMlXB2Mzw3yphAx7OG2clYcHwBZyLPMLjWYHI7Z30ficyQwpBGrk6uDA0cSmhUKHMOzzE7jhDiUeHHYE4zWDUISteC3juhbl9wdjE7GeF3w5l2YBr1verToHQDs+OkSgpDOgR6BtKsbDNmH5pNaFSo2XGEEABx0fDXGPj2Gbh2Gl6cAW/8BoXLmp3sX18Gf0lsQixDav5nHjObJIUhnR4MdPXFni/MjiKEuLDLKAibxxpNT/vuhSfbZ3sT1JQEXQ5i5ZmVvF3tbUoXLJ36BjZACkM6PZbvMXrU6MHGixvZFrrN7DhCOKb7t2DlQOPSUexdeH0xtJsF+cyd+exRcQlxjN4zGs98ntk27IU1SGHIgE4+nShbsCxj94wlJj7G7DhCOJbjq4yby0FzoHZv6L0LKjY1O1WSfjnxC6dunOLDmh/axAQ8aSWFIQNcnV0ZWmso52+dlxvRQmSX21dg0VuwsIPROa3remg+BnLnNztZksLvhjNl/xTqlqxLY+/GZsdJFykMGVS3VF2al23OzIMzORN5xuw4QuRcWsO+H2BqTTixGhr9D97ZDF7+ZidL0ejdo4lNiGVY4DC7G51ZCkMmDK41mDwueRi5Y6TMES1EVrj2D8xrDcv6QYnqRs/l+oPA2dXsZClaf349Gy5soNeTvShTsIzZcdJNCkMmFM1TlA9rfsi+8H38euJXs+MIkXPEx8LWL2FaHQg7CK2/MsY4KlrB7GSpioyOZNTuUVTxqEKnqp3MjpMhVikMSqnmSqkTSqnTSqn/NNRVSnVWSkUopQ5Yvroleq6TUuqU5cvufottHm9Dbc/aTNo3ict3LpsdRwj7FxpsDIu9YSQ88ZwxCqp/52ybKyGzJgVP4sb9G3xS9xNcnMzvXJcRmf5NK6WcgalAC8AH6KCU8kli1V+01r6Wr+8s23oAI4BAoBYwQillzpCHGaSUYnid4SToBEbtGoU9zognhE2IuQNrPoLvmsDdq/DafHjtJyjwmNnJ0mxP2B6WnFrCW1XfwqdIUh+D9sEaJbgWcFprfUZrHQMsBNqmcdtmwDqt9XWt9Q1gHdDcCpmyVekCpenj24dNIZtYc36N2XGEsD+n18O02rBzinF20Gc3VGlldqp0uR93n5E7R1K6QGl6PdnL7DiZYo3CUAq4mOhxiGXZo9oppQ4qpRYrpR50/0vrtiileiilgpRSQREREVaIbV2vV3mdqkWqMmb3GCKjI82OI4R9uHMNfnsHfmoHLm7w9p/QahK4uZudLN2m/z2dC7cvMKLOCLvqs5AUaxSGpNphPXo9ZTlQVmtdA1gPzEvHtsZCrWdqrQO01gHFihXLcNis4uLkwsi6I7kVfYsJQRPMjiOEbdMaDi4ymqAeXgL1P4R3tkKZumYny5Bj144x78g8Xqr4EoGegWbHyTRrFIYQIPEAIF7ApcQraK2vaa0fjFU9C/BP67b2pJJHJd6u9jZLTy9l56WdZscRwjbdOG+cIfzWHTzKwztboNFH4OpmdrIMiUuIY8SOERR2K8wA/wFmx7EKaxSGvUBFpVQ5pVQuoD2wLPEKSinPRA/bAMcsP68BnlNKFbbcdH7OssxuvfPkO5QtWJaRO0dyL+6e2XGEsB0J8bBzqnEv4eJuaDEeuqyBEvZ7kxbgx6M/cuz6MYYFDsM9t/1dAktKpguD1joO6IvxgX4MWKS1PqKU+lQp1cay2rtKqSNKqb+Bd4HOlm2vA59hFJe9wKeWZXYrt3NuRtQZQWhUKNMOTDM7jhC24fJho7XRmmFQ9hljfKPAHuDkbHayTLlw6wJTD0ylUelGNPFuYnYcq1H22LwyICBABwUFmR0jRZ/u/JQlp5awoOUCqhatanYcIcwRew82fwE7vga3QtDyC6j6kk0Ni51RWmu6r+3OkWtHWNp2KSXylTA7UqqUUsFa64DU1rOPHiN26H3/9yniVoQRO0YQmxBrdhwhst/ZrTD9adj2JdR4zZgroVq7HFEUAJaeXsruy7sZEDDALopCekhhyCIFchXgo9ofceLGCeYdmZf6BkLkFPduGGMbzWsFOh7eXAovTIO8HmYns5qr964yPmg8/iX8aVexndlxrE4KQxZq7N2YpmWaMv3AdBmBVeR8WsORpcZcCfvnQ913oddOeLyh2cmsSmvN6N2jiY6L5pM6n+Ckct7HaM57RTZmWOAw8rrm5aOtH8klJZFz3boEC1+HXztB/hLQfSM89xnkymt2MqtbcWYF686vo7dvb8q6lzU7TpaQwpDFiuYpyvA6wzl87TCzDs4yO44Q1pWQAHu/M84S/tkITT+F7n9BSV+zk2WJsKgwRu8ejV9xPzpX7Wx2nCxjn0P/2ZmmZZrS5vE2zDw4k3ql6lGjWA2zIwmReREnYNm7cHEXlHsWWk82OqzlUAk6gY+2f2QMmFlvFM523tQ2JXLGkE2G1BpC8bzFGbZtGHdj75odR4iMi4uBTePg23oQcRzaToO3/sjRRQGMjmx7L+9lSK0heBXwMjtOlpLCkE0K5CrAqHqjuHDrAhODJpodR4iMubgHZtSHTaOhShvoGwRPvZ5jmqAm5+SNk3y17ysalW7ECxVeMDtOlpPCkI1qPlaTTlU7sejkIraEbDE7jhBpF30bVn0As58zfu64CF6eDfltb0BLa4uJj2Ho1qEUyFWAEXVH2N38zRkhhSGb9XuqHxULV2T49uHcuH/D7DhCpO7kGuPm8p5ZEPgO9NkFTzQzO1W2mXJgCidvnOTTup/i4ZZz+mKkRApDNsvlnIsx9cZwK+YWI3eOlBnfhO2KCodf34YFr0LugtB1HbQYB7kLmJ0s2wRdDuL7w9/z8hMv82zpZ82Ok22kMJigkkcl+j3Vjw0XNrDsn2WpbyBEdtIa9v8EU2rC8RXQ8GNjaOzSNc1Olq2iYqL4aNtHeBXw4oOAD8yOk62kuapJ3vJ5i80hmxmzZwwBjwVQKn+SE9cJkb2un4Hl78HZzeBdF1p/BcWeMDuVKcbuGcvlu5eZ13weeV1zXke9lMgZg0mcnZwZVW8UAMO2DiM+Id7kRMKhxcfBtskwrQ5c2m9Mr9l5pcMWhfXn1/PHP3/QrXo3fIvnzM56KZHCYKJS+UsxLHAY+8L38cPRH8yOIxzVpQMwqyGsHwEVmkCf3RDQBZwc8+Ph6r2rjNw5Ep8iPvR8sqfZcUzhmO+8DWldvjVNyzTl6/1fc+L6CbPjCEcScxfWfmwUhagr8OqP0H4+FCxpdjLTaK0Zvn049+LuMabeGFydXM2OZAqrFAalVHOl1Aml1Gml1JAknh+glDqqlDqolNqglCqT6Ll4pdQBy5fD3YlVSvG/2v+jUO5CDNk6hPtx982OJBzBPxuNKTZ3fANPvQl99oBPm9S3y+F+PfkrW0O38r7/+5QvlLN7cqck04VBKeUMTAVaAD5AB6XUo5O47gcCtNY1gMXAF4meu6e19rV8OeS/zMJuhfns6c84ffM04/aOMzuOyMnuXoffe8GPL4KTi3Efoc3XkKeQ2clMd/z6ccbtGUfdknXpULmD2XFMZY0zhlrAaa31Ga11DLAQaJt4Ba31X1rrBwME7QJy9kAjGVCvVD26Ve/G4pOLWf7PcrPjiJxGazi02GiCemgRPDMQeu2AsvXMTmYTbsfcZuCmgRRyK8SYZ8bkyDkW0sMazVVLARcTPQ4BAlNYvyvwZ6LHbkqpICAOGKu1XmqFTHapj28f9ofv57Ndn+FTxIfHCz1udiSRE9y8CCsHwKm1UNIP2vwBj1UzO5XN0FozYscIQqNCmdNsjsP0bk6JNcpiUgOHJNmdVyn1BhAAjE+02NsyOXVHYLJSKslPQ6VUD6VUkFIqKCIiIrOZbZKLkwtf1P+CPC55GLBpgIzCKjInIR52fWsMZ3FuOzQbA93WS1F4xILjC1h3fh39/frjV8LP7Dg2wRqFIQQoneixF3Dp0ZWUUk2Aj4A2WuvoB8u11pcs388Am4CnkjqI1nqm1jpAax1QrFjOHbireN7ijKs/jrORZ/l81+cyZIbImCtHjAHvVg+GMnWN8Y3q9IYcPIdARhyMOMiEoAk08GqQoyfeSS9rFIa9QEWlVDmlVC6gPfBQ6yKl1FPADIyiEJ5oeWGlVG7Lz0WBp4GjVshk12p71qaXby+Wn1nOb6d+MzuOsCex92Hj58bQ2DfOwkvfweu/QiFvs5PZnMjoSAZtHkSJvCX4vN7nDjFqalpl+h6D1jpOKdUXWAM4A3O01keUUp8CQVrrZRiXjvIDv1p++RcsLZCqADOUUgkYRWqs1trhCwPAOzXe4UD4AUbvHk3VolWp7FHZ7EjC1p3bDsv7w7VTUKM9NBsN+YqYncomJegEPtr2ERH3IvixxY+453Y3O5JNUfZ4qSIgIEAHBQWZHSPLXb9/nVeWv4KbsxsLWy2kQC7HGdVSpMP9SFg3AoLnGmcGrSZDhcZmp7Jpcw7PYVLwJIbWGkrHKh3NjpNtlFLBlnu6KXLsNlk2zsPNg/H1xxMaFcqIHSPkfoP4r2PLYUot2DcP6vSF3rukKKQi+EowX+/7mmZlmzl8f4XkSGGwcX4l/Ojv159159ex4PgCs+MIW3ErDH55w/jKVwy6bYBmoyBXPrOT2bRr967x4eYP8SrgxSd1PpH7CsmQYbftQKeqndh3ZR8TgiZQvWh1ahSrYXYkYZaEBOPsYN0IiI+GJp8YZwrOjjmmT3rEJ8QzZOsQImMimdZkGvlz5Tc7ks2SMwY74KSc+Lze5xTPU5xBmwcRGR1pdiRhhqunYF4rWPEeeNYwei7Xe1+KQhrNPDiTXWG7GBY4jEoelcyOY9OkMNgJ99zuTGwwkYh7EQzdOlTmb3AkcTGwZTxMrwtXDkObb6DTcigiPePTakfoDqb/PZ02j7fhxQovmh3H5klhsCPVilZjSM0hbA3dylf7vzI7jsgOIUEw81mjb0KlltBnL/i9BXJtPM3ORJ5h0OZBVChcgY8CP5L7Cmkg9xjszGuVX+PUzVPMPTyX8u7leaHCC2ZHElkhOgo2fga7ZxjzI3RYCJVamJ3K7kRGR9JvQz9cnV2Z0miKw03RmVFSGOzQ4FqDOXfrHCN3jsS7gLeM75LTnFxrDHoXGQI1u0Hj4eBW0OxUdic2IZaBmwYSdieMOc3mUDK/405AlF5yKckOuTq5MvHZiXjl9+K9v94jNCrU7EjCGqIiYHFXWPAKuOaFLmvg+QlSFDJAa82Y3WPYfXk3I+uOdMh5mzNDCoOdcs/tzjeNviFOx9F3Q1/uxN4xO5LIKK3hwM8wtSYc/QMaDIWeW8E7pdHrRUoWHF/Aryd/pWu1rrR+vLXZceyOFAY7Vta9LBOfncjZyLMM3jJYWirZo+tn4ccXYGlPKPoE9NwGDYaAS26zk9mt7aHb+WLvFzQs3ZB3/d41O45dksJg5+qUrMPQWkPZHLKZyfsmmx1HpFV8HGz/GqbVgZBgeH4ivL0aistgiZlx5uYZPtj8ARULVWTsM2Mdfia2jJKbzznAa5Vf4/TN03x/5HvKu5fnxYrSTtumhf0Ny/oZ359oYRQF91Jmp7J7N+/fpO/Gvrg6u/JNo2+kBVImSGHIIR60VPp016d4F/TGv4S/2ZHEo2LuwuaxsGMK5C0Cr3wPPi9InwQriI2PZcDmAVy+c5k5zebgmd/T7Eh2Tc6zcggXJxcmPDsBr/xevP/X+4TcDjE7kkjszGaj5/L2r8C3I/TdA1VflKJgBVprRu0exd7Le6UFkpVIYchB3HO7M6XxFOJ1PP029iMqJsrsSOLudVjaB35oYxSBTsuh7RTIU9jsZDnG/GPzWXJqCd2qd5MWSFZilcKglGqulDqhlDqtlBqSxPO5lVK/WJ7frZQqm+i5oZblJ5RSzayRx5GVKViGLxt8ybnIcwzcPJDY+FizIzkmreHwbzC1Fvz9szHYXa8dUK6+2clylI0XNjI+aDyNSjei31P9zI6TY2S6MCilnIGpQAvAB+iglPJ5ZLWuwA2tdQVgEjDOsq0PxhzRVYHmwDTL/kQmBHoGMrzOcHZc2sGwbcOkGWt2iwyBn9vD4rehYCnosckYHts1j7m5cpi9l/fyweYPqFqkKmOeGSMtkKzIGjefawGntdZnAJRSC4G2QOK5m9sCn1h+XgxMUcZIVm2BhVrraOCsUuq0ZX87rZDLob1Y8UUioyOZGDyRgrkK8nHtj2XwsKyWkABBs2H9J6AT4LlRENgTnKWNh7UduXaEfhv7UbpAaaY1niYtkKzMGv9iSwEXEz0OAR7tsvnvOlrrOKVUJFDEsnzXI9tKuz0r6VytMzejbzL78Gzcc7tLZ5+sFH4Mlr0LIXvg8UbQahIULmt2qhzpTOQZeq3rRaHchZjRdAaF3AqZHSnHsUZhSOrP0EcnJ05unbRsa+xAqR5ADwBvb+/05HNo/f36ExkTyaxDs3DP7U6nqp3MjpSzxEXD1omw9UvIXQBenAk1XpXWRlkkLCqMHmt74KScmNl0JiXylTA7Uo5kjcIQApRO9NgLuJTMOiFKKRfAHbiexm0B0FrPBGYCBAQEJFk8xH8ppfg48GNuRd9iQtAECuYqKB3grOXCLuMs4eoJqPEaNBsN+YqanSrHun7/Oj3W9eBu7F3mNp+Ld0H5AzGrWONuzV6golKqnFIqF8bN5GWPrLMMePCn6svARq21tixvb2m1VA6oCOyxQiaRiLOTM2OfGUvdknX5ZOcnbLiwwexI9u3+LVgxAOY0g9h78PoSeGmmFIUsFBUTRc91Pbl85zJTGk+RqTmzWKYLg9Y6DugLrAGOAYu01keUUp8qpdpYVpsNFLHcXB4ADLFsewRYhHGjejXQR2stTWiygKuzK5MaTKJa0Wp8sPkDdoftNjuSfTq+EqYGQvBcqN0beu+Eik3MTpWj3Y+7T7+N/Th14xQTG0yU+UeygTL+cLcvAQEBOigoyOwYdikyOpLOqztzKeoSs5vNplrRamZHsg+3L8OfHxrDYpeoBm2+hlIy7EhWi0uI4/1N77P54mbGPjOWluVbmh3JrimlgrXWAamtJw1/HYx7bndmNJ1BYbfC9FrfizM3z5gdybZpDcHzYEotOLHamE2txyYpCtkgQScwYscINl3cxLDAYVIUspEUBgdUPG9xZjWdhYuTC93XdZcZ4JJz7R+Y1xqWvwuPVTd6Lj8zEJxdzU6W42mtGb93PMv+WUZf3760r9ze7EgORQqDgypdsDTfNvmWe3H3eHv121y4dcHsSLYjPtZogjqtDoQdhNZfGWMcFa1gdjKHoLVbG6ktAAAgAElEQVRm3N5x/HTsJ96o8gY9avQwO5LDkcLgwCp5VOK7577jXtw9Oq/uzD83/zE7kvlCg2FmA9jwKTzRzBgF1b8zOMl/lewQnxDPyJ0jmX9sPm/6vMmHNT+UHvsmkH/tDs6niA9zm81Fo3l79dscv37c7EjmiLkDq4fBd03g7jV4bT689iMUeMzsZA4jLiGOYduGseTUEnrU6MEHAR9IUTCJFAZBhcIV+L759+R2yU2XNV04GHHQ7EjZ6/R6mFYbdk01zg767IYqrcxO5VBi42P5YPMHrDq7iv5+/en3VD8pCiaSwiAAY7juec3n4Z7Lne5ruxN02QGaA9+5Br/1gJ/agYsbvP2nMcaRm7vZyRzK/bj7vPvXu6y/sJ7BNQfTrXo3syM5PCkM4l8l85fk++bfUyJfCXqt78WO0B1mR8oaWsPBRTC1pjFnwrODoec2KFPX7GQO527sXfps6MP20O2MqDOCN3zeMDuSQAqDeESJfCWY22wuZQqWoe/Gvvx14S+zI1nXjfPGGcJv3cGjPPTcCg2HgUtus5M5nNsxt3ln3TsEXwlmVL1RvPzEy2ZHEhZSGMR/FMlThNnNZlOpcCUGbBrA6rOrzY6UeQnxsHOacS/h4m5oMR66rIHiVcxO5pBu3r9Jt7XdOHztMOOfHS9TctoYmUFEJMk9tzuznptFnw19GLx1MNHx0bSt0NbsWBlz+ZAxCuqlfVCxGbT6Ety9zE7lsK7eu0r3td25cOsCXzX8ivpeMt2prZEzBpGs/LnyM73JdAIfC+Tj7R8z5/Ac7Gpsrdh7sH6k0S8h8iK8PAc6/iJFwURnI8/S6c9OhEaFMrXJVCkKNkoKg0hRXte8fNP4G5qVbcak4EkM3zGc2PhYs2Ol7uxWmP40bPsSarSHPnugWjuZQMdEOy/t5PVVrxMVG8XMpjOp7Vnb7EgiGXIpSaQqt3Nuvqj/BeXcy/Ht399y4dYFJjecTGG3wmZH+697N2DdcNj3gzG15lt/QPkGJocSi04sYvTu0ZRzL8eUxlMolV9m8LVlcsYg0sRJOdHHtw/jnhnH4auH6biyo20NoaE1HFlqjIK6fz483R967ZSiYLK4hDjG7hnLZ7s+o27JuvzY4kcpCnZACoNIl5blWzK3+Vzuxd3jjVVvsD10u9mRIDIUFnaEXzsZQ1h03whNP4Vcec1O5tBux9ym78a+/4579E2jb8ifK7/ZsUQaZKowKKU8lFLrlFKnLN//c21BKeWrlNqplDqilDqolHot0XPfK6XOKqUOWL58M5NHZI8axWrw8/M/Uyp/KXpv6M2CYwvMCZKQAHu/M2ZU++cvaPoZdP8LSso/I7NdvH2RN1e9ye5LuxlRZwQf1vwQZydns2OJNMrsGcMQYIPWuiKwwfL4UXeBt7TWVYHmwGSlVKFEz3+gtfa1fB3IZB6RTTzze/JDix+o71WfMXvG8Pmuz4lNyMab0hEnYG4LWDkQSvlB7x3w9LvgLLfNzBZ8JZjXV75OxL0IZjSdIR3X7FBmC0NbYJ7l53nAC4+uoLU+qbU+Zfn5EhAOFMvkcYUNyOual8kNJvN21bf55cQv9F7fm8joyKw9aFw0bBoL39aDiOPwwnTjBrNH+aw9rkiTpaeX0m1tN9xzu7Pg+QXU8qxldiSRAZktDCW01mEAlu/FU1pZKVULyAUkvms5ynKJaZJSSsYlsDPOTs4MCBjAp3U/JehKEG+seoOTN05mzcEu7IYZ9WHTGKjSBvoGgW9HaYJqA2LjYxm/dzz/2/4//Ev481PLnyhTsIzZsUQGpVoYlFLrlVKHk/hKVzdYpZQn8CPwttY6wbJ4KFAZqAl4AINT2L6HUipIKRUUERGRnkOLbPBixRf57rnvuB1zmw4rOrDg2ALrdYa7fwtWDoI5zSA6Cjr+Ci/Phvxy4mkLzkWe4/VVr/PD0R/oULkD05tMxz23jFBrz1Rm/vMqpU4ADbTWYZYP/k1a60pJrFcQ2ASM0Vr/msy+GgCDtNapDoQfEBCgg4IcYFhoO3Tt3jWG7xjOlpAt1Peqz2dPf4aHm0fGd3jiT+M+wq1LEPgONPoYchewXmCRYVprlp5eypg9Y8jlnItP635KI+9GZscSKVBKBWutA1JbL7OXkpYBnSw/dwL+SCJILuB34IdHi4KlmKCMGTleAA5nMo8wWZE8RZjSaApDaw1l16VdtFvWLmPDd0eFw6+d4ef2xvwI3dZDi3FSFGxEZHQkgzYPYviO4dQoWoMlrZdIUchBMnvGUARYBHgDF4BXtNbXlVIBQE+tdTel1BvAXOBIok07a60PKKU2YtyIVsAByzZRqR1Xzhjsw8kbJxm8ZTCnb56mk08n3vV7l1zOuVLeSGvY/xOs/Rhi70L9D43Oai6pbCeyTfCVYIZsHcLVu1fp59ePzlU746SkS5Q9SOsZQ6YKg1mkMNiP+3H3mRA0gV9O/EIVjyqMqz+Ocu7lkl752j+w4j04uwW860Lrr6DYE9kbWCQrLiGOb//+llmHZuGV34tx9cdRrWg1s2OJdMiuS0lCpMjNxY2Pa3/M1w2/JuxOGK+teI0lJ5c8fGM6Pha2TYLpdeHSAWN6zc4rpSjYkJDbIXRe3ZkZB2fQunxrFrVeJEUhB5MzBpFtwu+GM2zbMHaH7aaxd2OG1hpKichLsKyfMWdC5VbQcgIU9DQ7qrBI0AksPb2U8XvHAzC8znBalGthciqRUXIpSdikBJ3AvCPzmHpgKs4JcfS8do034vPi2nI8+LQxO55I5Mi1I4zeNZqDVw/iX8KfUfVGyQB4di6thUHGDxDZykk58XaeMjS5Hs0XLnf50qMQvxcow7DCxZHR+W1DZHQkX+/7ml9P/oqHmwej642mVflWKOlI6DCkMIjsc/c6rBkGf/9MaY/H+ablPDY7xzN2z1i6r+1Os7LNGBQwiMfyPWZ2UoeUoBP4/dTvTN43mdsxt3m9yuv09u1NgVzSRNjRSGEQWU9rOLQYVg+B+zfhmUFQ/wNwdeNZoHbJ2sw5PIfZh2azJWQLPZ/syZtV3sTV2dXs5A7jyNUjjNo9ikNXD+FX3I9hgcOo5PGfvqrCQcg9BpG1bl6AFQPg9Doo5Q+tv4bHkm7NEnI7hHF7x7Hp4ibKuZdjaK2h1ClZJ5sDO5ab92/y9f6vWXxyMUXyFGGA/wC5bJSDyc1nYa6EeNgzEzZ8Zjxu/D+o1QPSMCb/lpAtjNk9hpCoEBqVbkSPJ3tQtUjVLA7sWKJiolh0chFzD8/ldsxtOlbpSO8ne8tEOjmcFAZhnitHjCaoocFQoSm0+hIKeadrF9Hx0Xx/+HvmHZ3H7ZjbPF3yabrX6I5/Cf8sCu0Ybt6/yfzj85l/bD63Y25Tt2RdBgYM5InC0mfEEUhhENkv9j5sGQ/bJxvjGzUfB9VfztSw2FExUSw8sZAfj/7I9fvX8SvuR48aPahbsq5c7kiHiLsR/HD0B3458Qv34u7R2Lsx3ap3k05qDkYKg8he57bB8v5w7TQ82QGeGwX5ilht9/fi7vHbqd+Ye3guV+5ewaeID92rd6eRdyMZpycFoVGhzD08l99P/U6cjqNFuRZ0q9aNCoUrmB1NmEAKg8ge927C+hEQ/D0UKmMMZ1GhcZYdLjY+luVnljP70Gwu3L5AeffydKvejeblmuPqJK2YHjgTeYbZh2az8sxKnJQTbSu0pUvVLpQuWNrsaMJEUhhE1ju6DFZ9AHfCoXZvaDgMcuXLlkPHJ8Sz9vxaZh2axakbp/Bw86BluZa0Kt8KnyI+DnmZKTI6kjXn1rDyzEr2he/DzdmNl594mc5VO1MiXwmz46UqNjaWkJAQ7t+/b3YUu+fm5oaXlxeurg//sSSFQWSdW2GwahAcXwGPVTeaoJbyMyVKgk5gW+g2lp5eyqaLm4hNiKVswbK0Kt+KluVbUrpAzv4LOTo+mi0hW1jxzwq2hG4hLiGO8u7laVW+Fe2eaJe5SZKy2dmzZylQoABFihRxyMJuLVprrl27xu3btylX7uGRjKUwCOtLSIB938O6ERAfAw2GQJ2+YCMd0W7F3GLduXWsOLOCoCvGvw/fYr60Kt+KZmWbUcitkMkJrSNBJxB8JZiVZ1ay9txabsfepmieov+eMVX2qGyXH6zHjh2jcmX7zG5rtNYcP36cKlWqPLRcCoOwroiTxs3lCzug7DPGXAlFHjc7VbLCosJYeXYlK/5ZwT+R/+Di5EK9UvV4ptQz+Jfwp7x7ebv6AIqMjuRA+AH2Xt7LmvNruHznMnld8tKkTBOeL/88gY8F4pyGPiK27NixY//5IBMZl9TvM1sG0VNKeQC/AGWBc8CrWusbSawXDxyyPLygtW5jWV4OWAh4APuAN7XWMZnJJKwsLga2fwVbvgDXPNBmCjz1RqaaoGYHz/yedKveja7VunLixglW/LOCP8/9yaaLmwAolLsQfsX98CvhR0CJACp5VMLFyXZGiIm4G0FweDDBl4PZF76PUzdOodG4OLlQ27M27/u9T4PSDcjrmtfsqCIHyuzUnl8A17XWY5VSQ4DCWuvBSawXpbX+T5dKpdQi4Det9UKl1LfA31rr6akdV84YssnFvbD8XQg/ClVfNPolFLD9m5jJ0Vpz8fZFgq8EE3zF+MC9ePsiAHld8vJksSfxL+HPU8WfwrugN8XyFMuWv8KjYqK4dOcSR64eYV/4PoKvBP+bK49LHnyL+eJXwg//Ev5UL1odNxe3LM9kBjljsC7TzhiAtkADy8/zgE3AfwpDUpRxHt8I6Jho+0+AVAuDyGLRt2Hj57B7BhQsCR0WQiX7n5xFKYV3QW+8C3rzYsUXAWPyoH1X9hF0JYh94fuYcmDKv+u7KBdK5CuBZz5P4yu/JyXzlfz356J5iuKsUi8cd2LvcOnOJcKiwgi7E8alqEuE3TF+DosK43bs7X/Xdc/tjl9xP16r9Br+Jfyp5FFJmuFmM2dnZ6pXr05sbCwuLi506tSJ9957Dycn8/vLrF69mv79+xMfH0+3bt0YMmRIlhwns4WhhNY6DEBrHaaUKp7Mem5KqSAgDhirtV4KFAFuaq3jLOuEADILiNlOrjEGvbsVCjW7QePh4FbQ7FRZpnje4jQv15zm5ZoDxrX8w1cPExoV+u+H+OU7l9l7ZS/hZ8NJ0AmZPmYB1wL/Fhm/4n6UzF8Sz/yeVHCvQPlC5aXDnsny5MnDgQMHAAgPD6djx45ERkYycuTITO1Xa43WOsMFJj4+nj59+rBu3Tq8vLyoWbMmbdq0wcfHJ1O5kpJqYVBKrQeSGiD/o3Qcx1trfUkpVR7YqJQ6BNxKYr1kr2sppXoAPQC8vdM37o5Ig6gIWD0YDi+BYpWhyxrwDjQ7VbZzz+3O06WeTvK52IRYIu5G/PsX//X710nLpVg3F7d/zzI883nK/AZ2pHjx4sycOZOaNWvyySefMH/+fL7++mtiYmIIDAxk2rRpODs789lnnzF//nxKly5N0aJF8ff3Z9CgQZw7d44WLVrQsGFDdu7cydKlS9m6dWuS+/jpp5+SXP7Anj17qFChAuXLlwegffv2/PHHH+YUBq11k+SeU0pdUUp5Ws4WPIHwZPZxyfL9jFJqE/AUsAQopJRysZw1eAGXUsgxE5gJxj2G1HKLNNIa/v7ZmEAnOgoaDIN674FLbrOT2RxXJ1dK5i9JyfwlzY6S441cfoSjl5L62zHjfEoWZETr9I/SW758eRISEtiyZQu//PIL27dvx9XVld69ezN//nx8fHxYsmQJ+/fvJy4uDj8/P/z9/3+wxxMnTjB37lymTZvGsWPHktxHzZo1k1z+1ltv/buf0NBQSpf+/345Xl5e7N69O3O/lGRk9lLSMqATMNby/Y9HV1BKFQbuaq2jlVJFgaeBL7TWWin1F/AyRsukJLcXWej6GVjxPpzZBKUDjY5qxSubnUoIm6O1ZtOmTQQHB1OzZk0A7t27R/Hixbl+/Tpt27YlT548ALRu3fqhbcuUKUPt2sbEtRs2bEhyH7du3Upy+aMZHpVVTa4zWxjGAouUUl2BC8ArAEqpAKCn1robUAWYoZRKAJww7jEctWw/GFiolPoc2A/MzmQekRbxcbBrGvw1GpxcoOUECOgKNnBzTQggQ3/ZZ5UzZ87g7OyMh4cHnTp1YsyYMQ89P2nSpBS3z5fv/4eJ0VonuY9vvvkmyeWJeXl5cfHixX8fh4SEULJk1py9ZuqTQGt9TWvdWGtd0fL9umV5kKUooLXeobWurrV+0vJ9dqLtz2ita2mtK2itX9FaR2fu5YhUXToA3zWCdf+DxxtCn91Qq7sUBSGSEBERQc+ePenbty+NGzdm8eLFhIcbV8yvX7/O+fPnqVevHsuXL+f+/ftERUWxcuXKZPeX3D6SW55YzZo1OXXqFGfPniUmJoaFCxfSpk2bLHndttOjR2StmLuwaQzsnAr5isIr88Cnrc13VBMiu927dw9fX99/m6u++eabDBgwACcnJz7//HOee+45EhIScHV1ZerUqdSuXZs2bdrw5JNPUqZMGQICAnB3d09y3z4+PsnuI6nlZcqU+XdbFxcXpkyZQrNmzYiPj6dLly5UrZo1Z1YyJIYjOLPJGM7ixjl46k147jPIU9jsVEI8xJ47uEVFRZE/f37u3r1L/fr1mTlzJn5+5gws+YCZHdyELbt7Hdb+Dw78BB6PQ6cVUO4Zs1MJkeP06NGDo0ePcv/+fTp16mR6UcgsKQw5kdZw5Df4c7BRHOoNgGc/NMY6EkJY3YIFC8yOYFVSGHKamxdh5UA4tQZKPgVv/m7MmSCEEGkkhSGnSIiHvbNhw0jQCdBsNAT2BDsfilkIkf2kMOQE4cdg2bsQsgcebwytvoTCZc1OJYSwU1IY7FlcNGyZANsmQe4C8OJMqPGqNEEVQmSKFAZ7dX6nMVfC1ZNQ4zXj0lG+omanEkLkAFIY7M39SFj/CQTNAXdveGMJVEh2nEMhhEg3KQz25NgKWDUIoq5A7T7QcBjk/s/EeEIIkSlSGOzB7cuw6gM4tgyKV4XX5oOXf+rbCSFEBsjIabYsIQGCv4cptYyZ1RoPh3c2S1EQIgs5Ozvj6+tL1apVefLJJ/nyyy9JSMj8zH3W0KVLF4oXL061atWy9DhSGGzV1dMwr7UxxpFnDei1A54ZCM4y/68QWenB1J5Hjhxh3bp1rFq1KtPTeoIx5HZmC0znzp1ZvXp1prOkRgqDrYmPNZqgTq8LVw5Bm2+g03IoWsHsZEI4nAdTe06ZMgWtNT/99BO1atXC19eXd955h/j4eAA+++wzKleuTNOmTenQoQMTJkwA4Ny5c1SpUoXevXvj5+fHxYsXk91HcssTq1+/Ph4eHln+uuUegy0JCTaaoF45DD4vQIsvoEAJs1MJkf3+HAKXD1l3n49VhxZj072ZrUztmZ0yVRiUUh7AL0BZ4Bzwqtb6xiPrNAQST3FUGWivtV6qlPoeeBaItDzXWWt9IDOZ7FJ0FPw1CnZ/C/kfg/YLoPLzZqcSQljYwtSe2SmzZwxDgA1a67FKqSGWx4MTr6C1/gvwhX8LyWlgbaJVPtBaL85kDvt1ar0x73LkBWN6zSafgFtBs1MJYa4M/GWfVWxlas/slNl7DG2BeZaf5wEvpLL+y8CfWuu7mTyu/btzFZZ0h/ntjOGwu6wxxjiSoiCEzbClqT2zU2bPGEporcMAtNZhSqnUzn3aA18+smyUUmo4sAEYkuPnfdYaDv4Cq4dC9G14dgg8MwBccpudTAiB7U7tCdChQwc2bdrE1atX8fLyYuTIkXTt2tXqv4NUp/ZUSq0HHkviqY+AeVrrQonWvaG1TnLOSKWUJ3AQKKm1jk207DKQC5gJ/KO1/jSZ7XsAPQC8vb39zaymGXbjnHHZ6J+N4FUL2nwNxe1zKkMhrE2m9rSuLJ3aU2ud7EA8SqkrSilPy9mCJxCewq5eBX5/UBQs+w6z/BitlJoLDEohx0yM4kFAQIB9TVQdH2fcWP5rFCgnaDnBuJ/gJK2FhcgJZGrPhy0DOgFjLd//SGHdDsDQxAsSFRWFcX/icCbz2J7Lh2BZP7i0H55oDs9PBHcvs1MJIaxIpvZ82FhgkVKqK3ABeAVAKRUA9NRad7M8LguUBjY/sv18pVQxQAEHgJ6ZzGM7Yu/B5nGw/WvI6wEvz4GqL8lcCUIIm5epwqC1vgY0TmJ5ENAt0eNzQKkk1muUmePbrLNbjKEsrp8B3zfguc+M4iCEEHZAej5b070bsPZ/sP9HKFwO3voDyjcwO5UQQqSLFAZr0BqO/mEMjX33Gjz9Hjw7GHLlNTuZEEKkmxSGzIoMNSbPObEKPH3hjcXg+aTZqYQQIsOkMGRUQgIEz4F1n0BCHDz3OQT2Amf5lQoh7Jt8imVE+HHj5vLFXVC+IbSaBB7lzE4lhBBWIYUhPeKiYdsk2DoRcuWDF76FJ9tLE1QhRI4iXW/T6sJumFEfNo0Bn7bQZy/4dpCiIEQOY6tTe168eJGGDRtSpUoVqlatyldffZVlx5IzhtTcvwUbPoW93xk9ll9fDBWbmp1KCJFFHkztCRAeHk7Hjh2JjIzM9PSeWmu01jhlcCgcFxcXJk6ciJ+fH7dv38bf35+mTZvi4+OTqVxJkTOGlBxfBVMDjaIQ2BN675KiIIQDsaWpPT09Pf8dg6lAgQJUqVKF0NDQLHndcsaQlNtX4M8P4ehSKO4Dr/0IXqkOSCiEsJJxe8Zx/Ppxq+6zskdlBtcanPqKj7DFqT3PnTvH/v37CQwMzPDvIyVSGBLTGvb/BGs/MsY6avQx1O0PLrnMTiaEMJEtTe0ZFRVFu3btmDx5MgULZs3EXlIYHrj2j9EE9dxW8K4Lrb+CYk+YnUoIh5SRv+yzii1N7RkbG0u7du14/fXXeemll9L5StJO7jHEx8LWL2F6XQj72+iT0HmlFAUhhE1N7am1pmvXrlSpUoUBAwZk0Ss2OPYZQ+g+WPYuXDkElVsZE+gU9DQ7lRDCRLY6tef27dv58ccfqV69Or6+vgCMHj2ali1bWv13kOrUnrYoICBABwUFZXwHMXfgr9GwaxrkKw4tx4NPG+sFFEKkm0ztaV1ZOrVnjnN6A6x4D25eAP/O0GQk5CmU6mZCCJEcmdozEaXUK8AnQBWglmWCnqTWaw58BTgD32mtx1qWlwMWAh7APuBNrXVMZjKlaHl/CP4eilSAzqug7NNZdighhOPIaVN7Zvbm82HgJWBLcisopZyBqUALwAfooJR60FVvHDBJa10RuAF0zWSelHmUh/ofQM/tUhSEECIZmZ3a8xiASnm8oFrAaa31Gcu6C4G2SqljQCOgo2W9eRhnH9MzkylFT/fPsl0LIUROkR3NVUsBFxM9DrEsKwLc1FrHPbJcCCGEiVI9Y1BKrQceS+Kpj7TWf6ThGEmdTugUlieXowfQA8Db2zsNhxVC2ButdWpXIEQaZLa1aaqFQWvdJFNHMM4ESid67AVcAq4ChZRSLpazhgfLk8sxE5gJRnPVTGYSQtgYNzc3rl27RpEiRaQ4ZILWmmvXruHm5pbhfWRHc9W9QEVLC6RQoD3QUWutlVJ/AS9jtEzqBKTlDEQIkQN5eXkREhJCRESE2VHsnpubG15eXhnePrPNVV8EvgGKASuVUge01s2UUiUxmqW21FrHKaX6AmswmqvO0VofsexiMLBQKfU5sB+YnZk8Qgj75erqSrlyMkWuLXDMns9CCOGA0trzWQbRE0II8RApDEIIIR5il5eSlFIRwPlUV0xaUYwWUfZMXoNtkNdgG3LCa4DseR1ltNbFUlvJLgtDZiilgtJyjc2WyWuwDfIabENOeA1gW69DLiUJIYR4iBQGIYQQD3HEwjDT7ABWIK/BNshrsA054TWADb0Oh7vHIIQQImWOeMYghBAiBQ5VGJRSzZVSJ5RSp5VSQ8zOkxFKqXNKqUNKqQNKKbvo/q2UmqOUCldKHU60zEMptU4pdcryvbCZGVOTzGv4RCkVankvDiilrD8ruxUppUorpf5SSh1TSh1RSvW3LLeb9yKF12A374VSyk0ptUcp9bflNYy0LC+nlNpteR9+UUrlMi2jo1xKsswkdxJoijHi616gg9b6qKnB0kkpdQ4I0FrbTbttpVR9IAr4QWtdzbLsC+C61nqspUgX1loPNjNnSpJ5DZ8AUVrrCWZmSyullCfgqbXep5QqAAQDLwCdsZP3IoXX8Cp28l4oY+jYfFrrKKWUK7AN6A8MAH7TWi9USn0L/K21zrqJy1LgSGcM/84kZ5lXeiHQ1uRMDkFrvQW4/sjithiz9mH5/kK2hkqnZF6DXdFah2mt91l+vg0cw5gcy27eixReg93QhijLQ1fLl8aY0XKxZbmp74MjFYbkZpKzNxpYq5QKtkxeZK9KaK3DwPjPDhQ3OU9G9VVKHbRcarLZSzCPUkqVBZ4CdmOn78UjrwHs6L1QSjkrpQ4A4cA64B9saEZLRyoM6ZoxzoY9rbX2A1oAfSyXOIQ5pgOPA75AGDDR3Dhpo5TKDywB3tNa3zI7T0Yk8Rrs6r3QWsdrrX0xJiirBVRJarXsTfX/HKkwJDeTnF3RWl+yfA8Hfsf4R2WPrliuFz+4bhxucp5001pfsfwHTwBmYQfvheWa9hJgvtb6N8tiu3ovknoN9vheAGitbwKbgNpYZrS0PGXq55MjFYZ/Z5Kz3O1vDywzOVO6KKXyWW64oZTKBzwHHE55K5u1DGPWPrDT2fsefJhavIiNvxeWm56zgWNa6y8TPWU370Vyr8Ge3gulVDGlVCHLz3mAJhj3Sh7MaAkmvw8O0yoJwNKEbTL/P5PcKJMjpYtSqjzGWQIYs+8tsIfXoJT6GWiAMXrkFWAEsBRYBHgDF4BXtNY2e3M3mdfQAOPShQbOAe88uFZvi5RS9YCt/9feHdsgDANhGP0tNmAQ9kFiEkqmQazDNDRcCqe5BuhQlPdKV4ms6LMuhaONUicAAABYSURBVJM8k7zX5WvmjH4Te/HhHc7ZyF6MMU6ZP5cPmYfzR1Xd1u/7nuSYeaPlpapef3nGPYUBgO/2NEoC4AfCAEAjDAA0wgBAIwwANMIAQCMMADTCAECzAFUHoE/a5CPHAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for i in range(t.ranks_tucker[0]):\n",
" plt.plot(t.Us[0][:, i].numpy(), label='Degree ${}$'.format(i))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are now ready to proceed with our usual optimization target, this time just with this \"smarter\" tensor:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"iter: 0 | loss: 0.998845 | total time: 0.0026\n",
"iter: 500 | loss: 0.816849 | total time: 1.5710\n",
"iter: 1000 | loss: 0.087476 | total time: 3.0443\n",
"iter: 1500 | loss: 0.011566 | total time: 4.4006\n",
"iter: 1737 | loss: 0.011055 | total time: 5.0428 <- converged (tol=0.0001)\n"
]
}
],
"source": [
"tn.optimize(t, loss)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Test relative error: tensor(0.1027, grad_fn=<DivBackward1>)\n",
"This model is more restrictive: it only has 48 degrees of freedom\n"
]
}
],
"source": [
"print('Test relative error:', tn.relative_error(y_test, t[X_test]))\n",
"print('This model is more restrictive: it only has {} degrees of freedom'.format(tn.dof(t)))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XucXXV97//Xm8kEJ9wmSLDJQAxQBOEXIJpCWigFLCTYAoHfUUHUIFg8Flq0EAhVm+DlQSQiwqlCoaBYkYsKMSiYRsFDRRMSSAIEyCFcGjIgxAOBSEaYJJ/zx1o72ZnsvWfPvl/ez8djHrP3d6+113cla81nfe+KCMzMzMqxQ70zYGZmzc/BxMzMyuZgYmZmZXMwMTOzsjmYmJlZ2RxMzMysbA4mlpOkDkl/kDS23nkxs8bnYNIi0j/8mZ/Nkvqy3p851O+LiE0RsXNErK5Gfs2GqtLXeNb3LpT0sUrmtR0Nq3cGrDIiYufMa0nPA5+KiF/k217SsIjYWIu8mVXCUK9xqy2XTNqEpK9Iul3SrZLWAx+T9OfpU9k6SS9JukZSZ7r9MEkhaVz6/vvp5/dKWi/pt5L2qeMpmW0jrZr9oqRnJf1e0i2SutPPdpJ0m6RX0+t9kaSRkq4E/gz497SEc2V9z6J5OZi0l1OBHwC7AbcDG4ELgD2AI4EpwKcL7P9R4IvA7sBq4MvVzKzZEE0HTgCOAvYC+oGr0s8+RVIT00NyvZ8PvB0RFwKLSUo5O6fvrQQOJu3l1xFxd0Rsjoi+iFgcEYsiYmNEPAtcD/xVgf1/FBFLIqIfuAU4rCa5NivOp4EZEfFiRPwRuAz4iCSRBJZRwH7p9b44It6sZ2ZbjdtM2ssL2W8kHQhcCbwfGEFyPSwqsP/vsl5vAHbOt6FZLaUBY2/gHknZs9fuALwTuBH4E+BHknYGvgd8MSI21TyzLcolk/YycIrofwMeB/40InYF/gVQzXNlVqZIpj/vBY6LiO6sn3dExO8j4q2I+JeIOBA4GvgQcHpm93rlu5U4mLS3XYDXgTclvZfC7SVmje46YLakvQEk7SnppPT1X0s6SNIOwBsk7YWZUsnLwL71yHArcTBpbxcC04D1JKWU2+ubHbOyXAH8Argv7bH4G+B96Wc9wE9IrvXHgXuAO9LPrgI+Iek1SVfUNsutQ14cy8zMyuWSiZmZlc3BxMzMyuZgYmZmZXMwMTOzsrXkoMU99tgjxo0bV+9sWAt7+OGHfx8Ro2p9XF/bVk3lXNctGUzGjRvHkiVL6p0Na2GS/rsex/W1bdVUznXtai4zMyubg4mZmZXNwcTMzMrmYGJmZmVzMDEzs7K1ZG8us3LMXdrLnPkreXFdH2O6u5g++QCmTuipd7bMylLt69rBxCzL3KW9XHrnY/T1J7OT967r49I7HwNwQLGGly9g1OK6djWXWZY581duueEy+vo3MWf+yjrlyKw4mYDRu66PYGvAyASYal/XDiZmWV5c1zekdLNGUShg1OK6djAxyzKmu2tI6WaNolDAqMV17WBilmX65APo6uzYJq2rs4Ppkw+oU47MilMoYNTiunYwMcsydUIPl582np7uLgT0dHdx+Wnj3fhuDa9QwKjFde3eXGYDTJ3Q4+BhTSdzzebr/lvt69rBxMysRdTzQcjVXGZmVjYHEzMzK5uDiZmZlc3BxMzMyuZgYmZmZataMJH0DkkPSVouaYWky9L0fSQtkvS0pNslDU/Td0zfr0o/H5f1XZem6SslTa5Wns3MrDTVLJm8BRwXEYcChwFTJE0CvgZcFRH7A68B56TbnwO8FhF/ClyVboekg4DTgYOBKcC3JW07MsfMzOqqasEkEn9I33amPwEcB/woTb8ZmJq+PiV9T/r5ByQpTb8tIt6KiOeAVcDh1cq3mZkNXVXbTCR1SFoGvAIsAJ4B1kXExnSTNUBmhE0P8AJA+vnrwDuz03Psk32scyUtkbRk7dq11TgdMzPLo6rBJCI2RcRhwF4kpYn35tos/a08n+VLH3is6yNiYkRMHDVqVKlZNjOzEtSkN1dErAN+BUwCuiVlpnHZC3gxfb0G2Bsg/Xw34NXs9Bz7mJlZA6hmb65RkrrT113AXwNPAvcD/yPdbBrwk/T1vPQ96ef3RUSk6aenvb32AfYHHqpWvs3MbOiqOdHjaODmtOfVDsAdEfFTSU8At0n6CrAUuDHd/kbgPyStIimRnA4QESsk3QE8AWwEzouITZiZWcOoWjCJiEeBCTnSnyVHb6yI+CPwoTzf9VXgq5XOo1kpXnjhBYD3SHoS2AxcHxFXS5oF/B2Q6QHyzxFxDyRjpUi6v28C/jEi5qfpU4CrgQ7g3yNidi3PxaxSPAW92RANGzYMYE1EHCRpF+BhSQvSj6+KiK9nbz9grNQY4BeS3pN+/C3geJK2wcWS5kXEE7U4D7NKcjAxG6LRo0cDbACIiPVpCaXQIhJbxkoBz6VVuZnS+aq0tI6k29JtHUys6XhuLrMypNP+TAAWpUnnS3pU0k2SRqZp+cZKeQyVtQwHE7MSSdoZ+DHw2Yh4A7gW2I9k+qCXgCszm+bY3WOorKW4msusNCIJJLdExJ0AEfHylg+lG4Cfpm8LjZXyGCprCS6ZmA1RMvyJdwNPRsQ3MumSRmdtdirwePo631ipxcD+6Uzaw0ka6efV4BTMKs4lE7MhevDBByGZN+64dO45gH8GzpB0GElV1fPAp6HwWClJ5wPzSboG3xQRK2p4KmYV42BiNkRHHXUUwMMRMXHAR/fk2yffWKl0HEre/cyahau5zMysbA4mZmZWNgcTMzMrm4OJmZmVzcHEzMzK5mBiZmZlczAxM7OyOZiYmVnZPGjRzKyJzF3ay5z5K3lxXR9juruYPvkApk4otAJCbTiYmJk1iblLe7n0zsfo609WLu9d18eldz4GUPeA4mouM7MmMWf+yi2BJKOvfxNz5q+sU462cjAxM2sSL67rG1J6LTmYmJk1iTHdXUNKryUHEzOzJjF98gF0dXZsk9bV2cH0yQfUKUdbuQHezKxJZBrZ3ZvLrIYatQulWTmmTuhpyOvYwcRaUiN3oTRrRW4zsZbUyF0ozVpR1YKJpL0l3S/pSUkrJF2Qps+S1CtpWfrzwax9LpW0StJKSZOz0qekaaskzahWnq11NHIXSrNWVM1qro3AhRHxiKRdgIclLUg/uyoivp69saSDgNOBg4ExwC8kvSf9+FvA8cAaYLGkeRHxRBXzbk1uTHcXvTkCRyN0oTRrRVULJhHxEvBS+nq9pCeBQpXVpwC3RcRbwHOSVgGHp5+tiohnASTdlm7rYNLmCjWwT598wDZtJtA4XSjNhqoZOpPUpM1E0jhgArAoTTpf0qOSbpI0Mk3rAV7I2m1NmpYvfeAxzpW0RNKStWvXVvgMrNFkGth71/URbG1gn7u0F0ga2S8/bTw93V0I6Onu4vLTxjfcDWg2mMGu9UZR9d5cknYGfgx8NiLekHQt8GUg0t9XAmcDyrF7kDvgxXYJEdcD1wNMnDhxu8+ttRRqYM8EjEbtQmk2FMVc642gqsFEUidJILklIu4EiIiXsz6/Afhp+nYNsHfW7nsBL6av86Vbm3IDuzWjUqqrmuVar2ZvLgE3Ak9GxDey0kdnbXYq8Hj6eh5wuqQdJe0D7A88BCwG9pe0j6ThJI3086qVb2sOjTxHkVkupVZXNcu1Xs02kyOBjwPHDegGfIWkxyQ9ChwLfA4gIlYAd5A0rP8cOC8iNkXERuB8YD7wJHBHuq21sUaeo8gsl1LHPjXLtV7N3ly/Jnc7yD0F9vkq8NUc6fcU2s/aTyPPUWSWS6nVVc1yrXs6FWtabmC3ZlLO2KdmuNY9nYqZWQ00S3VVqVwyMTOrgWapriqVg4mZWY00Q3VVqVzNZWZmZXMwMTOzsjmYmJlZ2RxMzMysbA4mZkP0wgsvALwnx8Jvu0taIOnp9PfINF2SrkkXd3tU0vsy3yVpWrr905Km1eeMzMrnYGI2RMOGDQNYExHvBSYB56WLu80AfhkR+wO/TN8DnEgy19z+wLnAtZAEH2AmcATJ2j0zs5ZkMGsqDiZmQzR69GiADZAs/EYyZ1wPyaJtN6eb3QxMTV+fAnwvEguB7nTC08nAgoh4NSJeAxYAU2p2ImYV5GBiVoYBC7+9K11hNLPS6J7pZl74zVqeg4lZiQYu/FZo0xxpUSB924SI6yNiYkRMHDVqVGmZNasyBxOz0ogBC78BL2fW60l/v5Km51v4rdCCcGZNxcHEbIgiAuDdDFj4jWTRtkyPrGnAT7LSP5H26poEvJ5Wg80HTpA0Mm14PyFNM2s6npvLbIgefPBBgHeSLvyWJv8zMBu4Q9I5wGrgQ+ln9wAfBFaRNNx/EiAiXpX0ZZLVRAG+FBGv1uQkrOpKWaK3mTmYWENrxBvyqKOOAng4Iibm+PgDAxMiKcqcl+u7IuIm4KaKZtDqLrNEb2ZlxcwSvUDdr99qcTWXNaxS18w2q7dSl+htZg4m1rDy3ZCX3b2iTjkyK06pS/Q2MwcTa1j5brzXNvS7dGINLd9SvMUs0dusHEysYRW68Vq5usCaX6sv0ZuLg4k1rEI3XitXF1jzmzqhh8tPG09PdxcCerq7uPy08S3b+A7uzWUNJFfPre6uTtb19W+3bStXF1hraOUlenNxycQaQr6eW3976Oi2qy4wa0YFSyaS/qnQ5wNG/5qVLF/PrfufWsvlp42v+FiTb3yj8KX7T/9U8NI3swEGq+baJf19APBnJNNCAJwEPFCtTFn7ydcG0ruujznzV1Z8sOL69esBWLlyJYsXL+bkk08G4O677+boo4+u2HHM2kXBYBIRlwFI+k/gfenaDUiaBfyw0L6S9ga+B/wJsBm4PiKuThcEuh0YBzwPfDgiXpMk4GqSaSc2AGdFxCPpd00DvpB+9Vci4maspYzp7qK3QECp9OjhmTNnAnDCCSfwyCOPsMsuyXPTrFmz+NCHPlRoVzPLodg2k7HA21nv3yYJBoVsBC70anRWjFxdKbNVa/Tw6tWrGT58+Jb3w4cP5/nnn6/4ccxaXbG9uf4DeEjSXSTrLZxKUurIK50VNbNQ0HpJ2avRHZNudjPwK+ASslajAxZKyqxGdwzpanQAkjKr0d1aZN6tCWRKHBfesZxNsd2SHkB1ugN//OMf5/DDD+fUU09FEnfddRef+MQnKn4cs1ZXVDCJiK9Kuhf4yzTpkxGxtNiDFFqNTlLFVqMjKdEwduzYYrNmDWTqhB4+d/uyvJ9Xozvw5z//eU488UT+67/+C4DvfOc7TJgwoeLHMWt1Q+kaPAJ4IyKuBtZI2qeYnbwanQ1FvoAhCg9iLMeGDRvYddddueCCC9hrr7147rnnqnIcs1ZWVDCRNJOkKurSNKkT+H4R+3Xi1ehsCHK1nQg4c9LYqgwAu+yyy/ja177G5ZdfDkB/fz8f+9jHKn4cax9zl/Zy5Oz72GfGzzhy9n1tM49csSWTU4GTgTcBIuJFtnYbzintnXUjXo2ubZVyU+WahuKqjxzGV6aOr0oe77rrLubNm8dOO+0EwJgxY7Z0GzYbqnZeNqHYBvi3IyIkBYCknYrY50jg48BjXo2u/RRaHAgoOAixltNQDB8+HEkkzz7w5ptv1uS41poKrWPS6lOrFBtM7pD0b0C3pL8Dzgb+vdAOEfFrcrd3gFeja3n5bqpZ81bw1sbNDbMC3Yc//GE+/elPs27dOm644QZuuukmPvWpT9U8H9Ya2nEdk4xie3N9XdLxwBsko+H/JSIWVDVn1tTy3Ty5Jm3MfnKr9TK9F110EQsWLGDXXXdl5cqVfOlLX+L444+v2vGsteUbfNsOE5MW2wD/tYhYEBHTI+KiiFgg6WvVzpw1r6HePC+u66tLffMll1zC8ccfz5w5c/j617/O8ccfzyWXXFK141lra8d1TDKKbYDP9ah2YiUzYq0l3001ckRnzu3HdHfVZd3sBQu2L2Dfe++9VTuetbZ2XMckY7BZgz8D/D2wn6RHsz7aBfhNNTNmzS1z8wyssgK2aZiHrU9u+QYsVqO++dprr+Xb3/42zzzzDIcccsiW9PXr1/MXf/EXFT+etY92W8ckY7A2kx8A9wKXs3UOLYD17lFlgyl0U+VqF5kzf2XN6ps/+tGPcuKJJ3LppZcye/bsLem77LILu+++e8WPZ9bqBps1+HXgdUlXA69mzRq8i6QjImJRLTJprSVfkJk++YC8pZZK22233dhtt9244IIL2H333bfMGrx+/XoWLVrEEUccUfFjmrWyYttMrgX+kPX+zTTNLK+hDlqsR33zZz7zGXbeeect73faaSc+85nPVO14Zq2q2HEmSseBABARmyV5/XjLa+7SXqb/cDn9m5PLpnddH9N/uBwoPJ6k1vXNEbFlwCLADjvswMaNG2t2fLNWUWzJ5FlJ/yipM/25AHi2mhmz5jZr3ootgSSjf3Nw6Z2PbrdtPecy2nfffbnmmmvo7++nv7+fq6++mn333bdmxzdrFcUGk/8J/AXQSzLx4hGk072bDTR3aW/OwYkAff2btwkW9Z7L6LrrruM3v/kNPT097LXXXixatIjrr7++Jsc2ayXFjoB/BTi9ynmxFpAJDoVkz1NU77mM9txzT2677baqH8es1Q02zuTiiLhC0v8i9xoi/1i1nFlTyhUcBsru/luvuYyuuOIKLr74Yv7hH/5hmzaTjGuuuaaqxzdrNYOVTJ5Mfy+pdkas+c1d2ptznMhAHVl/vOs1l9F73/teACZOnFjV45i1i8HGmdyd/r65NtmxZlVM9VZG9hrvtRxbku2kk04CYNq0aYNsaWbFGKya625yVG9lRMTJFc+RNZ25S3u58I7l2wSJQnqySh35pl2pdnvJSSedlLN6K2PevHlVPb5Zqxmsmuvr6e/TgD9h61K9ZwDPVylP1kQyJZJiA4lI2kyOnH3flqBRj7mMLrroIgDuvPNOfve7321ZqvfWW29l3LhxNc2LNaZaL4fQ7Aar5vrfAJK+HBFHZ310t6QHqpozawqDNbh3d3Wy047D6F3Xh9hazK33olh/9Vd/BcAXv/hFHnhg66V80kkncfTRR+fbzdpEoZVCHVByK3acyShJW0ZySdoHGFWdLFkzKdTrqnMHMevkg3lwxnH0dHdtV19a7enli7F27VqefXbr+NvnnnuOtWvX1jFH1gjqsRxCsys2mHwO+JWkX0n6FXA/8Nmq5cqaRqFeV5si+Nztyzhy9n15e3n1poti1ctVV13FMcccs+Xn2GOP5Zvf/GbBfc4++2yAQyU9nkmTNEtSr6Rl6c8Hsz67VNIqSSslTc5Kn5KmrZI0A2sY+R6Siumt2K6KHbT4c0n7AwemSU9FxFvVy5Y1i+mTD+CzedYhycymMrCKa6B6Vh9MmTKFp59+mqeeegqAAw88kB133LHgPmeddRbf+c53ns7x0VUR8fXsBEkHkQz4PRgYA/xC0nvSj79FsvDcGmCxpHkR8URZJ2QVka/LukiqwFzVtb1il+0dAUwHzo+I5cBYSX9b1ZxZU5g6oYfurtyrJ2YLkhsxl3pWH2zYsIE5c+bwr//6rxx66KGsXr2an/70pwX3SdtUip0N8hTgtoh4KyKeA1YBh6c/qyLi2Yh4G7gt3dbqJHuOuDffyv3fG+CqrjyKreb6DvA28Ofp+zXAV6qSI2s6s04+mM6O/N1sMwr196r2iPd8PvnJTzJ8+HB++9vfArDXXnvxhS98odSvO1/So5JukjQyTesBXsjaZk2ali99O5LOlbRE0hK351THwDni8s0tB/W7VhtdscFkv4i4AugHiIg+8j9oWpMraRbfInoG93R3bTPGJFu1R7zn88wzz3DxxRfT2ZmUrrq6uogiuzkPcC2wH3AY8BJwZZqe6z7JV1DLeeCIuD4iJkbExFGj3O+lGoqZBiijXtdqoyt2TZK3JXWRXuyS9gPcZtKCSukSOWf+yu2mmx8oe1R7PUa85zN8+HD6+vq2DGB85plnBm0zySUiXs68lnQDkKkrWwPsnbXpXsCL6et86VZjxZY26nmtNrpig8lM4OfA3pJuAY4EzqpWpqx+SpnFd7AeLh3SdismNspgsMsuu4wpU6bwwgsvcOaZZ/Lggw/y3e9+d8jfI2l0RLyUvj0VyPT0mgf8QNI3SBrg9wceIimZ7J92s+8laaT/aHlnY6XqHtHJaxu2r9oa0bkDI3fasSGu1UY3aDBR8sj2FMko+EkkN8EFEfH7KufN6qBQl8jsXiyZ0cGDBZKuzo7tAkk9RrznEhEceOCB3HnnnSxcuJCI4Oqrr2aPPfYouN8ZZ5wBSc9GSVpD8rB1jKTDSErvzwOfTo+xQtIdwBMkjfbnRcQmkp3PB+YDHcBNEbGiGudpg8tXszl8WAcPzjiutplpUiqmfljSwxHx/iF9sXQT8LfAKxHx/6Vps4C/AzKtiP8cEfekn10KnANsAv4xIuan6VOAq0luuH+PiNmDHXvixImxZIknOi5FoTEhme693V2dvPn2Rvo3Fb52urs6mXXywQ0ROPJ5//vfz8MPPzzk/dJ7ouZTDvvaro59ZvwsZ4OVgOdm/02ts1M35VzXxTbAL5T0Z0P87u8CU3KkXxURh6U/mUCS3Rd/CvBtSR2SOkj64p8IHASckW5rVTJ98gF0dXbk/Cxzs63r6x80kAAsm3lCQwcSgEmTJrF48eJ6Z8PqLF+juhvbi1dsMDmWJKA8k3Z9fEzS9ot5Z4mIB4BXi/x+98VvEFMn9HD5aePL/p7M4K5Gd//99zNp0iT2228/DjnkEMaPH88hhxxS72xZjeV6iHJj+9AU2wB/YgWPeb6kT5AsuHVhRLxG0r9+YdY22X3uB/bFPyLXl0o6l3Rd+rFjx1Ywu+1n6oSeotpDCskM7mr0ksm9995b7yxYlRUz+2+9lkJoJYOtZ/IO4H8Cfwo8BtwYEcWO/M3lWuDLJH9rvkzSF/9s8ve5z1VyytsXH7geknrlMvJoJE9q03+4fNAuv4U08uCuP/7xj1x33XWsWrWK8ePHc8455zBsWLHPVtYshtLVvVE6hjSrwe6em0kGKv4XW9stLij1YO6L3/iK7aVVjEaub542bRqdnZ385V/+Jffeey9PPPEEV199db2zZRVS6DoerKu7lWawYHJQRIwHkHQjSf/4krkvfmMb+BRXjkavb37iiSd47LHkCfWcc87h8MMPr3OOrFKKuY4budTcrAYLJltG8UTExkLLnA4k6VbgGGAP98VvDkOZUiIfQVPUN2emTwFcvdViirmOG7nU3KwGu4sOlfRG+lpAV/peQETErvl2jIgzciTfWGD7rwJfzZF+D3DPIPm0Cij3aa2nu6tpBngtX76cXXdNLt+IoK+vj1133ZWIQBJvvPHGIN9gjWqw67jRS83NarBle3MPOLCWlG8Nh2I02w26aVP5VXnWmApdx7mm9rHKKHacibWBYw8sfUZa36DWKAo91GyO8HVaJQ4mBiSNlrcsWl3Svj3dXb5BrWFMndDDyBG5F2xzW0n1OJjYlt4vpSzj0dmhpqresvYw86SDPaK9xtyNpY2VO6Zk5IhOZp7U2BM5WnvyiPbaczBpU+WMKWmmXlvWOoqZFiWbR7TXloNJmypnTEklRsebDUWhaVHAJZBG4GDSpsoZU9IxhMGrZpWQbwXQWfNW8NbGzUNaZtqqww3wLWru0l6OnH0f+8z4GUfOvm+76eC78/R2KcamUlrqzcqQ7+FnXV9/3mWmrbZcMmlBuaoEPnf7Mj57+zJ6urs49sBR/OGPpU/+3OPulVZjQx1Q67m3as8lkxaUq0ogU5boXdfH9xeuLnlqeXevtHrIt3iVx5M0DpdMWlAln8pGjuhkxPBhbty0usrX1RfYrlfiwAeeofYCs9I4mLSgcubYyibwOBJrGIW6+uYLFkNZHMvK42DSgqZPPqDsdUkEnDlprG84a3iDBZl8DfS+tivLwaQFZW6SC+9YXlLPqw6JKz98qG82a1iDVV0NNruDG+grzw3wLWrqhB42l9iF94wj9nYgsYb1hbmP8bnbl9G7ro9ga9VVpvt7pmqrUFWvG+grz8GkhZV6w9z/1NoK58SsMuYu7eWWhasZ+JiUPbZksNkd3COxOlzN1aLmLu1lw9uljSXxdCnWqObMX7ldIMnIXLeFqrB63JurahxMWszcpb3MmreCdX399c6KWcUN1tYxd2lv3t6MnqC0ulzN1UIydcUOJNaqBqu6nTVvRd4Bjq7aqi4HkxZSzkzA2TyRozWqXIEi27q+fubMX8n///4eerq7EEmJxMtKV5+ruVpIpbo7nnHE3hX5HrNKywSEz96+LO82vev6+PHDvQ4gNeaSSQspZyZgAAk+NmksX5k6vkI5MivOYLNcZyu0xnuGZw6uPZdMWkg5M8MLeO7yv6lYXsyKVcqUJ39zyGi+v3B1we/1wMTacsmkhZTT8O5BXFYv+aY8ufCO5XlLKsWMhfI1XVsumTSh7KkkduvqRIJ1G0oPJO7pYvWUrwSRmQooV0llsFKHr+naq1rJRNJNkl6R9HhW2u6SFkh6Ov09Mk2XpGskrZL0qKT3Ze0zLd3+aUnTqpXfZpE9VUSQlEZe29CfdyDXYNzTxeqtmBLEwDaQQvv4mq6PalZzfReYMiBtBvDLiNgf+GX6HuBEYP/051zgWkiCDzATOAI4HJiZCUDtqlLdf7s6O/jmRw7jwRnH+aYborPPPhvgUD8oVcZg3X0zsgci5htL4mu6fqoWTCLiAeDVAcmnADenr28Gpmalfy8SC4FuSaOBycCCiHg1Il4DFrB9gGorlWpU9JNb6c466yyApwck+0GpRFMn9HD5aeMHXQ5asKXtJHsfjyVpDLVuM3lXRLwEEBEvSdozTe8BXsjabk2ali99O5LOJblZGTt2bIWz3Ti6R3TyWhntI5DceL7pSnf00UcDbGTbh7FTgGPS1zcDvwIuIetBCVgoKfOgdAzpgxKApMyD0q3VP4PGk1mTZJ8ZP8tbZRuwzTokhdYxsdprlN5cuYZcR4H07RMjro+IiRExcdSoURXNXCP5wx/LnyrFDZNVsc2DElDRByVJSyQtWbu2tWd0Hqz9xN19G1etSyYvSxqdlkpGA686k6f2AAAOm0lEQVSk6WuA7GHXewEvpunHDEj/VQ3y2RDmLu3lsrtXbCmJjOjcgf7N5X1nd1enn+ZqqyIPSsD1ABMnTixjNFHjG2yVUHf3bVy1LpnMAzINjdOAn2SlfyJtrJwEvJ4+3c0HTpA0Mq1PPiFNa3lzl/Yy/UfLt6nS2lBmJOnq7GDWyQeXmzXL7eX0AYkhPCjlSm9rmbaQ7q7tR7i7u29jq2bX4FuB3wIHSFoj6RxgNnC8pKeB49P3APcAzwKrgBuAvwdI65O/DCxOf76UqWNudXPmr6R/U/kPoR2SGyhrww9KRShm2pSpE3pYNvMEvvmRw9zA3kSqVs0VEWfk+egDObYN4Lw833MTcFMFs9YUKrFAVVdnh2/AKjjjjDMADiTp+buGpFfWbOCO9KFpNfChdPN7gA+SPChtAD4JyYOSpMyDErTBg9JQp01xA3tz8Qj4BtUhbRkBPBQ7DtuBtzduZoxXlKuaW2+9ldtuu+3RiJg44CM/KBWQb9qU7B5a1rwcTBpUKYEEYI+dd/RqclZ32VP+ZB5s8vXEcg+t1uBg0kCyb0CptFmAfWNaveWqzvrc7cvyjh9xD63W4GDSIOYu7eWfbl/Glv5aJba9+8a0estVnZXvcu7skHtotQgHkzrKLolUYvCAu05aIxhK6Xin4cPcXtIiHEzqZGBVQLl63OBuDWJMd1fRvRFfL2MNHmssjTKdStup1Oy/kAQSz5RqjaLYWYDB1bKtxCWTOqlUQ7mrtqzRZB5q5sxfWbCE4mu3tbhkUgdzl/ayg3JNzTQ0HhVsjWrqhB4enHEc3/zIYXTusP21PnJEp6/dFuOSSY19Ye5j3LJwdVkN7jsN72DFl9p6WRdrEnPmr6R/8/ZX+wg3vLccB5Mamru0t+xAAvDVU8dXJD9mlZBrgGImUOSr5vJ4qNbjYFJFmZusd11fydOjDCTlnsfIrB4KzbcFyTz7ua763bo6OXL2fTkDkDUnB5MqGXiTVSKQQGmj4s2qpdB8W5B/sOKbb29kXdoteLAJH605uAG+SirZ9TfbYOtkm9VSvuqq3nV9BXtyDVxeITsAWXNyyaRKqlEn7K6UVi/52kWGMkBxMG5HaW4umVRJJQZjdXd1MnJEpxcHsrrKVNn2ptP+ZKql5i7tHdIARUjaUEaO2H4VRfAAxmbnkkkFZT+9dXaUN45EwLKZJ1QmY2ZlyNcuMmveCpbNPIEfLlnNg88Ut65XADNPOni7qYRc6m5+DiYVMHdpL5fdvWKb9drfLnPJXT+lWaPIV/20rq+fL8x9jN8UGUggKWFnj5B3b67W4WBSpkpP2Ah+SrPGUqhd5NZFLxQ9bir7uvaSvK3HbSZlmLu0lwvvWF7RQOK2EWs0hR5shtLl3dd1a3PJpESZEkmlxo90dXb4ZrO6y9dra2A1bkaxg3Gzq7esNblkUqJKjCPpSCd7dGnEGkGhXlszTzp4u15bXZ0dnHHE3oP25hJw7IGjqpdxawgumZSonL71LoVYIyo0mv3BGcdt2WZgqWXiu3ffJn3cO7v4zTOvbmlLCeDHD/cy8d27+5pvYQ4mJRrqXFsdEpsj3HPFGkKu6qx8vbYy6fkazQemHzn7vu0a5TNBydd963IwKcHcpb1DCiQCrvzwob6RrCHkm5yxe0RnznaRoXZTHywoWWtym8kQZW7EYgk4c9JYBxJrGPmqsyLI2S4y1G7q+YKPx061troEE0nPS3pM0jJJS9K03SUtkPR0+ntkmi5J10haJelRSe+rR54zhtLw3tPdxVUfOYyvTPX6I9Y48pUQXu/r5/LTxtPT3VXWFD65pljx2KnWV89qrmMj4vdZ72cAv4yI2ZJmpO8vAU4E9k9/jgCuTX9XTaHFfootqo8c0bml0dKskeQbhDgm7b5bbinaI9zbUyO1mZwCHJO+vhn4FUkwOQX4XkQEsFBSt6TREfFSJQ46MHAce+Aofvxwb87FfoYyS+q6HHXPZo1g+uQDBp0bq9ADVTE8wr391KvNJID/lPSwpHPTtHdlAkT6e880vQd4IWvfNWnaNiSdK2mJpCVr164tKhO5+tXfsnB1wcV+ip0l1fXD1qimTugpWJ1VaLyJWT71KpkcGREvStoTWCDpqQLb5pp+d7uuVBFxPXA9wMSJE4vqapWr/SPfjtndIzP7vriuj+4Rnfzhjxvp37x1T9cPW6MrVHIoNN7EpQ3Lpy7BJCJeTH+/Iuku4HDg5Uz1laTRwCvp5muAvbN23wt4sRL5GEpXxeySxsAbsdwqAbNqKeXadNdeK0XNg4mknYAdImJ9+voE4EvAPGAaMDv9/ZN0l3nA+ZJuI2l4f71S7SXFtn8MVtJw/bA1onzjSaDwWuuFGujN8qlHm8m7gF9LWg48BPwsIn5OEkSOl/Q0cHz6HuAe4FlgFXAD8PeVykgx7R8jR3R66hNrSoWqqwpx114rRc1LJhHxLHBojvT/C3wgR3oA55V73Fy9tu5/au2gY0ZGDB/mQGJNqdTqKnfttVI0UtfgqslV3P/+wtVF7et6Ymt0+dpFyqmuctWtDVVbTKdSznTxrie2RlaoG6+rq6yW2iKYlFO68I1nQ1HrqYIG68ZbielRzIrRFtVcxfbaysU3npWgZlMFlTptvFmltUXJpNhR6wP1uIrLKuMUkimCSH9PzUr/XiQWAt3pGKuieYZeaxRtEUxyFfc/NmnslvcjR3TSucO2A+1dt2wlqulUQW4XsUbRFtVcMHhx36PYrUJqOlWQu/Fao2ibYDIY1y1bJdRjqiBfu9YI2qKay6wWJO0kaZfMa5Kpgh5n61RBsP1UQZ9Ie3VNooJTBZnVmksmZpXzLuAuSZDcWz+IiJ9LWgzcIekcYDXwoXT7e4APkkwVtAH4ZO2zbFYZDiZmFVKvqYLMGoGruczMrGwOJmZmVjYlJe3WImkt8N95Pt4D+H2ez5qVz6n23h0Ro2p90EGu7XLU+9+7nsf3uW9V8nXdksGkEElLImJivfNRST4nK1e9/73reXyfe2WO7WouMzMrm4OJmZmVrR2DyfX1zkAV+JysXPX+967n8X3uFdB2bSZmZlZ57VgyMTOzCnMwMTOzsrVVMJE0RdLKdJnUGfXOTyGSbpL0iqTHs9KGvPyrpGnp9k9LmpbrWLUgaW9J90t6UtIKSRc0+zk1m1ouKVzv6zfP8WdJ6k3Pf5mkD2Z9dml6/JWSJmelD/lvRj2v9QLHrv65R0Rb/AAdwDPAvsBwYDlwUL3zVSC/RwPvAx7PSrsCmJG+ngF8LX39QeBekvUxJgGL0vTdgWfT3yPT1yPrdD6jgfelr3cB/g9wUDOfU7P9AM8DewxIG9K/f7Ncv3mOPwu4KMe2B6V/D3YE9kn/TnSU+jejntd6gWNX/dzbqWRyOLAqIp6NiLeB20iWTW1IEfEA8OqA5KEu/zoZWBARr0bEa8ACYEr1c7+9iHgpIh5JX68HniRZVbBpz6lFVGVJ4Xpfv3mOn88pwG0R8VZEPEcyi/PhlPg3o57XeoFjV/3c2ymYFLVEaoMb6vKvDXnOksYBE4BFtMg5NYmKLyk8RI3wf31+WpV0U6aaqZrHr+e1PuDYUOVzb6dgUtQSqU0q37k13DlL2hn4MfDZiHij0KY50hrynJrIkRHxPuBE4DxJRxfYtpb/zrX6v74W2A84DHgJuLKax6/ntZ7j2FU/93YKJhVZIrXOXs5UNai45V8b6pwldZJc4LdExJ1pclOfUzOJrCWFgW2WFIai//3LUdf/64h4OSI2RcRm4AaS86/K8et5rec6di3OvZ2CyWJgf0n7SBoOnE6ybGozGeryr/OBEySNTIu1J6RpNSdJwI3AkxHxjayPmvacmokaY0nhuv5fD2jzOZXk/DPHP13SjpL2AfYHHqLEvxn1vNbzHbsm5z5Yz4RW+iHpNfF/SHopfL7e+Rkkr7eSFEf7SZ4SzgHeCfwSeDr9vXu6rYBvpef1GDAx63vOJmlUWwV8so7ncxRJMflRYFn688FmPqdm+iHplbM8/VmRuf5L+fdvhus3z/H/I/3+R9M/jKOztv98evyVwIlZ6UP+m1HPa73Asat+7p5OxczMytZO1VxmZlYlDiZmZlY2BxMzMyubg4mZmZXNwcTMzMrmYNJgJL1L0g8kPZtOe/FbSafWOA/PS9pDUrekvy/xO86SNGbgd1Yul9ZsfG23NgeTBpIOOJoLPBAR+0bE+0kGC+01YLthNcpSN5DzhpPUMci+ZwFjBtnG2oSv7dbnYNJYjgPejojrMgkR8d8R8b/Sp6EfSrqbZLI+SZoj6XEla1R8BEDSMZJ+mtlf0r9KOit9/bykyyQ9ku5zYJr+Tkn/KWmppH9j67w8s4H9lKx/MCf97vsl/QB4TNI4bbtexEVK1k34H8BE4JZ03650k38YeGxrG762W5yDSWM5GHikwOd/DkyLiOOA00gmbTsU+GtgjoqbIvz3kUz2dy1wUZo2E/h1REwgGR07Nk2fATwTEYdFxPQ07XCS0bAH5TtARPwIWAKcme7bV+DY1h58bbc4B5MGJulbkpZLWpwmLYiIzBoNRwG3RjJ528vA/wb+rIivzUw69zAwLn19NPB9gIj4GfBagf0fimTdg1LkOra1IV/brcfBpLGsIFkdDoCIOA/4ADAqTXoza9tcU0QDbGTb/9d3DPj8rfT3JiC7frrYeXWy8zDYsQbKd2xrfb62W5yDSWO5D3iHpM9kpY3Is+0DwEckdUgaRfIE9hDw38BBSmYB3Y3khh3MA8CZAJJOJFkiFGA9ydKf+bwM7JnWS+8I/G3WZ4Pta+3F13aLa8sI2qgiIiRNBa6SdDGwluRp6RKga8Dmd5HUMy8nefK6OCJ+ByDpDpLZQZ8GlhZx6MuAWyU9QlKlsDrNz/+V9GDaEHkv8LMB+e2X9CWSldyeA57K+vi7wHWS+tJ8Whvztd36PGuwmZmVzdVcZmZWNgcTMzMrm4OJmZmVzcHEzMzK5mBiZmZlczAxM7OyOZiYmVnZ/h/zIL2pG414NwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"show()"
]
}
],
"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.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}