Py_lectures_06_Problems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem: generate the first n-th element of the Fibonacci sequence\n",
    "\n",
    "https://en.wikipedia.org/wiki/Fibonacci_number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can try yourself. It is a simple but not so easy problem to solve. Here below some of the solutions suggested by various groups of students. \n",
    "\n",
    "Let's take n = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500: 139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125\n",
      "0.45919 ms\n"
     ]
    }
   ],
   "source": [
    "# Solution 1: n = 500\n",
    "n = 500\n",
    "from time import time\n",
    "t0 = time()\n",
    "serie = []\n",
    "serie.append(0)\n",
    "serie.append(1)\n",
    "for i in range(0, n-2):\n",
    "    serie.append(serie[-1]+serie[-2])\n",
    "print('{x:3}: {f:3}'.format(x = i+3, f = serie[-1]+serie[-2]))\n",
    "ts = (time() - t0)\n",
    "print(\"{0:.5f} ms\".format(ts*1000.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Used Memory: 4.272 kbytes\n"
     ]
    }
   ],
   "source": [
    "print(\"Used Memory: %.3f kbytes\" % (sys.getsizeof(serie)/1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comment\n",
    "It is a good solution, it stricty apply the definition to solve the problem. Using range can be a little problem with very large N (why generate N numbers and make a loop over it?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125\n"
     ]
    }
   ],
   "source": [
    "# Solution 2:\n",
    "l = [0,1]\n",
    "count = 0\n",
    "while count <= n:\n",
    "    f = l[-2] + l[-1]\n",
    "    l.append(f)\n",
    "    count+=1\n",
    "print(l[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comment\n",
    "Another good solution, this time using a \"while\" which solves the problem to make a list of n numbers using range.\n",
    "The problem is now: do we need to save all the other n-2 elements of the sequence in the list \"l\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Used Memory: 4.392 kbytes\n"
     ]
    }
   ],
   "source": [
    "print(\"Used Memory: %.3f kbytes\" % (sys.getsizeof(l)/1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125\n"
     ]
    }
   ],
   "source": [
    "# Solution 3:\n",
    "a, b = 0, 1\n",
    "count = 2\n",
    "while count <= n:\n",
    "    c = a + b\n",
    "    a, b = b, c\n",
    "    count += 1\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comment\n",
    "Excellent solution, no memory used!\n",
    "\n",
    "A little improvement, by the way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125\n"
     ]
    }
   ],
   "source": [
    "# Solution 3b:\n",
    "a, b = 0, 1\n",
    "count = 2\n",
    "while count <= n:\n",
    "    a, b = b, a + b\n",
    "    count += 1\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comment\n",
    "It would be nice to make a loop without setting a counter. Something giving back a number of the sequence, everytime is asked. By the way, with the property of holding the previous values... is it too much to ask? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "139423224561697880139724382870407283950070256587697307264108962948325571622863290691557658876222521294125"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution 4\n",
    "# http://stackoverflow.com/questions/494594/how-to-write-the-fibonacci-sequence-in-python\n",
    "def fib_eff():\n",
    "    a, b = 0, 1\n",
    "    while True:            # First iteration:\n",
    "        yield a            # yield 0 to start with and then\n",
    "        a, b = b, a + b\n",
    "\n",
    "def fib(n):\n",
    "    for i, f in enumerate(fib_eff()):\n",
    "        if i == n:\n",
    "            return f\n",
    "        \n",
    "fib(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comment\n",
    "It is also possible to use another method, calling the function itself as many times as need. \n",
    "\n",
    "Unfortunately this is pretty inefficient! Do not try with n > 35!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "55\n"
     ]
    }
   ],
   "source": [
    "# Solution 5\n",
    "def fib_rec(n):\n",
    "    '''inefficient recursive function as defined, returns Fibonacci number'''\n",
    "    if n > 1:\n",
    "        return fib_rec(n-1) + fib_rec(n-2)\n",
    "    return n\n",
    "print(fib_rec(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9227465\n",
      "7399.85967 ms\n"
     ]
    }
   ],
   "source": [
    "t0 = time()\n",
    "print(fib_rec(35))\n",
    "ts = (time() - t0)\n",
    "print(\"{0:.5f} ms\".format(ts*1000.))\n",
    "# n = 500 takes too much time!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "610\n",
      "0.99254 ms\n"
     ]
    }
   ],
   "source": [
    "# Solution 5b: use an internal cache!\n",
    "def fib_mem(n, _cache={}):\n",
    "    '''efficiently memorized recursive function, returns a Fibonacci number'''\n",
    "    if n in _cache:\n",
    "        #print(_cache)\n",
    "        return _cache[n]\n",
    "    elif n > 1:\n",
    "        # https://docs.python.org/3.6/library/stdtypes.html?highlight=setdefault#dict.setdefault\n",
    "        return _cache.setdefault(n, fib_mem(n-1) + fib_mem(n-2))\n",
    "    return n\n",
    "\n",
    "t0 = time()\n",
    "print(fib_mem(15))\n",
    "ts = (time() - t0)\n",
    "print(\"{0:.5f} ms\".format(ts*1000.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "RecursionError",
     "evalue": "maximum recursion depth exceeded in comparison",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRecursionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-14-c25d458c9a8f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      9\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 11\u001b[1;33m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-14-c25d458c9a8f>\u001b[0m in \u001b[0;36mfib\u001b[1;34m(n)\u001b[0m\n\u001b[0;32m      7\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mn\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "... last 1 frames repeated, from the frame below ...\n",
      "\u001b[1;32m<ipython-input-14-c25d458c9a8f>\u001b[0m in \u001b[0;36mfib\u001b[1;34m(n)\u001b[0m\n\u001b[0;32m      7\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mn\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mn\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[0mfib\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mRecursionError\u001b[0m: maximum recursion depth exceeded in comparison"
     ]
    }
   ],
   "source": [
    "# Now in Python 3\n",
    "# https://stackoverflow.com/questions/18172257/efficient-calculation-of-fibonacci-series\n",
    "import functools\n",
    "\n",
    "@functools.lru_cache(None)\n",
    "def fib(n):\n",
    "    if n < 2:\n",
    "        return n\n",
    "    return fib(n-1) + fib(n-2)\n",
    "\n",
    "fib(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem:\n",
    "\n",
    "Show that using the [Box-Muller transform](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) is possible to general a random sequence of gaussian distributed number taking two random sequences of uniformly distributed numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAHVCAYAAADl4K3UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X28VXWZ///3JYSWaVocpwYwqHCKstSOaHdmaYU5A81kPbT8dTNOpr+o5tvMfLOxtAc1k9U0ZUkmTaaTGaLdISAogiAqcA73AkKHG+MAyRGQezicc67vH2uf49pr77POvl/75vV8PHhwPp+91l7X2mefa1/7sz5rLXN3AQAAAMjuhKQDAAAAAKoZBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgxuCkA4gaOnSojxw5MukwAKAgy5Yte97dm5KOo5LI2wBqVa45u+oK5pEjR6q1tTXpMACgIGb2bNIxVBp5G0CtyjVnMyUDAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgRtVdhxkNwuzFn92TiwNA1SJNAKgWFMyofnxqAmgwpD1gYJX8OylqSoaZ3Wlmu8zs6X4eNzP7sZm1mdlqMzuvmO0hR2bp76JKbzeJbVeDRt//elDnv0NydjLq+C0F5K1W02yxc5jvkjQu5vHLJI1O/btW0u1Fbi9ZA/2Ww4/X4jsin3gHWrbY/a/F16+e1fJ7u9biLa+71EA5O8lffT4fF8U+VymUIp5a+zOrl/1AZRQ1JcPdF5rZyJhFJkj6X3d3SYvN7DQze4277yxmu3kpdLy+d71c1qnWv6z+9n2g16S//cln2XxUw7HHgV6r/uIaaP/z2bdSvA75/O6q2UCvaz6v5UDL1NprU4Rqy9nV9Cso9qMi3/VqRbEpcKDnDT93qVNguX4f+XyEliKGQv5OyvXRU+hHTNwXx1yfI2nlnsM8TNK2ULs91ZeWfM3sWgWjGTrzzDML21Ip3h3ZfqOVKIZzeQdm+4vJJ7ZSF7ZJrB99jmL/uir9RacUmaaUnywDyafo7E8+7+FKfEIm9TdeO3LK2Uko9u2TT9rP9zkqKZ/9KOX37mr4CCl0e0kWYqX8veSyjWIL+HIV+4UuW+ovasUo92Xlsu1Sxu67+xR3b3b35qampjJGE3PsJd9XvxTHyJLabr0o5THPWpbLfub6eKlfq2o77oyB5JSzpWCgw8xazay1o6OjzGHlp5Rv92qbDpHPeklNOSj1a5bt+aotXZSyJCjmfVmu935Sr3U1/Z7LPcLcLmlEqD1c0o4yb7M6Xtn+VHNstWyg17XSo+PVMBo/0PNW+/GvQvD3Vaycc7a7T5E0RZKam5uLejPV+9syH434Fi52FLIU6b/aDj5VYtu19l5LOt5yjzBPl/Sp1JnXF0raV9H5y9Wimr4iJSGffa+G16nSv69q++pea8NPKKXEc3aSB/2KxUGv0qnn1yrbyG21HO0oVDXHVipFjTCb2W8kXSxpqJm1S7pZ0kskyd1/JmmWpA9LapN0WNJni9lew+ErJoASImc3rmobQQVqTbFXybhqgMdd0heK2QaAED7hArwOBamHnF0t4wi8BXNT769Tve9fIer1NSn3lAwAAACgplEwAwAAADEomAEAAIAYFMwAAABADApmAAAAIAYFMwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEoGAGAAAAYlAwAwAAADEomAEAAIAYFMwAAABAjKIKZjMbZ2YbzKzNzG7I8viZZjbfzFaY2Woz+3Ax2wMAFIe8DQD5K7hgNrNBkiZLukzSGElXmdmYyGJflzTN3c+VdKWknxa6PQBAccjbAFCYYkaYx0pqc/fN7t4paaqkCZFlXNKpqZ9fIWlHEdsDABSHvA0ABSimYB4maVuo3Z7qC/umpKvNrF3SLElfzPZEZnatmbWaWWtHR0cRIQEAYpC3AaAAxRTMlqXPI+2rJN3l7sMlfVjSr8wsY5vuPsXdm929uampqYiQAAAxyNsAUIBiCuZ2SSNC7eHKPHR3jaRpkuTuT0k6SdLQIrYJACgceRsAClBMwdwiabSZjTKzIQpODpkeWebPki6RJDN7k4LEy7E7AEgGeRsAClBwwezuXZImSpojab2Cs6rXmtkkMxufWuxfJH3OzFZJ+o2kz7h79PAfAKACyNsAUJjBxazs7rMUnBQS7rsp9PM6Se8qZhsAgNIhbwNA/rjTHwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEoGAGAAAAYlAwAwAAADEomAEAAIAYFMwAAABADApmAAAAIAYFMwAAABCDghkAAACIQcEMAAAAxCiqYDazcWa2wczazOyGfpb5uJmtM7O1ZnZvMdsDAAAAKq3ggtnMBkmaLOkySWMkXWVmYyLLjJb0NUnvcvc3S/rnImIFABSBQQ4AKMzgItYdK6nN3TdLkplNlTRB0rrQMp+TNNnd90qSu+8qYnsAgAKFBjk+IKldUouZTXf3daFlwoMce83sjGSiBYDqUsyUjGGStoXa7am+sLMknWVmT5jZYjMbV8T2AACF6xvkcPdOSb2DHGEMcgBAFsUUzJalzyPtwZJGS7pY0lWS/sfMTst4IrNrzazVzFo7OjqKCAkA0I+SDnKQtwE0kmIK5nZJI0Lt4ZJ2ZFnmj+5+3N23SNqgoIBO4+5T3L3Z3ZubmpqKCAkA0I+SDXJI5G0AjaWYgrlF0mgzG2VmQyRdKWl6ZJk/SHqfJJnZUAWjF5uL2CYAoDAlG+QAgEZTcMHs7l2SJkqaI2m9pGnuvtbMJpnZ+NRicyTtNrN1kuZL+jd3311s0ACAvDHIAQAFKuYqGXL3WZJmRfpuCv3skr6S+gcASIi7d5lZ7yDHIEl39g5ySGp19+mpxz6YGuToFoMcACCpyIIZAFA7GOQAgMJwa2wAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEoGAGAAAAYlAwAwAAADEomAEAAIAYFMwAAABADApmAAAAIAYFMwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhRVMFsZuPMbIOZtZnZDTHLXWFmbmbNxWwPAFAc8jYA5K/ggtnMBkmaLOkySWMkXWVmY7Isd4qkL0laUui2AADFI28DQGGKGWEeK6nN3Te7e6ekqZImZFnuW5K+J+loEdsCABSPvA0ABSimYB4maVuo3Z7q62Nm50oa4e4z4p7IzK41s1Yza+3o6CgiJABADPI2ABSgmILZsvR534NmJ0j6oaR/GeiJ3H2Kuze7e3NTU1MRIQEAYpC3AaAAxRTM7ZJGhNrDJe0ItU+R9BZJj5nZVkkXSprOCSQAkBjyNgAUoJiCuUXSaDMbZWZDJF0paXrvg+6+z92HuvtIdx8pabGk8e7eWlTEAIBCkbcBoAAFF8zu3iVpoqQ5ktZLmubua81skpmNL1WAAIDSIG8DQGEGF7Oyu8+SNCvSd1M/y15czLYAAMUjbwNA/rjTHwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEoGAGAAAAYlAwAwAAADEomAEAAIAYFMwAAABADApmAAAAIAYFMwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGIUVTCb2Tgz22BmbWZ2Q5bHv2Jm68xstZk9amavLWZ7AIDCkbMBoDAFF8xmNkjSZEmXSRoj6SozGxNZbIWkZnd/q6QHJH2v0O0BAApHzgaAwhUzwjxWUpu7b3b3TklTJU0IL+Du8939cKq5WNLwIrYHACgcORsAClRMwTxM0rZQuz3V159rJD2U7QEzu9bMWs2staOjo4iQAAD9KFnOlsjbABpLMQWzZenzrAuaXS2pWdL3sz3u7lPcvdndm5uamooICQDQj5LlbIm8DaCxDC5i3XZJI0Lt4ZJ2RBcys0sl3Sjpve5+rIjtAQAKR84GgAIVM8LcImm0mY0ysyGSrpQ0PbyAmZ0r6Q5J4919VxHbAgAUh5wNAAUquGB29y5JEyXNkbRe0jR3X2tmk8xsfGqx70t6uaT7zWylmU3v5+kAAGVEzgaAwhUzJUPuPkvSrEjfTaGfLy3m+QEApUPOBoDCcKc/AAAAIAYFMwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEoGAGAAAAYlAwAwAAADEomAEAAIAYFMwAAABADApmAAAAIAYFMwAAABCDghkAAACIUVTBbGbjzGyDmbWZ2Q1ZHj/RzO5LPb7EzEYWsz0AQHHI2wCQv4ILZjMbJGmypMskjZF0lZmNiSx2jaS97v4GST+U9N1CtwcAKA55GwAKU8wI81hJbe6+2d07JU2VNCGyzARJd6d+fkDSJWZmRWwTAFA48jYAFKCYgnmYpG2hdnuqL+sy7t4laZ+kVxWxTQBA4cjbAFCAwUWsm23EwQtYRmZ2raRrU82DZrahgHiGSnq+gPVqRT3vH/tWu+p3/8wK3bfXljqUEiJvVw77Vrvqef/qdt/MCt63nHJ2MQVzu6QRofZwSTv6WabdzAZLeoWkPdEncvcpkqYUEYvMrNXdm4t5jmpWz/vHvtWuet6/Ot038naFsG+1q573j30rXDFTMlokjTazUWY2RNKVkqZHlpku6dOpn6+QNM/dM0YqAAAVQd4GgAIUPMLs7l1mNlHSHEmDJN3p7mvNbJKkVnefLukXkn5lZm0KRiiuLEXQAID8kbcBoDDFTMmQu8+SNCvSd1Po56OSPlbMNvJQ1KHBGlDP+8e+1a563r+63DfydsWwb7WrnvePfSuQcaQNAAAA6B+3xgYAAABi1FzBXM+3dc1h375iZuvMbLWZPWpm1Xz5qgwD7V9ouSvMzM2sZs7kzWXfzOzjqd/fWjO7t9IxFiqH9+WZZjbfzFak3psfTiLOQpjZnWa2y8ye7udxM7Mfp/Z9tZmdV+kYa10952ypvvM2Obs2c7ZUv3k70Zzt7jXzT8FJKpskvU7SEEmrJI2JLPP/S/pZ6ucrJd2XdNwl3Lf3SXpZ6ufra2Xfct2/1HKnSFooabGk5qTjLuHvbrSkFZJOT7XPSDruEu7bFEnXp34eI2lr0nHnsX8XSTpP0tP9PP5hSQ8puDbxhZKWJB1zLf2r55ydx/7VZN4mZ9dmzs5j/2oybyeZs2tthLmeb+s64L65+3x3P5xqLlZwDdVakcvvTpK+Jel7ko5WMrgi5bJvn5M02d33SpK776pwjIXKZd9c0qmpn1+hzOv6Vi13X6gs1xgOmSDpfz2wWNJpZvaaykRXF+o5Z0v1nbfJ2bWZs6U6zttJ5uxaK5jr+bauuexb2DUKvkXVigH3z8zOlTTC3WdUMrASyOV3d5aks8zsCTNbbGbjKhZdcXLZt29KutrM2hVcfeGLlQmtIvL9u0S6es7ZUn3nbXJ2beZsqbHzdtlydlGXlUtAyW7rWoVyjtvMrpbULOm9ZY2otGL3z8xOkPRDSZ+pVEAllMvvbrCCQ3wXKxhhetzM3uLuL5Q5tmLlsm9XSbrL3X9gZu9QcA3ft7h7T/nDK7tazSfVop5ztlTfeZucXZs5W2rsvF22fFJrI8z53NZVFnNb1yqUy77JzC6VdKOk8e5+rEKxlcJA+3eKpLdIeszMtiqYezS9Rk4iyfV9+Ud3P+7uWyRtUJCMq10u+3aNpGmS5O5PSTpJ0tCKRFd+Of1dol/1nLOl+s7b5OzazNlSY+ftsuXsWiuY6/m2rgPuW+rw1x0Kkm4tzaeSBtg/d9/n7kPdfaS7j1Qw12+8u7cmE25ecnlf/kHByT8ys6EKDvdtrmiUhcll3/4s6RJJMrM3KUi8HRWNsnymS/pU6szrCyXtc/edSQdVQ+o5Z0v1nbfJ2bWZs6XGztvly9mVPLuxFP8UnAG5UcEZoDem+iYp+EOVgl/6/ZLaJC2V9LqkYy7hvs2V9Jyklal/05OOuZT7F1n2MdXIGdc5/u5M0n9LWidpjaQrk465hPs2RtITCs7EXinpg0nHnMe+/UbSTknHFYxMXCPpOknXhX5vk1P7vqaW3pPV8q+ec3aO+1ezeZucXZs5O8f9q8m8nWTO5k5/AAAAQIxam5IBAAAAVBQFMwAAABCDghkAAACIQcEMAAAAxKBgBgAAAGJQMAMAAAAxKJgBAACAGBTMAAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEyKlgNrNxZrbBzNrM7IaY5a4wMzez5lR7pJkdMbOVqX8/K1XgAIDsyNkAUFqDB1rAzAZJmizpA5LaJbWY2XR3XxdZ7hRJX5K0JPIUm9z9nBLFCwCIQc4GgNLLZYR5rKQ2d9/s7p2SpkqakGW5b0n6nqSjJYwPAJAfcjYAlNiAI8yShknaFmq3S7ogvICZnStphLvPMLN/jaw/ysxWSNov6evu/nh0A2Z2raRrJenkk09++xvf+MY8dgEAqseyZcued/emBEMoe85OPQd5G0DNyzVn51IwW5Y+73vQ7ARJP5T0mSzL7ZR0prvvNrO3S/qDmb3Z3fenPZn7FElTJKm5udlbW1tzCAsAqo+ZPZt0CFn6SpqzJfI2gPqQa87OZUpGu6QRofZwSTtC7VMkvUXSY2a2VdKFkqabWbO7H3P33ZLk7sskbZJ0Vi6BAQAKQs4GgBLLpWBukTTazEaZ2RBJV0qa3vugu+9z96HuPtLdR0paLGm8u7eaWVPqBBSZ2eskjZa0ueR7AQDoRc4GgBIbcEqGu3eZ2URJcyQNknSnu681s0mSWt19eszqF0maZGZdkrolXefue0oROAAgEzkbAErP3H3gpSqIuXAAapmZLXP35qTjqCTyNoBalWvO5k5/AAAAQAwKZgAAACAGBTMAAAAQg4IZAAAAiEHBDAAAAMSgYAYAAABiUDADAAAAMSiYAQAAgBgUzAAAAEAMCmYAAAAgBgUzAAAAEIOCGQAAAIhBwQwAAADEyKlgNrNxZrbBzNrM7IaY5a4wMzez5lDf11LrbTCzD5UiaABA/8jZAFBagwdawMwGSZos6QOS2iW1mNl0d18XWe4USV+StCTUN0bSlZLeLOmvJc01s7Pcvbt0uwAM7PDhw1qyZIkOHjyoQYMGqbm5WWeccUbSYQElR85GPeju7tayZcv03HPPSZLe+MY3avTo0QlHhUaWywjzWElt7r7Z3TslTZU0Icty35L0PUlHQ30TJE1192PuvkVSW+r5gIo4fvy4brvtNr32ta/V+9//fo0fP16XX365hg0bpi9/+ct6/vnnkw4RKDVyNmqWu2vWrFl661vfqgsuuEDjx4/X+PHjddZZZ2n8+PFav3590iGiQeVSMA+TtC3Ubk/19TGzcyWNcPcZ+a6bWv9aM2s1s9aOjo6cAgcG0tHRoQsuuEBf/OIXMwrjrq4u/fjHP9bf/M3f6Mknn0woQqAsyp6zU89B3kZJdXV16VOf+pQuv/xyrVu3LuPxBx98UGeffbbuuOOOBKJDo8ulYLYsfd73oNkJkn4o6V/yXbevw32Kuze7e3NTU1MOIQHxDhw4oA+fcYZWrFgRu9yePXt0+eWX6+mnn65QZEDZlT1nS+RtlJa76/rrr9c999wTu1x3d7euv/56TZ06tUKRAYFcCuZ2SSNC7eGSdoTap0h6i6THzGyrpAslTU+dRDLQukDJHT9+XB/96EfVGukfIuldkl4Z6X/hhRf0oQ99SH/+858rEyBQXuRs1Jybb75Z//M//5PlkfMkjUzrcXd96lOf0ty5cysRGiApt4K5RdJoMxtlZkMUnBAyvfdBd9/n7kPdfaS7j5S0WNJ4d29NLXelmZ1oZqMkjZa0tOR7AYT813/9lx555JG0vkskPStpkYJP/3+LrLNjxw599rOflXvWwTSglpCzUVMWLFigb33rW5HeEZJWSFomabOkXyt8nYLjx4/rk5/8pPbu3VuxONHYBiyY3b1L0kRJcyStlzTN3dea2SQzGz/AumslTZO0TtJsSV/gbGuU0+bNmzVp0qS0vmZJv5f06lT7REnflfSFyLrz5s3Tb37zm7LHCJQTORu1pLOzU9dff32k91WSHpZ0Tqptkj4h6e60pXbt2qUbb7yx7DECkmTVNqLW3Nzsra3Rg+nAwNxdl19+uR566KG+vlcp+OTPdgG5bkkflDQv1HfGGWfomWee0emnn17WWFG/zGyZuzcPvGT9IG+jUN/5znf07//+75He+ZIu7meNmxRc4CVgZnrqqad0wQUXlCdA1L1cczZ3+kPdmDFjRlqxLEn/pezFsiQNknS7grnNvXbt2qVvfvObZYkPAPCi7du3Z5mK8Y/qv1iWpH9XMFMo4O76whe+wHQ6lB0FM+qCu2ck3oskfXqA9c6S9LVI35QpU/oulg8AKI8f/OAHOnLkSKjnVQomzMU5SdJP03qWLVum2bNnlzg6IB0FM+rCo48+qpaWlrS+nyj7NbKiblBwKYBeR48e1a233lrC6AAAYbt379aUKVMivZMkDc1h7Usl/UNaz3/+53+WKDIgOwpm1IXvfOc7ae3xkt6a47onSfrXSN/kyZO1b9++EkQGAIj6yU9+okOHDoV6Xq1gOkau0k/2W7RokRYtWlSK0ICsKJhR85YuXap58+al9UWnWQzknxQcDOy1f/9+3X777cWGBgCIOHjwoH7yk59Eer+iYPgiV+cpOG37RdGBE6CUKJhR82677ba09sUK7sSQj5MlfTnS99Of/lTd3VxRCwBK6Te/+Y327NkT6jlN0nUFPFP60MisWbO0adOmYkID+kXBjJr2wgsv6P7770/r+78FPtdESS8Ltbdt25ZxAxQAQHEy7+h3vYIbUObrvZLGpvXceeedBUYFxKNgRk279957dfTo0b72ayV9qMDnOl3SxyN92W/VCgAoxOrVq7V0afjmkSbp2gKfzSR9Pq3nl7/8pbq6ugp8PqB/FMyoadGC9h9V3Jv6nyLtP/7xj9q1a1cRzwgA6PWLX/wi0nOppJFFPOPHJb28r7Vz586M6/EDpUDBjJq1fPlyrVixoq9tkj5b5HO+U9IbQ+2uri797//+b5HPCgA4evSofvWrX0V6o8MU+Xq5pCvTejgyiHKgYEbNuuuuu9La4ySNKPI5TZnpO7odAED+HnzwQe3duzfU8ypJE0rwzOlZe+bMmRwZRMlRMKMmdXd364EHHkjry+cKnnE+JWlwqL127VqtXbu2RM8OAI3pvvvui/RcLenEEjzzWElv7mt1d3fr97//fQmeF3gRBTNq0qJFi7Rz586+9ssl/W2JnrtJway6sGnTppXo2QGg8Rw8eFAzZ86M9H6iRM9ukq5K68kszoHi5FQwm9k4M9tgZm1mdkOWx68zszVmttLMFpnZmFT/SDM7kupfaWY/K/UOoDFFC9gJyu+S9wOJXi1j2rRpcvcSbgEoH3I2qs2MGTPSrmgUXNPo/BJu4WNprQULFugvf/lLCZ8fjW7AgtnMBkmaLOkySWMkXdWbXEPudfez3f0cSd+T9N+hxza5+zmpf4VcmRxIk206RrTALdZHJL0k1H7mmWe0Zs2aEm8FKD1yNqpR5ojvxxWMDJfKWZLO6Wv19PTod7/7XQmfH40ulxHmsZLa3H2zu3dKmqrILH133x9qniyJoTiUzcKFC9NO6DhV0RukFu90SR+I9DEtAzWCnI2qsn///iyXeiv1MEfmc5KzUUq5FMzDJG0LtdtTfWnM7AtmtknBaMWXQg+NMrMVZrbAzN6TbQNmdq2ZtZpZa0dHRx7hoxFF7+xX6ukYvaLpPLpdoEqVPWen1idvIyczZszQsWPHQj2vk/T2MmwpPWsvXLiQaRkomVwK5mzHTDJGI9x9sru/XtJXJX091b1T0pnufq6kr0i618xOzbLuFHdvdvfmpqam3KNHw3F3Pfjgg2l9H+tn2WJNUPq0jI0bN2rDhg1l2hpQMmXP2an1ydvIyfTp0yM9V6i00zF6vV7SeX0td89yoiFQmFwK5nalX952uKQdMctPVTAFVO5+zN13p35eJmmTgolGQEFWrlyp9vb2vvZLlXlFi1I5TdJ7I33RYh2oQuRsVI3jx49r9uzZkd5SXHu5P+PTWuRslEouBXOLpNFmNsrMhii4pU7a10UzGx1qXi7pT6n+ptQJKDKz10kaLWlzKQJHY4omvw8oKJrL5e8G2D5QhcjZqBqPP/649u3bF+oZKumCMm4xPWs/8sgjkatzAIUZsGB29y5JEyXNkbRe0jR3X2tmk8ys96vcRDNba2YrFRzG+3Sq/yJJq81slaQHJF3n7ntKvhdoGNGCdXw/y5VKtGB+4okntGcPb2FUL3I2qknmIMPfShpUxi2eq/CU/cOHD2vevHll3B4ahVXbtWWbm5u9tbUeg7EtAAAgAElEQVQ16TBQhXbs2KFhw9LPXdop6dVl3u7Zkp4Ote+55x598pOfLPNWUavMbJm7NycdRyWRt5GNu2v06NHatGlTqPe3kv6hzFu+XtKLlxC/7rrrdPvtt5d5m6hVueZs7vSHmhE9eWOsyl8sS0zLAIBCPPPMM5FieYhKfxHQbNKz9owZM7jxFIpGwYyaES2Yo4VsuUS3M3v2bHV1dVVo6wBQmzKvUPE+SS+vwJbfL+llfa329natWrWqAttFPaNgRk04fvy4Hn300bS+yyu07bEKTlPptW/fPi1durRCWweA2jRnzpxIT6Wy9kmKXj/p4YcfrtC2Ua8omFETnnrqKR08eLCv/VeS3lahbQ9S5l3/Mj8IAAC9Dh8+rMcffzzSO66CEXworUXORrEomFETosnug6rsm/dDkTbJFwD6t2DBgsjd/UZJekMFI0jP2osWLdKhQ4cquH3UGwpm1IRogRotYMstOsLc0tLC5eUAoB+ZUyA+qPLc3a8/r1dwC+5AZ2enHnvssQpuH/WGghlVr6OjQ8uXL0/rixaw5fbXCi4v16unp0dz586tcBQAUBsyj8JVepgjc5scGUQxKJhR9R555JG0SwKdK+mMBOJgWgYADGzbtm1av359qGeQgitXVBoFM0qHghlVL3poL4lximzbffjhh7m2JwBEZE7HeIekVyQQyfskDe5rbdy4UVu3bk0gDtQDCmZUNXfPmPqQVMH8bkkvDbXb29vV1taWUDQAUJ0yp6tV4mYl2Zwq6Z1pPfPnz08mFNQ8CmZUtT/96U/avn17X/ulCsYqknCSpHdF+ki+APAid9e8efMivZdmXbYy0qeCkLNRKApmVLXozUreLenEZEKRFBzgCyP5AsCL1q5dq127doV6Xi6pOalwFM3a8+bNYyodCkLBjKoWHam4JKE4emUrmEm+ABDIHF1+r6SXJBFKygUKjg8Gtm/fzlQ6FCSngtnMxpnZBjNrM7Mbsjx+nZmtMbOVZrbIzMaEHvtaar0NZpbU9FPUoJ6enowR3CTOsw5rlnRyqP3cc8/pmWeeSSocICtyNpISPSqYfNY+UdHJdBwZRCEGLJjNbJCkyZIukzRG0lXh5Jpyr7uf7e7nSPqepP9OrTtG0pWS3qzgnpg/TT0fMKDVq1dr9+7dfe1XSDovuXAkBeMk74n0kXxRTcjZSEpXV1eWm4MkfVxQih4bJGejELmMMI+V1Obum929U9JUSRPCC7j7/lDzZEm9x6gnSJrq7sfcfYukttTzAQOKjlRcrOBqnkljHjOqHDkbiVi+fLn27w+/tYYq/ZZPScksmJlKh3zlUjAPk7Qt1G5P9aUxsy+Y2SYFoxVfynPda82s1cxaOzo6co0dda7a5i/3ihbMjz32mHp6ehKJBcii7Dk7tT55G2ky5y+/T9VxqtT5Ck+mYyodCpHLOznbzd8zvpq5+2R3f72kr0r6ep7rTnH3ZndvbmpqyiEk1Lvu7m4tWrQorS9aqCblXAVX9+z1/PPPa+3atUmFA0SVPWen1idvI83ChQsjPdWStV+i4BpLL+LIIPKVS8HcLmlEqD1c0o6Y5adK+kiB6wKSgvnL4UN7QxVMqqwGgyVdFOkj+aKKkLNRcdkGOYKJdNWCecwoTi4Fc4uk0WY2ysyGKDghZHp4ATMbHWpeLulPqZ+nS7rSzE40s1GSRktaWnzYqHfRkYr3KPvQV1KYx4wqRs5Gxa1atUoHDhwI9QyV9MakwskiPWszlQ75GjzQAu7eZWYTJc1RcM7Vne6+1swmSWp19+mSJprZpZKOS9or6dOpddea2TRJ6yR1SfqCu3eXaV9QR6IFc3REN2nRgnnBggXq6enRCSdUw3w9NDJyNpKQOR3jIlXXMMd5kk6RFBT1vVPpzj67Gk5KRC0YsGCWJHefJWlWpO+m0M9fjln3PyT9R6EBovG4e9UXzG+TdLqCSkOS9u7dq1WrVuncc89NMCogQM5GpWUvmKvJYAU3UZnR1zN//nwKZuSM4TBUnWeeeUbPP/98X/sUBQVqNTlBmbPzmJYBoBFlG+SovoJZYh4zikHBjKoTTbzvVnVcfzmKecwAIK1fvz7tJlPBdYTemlQ4MdKzdu9UOiAXFMyoOtU+HaNXtGBeuHChurq6EokFAJKSObpcrcMcvZPpAr1T6YBcUDCjqri7FixYkNZXrQXzmyWFrz67f/9+ki+AhlMb0zGkoOR5b1pP9PMG6A8FM6rK1q1btX379r72SZKakwsnlinzYyHzgwMA6lftzF/ulV4wk7ORKwpmVJVo8nqHpCHJhJITCmYAjWzLli1pgxzSSyW9PalwcpCetRcuXCj3rDezBNJQMKOqZLthSTWLFsyPP/44J5EAaBiZUxqqfZjjbQquvRTYvXu31q9fn1w4qBkUzKgqtXLCX6+zJb0i1Cb5AmgktTUdQwpORnx3Wg9HBpELCmZUjR07dqitra2vPVjShcmFk5PM1EvyBdA4aq9glrJNywAGQsGMqvH444+ntZslnZxMKHlhHjOARtTe3q7NmzeHel4i6YKkwskD85iRPwpmVI1am47RK1vBTPIFUO+igxzS+ZJelkQoeWpWcA2mwPbt27Vly5bkwkFNoGBG1ajVgvk8pX9E7NixIzLqAgD1pzanY0jBSYnvSOvhyCAGQsGMqrB79249/fTTfW2T9K7kwslLZuol+QKof7VbMEvMY0a+ciqYzWycmW0wszYzuyHL418xs3VmttrMHjWz14Ye6zazlal/00sZPOrHokWL0tpvk3RaMqEUhHnMqCbkbJRbR0eH1q1bF+o5QdI7kwqnABTMyM+ABbOZDZI0WdJlksZIusrMxkQWWyGp2d3fKukBSd8LPXbE3c9J/RtforhRZ2p1OkYvCmZUC3I2KiE6yCGdo/SLbFa7CxVciymwadOmyA1YgHS5jDCPldTm7pvdvVPSVEkTwgu4+3x3P5xqLpY0vLRhot7VesF8gYLzw3tt3rxZ7e3tSYWDxkbORtnV9nQMKTjz5Py0nsyTGIEX5VIwD5O0LdRuT/X15xpJD4XaJ5lZq5ktNrOPZFvBzK5NLdPa0dGRQ0ioJwcOHNDy5cvT+qr9Dn9RL1VQpYSRfJGQsudsibzd6Gq/YJaYloF85FIwW5a+rNfMMrOrFVyv5fuh7jPdvVnSJyT9yMxen/Fk7lPcvdndm5uamnIICfXkySefTLud9BslnZFcOAVjWgaqRNlztkTebmT79u3TypUrI73RWzjVAgpm5C6Xgrld0ohQe7ikHdGFzOxSSTdKGu/ux3r73X1H6v/Nkh6TdG4R8aIO1fp0jF4UzKgS5GyUVXSQI5gqX4tfmt6l8PfLtWvX6vnnn08uHFS1XArmFkmjzWyUmQ2RdKWktDOnzexcSXcoSLy7Qv2nm9mJqZ+HKnh3hk+rBeqmYH6n0v+g1q1bJw5VIwHkbJRVfUzHkIKTFM9J68k8mREIDFgwu3uXpImS5khaL2mau681s0lm1nsG9fclvVzS/ZFLEb1JUquZrZI0X9It7k7yRZ8jR45o6dKlaX21mnpPVeZQHMkXlUbORrnVT8EsMS0DuRo88CKSu8+SNCvSd1Po50v7We9JSWcXEyDq29KlS9XZ2dnXHqn0Y8m15iJJy0LthQsX6u///u+TCgcNipyNcjl8+LBaWloivbV2mnbYRZJu7WtRMKM/3OkPiaqX6Ri9mMcMoJ4tWbJEx48fD/W8TrV9VcL0Yn/FihXav39/QrGgmlEwI1H1VjBHzxNfuXKl9u3bl0gsAFBq9TUdQwpOVnxTX6unp0dPPvlkcuGgalEwIzHHjx/PSEy1nnqHSnpzqE3yBVBP6q9glpjHjFxQMCMxy5cv1+HDh/var5b0huTCKRmmZQCoR52dnXrqqacivRTMaAwUzEhMtukY2e64UGuiHx8LFixIJA4AKKVly5bpyJEjoZ6/VjCHudalz2NeunRpZD8BCmYkqN7mL/eKni/e0tKSNpIOALUo+3SMehjmGCFpVF/r+PHjWrJkSXLhoCpRMCMR3d3dGdcorpeCeZik8L2Eu7q6tHjx4qTCAYCSqM/5y72YloF4FMxIxJo1a/TCCy/0tU9X+slytY55zADqSbZBDgpmNBIKZiQiOq/3ItXXm5GCGUA9WblyZeT6xEMljUkqnDJIz9pPPvlk2k21gHqqUVBDogXzexOKo1yiBfNTTz1F8gVQszJPXn6P6mP+cq/XS3pNX+vIkSNavnx5cuGg6lAwo+J6enoyRlzrrWAepWAuc6+jR4+qtbU1qXAAoCiZBXO9ZW0T0zIQh4IZFbdu3Trt3r27r32qpLclF05ZZKZeki+A2tTT06PHH3880ltvBbNEwYw4FMyouOhIxbslDUomlLKiYAZQD55++mnt3bs31HOapLOTCqeM0rP2okWL1N3dnVAsqDY5FcxmNs7MNphZm5ndkOXxr5jZOjNbbWaPmtlrQ4992sz+lPr36VIGj9pU7/OXe0ULZpIvKoWcjVLKPn+5Hoc5xkh6ZV9r3759WrNmTXLhoKoMWDCb2SBJkyVdpuDddJWZRU+NXSGp2d3fKukBSd9LrftKSTdLukDSWEk3m9nppQsftcbd637+cq83KTiPvNeBAwe0atWqpMJBgyBno9Tqf/5yrxMUvfUURwbRK5cR5rGS2tx9s7t3SpoqaUJ4AXef7+69tzJbLGl46ucPSXrE3fe4+15Jj0gaV5rQUYs2btyo5557rq99sqTzkgunrEyZd/0j+aICyNkomWyDHPVbMEvMY0Z/cimYh0naFmq3K/0CAFHXSHoon3XN7FozazWz1o6OjhxCQq2KjlS8S9JLkgmlIpjHjASUPWdL5O1GsX79eqX/fk+RdE5S4VRAZsHs7gnFgmqSS8Gc7UKLWd89Zna1pGZJ389nXXef4u7N7t7c1NSUQ0ioVY0yf7lXtoKZ5IsyK3vOlsjbjSLzS/67JA1OIpQKOUfSy/taHR0d2rBhQ3LhoGrkUjC3SxoRag+XtCO6kJldKulGSePd/Vg+66IxuHvDFcxvUzAe02v37t1av359UuGgMZCzUTKNM3+512AFXwpexJFBSLkVzC2SRpvZKDMbIulKSdPDC5jZuZLuUJB4d4UemiPpg2Z2eurEkQ+m+tCANm/erO3bt/e1Xyrp/OTCqYhBCi6bF0byRZmRs1ES2QY56r9glpjHjGwGLJjdvUvSRAVJc72kae6+1swmmdn41GLfV3AM434zW2lm01Pr7pH0LQUJvEXSpFQfGlA08b5D0pBkQqko5jGjksjZKJW2tjbt3Lkz1PMyBTN46l161l6wYAFT6ZDbRCR3nyVpVqTvptDPl8ase6ekOwsNEPWj0aZj9OpvHrNZtumiQPHI2SiFzNHld6q+T9Pudb6kEyUFM5Xa29v17LPPauTIkUkGhYRxpz9UTKNcfzmqWdJJofb27du1efPmpMIBgJw01uXkwk6UdGFaT+aXBzQaCmZUxJYtW7R169a+9hAFF4ttBEMUTD8Jmz9/fhKhAEBO3F3z5s2L9EaPl9Wz9H0lZ4OCGRXx6KOPprXfqeCkv0bxvkg7+noAQDXZuHFj2knaQca+IKlwEpCetR999FHmMTc4CmZUxNy5c9Pa/U6grFPR/X300UfV09OTSCwAMJBozg5GXE9MIpSEvEPhYZ329nZt3LgxuXCQOApmlF1PT0/Gob1LEoolKecr/XrMHR0devrpp5MKBwBiZR4Fa7SsfZKk96T1cGSwsVEwo+zWrFmTdmvVU9UYFyYKGyzp4khf5ggOACSvu7s7y/zlRjsuKEW/JJCzGxsFM8ou+q38YtX3jVX7Ex2fYbQCQDVatmyZ9u3bF+p5lYL7ljaa9C8J8+fPV3d3d0KxIGkUzCi7Rp+/3Cu63wsWLFBnZ2cisQBAfzK/zL9fjVkunCPplX2tF154QcuXL08uHCSqEf8CUEGdnZ0Z1/JstJlwvcZIenWofejQIS1dujSpcAAgq8ypB42atU9Q8GXhRRwZbFwUzCirJUuW6NChQ33t10h6U3LhJMqU+bHDnDgA1eTIkSN64oknIr2NelxQYh4zelEwo6yi38YvUVA4NirmMQOoZk888YSOHTsW6nmtpNclFU4VSP+ysGjRIh05ciShWJAkCmaUFfOX00UL5sWLF+vgwYOJxAIAUZlf4i9VYw9zvF7SmX2tY8eO6cknn0wuHCSGghllc+DAAS1ZsiStr1FnwvU6U9LoULurqytjjjcAJIX5y1Gm6FAPRwYbU04Fs5mNM7MNZtZmZjdkefwiM1tuZl1mdkXksW4zW5n6N71UgaP6LVy4UF1dXX3tv5E0PLlwqka2u/4BpUTORiH27t2rZcuWRXrfn3XZxpL+pYGc3ZgGvByumQ2SNFnSByS1S2oxs+nuvi602J8lfUbSv2Z5iiPufk4JYkWNyTZ/GcHrcHuozUkkKCVyNgo1f/58uXuo52xJf5VUOFUk/dOrtbVVL7zwgk477bSE4kESchlhHiupzd03u3unpKmSJoQXcPet7r5aUk8ZYkSNYv5ydu9T+ozA1atXa9euXUmFg/pDzkZBss9fRvCl4S19rZ6eHj322GOJRYNk5FIwD5O0LdRuT/Xl6iQzazWzxWb2kWwLmNm1qWVaw7dQRu167rnntGbNmr72Ccq8NXSjeqWk8yJ9jDKjhMqesyXydj165JFHIj0cF3xR+peHhx9+OKE4kJRcCuZsp8d6lr7+nOnuzZI+IelHZvb6jCdzn+Luze7e3NTUlMdTo1rNnj07rf12SacnE0pV+kCkHX29gCKUPWdL5O16s2nTJv3pT38K9bxE0kVJhVOF0rP27NmzI9NXUO9yKZjbJY0ItYdL2pHrBtx9R+r/zZIek3RuHvGhRj300ENp7csSiqNajYu0Z8+erZ4ejo6jJMjZyFs0Z0vvlnRKEqFUqYslndjX2rJlizZu3JhYNKi8XArmFkmjzWyUmQ2RdKWknM6cNrPTzezE1M9DJb1L0rr4tVDrurq6NGfOnLQ+CuZ075R0aqjd0dGh5cuXJxUO6gs5G3nLLJjJ2uleJum9aT2zZs1KJhQkYsCC2d27JE2UNEfSeknT3H2tmU0ys/GSZGbnm1m7pI9JusPM1qZWf5OkVjNbJWm+pFsiZ2qjDi1ZskQvvPBCX/tVks5PLpyq9BJlnk6T+YEF5I+cjXwdPXpU8+fPj/RSMGdKf03I2Y1lwMvKSZK7z5I0K9J3U+jnFmW5xK67P6ngujRoINEk8kFJg5IJpap9WNLvQu2HHnpI3/jGN5IKB3WEnI18LFiwIHK75xGS3pxUOFXsMkn/p6+1YMECHTp0SCeffHJyIaFiuNMfSo75y7mJzmNevHixdu/enUgsABpX5tSCy9TYt8Puz1mSRvW1Ojs7s4zMo15RMKOkdu7cmTEX90MJxVLthkl6a6jt7hlzvwGg3Ji/nCtT9LWZOXNmMqGg4iiYUVLR5NEs6YxkQqkJH460Z8yYkUgcABrThg0bslxOjusv9y89a8+YMYPLyzUICmaU1B//+Me09viE4qgVfxdpz5o1S8ePH08kFgCNJ5qzg3uRcjm5/r1f0kv7Wu3t7VqxYkVy4aBiKJhRMocOHcq4Yx0Fc7wLJIVv+bBv3z4tXLgwqXAANJjp06NXHCRrx3upglPZX5T5GqIeUTCjZObOnaujR4/2tc9U+hxdZBok6W8jfSRfAJXQ0dGhJ598MtIbPe6FTOlfKsjZjYGCGSUTTRrjxXnWuYiO50yfPp05cQDKbubMmZFcc46CoQ7Eu1zhT7cVK1Zo27ZtyYWDiqBgRkl0d3frwQcfTOvjwF5uPiDppFB769atevrpp5MKB0CDYDpGof5K0oVpPdHPP9QfCmaUxOLFi9XR0dHXPlXRm4iiPycr865/v//975MIBUCDOHz4cJbLWE5IJJbalP5akbPrHwUzSuKBBx5Ia18maUgyodSk6MfUb3/720TiANAYZs+ercOHD4d6hks6N6lwalB61p4/fz43nqpzFMwoWk9PT0bB/LGEYqlVH1H67cNXr16tjRs3JhUOgDoXzdnSFeKsk3y8UeHbh3d3d2e5RB/qCQUzirZ06VK1t7f3tV8m7hOVr6GSLo70ZX6gAUDxjhw5kmXO7RWJxFLb0l+z+++/P6E4UAkUzChatLD7sIKiGfmJflxRMAMoh4cfflgHDx4M9bxG0juSCqeGpWftuXPnau/evQnFgnLLqWA2s3FmtsHM2szshiyPX2Rmy82sy8yuiDz2aTP7U+rfp0sVOKqDu2cUdoxTFObvlf4HuWLFCm3atCmpcFDDyNmIk/ll/KNi/KwQb5b0N32trq4urslcxwb8CzGzQZImKzjKPkbSVWY2JrLYnyV9RtK9kXVfKelmBTc0GyvpZjM7vfiwUS1aWlr07LPP9rVPUnCFSuTvryRdFOnjEB/yRc5GnKNHj2Yp6jjrpDCm6Gs3bdq0ZEJB2eXylXKspDZ33+zunZKmKnJ6qLtvdffVknoi635I0iPuvsfd90p6RNK4EsSNKvHrX/86rX2ZpJcnE0pdiI7O33PPPdzEBPkiZ6NfM2bM0P79+0M9fyXpXUmFUwfSs/acOXPSLrGK+pFLwTxMUvgWNu2pvlzktK6ZXWtmrWbWyhutdnR1dWnq1KlpfZ9MKJZ68TGlXy1j7dq1Wr16dVLhoDaVPWdL5O1adc8990R6rlR61kF+3qro1TLuu+++5MJB2eRSMGe7zkyuQ145revuU9y92d2bm5qacnxqJG3u3LnatWtXX/sVYjpGsc5QMMQXFh3FBwZQ9pwtkbdr0Z49ezRr1qxI79WJxFI/TNGhInJ2fcqlYG6XNCLUHi5pR47PX8y6qHLRkYorlH6LZxQmOkp/7733qru7O5FYUJPI2cjq/vvv1/Hjx0M9Z0l6e1Lh1JFPpLUWL16stra2hGJBueRSMLdIGm1mo8xsiILjN7meBjpH0gfN7PTUiSMfTPWhxh06dEh/+MMf0vqYjlEaExTcLrvX9u3btWDBgqTCQe0hZyOrzJHPq8XNSkrhtZLek9bDKHP9GbBgdvcuSRMVJM31kqa5+1ozm2Rm4yXJzM43s3YFUzDvMLO1qXX3SPqWggTeImlSqg817v7779ehQ4f62sMlvTe5cOrKyQouMRd21113JRAJahE5G9m0tbXp8ccfj/R+IuuyKET61Ja7775bPT3Rc2pRy6zazsBvbm721tbWpMPAAN797nfriSee6Gt/VdItyYVTdx5W+lzmk046STt37tRpp52WVEjIkZktc/fmpOOoJPJ29fva176mW24JZ+l3Snqiv8WRt70KbgBzrK9n7ty5uuSSSxKLCLnJNWdzpXLkbd26dWnFsiRdk1As9epSBQf5eh09epRDfAAKcvz4cf3yl7+M9P5TIrHUr9MV3ADmRT//+c+TCQVlQcGMvP3iF79Ia18saXQikdSvE5T5JeTnP/8512QGkLeZM2fqueeeC/WcIunjSYVTxz6X1vr973+v559/PqFYUGoUzMjLsWPHdPfdd6f1MU5RHp9R+h/oqlWrtGzZsoSiAVCrMkc6P6H0U4tRGu+V9Ia+Vmdnp371q18lFw5KioIZeXnggQe0e/fuvvZpkv4huXDq2ghl3mLt9ttvTyIUADVqy5Ytmj17dqSXYY7yMEWPDf7sZz/j5L86QcGMnLm7br311rS+/0/SS5MJpyFcG2n/+te/5rarAHJ22223RQq2c8S1l8vpM5IG97U2btyoOXO4MmM9oGBGzhYvXqyWlpa0vokJxdIo/lbSyFD72LFjmjJlSkLRAKglBw4cyDjnRPqSuPZyOb1a0fnh0YEm1CYKZuTsRz/6UVr7wwruE4XyGaTMLyWTJ09WZ2dnEuEAqCF333239u3bF+ppknRVUuE0kC+ntebMmaP169cnFAtKhYIZOdm2bZt++9vfpvV9uZ9lUVrXKP30nJ07d+r+++9PKhwANaCnp0c/+clPIr2fl3RSEuE0mLGSLkzrYZS59lEwIyc/+MEP1N3d3dd+k6QPJBdOQzlN0qcjfd/97nc5kQRAv373u99p48aNoZ7Bkq5PKpwGlD6kdNddd+kvf/lLQrGgFCiYMaDnnnsuY97sP4tZcJX0ZaW/3mvWrNGDDz6YVDgAqpi769vf/nak90pJf51EOA3qowqudRQ4duyYfvCDHyQXDopGwYwB/fCHP9SRI0f62sOUOeKJ8jpLmbcZ+Pa3v82NTABkmDlzplatWhXqMUlfSyqcBvUSSV9N67n99tu5kUkNo2BGrD179mjy5Mlpff9X0onJhNPQboy0W1tbuVwRgDTZR5c/KmlMEuE0uH9UcNWMwKFDhzJOnkftoGBGrFtuuUUHDx7sa5+h6M0/USlnS/pIpO/rX/86c5kB9HnwwQe1ZMmSSO/XE4kFL5X0b2k9t956a+Q25agVORXMZjbOzDaYWZuZ3ZDl8RPN7L7U40vMbGSqf6SZHTGzlal/Pytt+CinZ599Vj/+8Y/T+v5V3KgkSdGPvWXLlum+++5LJBZUL3J2Y+rq6tJXv/rVSO/fSXpbEuFAUnBlkqF9rYMHD2rSpEnJhYOCDVgwm9kgSZMlXabgmM5VZhY9tnONpL3u/gZJP5T03dBjm9z9nNS/60oUNyrgG9/4ho4dO9bXHibpC8mFAwX357oi0nfjjTem/Z7Q2MjZjeuXv/ylnnnmmVCPSfqPpMKBpOCioN9I67njjjsiVzBBLchlhHmspDZ33+zunZKmSpoQWWaCpLtTPz8g6RIz4yIKNWz58uW655570vomSXpZMuEg5D8VvvGqtGXLlox55mho5OwGdODAAd18882R3s8omMyFZF0n6XV9re7u7ixHAlDtcimYh0naFmq3p/qyLuPuXZL2SXpV6rFRZrbCzBaY2XuybcDMrjWzVjNr7ejoyGsHUHrd3d267rrr0q7A8JRmd/EAAA9fSURBVGZxZYxqMVpB+g27+eab1d7enkQ4qD5lz9kSebva3HTTTdq5c2eo5yQFwxxI3hAFQx0v+sMf/qCZM2cmEw4KkkvBnG3UIXotq/6W2SnpTHc/V9JXJN1rZqdmLOg+xd2b3b25qakph5BQTlOmTFFLS0ta3/cV3KYZ1eEbkl4Rah88eFD//M//nFQ4qC5lz9kSebuarFixIuN8k+DXNzyJcJDVxxQc/HnRxIkTdfjw4WTCQd5yKZjbFb76dvAXuKO/ZcxssILP8j3ufszdd0uSuy+TtEnBJWVRpXbu3KmvfS39ep3/oGAyJKrHGYqOV0i//e1vNWPGjCTCQXUhZzeQrq4uff7zn49cLWekMi9EiWSdIOl2hcuurVu3cgJgDcmlYG6RNNrMRpnZEAW3C5oeWWa6Xjxif4Wkee7uZtaUOgFFZvY6BUeTN5cmdJRaT0+PPvvZz2rfvn19fS+XdGtyISHG5yWdH+n73Oc+Jw6PNzxydgO55ZZbMo4ISreJM06q0XmSvpTW8/3vf19PPPFEMuEgLwMWzKn5bRMlzZG0XtI0d19rZpPMbHxqsV9IepWZtSk4DtR7GaOLJK02s1UKTiy5zt33lHonUBq33XZbxo0wviUO6lWrQZLuUPof8V/+8hd97nOf4w6ADYyc3TiWLl2qb37zm5Hef5B0eQLRIDeTFD6loKenR1dffXXaQBWqk1XbB2tzc7O3trYmHUbDWfX/2rv/6CjLK4Hj3xtiwNqFrIsUJDkLWKiypFiayEqXX7ZEqhE8brGIuyr0HA5QhXbr8UcRt6euVm219Viwx1Lh0EWx0l8KSrGnKxpbpFqEKUpjStGgXQryI8QS8uv2j2dCM5PJZDKZzPPOO/dzznsyeead5D6QuXNn3vd97q5dTJw4MWZ5sinAr7Bzl4PuTtwbm45WrVrF4sWLfYST90TkNVUt9x1HNlnezr7jx49TXl5ObW1th9FzgAjwEU9RmdRsBS6NGbnmmmtYv349tlhN9qWas63Tn+HQoUPMnj07plgeBKzDiuVcsIL4S0lg2bJlVFdX+wjHGNPHWltbmTdvXlyxDPAYVizngkpgWczIE088wYMPPugnHJMSK5jzXHNzM3PmzOHtt9+OGV8F/LOfkEwPnQGsx51v3q65uZmrrrqKd955x1NUxpi+snz5cp599tm40UVAlY9wTFruBcbFjNxyyy1s2bLFTzimW1Yw57HW1lauu+46tm3bFjO+BJjnJySTpo8Ca+PGDh06xKWXXmoXARoTIg8//DD33Xdf3OhFuIaNJncMAH5MxwVC29ramDNnDjt27PAWlemaFcx5SlVZtGgRGzZsiBmfCnzHT0iml/4ddz5zR3v37qWyspJjx475CMkYk0Fr165l6dKlcaPDgJ/iCjCTW8bgGnH+vRRraGhg5syZ7N6921tUJjErmPNQc3MzN9xwA6tXr44ZPw94CneI3+Sm/8YVzh29/vrrTJs2La4LmDEml6xcuZIFCxbEjZ6JK5bP9RCRyYyZwLdiRo4ePcr06dPZvn27n5BMQlYw55n6+nquvPJK1q1bFzNeAvwSd421yV0FuPOZZ8SN79q1i0mTJrF3714PURlj0tXW1sby5cu58cYb45aLLAJ+Bkz0FJnJnC8T32jmyJEjXHLJJTzzzDN+QjKdWMGcRyKRCOXl5Z0uFhmCK5ZH+AjKZFx/3GdO/xY3vn//fioqKnjyySc9RGWM6anDhw9z+eWXc8898X09+wFP4FZbMOFwF/ErZ5w8eZJZs2Zxxx130Nra6icsc5oVzHmgpaWFBx54gIkTJ/LWW2/F3DcSqAY+5iUy01fOArYAl8WNNzQ0MHfuXBYsWMDRo0c9RGaMScWmTZu48MILE6ya8CHgGVyDEhMegrtwc0Wne+6++26mTp1qRwg9s4I55F566SUuvvhibr75Zk6ePBlz3wTg17jetyZ8zsIdsJ2f4L41a9ZwwQUXsHbtWvvkwpgA2bdvH1dffTVXXHEF7777bty9g3HtpD7rITLT9wTXCXAl8eXZyy+/zPjx47nzzjs5ceKEj+DynhXMIaSqVFdXU1VVxZQpU0jUgesLuE+Wh2Y9OpNNZ+B6ID+CO+Oxo4MHDzJ//nzGjx/P448/TnNzc/YDNMYAUFtby5IlSzj//PN56qmnEuwxCdiJnbOcD5bg3hjFvkI3NTVx1113cd5553H//ffb6kdZZgVziLz//vusWrWKiooKJk+ezObNmzvtU4zr4Lcad321CT/BtTT4NTA2wf179uzh2muvZeTIkaxYsYKamprsBmhMnmpsbGTjxo1UVVUxZswYHnnkkQRvXPsBtwEv4C7PNvlhKu4NUvyJdW6N/VtvvZWSkhIWLlxIdXU1bW1tWY8w30jsVbf+lZeXa6JPRE1nLS0t7Ny5kxdeeIHNmzdTXV2d9PD654CHsU+V89kp4BvRrSnJfqNHj6ayspLKykqmTZvGwIEDsxNgCIjIa6pa7juObLK8nRpVpaamhhdffJHnnnuOrVu38sEHHyR5xATg+9GvJj8pbq3mZUDXTaiGDBnCjBkzTm/nnmtLDaYq1ZxtBXMOqK+vp66ujrq6Ot544w0ikQiRSIQ9e/bQ2NjY7eMnAffRedUEk7/exjU5+SEuHSdTUFDA6NGjKSsro6ysjHHjxjFixAhKSko455xzEJG+DziHWMFsGhsbOXDgAHV1ddTU1BCJRNi9ezeRSCTFw+gjcKsmzMMOBBvnBPAAbs3mZG+ynNLS0tM5u6ysjFGjRlFaWsrQoUMpLCzs62BzSkYLZhGZCTyEOza0WlXvjbu/P+5I/yeB94HPq+r+6H23406ZbQWWquovkv2usWPH6rp1606vN9kxvnTHgvAzmpqaTm+nTp3q9H19fT3Hjh3j+PHjMV/fe++9tE/wvwL4L9yBHStpTCJv4Z7Ya4C/pvH4oqIihg8fztlnn82gQYMoLi6O+TpgwACKioooKiqif//+nW4XFBQgIqe/dtzixxLt075lQqZ+TkVFhfeCOZs5G3qWt/s692Zq/5aWloS5uv12Q0NDTK5uv33w4EEOHz7c3T9ZFz4BfAWYQ+erDowB93R9FHe8uOfNqAoKChg2bBiDBw/ulK+Li4s588wzE+bq9q2wsLDHeTrRWCZkPWeratINl3D/CIzCPYN3AWPj9lkCfC96ey7wZPT22Oj+/XErmP0R6NfN71Pb0t/Ggv4P6H5Qtc22FLfjoD8AnQ5aiP+/4xzfXu0ur/blRpZztuXt3m7DFJYqvKrQpvhPB7blxNaksEnhaoV/UP9/xzm9pZSzU/lc/iKgVlX3AYjIBmA28EaHfWYDX4ve3gh8V1zpPxvYoKqngD+JSG305/0mhd9rUvARYDLuU+TLcK+QxvTUQGBBdDsBbAO2As8DtvJnzrGcHWgfxp0oNxnXFnkCdtqF6bkzgMujWzPwCi5jbwVeBVr8hRZSqRTMw4G6Dt8foPO6Nqf3UdUWETkO/FN0fHvcY4fH/wIRWQgsjH57Cvh9KsHngMFAusfmUnIQ92q3EbipL39RFuaSRWGZS1jmAeGai+8+QH2esyG0eTsLf4cNuKJmK4maVGRQmJ5TNpfgCcs8IMWcnUrBnOgkEU1xn1Qei6o+ijspBxF51ff5f5licwmmsMwlLPOA8M3FdwgJxjKasyGceTss8wCbS1CFZS5hmQeknrNTOQ50ACjt8H0J8F5X+4hIITAIOJLiY40xxmSO5WxjjMmwVArm3wKjRWSkiBThLhB5Om6fp4Hro7c/B/xK3ZUgTwNzRaS/iIzEdWHekZnQjTHGJGA52xhjMqzbUzKi57fdCPwCd/X1Y6q6R0S+jruy8Glc990fRi8QOYJL0ET3+xHuYpMW4Iuq2nVnDefR9KcTODaXYArLXMIyD7C5ZIyHnA3h+f8LyzzA5hJUYZlLWOYBKc4lcI1LjDHGGGOMCRJby8YYY4wxxpgkrGA2xhhjjDEmicAWzCJyk4j8QUT2iMj9vuPpLRG5WURURAb7jiVdIvJNEdkrIrtF5KciUuw7pp4QkZnRv6laEbnNdzzpEpFSEfk/EXkz+vxY5jum3hCRfiKyU0Q2+Y6lN0SkWEQ2Rp8jb4rIxb5jyibL2cFjOTsYwpazIT/zdiALZhGZjus49XFV/RfgW55D6hURKQVmAO/4jqWXngfGqerHgRrgds/xpExE+gErgc/i2v9eIyJj/UaVthbgK6p6AfCvwBdzeC4Ay4A3fQeRAQ8BW1T1fGA84ZhTSixnB5bl7GAIW86GPMzbgSyYgcXAvdH2rKjqXzzH01vfBm6hiwYAuUJVt6pqe7/N7bg1WnPF6XbBqtoEtLcLzjmq+mdV/V309gncEzxhN7agE5ESXG/X1b5j6Q0RGQhMwa0+gao2qeoxv1FlleXsALKcHQxhytmQv3k7qAXzGGCyiLwiIttEpMJ3QOkSkVnAu6q6y3csGbYAeM53ED2QqF1wziasdiIyAvgE8IrfSNL2HVxh0uY7kF4aBRwC1kQPU64WkbN8B5VFlrODz3J2AIQgZ0Oe5u1UWmP3CRH5JTA0wV3LcXH9I+7QRQXwIxEZpQFdA6+buXwVqMxuROlLNhdV/Xl0n+W4Q0zrsxlbL6Xc8jdXiMiHgR8DX1LVet/x9JSIVAF/UdXXRGSa73h6qRCYANykqq+IyEPAbcAKv2FljuXsYLKcnTtyPWdDfudtbwWzqn6mq/tEZDHwk2iy3SEibcBg3DuBwOlqLiJSBowEdokIuMNhvxORi1T1/7MYYsqS/b8AiMj1QBXw6aC+GHYhVC1/ReQMXOJdr6o/8R1Pmj4FzBKRy4ABwEAR+V9V/Q/PcaXjAHBAVds/NdqIS7yhYTnbcnaWWc4OprzN20E9JeNnwCUAIjIGKAIOe40oDaoaUdUhqjpCVUfg/nMmBDXxdkdEZgK3ArNU9a++4+mhVNoF5wRxr+Q/AN5U1Qd9x5MuVb1dVUuiz425uPbMuZh0iT6n60TkY9GhT+O65eULy9kBZDk7GMKSsyG/87a3T5i78RjwmIj8HmgCrs+xd8Zh9V2gP/B89NOX7aq6yG9IqemqXbDnsNL1KeA/gYiIvB4d+6qqPusxJgM3AeujL+77gPme48kmy9nBZDk7GCxnB1fKedtaYxtjjDHGGJNEUE/JMMYYY4wxJhCsYDbGGGOMMSYJK5iNMcYYY4xJwgpmY4wxxhhjkrCC2RhjjDHGmCSsYDbGGGOMMSYJK5iNMcYYY4xJ4m9S+kNsOYTj8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pylab as plt\n",
    "from scipy.stats import norm\n",
    "%matplotlib inline\n",
    "\n",
    "n, bins = 10**6, 100\n",
    "\n",
    "# transformation function\n",
    "def box_miller(u1,u2):\n",
    "    zz = np.sqrt(-2*np.log(u1))\n",
    "    z1 = zz * np.cos(2*np.pi*u2)\n",
    "    z2 = zz * np.sin(2*np.pi*u2)\n",
    "    return z1,z2\n",
    "\n",
    "def gaussian_th(x, mu=0, sig=1):\n",
    "    norm = ((2 * np.pi)**0.5 * sig)\n",
    "    g = np.exp(-((x - mu) * (x - mu)) / (2 * sig * sig))\n",
    "    return g/norm\n",
    "\n",
    "# uniformly distributed values between 0 and 1\n",
    "u1 = np.random.rand(n)\n",
    "u2 = np.random.rand(n)\n",
    "\n",
    "# run the transformation\n",
    "z1,z2 = box_miller(u1,u2)\n",
    "\n",
    "# plotting the values before and after the transformation\n",
    "fig, axs = plt.subplots(2,2, figsize=(12,8))\n",
    "res = [u1, u2, z1, z2]\n",
    "colors = \"rb\"\n",
    "for i in range(4):\n",
    "    row, col = i//2, i%2\n",
    "    axs[row, col].hist(res[i], bins=bins, color=colors[col], density=True)\n",
    "    #np.histogram(v, bins=bins)\n",
    "\n",
    "x = np.linspace(-6,6,150) \n",
    "for i in range(2):\n",
    "    axs[1,i].plot(x, norm.pdf(x,0,1), lw=4, c='k')\n",
    "    axs[1,i].axis((-6,6,0,0.45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pippo(u1,u2):\n",
    "    Z0 = (-2*np.log(u1))**0.5*np.cos(2*np.pi*u2)\n",
    "    return Z0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(np.e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log10(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}