{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Refactoring Trees: An exercise in Research Software Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise, you will convert badly written code, provided here, into better-written code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will do this not through simply writing better code, but by taking a refactoring approach, as discussed in the lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As such, your use of `git` version control, to make a commit after each step of the refactoring, with a commit message which indicates the refactoring you took, will be critical to success.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will also be asked to look at the performance of your code, and to make changes which improve the speed of the code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The script as supplied has its parameters hand-coded within the code. You will be expected, in your refactoring, to make these available as command line parameters to be supplied when the code is invoked."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some terrible code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's our terrible code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABFXklEQVR4nO3dd1yVdf/H8df3sEFZAoKIgDhx4MBtOVKzZcOGDW9tmQ3b427f1e3dsGllplmapWZlZablzK05cQ9kiIACsjeH8/39AfVDRTkHDxw4fp6Px3l4uK7rfK+3h8PnXON7XV+ltUYIIUTjZ7B1ACGEENYhBV0IIeyEFHQhhLATUtCFEMJOSEEXQgg74WirFfv5+emwsDBbrV4IIRqlHTt2ZGit/aubZ7OCHhYWxvbt2221eiGEaJSUUonnmyeHXIQQwk5IQRdCCDshBV0IIeyEFHQhhLATUtCFEMJOSEEXQgg7IQVdCCHshBR0YR356bDs31CUfd5F9p/ez4w9MygpL7G4eVNxMenTplF88OB5lykuKGPDwqMU5ZVa3L41FZeb+CjhFDF5heddJqeojNd+3U96nuXvBVrD7vmw78cLLKJ5b/t7xGXHWd6+aLSkoIszHEvP5/3lh7H4Pvl5KbB1OqybUu1srTXv/PUO3x78lrLyMotz6ZISsr6ey6m33j5vtm1L4tmzJonCXMsKutaajd/NJTMl+YzppaWZxMa+jdFYYFF7ZVoz80Q6rxxNPm/WT1YfZfamBE7lFlvU9j92zIZlz0FxbrWzl8UvY/b+2ezJ2GNx05nz5lGw9a/a5RI2JQX9EpNWUsbGrLzzzl93JJ2pq2NZdTDNsoaDoqD7nbD1czh97JzZKxJXsDNtJ5O6T6KJcxNLY+Pg5YXfpEco3LqV/NWrz5mfdbKAfWuTiRzYgmbB5rVvMpk4cOAAR7ZuYsui7zi+L+aM+UVFiSQen0Fi4nSKY7Mozzfvi6KpowPPtQ5ka04BS9JzzpmfkFHA7E0J3NKzJZ2Dvcxq8wxKwcj/QUE6bHj/nNnFxmI+2PkBHX07MipilEVNlyYlkfbmW+T89NN5l8lMLSDjxPk/Q8J2pKA3Elpr1q9fz9q1ay+wjIm9+yaRcfrP8y7zxYl0Ru8+xq27Y9mZe+6W5119Q2nt78H/lh6k1GiyLOTQl8HBGVa8csbkkvIS3t/xPu182nFjmxsta7MKn9tuw7lNBKfeeQddemZx3fhjLI7OBnpf19rs9o4cOcLChQv5fslSXMPa0mXoiDPme3l1x7/JVSTGf0HqN3+SvynF7LbvCGpGpIcrrx9Lobj8zPfxf0sP4uxg4OkR7c1u7xzBPSHqdtg8DbISzpg1Z/8cThac5Nlez2JQlv2Jp015Fxwd8X/iiXPm5WYUsXL2ARa8vpVNi8790v5b+vEEFr35KkX55y/6hXvTyfo5FlNpuUX5xIVJQa9nphIjWb/EUnw406LXKaXIyMhg3bp1ZGZW/9qTJ38iLW0pZWXZ523nybBAXmvTgn35RVy94yjj98ZxML/on/lODgZeuqYjcRkFzN1y3ltGVK9pIFz2JBxaAvHr/pn8zYFvSM5P5tlez+JgcLCszSqUoyPNn3uOssTjZH4775/pxw+cJnHvaXpeHYa7p7PZ7bVr144e4a0wlRtJd/Pii1mzOHr0KFprSlPyyZi9n6ZLh4HWZA9ehueQVma37aAUr7UJJqm4lJkn0v+ZvulYBssPnOKhIW0I8HQ1u71qXfEKGBxgxav/TEorTGPWvlkMDx1OdGC0Rc0VbttG3vLlNLv/PpyaB/wzvSC7hLXzDvPtq1uI3ZFG1LBWDL8n8oJtJcTsYvMP86qdp8tM5PwWT2liLsrRshKUnZ3NwoULyc7Otuh1lwop6HVg9x+/kXWy+q055WCg5Gg22b/Focst2wK+4oorMBgMrFy58px5RmMBx469h6dnFIHNz7+b7epg4IGQAP7qG8lz4YFsys5n6LbDPHQgkfjCihN0Q9oHcFlbPz5aeYSsAgtPMPZ7GLxC4I8XwFRORlEGM/fOZHDIYPoE9bGsrWo0uewyPC6/jIxp0zBmZmIqN7Hxh1g8/VyJGhJiUVvFebkkrvmdzl5u3HDDDRQVFfHtt98y851p7Pl4NSUJuTQb1JNWoRPINK0ip3CXRe1f5tuUK/08+SjxFOmlZZSbNG8sOUiwtxv3Dgy3qK1qebaAAY/DgZ8hcTMAU3dOxWgy8kTPc7ewL0SbTJx68y0cg4JodvfdABTll7Lxx1jmvryZAxtSiBzQgrFv9GPA6Da4NTn/F6d/qzC6XDGCmOVLOZ2cdM78vA3JlGeX4HVta5RBWZRzxYoVHDlyBKWqf11ZWRYJCdMxmYwWtWsvpKBbWWFuDuu+/YqvnpjI8hkfk5uRfsZ85WjA6+pwjGlFFGw9aVHbnp6eDBgwgAMHDpCQkHDGvMTjMygpPUXbti+izNjNbuLowBNhgWztG8kjrQJYlp7DwL8O8vShJFJKynjpmkjyS4x8uPKIRRlxcoPhr8HJvbB7Hp/s+oSS8hKejn7asnYuoPlzz2EqLCTjk084sDGVzJQC+o9ug4OTZR/njQu/wVhawqCx99I5tAN3triS/sb2ZBXmsMRlJ6vDY8lv70BY6wdwdg7g6NHJaG3Zl/ArES0oNpl4O+4kP+xI4mBqLv++qgOuTrXfUzlD/0ngGQx/PM/+9H0sPraYuyLvIqSpZV9uOT//QvGBAwQ8+SRGnPjr1zjmvrSZmJXHadMzgDte68ugO9rj4e1iVnsDbr0LR2cX1n3z5RnTy/NKyVuThGtkM1wjvC3KePz4cfbv38+AAQPw8jrz3IPRmEdc/FQ2bhrMsbh3ycnZaVHb9kIK+nlok+bYrjS0ybLeHu6eXtzz0Qyihl/N/j9X8eXjE1gzZyaFOdn/LOPa0ReXCC9yVyZiKrSsx0f//v3x9PTkjz/+wGSqKC7FxSkcPz6T5gHX4u3V06L2fJwceTGiBVv7dmR8Cz8Wnsyk35aDfJOXw43RLflm63Fi0yw8AdbpJgjpw+G1b7Do6CJu73A7oZ6hlrVxAS4REfiMGUPa94vZ+tNRWrT1pnW3am8PfV7pxxPYu2o5PYZeh8OOUk6+u52SmNP07tObxx5/jOHDh5NyMoUZM2bw449LCPCfQG5uDKdO/WrReiLcXbk32J9vk9J58/fD9Az14dquQRa1cUHO7jDsP+iUXbyz9ll8XH2Y0GWCRU2YCgpI++B9nKN6cMy5M1+/tIltvyXQqqMvY17uw7DxkXj5u1nUpruXN31vuo24ndtIiPn/4przRwK63ITX1ZbtoZhMJn7//XeaNm3KgAED/pleXl5M4vGZbNo8hPj4j/D1HUCf3kvx8eltUfsAKzJyKLRwr7mhuWQLutaaspTzn+SKj8ng98/3cWTbKYvbbuLjyxX3TOSeDz+nw4BB7Fr2K19Muo8NC+ZSXJCPUgqvayMwFRnJXXXcoradnZ0ZNmwYqampxMRU9MqIPVbRVTAi4lmLs/4twMWJye1asrFvR0YH+jDrRAaLPE04OBp45dcDljWmFHrE/3jHTeOuHXig6wO1znU+fo88zPHwkRQXljPgljbn3QWvjtaa9bNn081/KBGJHcnfkoJHz+YEPtML7+sicPXxYMCAATz22GMMGjSIY8eOsWBBIuXGlhw5+hbl5UU1r6SKJ8Ka43Ysj+yCUl6+pqNFWc3S+WYWe4awsyCJh7vcb3EvorTPZ5HoFMna4PvY/FMczcM8ueX5aEY+0AXfFh61jtX9qlF4NQ/kz6+/wFReTmlKPoU7TtGkXwuc/Cz7gti7dy8pKSkMGzYMZ2dnTKZSTpz4lk2bhxAb+xZNm3amV/RPdO0yjSZN2lmc9WB+EeP2xjM18fx/75kFpRQ18JO4l2xBL9i4idjhI0h56aVqC3t4lB8BoU3Z/NMxykpq90v0CmjOyAcfZ9x7nxLeoxdbf/qOLybdy9afFqJ8HPDoFUj+5lTK0s9/AUp1OnfuTHBwMKtWrSI94y9OnVpMq5B7cXMLrlXOqkJcnfmgQyvW9enAiBY+FIR7sOloBvestqyoG4OjcDWF0LugE14uteiaVwNHHx/8WhbS1rgR/xbuFr121exv6ZY/kLYePXDr1IzmT0bjc1NbHM86nODq6sqQIUN47LHH6NevPwcPdqGsLI3l3/wLbTJ/S87byZGuRQaC/DyIbGn99wKDgaSydvQrUNzUcphFL01cto/f9gRxpN1teLfw5ManenDdpG4EhHpedCxHJycG3XkPp08cZ++qP8hZEofBzRHPK8w/uQxQWlrKypUradGiBZ07R5KauojNW0Zw+MgruLm1okf3+XTvNhtPz661yqm15pXYZDwdHXgg5Nw9vdziMt5ffpjL3l7NV5via7WO+nLJFnTXDu3xufMOcn9ZzLErR3Jy8v8wZmT8M18ZFANuaUtBdgm7Vli2FX22ZsEhXPf4c4x9eyrB7SPZsOBrvnj0PhIcDqGcDOQstexDYjAYGDlyJPn5eezZ8yLOzv6Ehk68qIxna+Puyuedwvj22q4Y/Fz5/XQu78afNPuCIyeDExkOkznl+LBVc1XVt9dhRgzcinI0b+AtkzYxdedUfjm5mFSScJ8QQbMxHWrcWvTw8GDEiBHcPfZ1ilPboI+kkPzY45iKzN9SH2fy4FnlibOhbv7k/IvbcUX61Th6Bpr9msI96ZQuPYSbqYgR1wZw41M9aNHW26q52vTuR8vIzhz7eQMlcTl4Dg/F4GbZQGkbN24kLy+XAQPd+GvbtRw4+AxOTp50i/qSnj0W1OrwSlUrTueyPiufp8MD8XH6/2xFpeV89ucxLnt7DVNXxzK4fQAjIptf1Lrqmt0UdK01hzMPm728o58fgS+8QMQfv+N1w/VkzZtH7PARpL3/AeU5FReDtGjjTUSPAHYtTyQ/q5ZX9FURENaaG597ldvfmEKz4BBWz/ucg7lbKD6YSc7eVIvaCgkJoWPHRCCOoKCHcHSs/a7xhQz29+Lgk0O4uXMQ7yac5KEDief0q24MioxFPLP2GWbunYnvoI70fvMumrW2bI/G09eXq+/4nT5t7ydv5UoS7xpL2SkLL8BqALTW5K46Tua8Q3hFtuW2aWNoe21n6x8KoqK77aA77qWjW1/yVT4evS07f5Cdncn+/d/Rr/8aTp78DwBdOn9Kr+hfaNZs0EVnLjWZeC02hTbuLoxr4QdAibGcOZsSuHzKGt7+/RA9WnmzZNJAPr2zB20CmprfeGY8lNTvBVh2U9DXJ6/n5l9vZtKqSRYVdqcWLQh64w0ifltC06FDOT1jBrHDhpMxfTqmggL63xSBNsGWn613T4wW7Tpyyyv/4+YX/0uaWwr5Zdkkzd5C+knLjtcHtcgnPz+MFkGjrZatOi4GAx91aMWLrYP4KS2b0btjSS+1/PJ9W0kvTOee3+9hReIKno5+mv/0+w9ODk61akspRbPx42n56aeUxMeTcOutFO3fb+XEdUeXmcj67jC5KxJx7x6A//1dcLhAF0RraN46gkSXFPY3TUQ5mF+AjcZ81q67gY6RK3B3NxHZcQp9+ywlIGCk1b58ZidncKyohFcjWqC0ZuH2JIa+u5ZXF+8n3M+DHyb246u7e1t2RW9uCix5Aj6Jhi2fWSWnuRpfQTeZKm5OdJbo5tFM6j6JHad2cMuvt/Ds2mdJyEkwu1nnsDCC33uX8F9+xr1XL9I//IjY4SMw/raQroODOLz1JKcSqr9vRm0opQjt2o3bJ0+hsKsDcSQxe87XHD9u/uEdDw9nWrZsjrt73WydV6WUYlJoc77oFMaB/CJGbj9yxgVJDdWhzEPc/tvtHMs5xkdDPmJcp3FWKQZNhw4hbN63YDCQeNdY8qq5NqChKc8vJX3mHgp3p+M5IhSfW9tZfGFPbSgHA6ltXcgNMq/LI0BRUTI7dtyKNpVSbryCAf1XExR0E0pZqbsnkFlm5L2EUwzybkJxcgEjPlzHsz/soVkTZ76+pzffTehLdJiv+Q0WnIY/XoSp3WHnXOg5HrqPrXZRk4XdX83V+Ap63Bp4rwMs/Bds/hRO7ABjKe5O7kzoOoFlo5dxb5d7+fPEn9zwyw28uulVUvPNP5zh2r49IdM+JWzBfFzat+PUm2/hM+1RXJ1NbFh4xPKbVtVAKUWPcVfT//k7cHFzYc6cOezZY/kNlerLtQHe/NyjLUatuXbnUVZknHuvkoZizfE1/GvZvwD4+qqvGdJqiFXbd+3QgfCF3+HSti0nJj3K6S++sPrnw1rKThaQ9uluSlMK8L2zA55DW9XJIRZryMnZxbbtN1JcksLll7/HiBEzMBisvxcxJS6V/NQC0v9MZtL8XTgaFNPv6skvDw/g8nb+5r8/xTmwejJ81BW2TKvotjtpO1zzHngGYdImYrNi+eHID7y44UWuWXQNX+37yur/H4Aaz04opVyBdYBL5fI/aK1fPWsZF+BroCdwGrhNa51g9bQArt4QfjkkbYEDv1RMc3SD4B4Q0gevkD481mEsd3a8ky/2fsHCwwv59div3Nr+Vu7rch9+bn5mrcatWzdCv/qKgi1bSP/gQ8L2zudQ6Z3sevdXejxj2Q2PzOHn58d9993Hd999x6JFi8jIyGDw4MEY6ugk2sWIaurOsp7tGLc3nnF743mtTTD3tfRrMAVCa82c/XN4f8f7dGrWialDp+Lvblk/dXM5+vsT+vUcUp5/nrR336MkLp6g/7yKcq7bwxiWKDqcSea8QyhnAwEPdMU5xILjwPXs5KlfOXjwWVxcAonqOh8Pj4g6Wc9ne5L4dukhnLJLMfq688FtUYyKCsbBkitXSwvhr89hw4dQnA2R18OQFyn0DmH/6f3s2rOc3Wm7iUmPIbe0Yu/e19WXKP8owjzD6uK/VXNBB0qAoVrrfKWUE7BBKbVMa72lyjL3Alla6zZKqTHA28BtdZAXWvaEljMrnuemQtLWisfxLbBpKlRe8uvn155/h/RmXLv7mF5whAWHFrDo6CLu7Hgn4zuNN7srnUffvrgvmI/38tUkzztFakw+eRuSaTKghdULmLu7O2PHjuW3335j3bp1nD59mhtuuAEnp9od761LLVyd+blHGx45cJyXY5M5WljM5LYtcbLwUm5rKysvY/LWyfx49EeGhw5n8sDJuDla1ufZUgZXV4Lfe4+M8NZkTJtG2fHjBH88FUcfnzpdrznyNyaTvSQOp0APmo3rdE7XzIZCa018wsfEx3+Et1cvunSZhrOzBYc7LLD6dC5vHUhGFRl5YGQ7nr4sAicHCzacjCWwYw6sfxfyT3EqYjC7Ol1NjDGbXVtf5XDmYYy6og5FeEUwPHQ43QK60T2gO62a1u2eUY0FXVfsQ+ZX/uhU+Th7v/J64D+Vz38APlFKKV3X+5+eQdDphooHVHxjpuysKO5Jf8HBXwnaNZfXgLub+jEtIIgv9n7Bdwe+ZXyn8dzVeTzuTjX3YVZK4X3lFdzar5i8Hw6RsyQOY3oh3qMiUJZ8EMzg6OjIqFGj8PPzY8WKFWRnZzNmzBiaNm14W1UeDg7M6hzG5LhUPj2eRkJRCTM7heHlZFm3NGvJKcnhyT+f5K+TfzGh6wQe7vawxXcbrC1lMOD/6CScw8NIfeFFEm4bQ8j0z3Bpbf7dH61Jl2uyfz1GwZZUXDv64jumAwYX6x1/tqby8hIOHnqOU6d+JSjwJjp0+C8GQ9188XxReZ/6jq29+fiqLkR6WXD+qdyIMWYeRzdMYXdZJrv8WrK7ZRCpJXGw9xNcHVzp7NeZuzvfTbeAbkT5R9XJNRgXYtZfnqo4E7EDaAN8qrXeetYiwUASgNbaqJTKAZoBGWe1MwGYANCqlWUXF5jF2R3CBlY8oOIEasYRSNpKWNJW3jm+hXvzU/nE24uP93zG3F3TuMMUwW1Xf4pvi5Y1Nu/q6YrL+Chy/0ggb+0JjKeLaXZHBwzu1t2CVkoxYMAAmjVrxo8//sjMmTO54447CAw0v49xfTEoxcsRLWjj7sKzh09w7c6jzO3amjC3+t0STMhJ4JHVj5CSn8L/Bv6P6yKuq9f1/83ruutwCm7JiUceIWHM7bT86EM8+vWr1wymIiOn5x2k5Gg2TS4PxmtkuMU3waovJaUZ7NkzkdzcXUREPEtoqwl1sgVrNGleik1mdnIGV/p5Mq1jKB6O5n3BpeUVM2/5O+zMWMhhZ02hlwHwxd/Ng24B3Rgb0J3uAd1p79seJ4Nt96bNKuha63Kgm1LKG/hJKdVZa73P0pVprWcAMwCio6Pr/uyRwQABHSoePccB0L4gg4+TtrL76G+8n7yWE5lxfPXERALCIugw4HLa978MT7+A8zapDAqvq8Jx9Hcj66dY0j6LwW9cJxwtvJTZHB06dOCee+5h3rx5fPnll4wePZr27S/iHtp16PagZoS6unDvvniu3nGEWZ2tcDdBM/2V+hdP/PkEDsqBWVfOontA93pbd3Xce3QnbOFCTjw4keP33U/gyy8Dbetl3cbTRWTM2Y8xoxif0W3x6NXwNgL+lp9/mJiY+ygty6RL52kEBFxZJ+vJKTPywP5E/szK46GQAF6MCMKhhi+NnMIylu1LZXFMClviTtO16QlM/jC8WW/6dhpN9+bdaeFh/cOuF8uifWOtdbZSag0wEqha0JOBEOCEUsoR8KLi5GjD4+EHHa6hW4dr+BrIPn2KY1u2cGjTWtZ9+xXrvv2KFu0jK4p734G4e3lX30x0II6+rpz+5iBp03bT7K5IXFpbf/cqKCiI+++/nwULFjB//nxGjBhBv379GtwHCaC/TxOW9mzH2L1x3Lr7GG1KywhwrrstFo1mEQX8d8UDhHqG8vEVH1t8l8G64twymND580l+8klO/uc/lF75Js4hdZhNQ1OTD2mf7kZr8Lu3s8V3M6xPGRlr2Lf/MRwdmtCzxwI8PbvUyXoSi0q4a08c8UUlvN8+hDtaNDvvsgUlRlYePMXi3SmsO5pOWbkmrJk7jwxpw3Vd36VtgEfF/ecbMHN6ufgDZZXF3A0YTsVJz6oWA+OAzcDNwOo6P35uJd7NmtPzmuvpec31ZJ9M5fDm9RzauJbVX05nzewZtOocRYf+l9Omdz9cPc686ZFLa28CHupGxpz9pM/ai8+NbfGItv6lwZ6enowfP56ff/6Z5cuXc/r0aa6++mqrr8cawt1dWNKjLfftS+CvwjTSSuvuvtSvF2fwg3sR/QP78+7gd2nq3LDOMzg0aULItGmcevsdjHtPonMzMBV3wuB6kQNbVMOt2ItOujemJorm90RZfPOr+qK1JunEbI4e/R9Nm3aka9cZuLrUzV7E1ux87t4Xj0nDd1ERDPA59/NRYixn7eF0FseksOpgGkVl5QR5uTK+fxijooLpHOzZIDeezsecLfQgYE7lcXQDsFBrvUQp9TqwXWu9GJgFzFVKxQKZwJg6S1yHvAOD6HPjrfS58VYyjidwaNM6Dm1axx/TP2LlF58S3j2a9v0vJ6Jnb5xcKv4oHf3cCHgwitPzDpH1wxGM6YV4Xhlm9WOWzs7O3HzzzaxZs4b169eTmZlJZKdyDA1wi8HbyZH5URFEbUslr7yc5OJSgl2t240vKbOQzAIvLndoykfDPsXRYJsTsTVRjo4EvvgCAU9PpzwlFpW4Ftpb99BCWXkZhx0S8DG1oM9Dt+LUpGEWcyjn8JFXSE6eh7//CDpFvoeDg2U3VjPXwpOZPH0oiRBXZ+Z2bU1r9/8/p2MsN7E57jSLd6fw+/6T5BUb8fVwZnTPYEZFBRMd6oOhgZ5zqImy1YZ0dHS03r59u03WbQmtNSePHeHQxnUc3ryegqxMnFxciYjuQ4cBlxMW1QMHRyd0uYnsxcco2HoS107N8L2tPQbnuim2MTExLF68mKhuKwjw96Nv30V1sp6LlVRcysCtB7nW35tPI613P3SAh+ftZNXBU6x5ejBBXg21gFVhLIVpfSt22R/cBLW89UB1vt7/NVO2T+GzYZ8xMHig1dq1ptlffoJf859xc4snNHQiEa2fMmsgFkuZtObt+JN8lHiKAd5N+KJzGD5OjphMmp3Hs1gck8LSvalk5JfS1MWREZ0CGdWtBf0jmlnWddGGlFI7tNbVji/YMDdrGhClFEFt2hPUpj2Dxt5D8sH9HNq0jiNbN3Fo41pcPZrQtk9/OgwYRPCoTjj6u5PzWxzp0ytOljp4Wb+3R1RUFN7e3uzctYyU1BROJB+gZfCFx3i0hRBXZyaGBPBR4inuDfajhyVdxC5gW0Imv+1J5bEr2jaOYg7g6Awj/gsLboftX0If69wjPqs4i+kx0xkQPKDBFvOyoiwCms3E2aWQjh3fpkXQzXWynsJyE5MOJvJbeg53BvnyZtuWHDmZx/SYFJbsSSU5uwgXRwPDOjbnuqgWDG7vb72RoxoI2UKvpXKjkcS9uzi8cR1Ht22hrLgID28f2vUbSPuQvrC2AOXqgN+/InFuWTfHdnftepy4uG20b9+fzp2n1Mk6Lla+sZz+Ww8S4urMkh5tL/p4pMmkuWHaRtJyS1j99CDcnRvRNonW8PX1cHIPTNoJ7hd/4czkLZP5/sj3/DjqRyK86+aqyouV+fVc1p+aiUPgIK4d+0adrONkSRn/2hvH3rwiHvH1wT2tmF/3pBKXXoCjQXF5O39GRbVgWGRzmrg0os9MNS60hS4F3QrKSkuI37mNQ5vWEbdzG+VlZQQ370Bfr6tx1M743tYe9y51c+l5bOzbJB6fQa/on2p9g/+6Ni/1NE8eSuKzyFBubH5xV0/+uOMET30fw/u3RnFTj5qvHWhwTu6Dzy+D3g/AVW9dVFPHso8xevFobml3Cy/2fdFKAa2rPDub2CtH4tYpkpBZs+rkBOPevELu2nqU3KQ8WmYZSUorQCnoG96MUd1aMLJTID4eDedWDBdLCno9KiksJHbbZg5vWsfJfUcY4H8DzVxbkNAkkeinb8bVyj0cjMY8Nm2+Anf3cHr2WNAgz8iXa83I7UfILDOyoU9H3Gp5rLKw1MiQd/8k0NOVnx4a0GhPXPHrY7DrG3hoC/jVvn/6xJUT2ZO+h99u/A0fV9vfZqA6Jyf/j6xvvyX8p59wbW/50HAXYjJpnvg+hp/j01HZpQB0C/FmVFQLrukaRHNP6/cmagguVNAbx1mARsTF3Z1Og67gpudfY/xn02FkU04Y43HLdOSTqZ+wdetWjEbrdeVzdGxKROsnyMnZTlr6Mqu1a00OSvFam2CSS8qYnlT7ASGmr43jVG4JL18b2XiLOcCQFytuKLf8pVo3sf7EejYmb2Ri14kNtpiXxMWRNW8e3rfcYtVirrXm2M40Fry+lQ3HTuNYamLi0AjWPzuEnx8ewD0Dw+22mNdEttDrgdaaxCPHWLNpPYmJiXh5eTFo0CCioqJwcLj4kzJal/PXtlEYjfn07bMcB4eGeQOme/fFsyYzj019OhLoYlkvj5TsIoa+9yfDOjbnkzt61FHCerThQ1j5Koz9CSKGWvTSMlMZNy++mXJdzk+jfqr1YB11LemBiRTu2EHEH7/j2Oz8F/SYS2vN8QOZbP0ljvTjefgEutPl6lAiugfgbmcnNy9EttBtTClFWPs2jB8/nrFjx+Lh4cHixYuZNm0a+/btw2TBgMPVt+9A2zYvUFx8gqSkurnPsjW8HNECo0nzZpxlw+0BvPP7IUwa/n1VhzpIZgN9HwTv0IoBEcot22P7/vD3xOXE8VTPpxpsMc/fsJH8tWvxe3CiVYp5ytFsfnpvJ0s+jqG4oIwrxnVkzMu96dI76JIq5jVp3Kd7GxmlFBEREbRu3ZpDhw6xevVqfvjhB5o3b87QoUNp165drY+B+/oOwM9vGAmJnxEUNBoXl7o5CXsxwtxcuK+lP58lpXFPSz+impp3Ucmu41n8vDuFh4dE0NKnbi5EqXeOLjDijYqBWnZ9DdH3mPWynJIcpsVMo09QHwaHDK7bjLWkjUbS3n4Lp5AQfMZWP2KPudISc9n6SxzHD2Ti7unM5WPaETmwBQ71MNJSYyTvig0opejYsSMPPvggN910E6WlpcyfP59Zs2YRHx9f63bbtvk3JlMxcXHvWzGtdT0e1hxfJ0dePZps1ug+WmveWHIA/6YuPDi4TT0krEcdR0HogIrRborNG/lpesx08krzeCb6mQZ5Ahwg+/vvKTkaS8AzT2Oo5UAfmSkFLPt8L9+/uZ1Tibn0uymCu/7bjy6DW0oxvwB5Z2zIYDDQtWtXHnnkEa699lpycnKYM2cOc+bM4cSJExa35+4eTsuW/yIl9Xvy8g7UQeKL5+nowHPhgWzJKeC39JqL2K97Utl5PJtnRrRv9P2Hz6EUXDkZCk/DundrXDwhJ4EFhxZwU9ubaO/bMO+6WZ6bS/rUj3Hv1Yumw4db/Pqc9CJWfnWA+W9sJelgJr2uCWPsf/vTY0QoTnV05bU9kZOiDUhZWRnbt29n/fr1FBYW0r59e4YMGWLRfdDLynLYvOUKmni0p3v3bxrkVpzRpBm2/TCF5SbW9e6A63m6MRaXlTP03T/xdnfm10kDLRserDH5+SHYsxAe+Qt8zz8gxqRVk9h2ahtLblxi9lCK9e3U2++QOXs24T/+gGuk+Vcv52eVsH1ZAgc3pKAcFF0Gt6THla1wa2I//cetRU6KNhJOTk7069ePxx57jCFDhpCQkMD06dP54YcfOH3avLsROzl5ER7+GFnZW8jIaJgj0TsaKroxHi8u5YsT6edd7ov1caTkFPPytZH2W8wBhr4MDs6w4pXzLrI5ZTN/nviTCV0nNNhiXpqYSOY33+B1041mF/Oi/FI2/nCUb17ZzMENKUQObMHYN/oxYHQbKea1IFvoDVhhYSGbNm36p+969+7dGTRoEF5eF77vuslkZOtf16B1GX37/F4nI6Zbw9g9cWzOzmdz3474n3Xf9FO5xQx5908ua+vH52Or3RixL2unwJr/wrglEH7ZGbPKTeXcsuQWCssKWXzDYpwdGubvM+mRRyjYtJmI35fhFHD+QWIASoqM7F55nJiVSRhLy2nXJ5Be14Tj5d9I7s1jQ7KF3ki5u7szbNgwHn30UXr16kVMTAxTp05l2bJl5Ofnn/d1BoMj7dq+QFFRIidOzK3HxJZ5tU0Lik0m3ok/ec68d/84TFm5iReu7miDZDbQ/xHwbAl/vACm8jNmLYpdxNGsozwV/VSDLeYFW7aSv3IVfhMmXLCYl5WWs/OPROa+tIntvyXQKtKXMS/3Ydj4SCnmViBb6I1IdnY2a9euZffu3Tg6OtK3b1/69++Pm1v1fwi7d99NTu4u+vVdhbPzxfcFrgsvHz3BrBMZrOzVnsjK+3jvS87huk82cP9lrS+dgg6w9wf48V4Y9Qn0qOjul1eax7U/XUu4VzhfXflVgzwnosvLiR99M+W5OUQsXVrtAB7lRhMHNqSwfWkChbmltOrkS59RrQkI9bRB4sZNttDthLe3N9dffz0PP/ww7dq1Y/369Xz00UesW7eOkpKSc5Zv0/YFyssLiYufaoO05nkyLBAvRwdeja3oxqi15vUlB/B1d+aRoXbWTbEmnUdDy16w+g0oyQNg5t6ZZBVn8UyvhttNMeennyg5dIjmTz99TjE3lZs4uCmVb1/ZwroFR/AKcOPGp3pw3aRuUszrgJ31A7s0+Pn5ccsttzBw4EDWrFnD6tWr2bp1K5dddhk9e/bEyanieHQTj7YEt7iDE8nf0jL4Tpo0se7NkazBx8mRp8IDeeloMstP52JKLeSv+Ez+e0NnPF0b5lWQdUYpGPkWfHEFbPiApN53882BbxgVMYpOzTrZOl21yvMLSPvwI9y6d6fpVVf9M12bNMd2pbN1cRzZpwrxb9WUwXe2JyTSt8F+MdkDOeRiB5KSkli1ahUJCQl4enoyaNAgunXrhoODA6WlmWzecgWenlF0i2qYu+xlJs3QbYcwGjUOG07i4ezIb48OxLGRjCBjdT/eDwd+4ck+N7EhfSdLblxCgPuFTzLaStr7H3B6xgzCFn6HW9euFfct2nearYvjyEjKxyfQnT7Xt6Z1N/8G+dlrjGTEIjsXEhLC+PHjiYuLY9WqVfz6669s3LiRIUOG0KlTJ8LDJnE0djKnT/+Jn98QW8c9h5NB8WqbYMb9EoNTVhFz7+196RZzgGGvsvXob6xI3cgj3R5psMW89EQymbNn4znqOty6diX5SBZbfo7jZFwOnn6uDBvfkba9Axv3nTEbGSnodqR169aEh4dz+PBhVq9ezY8//sj69esZMuQy3NzCOBr7Jr6+AzEYGt6hjKE+TWiaWozycWFgm4bZz7reeLVksVtLAo35jGt9g63TnFfae++CwQC3Pcjij3aRdDALDy9nBt3Rno79g+QSfRuQgm5nlFJ06NCBdu3asX//ftasWcN33/1I27Y9CQz6keTkeYSEjLN1zHMYDAZecPdFGbXsmgPdsy6jfU4prk2b2zpKtQp37CB1XQzJw18leWY8rh5O9B/dhi6DgnGUS/RtpsaCrpQKAb4GmgMamKG1/uisZQYDvwB/31lqkdb6dasmFRYxGAx06dKFyMhIYmJi+PPPNbi4BnJw/1uonA607NzH1hHP0dSkQHbPAShXTrjohrcnBVCQUcgf720hNfpFnMsc6X1dK6KGhuDsJtuHtmbOb8AIPKW13qmUagrsUEqt0Fqfffen9Vrra60fUVwMBwcHevToQZcuXdi02JGC8i9Ieftpmj39CW5RUbaOJxoZY3YJWTO2kOvkT7tmeVz2wjW4NmmYXzyXohoPcmmtU7XWOyuf5wEHgeC6Diasy8nJiUGjH2Nwm4X4prqS+K9x5C5rmEPWiYapNCmPtE93oUpcufmZwQz77ygp5g2MRWctlFJhQHdgazWz+ymlYpRSy5RS1XaaVUpNUEptV0ptT08//02ZRN1xb9+RsIXf4dq5M8lPPEn6tGlm3ZdcXNoK96ST9vkelKOBgAejaNI5AGWQk54Njdm/EaVUE+BH4HGtde5Zs3cCoVrrKOBj4Ofq2tBaz9BaR2uto/39G96IOpcKR19fWn31JV7XX0/G1I9JeeZZTNVcaSqE1prc1cfJnHcI5+AmBDzcDadAD1vHEudhVkFXSjlRUcy/1VovOnu+1jpXa51f+Xwp4KSUusT7njVsBmdngt56E/8nniB3yRKOjxuP0cxb9IpLgzaayFp4hNzlibh388f/vi44yC1tG7QaC7qq6EM2Cziota52bDOlVGDlciilele2K9WhgVNK4ffABII//JDiQ4dIuPU2io8csXUs0QCU55eSPnMvhbvS8Bweis9t7VFOcoiloTPnNzQAGAsMVUrtrnxcrZSaqJSaWLnMzcA+pVQMMBUYo+XAbKPhOfJKQufORZeWknj7HeSvW2frSMKGyk4VkDYthtLkfHzv6IDnFa3k2oBGosZui1rrDcAFf5ta60+AT6wVStQ/ty6dCft+IUkPPkTSxAdp/vzz+Nx1p/whX2KKj2Rx+tuDKCcDAQ90xTmkqa0jCQvIPpT4h1NgIGHfzKXJkCGcmjyZU2+8gTYabR1L1JP8zSlkfLUPRx9XAh7pJsW8EZKCLs5g8PCg5cdTaXbfvWTNm0/SAxMpzz27U5OwJ7pck/VLLNm/HMO1gy/+D3bF0fvcQSpEwycFXZxDGQwEPP00QZP/S8HWrSTcfgelSUm2jiXqgKnYSMac/RRsTqXJZcE0GxuJwUUu4W+spKCL8/IePZpWs2ZRnpFBwi23Urhjh60jCSsyZhaTNi2GkthsvG9qg/c1rVFyL51GTQq6uCCPPr0J+24BDt7eHB9/N9k//2zrSMIKShJySPt0F+W5pfjd05kmvYNsHUlYgRR0USPnsDDCvluAW8+epP77edI++BBtMtk6lqilgl1ppM/ci8HNiYCHo3Bt423rSMJKpKALszh4edFq5gy8b72V059/TvLjT2AqKrJ1LGEBbdLkLE8g67vDuIR6EvBQFE7+7raOJaxIzn4IsyknJwJf+w/OrcNJe/sdEpOTaTltGk7NG+YQaeL/mUrLyfr+CEV7M3CPbo7PDW1QMqKQ3ZHfqLCIUopm48fT8tNPKY2PJ+HWWyk+cPat8UVDUp5bSvqMPRTty8Dr6nB8RreVYm6n5LcqaqXp0CGEzp8HBgMJd95F3sqVto4kqlGakk/ap7swphXSbGwkTS9vKVf/2jEp6KLWXNu3J3zhd7i0bcuJSY9yetYsubd6A1J04DTp02MA8J8YhVtkMxsnEnVNCrq4KI7+/oR+PYemI68kbcq7pL70Erq01NaxLmlaa/LWneD03AM4BrgT8HB3nFs0sXUsUQ/kpKi4aAZXV4Lfe4+M8NZkTJtG2fEkgqd+hKOPj62jXXK00UT2L8co2HYSty5++NzSDoOzg61jiXoiW+jCKpTBgP+jk2gx5R2KYmJIGDOGkrh4W8e6pJgKy8j4ch8F207SdGgIvrd3kGJ+iZGCLqzK67rraDV7Nqa8fBLGjKFg82ZbR7oklKUXVlzGn5iLz63t8BoRJpfxX4KkoAurc+/RnbCFC3FqHsDx+yeQtXChrSPZteJj2aRNi8FUVIb//V3w6NHc1pGEjUhBF3XCuWUwofPn49G/HydfeZVTb72NLi+3dSy7U/DXSTJm7cOhqTMBD3XDJczL1pGEDUlBF3XGoUkTQqZNw2fsWDJnz+bEw49Qnl9g61h2QZs02UvjyFp0FJc23gQ8FIVjMzdbxxI2JgVd1Cnl6Ejgiy8Q+Oor5K9fT+Kdd1KWkmLrWI2aqaSc03MPkL8uGY9+QfiN64TBVTqsCSnoop743H47IZ9/TllyMvG33kZRTIytIzVKxuwS0qfHUHwoE+9REfhc3wblICc/RQUp6KLeNBk4gLAF8zG4uZH4r3HkLltm60iNSmlSXsVl/JnF+I3vRJP+LWwdSTQwNRZ0pVSIUmqNUuqAUmq/UuqxapZRSqmpSqlYpdQepVSPuokrGjuXNm0IW/gdrp07k/zEk6RPmya3CzBD4Z500j7fg3I0EPBQFK7tfW0dSTRA5myhG4GntNaRQF/gYaVU5FnLXAW0rXxMAD6zakphVxx9fGj11Zd4XX89GVM/JuWZZzGVlNg6VoOktSZ39XEy5x3CObgJAQ93w6m5h61jiQaqxjMpWutUILXyeZ5S6iAQDFS9Z+r1wNe6YlNri1LKWykVVPlaIc5hcHYm6K03cY6IIP399yk7cQLd/XGUo5OtozUYShvIWniEwl1puHcPkNveihpZdGpcKRUGdAe2njUrGKg6LPyJymlnFHSl1AQqtuBp1aqVhVGFvVFK4TfhfpxDQ0l57jmKnffjEhZs61gNgqlM08kYTeGuNDxHhNJ0SIjc9lbUyOyve6VUE+BH4HGtdW5tVqa1nqG1jtZaR/v7+9emCWGHPK8cQeg3c2maG4/rwQ3osmJbR7I5P1rTFG/cbgvFc2grKebCLGYVdKWUExXF/Fut9aJqFkkGQqr83LJymhBmcevcmZHPd2N4q09Rf31u6zg2tTllMy+1/oJtV6XQrLvsyQrzmdPLRQGzgINa6/fPs9hi4F+VvV36Ajly/FxYyqnPaJx7DIO1UyA/3dZxbKLcVM6U7VPw8WrG9QNvsXUc0ciYs4U+ABgLDFVK7a58XK2UmqiUmli5zFIgDogFZgIP1U1cYfdG/BeMRbBmsq2T2MSi2EUczTrKkz2fxMXBxdZxRCNjTi+XDcAFD+BV9m552FqhxCXMvx30ug/+mgG974fmnWydqN7klebxya5P6BHQg+Ghw20dRzRC0gdKNDyDngMXT/j9ebiELjqauXcmWcVZPNv7WTkJKmpFCrpoeNx9YfDzEL8Wjvxu6zT1IikviW8OfMN1EdfRqdmls1cirEsKumiYet0LzdrC8pfAaP+DTn+w4wMcDY481uOcO2sIYTYp6KJhcnCCKyfD6VjY9oWt09SpbSe3sSJxBfd0vocA9wBbxxGNmBR00XC1HQERQ2HtW1CYaes0daLcVM6UbVMI9AhkXKdxto4jGjkp6KLhUgpGTIaSPPjzTVunqROLjy3mYOZBHu/xOG6OMuKQuDhS0EXD1jwSet4N22ZB+mFbp7GqwrJCpu6aSlf/rlwdfrWt4wg7IAVdNHxDXgDnJhUnSO3IrH2zyCjK4Nle0k1RWIcUdNHwefjBoGfg6HKIXWnrNFaRmp/KnP1zuCr8KqL8o2wdR9gJKeiiceg9AXzC4Y8Xodxo6zQX7YOdHwDwRI8nbJxE2BMp6KJxcHSpuM9L+iHY8ZWt01yU3Wm7WRa/jHGdxhHUJMjWcYQdkYIuGo8O10DYZbDmf1CUbes0tWLSJqZsm4K/mz/3dr7X1nGEnZGCLhoPpeDK/0FRFqybYus0tbI0fil7MvbwaI9HcXdyt3UcYWekoIvGJagrdL8Ltn4Op4/ZOo1FioxFfLjjQzr6dmRUxChbxxF2SAq6aHyGvlxxTH35y7ZOYpHZ+2dzqvAUz/V+DoOSPz1hffKpEo1P0+Zw2ZNw+DeIW2vrNGY5VXCKr/Z9xfDQ4fRs3tPWcYSdkoIuGqe+D4NXK/jjBTCV2zpNjabumorRZOSJntJNUdQdKeiicXJyheGvwal9sOsbW6e5oP0Z+1l8bDFjI8cS0jSk5hcIUUtS0EXj1elGCOkLq9+A4lxbp6mW1pq3t72Nr6sv93e539ZxhJ2Tgi4aL6Vg5P+gIB02vG/rNNVanricXWm7mNR9Ek2cm9g6jrBzUtBF4xbcE6Juh82fQlaCrdOcoaS8hA92fEA7n3bc2OZGW8cRl4AaC7pS6kulVJpSat955g9WSuUopXZXPl6xfkwhLuCKV8DgCCtetXWSM8w9MJfk/GSe7fUsDgYHW8cRlwBzttBnAyNrWGa91rpb5eP1i48lhAU8W8CAx+HAz5C4ydZpAMgoymDmnpkMDhlMn6A+to4jLhE1FnSt9TrAPsf/Evaj/yTwDIbfnweTydZp+GTXJ5SaSnk6+mlbRxGXEGsdQ++nlIpRSi1TSnU630JKqQlKqe1Kqe3p6elWWrUQgLM7DPsPpO6GPQtsGuVQ5iEWHV3E7R1uJ9Qz1KZZxKXFGgV9JxCqtY4CPgZ+Pt+CWusZWutorXW0v7+/FVYtRBWdb644SbrqdSjJt0kErTXvbHsHLxcvHuj6gE0yiEvXRRd0rXWu1jq/8vlSwEkp5XfRyYSwlMEAI9+CvFTY+JFNIqxOWs22k9t4qNtDeLl42SSDuHRddEFXSgWqygERlVK9K9s8fbHtClErIb2h82jYNBWyk+p11aXlpby3/T0ivCK4pd0t9bpuIcC8bovzgc1Ae6XUCaXUvUqpiUqpiZWL3AzsU0rFAFOBMVprXXeRhajBsP9U/LvqtXpd7fxD80nKS+LpXk/jaHCs13ULAVDjp05rfXsN8z8BPrFaIiEulncr6PcIrH8Xej8AIb3qfJWZxZl8HvM5A4MHMjB4YJ2vT4jqyJWiwj4NfAKaNIc/nod62GGctnsahcZCnol+ps7XJcT5SEEX9smlScUVpCe2wb4f63RVsVmxfH/ke25tfyutvVvX6bqEuBAp6MJ+Rd0BgV0rbglQVlQnq9BaM2X7FDycPHgo6qE6WYcQ5pKCLuzX390Yc0/Apro5zbM+eT2bUjYxsetEvF2962QdQphLCrqwb2EDoOMo2PAB5KZatekyUxnvbn+XUM9Qbu9wwb4DQtQLKejC/g1/HUxlFQNhWNHCwwuJz4nnqZ5P4eTgZNW2hagNKejC/vmGQ98HYfc8SNlllSZzSnL4LOYz+gT1YXDIYKu0KcTFkoIuLg2XPQ3uzeD3F6zSjXF6zHTySvN4JvoZKi+UFsLmpKCLS4OrJwx9CY5vgoOLL6qp+Jx4FhxawE1tb6K9b3srBRTi4klBF5eOHv+CgE6w/GUoK651M+9tfw8XRxce7vawFcMJcfGkoItLh8GhYlDp7ETYOr1WTWxK2cTaE2uZ0HUCfm5yU1HRsEhBF5eW1oOh3VWw7l3IT7PopUaTkSnbphDcJJi7Ot5VN/mEuAhS0MWlZ8R/wVgEayZb9LJFRxcRmx3LU9FP4ezgXEfhhKg9Keji0uPXBnpPgJ1fw8l9Zr0krzSPT3d/Ss/mPRnWalgdBxSidqSgi0vToGfB1Qv+MK8b48w9M8kqzuKZXtJNUTRcUtDFpcnNBwa/APFr4fCyCy6alJvE3INzGRUxik7NzjsGuhA2JwVdXLqi7wa/9rD8JTCWnnex93e8j5PBiUd7PFqP4YSwnBR0celycIIrJ0PmMdg2s9pFtp3cxsrjK7m3870EuAfUc0AhLCMFXVza2g6HNsNg7dtQcObY5uWmcqZsm0KgRyDjOo2zUUAhzCcFXYgRk6EkH/5884zJi48t5mDmQZ7o8QSujq42CieE+aSgCxHQAaLvge1fQtohAArKCpi6aypd/btyVfhVNg4ohHlqLOhKqS+VUmlKqWo77KoKU5VSsUqpPUqpHtaPKUQdG/w8ODeB5S8CMGvvLDKKMniu13PSTVE0GuZsoc8GRl5g/lVA28rHBOCzi48lRD3zaAaDn4PYlaTs/Y45++dwTetr6Orf1dbJhDBbjQVda70OyLzAItcDX+sKWwBvpVSQtQIKUW963Q++EXyw9U0MysDjPR63dSIhLGKNY+jBQFKVn09UThOicXF0Zn//ifzuVE7/sjYEegTaOpEQFqnXk6JKqQlKqe1Kqe3p6en1uWohzNKh+z0Mze9Feol0UxSNjzUKejIQUuXnlpXTzqG1nqG1jtZaR/v7+1th1UJYl4ODIycd76dcedk6ihAWs0ZBXwz8q7K3S18gR2udaoV2hRBCWMCxpgWUUvOBwYCfUuoE8CrgBKC1ng4sBa4GYoFC4O66CiuEEOL8aizoWuvba5ivARlcUQghbEyuFBVCCDshBV0IIeyEFHQhhLATUtCFEMJOSEEXQgg7IQVdCCHshBR0IYSwE1LQhRDCTkhBF0IIOyEFXQgh7IQUdCGEsBNS0IUQwk5IQRdCCDshBV0IIeyEFHQhhLATUtCFEMJOSEEXQgg7IQVdCCHshBR0IYSwE1LQhRDCTkhBF0IIOyEFXQgh7IRZBV0pNVIpdVgpFauU+nc188crpdKVUrsrH/dZP6oQQogLcaxpAaWUA/ApMBw4AWxTSi3WWh84a9HvtNaP1EFGIYQQZjBnC703EKu1jtNalwILgOvrNpYQQghLmVPQg4GkKj+fqJx2ttFKqT1KqR+UUiHVNaSUmqCU2q6U2p6enl6LuEIIIc7HWidFfwXCtNZdgRXAnOoW0lrP0FpHa62j/f39rbRqIYQQYF5BTwaqbnG3rJz2D631aa11SeWPXwA9rRNPCCGEucwp6NuAtkqpcKWUMzAGWFx1AaVUUJUfRwEHrRdRCCGEOWrs5aK1NiqlHgH+AByAL7XW+5VSrwPbtdaLgUeVUqMAI5AJjK/DzEIIIapRY0EH0FovBZaeNe2VKs+fB563bjQhhBCWkCtFhRDCTkhBF0IIOyEFXQgh7IQUdCGEsBNS0IUQwk5IQRdCCDshBV0IIeyEFHQhhLATUtCFEMJOSEEXQgg7IQVdCCHshBR0IYSwE1LQhRDCTkhBF0IIOyEFXQgh7IQUdCGEsBNS0IUQwk5IQRdCCDshBV0IIeyEFHQhhLATUtCFEMJOmFXQlVIjlVKHlVKxSql/VzPfRSn1XeX8rUqpMKsnFUIIcUE1FnSllAPwKXAVEAncrpSKPGuxe4EsrXUb4APgbWsHFUIIcWHmbKH3BmK11nFa61JgAXD9WctcD8ypfP4DcIVSSlkvphBCiJo4mrFMMJBU5ecTQJ/zLaO1NiqlcoBmQEbVhZRSE4AJAK1ataplZCHqVmQLT1tHEKJWzCnoVqO1ngHMAIiOjtb1uW4hzPXqdZ1sHUGIWjHnkEsyEFLl55aV06pdRinlCHgBp60RUAghhHnMKejbgLZKqXCllDMwBlh81jKLgXGVz28GVmutZQtcCCHqUY2HXCqPiT8C/AE4AF9qrfcrpV4HtmutFwOzgLlKqVggk4qiL4QQoh6ZdQxda70UWHrWtFeqPC8GbrFuNCGEEJaQK0WFEMJOSEEXQgg7IQVdCCHshBR0IYSwE8pWvQuVUulAok1Wbj4/zrratYGSnNbXWLJKTutqDDlDtdb+1c2wWUFvDJRS27XW0bbOURPJaX2NJavktK7GkvN85JCLEELYCSnoQghhJ6SgX9gMWwcwk+S0vsaSVXJaV2PJWS05hi6EEHZCttCFEMJOSEEXQgg7cckXdKWUr1JqhVLqaOW/PtUsM0QptbvKo1gpdUPlvNlKqfgq87rZKmflcuVVsiyuMj28cgDv2MoBvZ1tlVMp1U0ptVkptV8ptUcpdVuVeXX6fl7MgOdKqecrpx9WSl1pzVy1yPmkUupA5fu3SikVWmVetZ8BG2Ydr5RKr5LpvirzxlV+Vo4qpcad/dp6zvlBlYxHlFLZVebV63taa1rrS/oBvAP8u/L5v4G3a1jel4pbBLtX/jwbuLmh5ATyzzN9ITCm8vl04EFb5QTaAW0rn7cAUgHvun4/qbj98zGgNeAMxACRZy3zEDC98vkY4LvK55GVy7sA4ZXtONgw55Aqn8EH/855oc+ADbOOBz6p5rW+QFzlvz6Vz31slfOs5SdRcavwen9PL+ZxyW+hc+YA13OAG2pY/mZgmda6sC5DVcPSnP+oHLB7KBUDeFv8egvVmFNrfURrfbTyeQqQBlR75ZuVXcyA59cDC7TWJVrreCC2sj2b5NRar6nyGdxCxUhitmDOe3o+VwIrtNaZWussYAUwsoHkvB2YX0dZ6owUdGiutU6tfH4SaF7D8mM49xc9uXLX9wOllIvVE1YwN6erUmq7UmrL34eFqBiwO1trbaz8+QQVA3vbMicASqneVGwxHasyua7ez+oGPD/7fThjwHPg7wHPzXltfeas6l5gWZWfq/sM1BVzs46u/J3+oJT6e0jLBvmeVh6+CgdWV5lcn+9prdXrING2opRaCQRWM+vFqj9orbVS6rz9OJVSQUAXKkZv+tvzVBQuZyr6sD4HvG7DnKFa62SlVGtgtVJqLxVFyWqs/H7OBcZprU2Vk632fl4KlFJ3AdHAoCqTz/kMaK2PVd9CvfgVmK+1LlFKPUDFHtBQG+apyRjgB611eZVpDe09rdYlUdC11sPON08pdUopFaS1Tq0sMGkXaOpW4CetdVmVtv/eGi1RSn0FPG3LnFrr5Mp/45RSfwLdgR8Bb6WUY+VWZ3UDfddrTqWUJ/Ab8KLWekuVtq32flbDkgHPT6gzBzw357X1mROl1DAqvkQHaa1L/p5+ns9AXRWfGrNqrasOGP8FFedZ/n7t4LNe+6fVE/7/usz9/Y0BHq46oZ7f01qTQy5nDnA9DvjlAsuec1ytsmj9fZz6BmCf9SMCZuRUSvn8fYhCKeUHDAAO6IqzOmuoOP5/3tfXY05n4Cfga631D2fNq8v382IGPF8MjKnsBRMOtAX+smI2i3IqpboDnwOjtNZpVaZX+xmoo5zmZg2q8uMo4GDl8z+AEZWZfYARnLn3W685K7N2oOIE7eYq0+r7Pa09W5+VtfWDiuOjq4CjwErAt3J6NPBFleXCqPhGN5z1+tXAXioKzzdAE1vlBPpXZomp/PfeKq9vTUUBigW+B1xsmPMuoAzYXeXRrT7eT+Bq4AgVW1cvVk57nYrCCOBa+f7EVr5frau89sXK1x0Grqrjz2VNOVcCp6q8f4tr+gzYMOubwP7KTGuADlVee0/lex0L3G3LnJU//wd466zX1ft7WtuHXPovhBB2Qg65CCGEnZCCLoQQdkIKuhBC2Akp6EIIYSekoAshhJ2Qgi6EEHZCCroQQtiJ/wN3qyBI6iA4JwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import sin, cos\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "s = 1\n",
    "d = [[0, 1, 0]]\n",
    "plt.plot([0, 0], [0, 1])\n",
    "for i in range(5):\n",
    "    n = []\n",
    "    for j in range(len(d)):\n",
    "        n.append(\n",
    "            [\n",
    "                d[j][0] + s * sin(d[j][2] - 0.2),\n",
    "                d[j][1] + s * cos(d[j][2] - 0.2),\n",
    "                d[j][2] - 0.2,\n",
    "            ]\n",
    "        )\n",
    "        n.append(\n",
    "            [\n",
    "                d[j][0] + s * sin(d[j][2] + 0.2),\n",
    "                d[j][1] + s * cos(d[j][2] + 0.2),\n",
    "                d[j][2] + 0.2,\n",
    "            ]\n",
    "        )\n",
    "        plt.plot([d[j][0], n[-2][0]], [d[j][1], n[-2][1]])\n",
    "        plt.plot([d[j][0], n[-1][0]], [d[j][1], n[-1][1]])\n",
    "    d = n\n",
    "    s *= 0.6\n",
    "plt.savefig(\"tree.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Suggested Marking Scheme\n",
    "If you want to self-assess your solution you can consider using the marking scheme below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part one: Refactoring (15 marks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Copy the code above into a file tree.py, invoke it with python tree.py, and verify it creates an image tree.png which looks like that above.\n",
    "* Initialise your git repository with the raw state of the code. (**1 mark**)\n",
    "* Identify a number of simple refactorings which can be used to improve the code, *reducing repetition* and *improving readability*. Implement these one by one, with a git commit each time.\n",
    "    * **1 mark** for each refactoring, **1 mark** for each git commit, at least five such: **10 marks total**.\n",
    "* Do NOT introduce NumPy or other performance improvements yet (see below.)\n",
    "* Identify which variables in the code would, more sensibly, be able to be input parameters, and use Argparse to manage these.\n",
    "   * **4 marks**: 1 for each of four arguments identified."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part two: performance programming (10 marks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* For the code as refactored, prepare a figure which plots the time to produce the tree, versus number of iteration steps completed. Your code to produce this figure should run as a script, which you should call `perf_plot.py`, invoking a function imported from `tree.py`. The script should produce a figure called `perf_plot.png`. Comment on your findings in a text file, called `comments.md`. For your performance measurements you should turn off the actual plotting, and run only the mathematical calculation using an appropriate flag. **5 marks**:\n",
    "   * Time to run code identified. (**1 mark**)\n",
    "   * Figure created. (**1 mark**)\n",
    "   * Figure correctly formatted. (**1 mark**)\n",
    "   * Figure auto-generated from script. (**1 mark**)\n",
    "   * Performance law identified. (**1 mark**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The code above makes use of `append()` which is not appropriate for `NumPy`. Create a new solution (in a file called `tree_np.py`) which makes use of `NumPy`. Compare the performance (again, excluding the plotting from your measurements), and discuss in comments.md. **5 marks**:\n",
    "    * Array-operations used to subtract the change angle from all angles in a single minus sign. (**1 mark**)\n",
    "    * Array-operations used to take the sine of all angles using np.sin. (**1 mark**)\n",
    "    * Array-operations used to move on all the positions with a single vector displacement addition. (**1 mark**)\n",
    "    * Numpy solution uses `hstack` or similar to create new arrays with twice the length, by composing the left-turned array with the right-turned array. (**1 mark**)\n",
    "    * Performance comparison recorded. (**1 mark**)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}