Enkel simulering_inte_enkel_gustav.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Information\n",
    "\n",
    "Toy simulation for the HGTD.\n",
    "\n",
    "#### Outline:\n",
    "* Imports &nbsp;--&nbsp; imports used for the simulations and functions\n",
    "* Simulation model &nbsp;--&nbsp; Code for the simulation of activations in HGTD\n",
    "* Simulation examples &nbsp;--&nbsp; Example usages of the simulation model\n",
    "   * Sequential / Single core\n",
    "   * Parallel / Multi core\n",
    "* Tools / data functions &nbsp;--&nbsp; help functions to extract and process simulation data\n",
    "* Run &nbsp;--&nbsp; Main code to run, i.e. what simulation parameters, what to process, what to plot\n",
    "* Testing &nbsp;--&nbsp; to test small code parts while developing (not relevant to the simulation itself)\n",
    "\n",
    "\n",
    "#### Words:\n",
    "* BC &nbsp;=&nbsp; bunch crossings -- entire group of protons cross each other's path.  &nbsp;(happens about every 25ns)\n",
    "* collisions_per_BC &nbsp;=&nbsp; number of proton-collisions in each bunch crossing event.  &nbsp;(an average)\n",
    "* particles_per_collision &nbsp;=&nbsp; number of particles created in each single proton-collision.  &nbsp;(an average)\n",
    "* threshold &nbsp;=&nbsp; the threshold of the minimum voltage detected by the detector.  &nbsp;(can be calibrated)\n",
    "* activations &nbsp;=&nbsp; number of channels reporting a signal above the threshold voltage.  &nbsp;(both noise and particles)\n",
    "* deviation &nbsp;=&nbsp; a relative error between [the actual number of particles detected] and [the calculated estimate of particles from the number of activations]. \\\n",
    "&nbsp; Corresponds to an error in our model.\n",
    "* guess_voltage &nbsp;=&nbsp; guessed_expected_voltage &nbsp;=&nbsp; Our guess of the expected voltage of the signals from detected particles. \\\n",
    "&nbsp; Note that the real expected value of the voltage from particles detected is unknown, only approximated. \\\n",
    "&nbsp; But different guesses can be made to test the robustness to error in the unknown approximate value.\n",
    "* expected_count &nbsp;=&nbsp; the average or expected number of particles in hitting the detector each bunch crossing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment these, if the enviorment doesn't have those installed\n",
    "# !pip install scipy\n",
    "# !pip install matplotlib\n",
    "\n",
    "import ROOT\n",
    "import time\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import multiprocessing as mp\n",
    "from scipy.ndimage import uniform_filter\n",
    "\n",
    "\n",
    "global_counter = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": [
    "class noise_generator:\n",
    "    \"\"\"Class for simulation of noise. It's using information from the larger time window.\"\"\"\n",
    "    def __init__(self, threshold, channel_count):\n",
    "        global global_counter\n",
    "        low = 0\n",
    "        high = 1000\n",
    "        noise = ROOT.TF1(\"noise_gauss\", \"TMath::Gaus(x, 0, 5)\", low, high)  # we know it is centered at 0 and sigma=5.\n",
    "        #  activations:  channel_count * probability mass     (the *2 is because low=0 not -1000, but it is symetric so just double it)\n",
    "        self._mean = channel_count * noise.Integral(threshold, high)/(noise.Integral(low, high)*2)\n",
    "        \n",
    "        global_counter += 1\n",
    "        self._prng = ROOT.TRandom3(int(time.time()+global_counter))\n",
    "\n",
    "    def generate(self):\n",
    "        prng = self._prng\n",
    "        total = prng.Poisson(self._mean)\n",
    "        narrow_exact = 0\n",
    "        if total <= 5:\n",
    "            for i in range(total):\n",
    "                if prng.Uniform(25) < 3.125:\n",
    "                    narrow_exact += 1\n",
    "        else:\n",
    "            # faster to use numpy\n",
    "            #narrow_exact = (np.random.uniform(low=0, high=25, size=total) < 3.125).sum()\n",
    "            narrow_exact = (np.random.rand(total) < 0.125).sum()  # 3.125/25 = 0.125\n",
    "        \n",
    "        # using information from the large time-window to estimate noise in small time-window\n",
    "        estimate = 3.125 * (total - narrow_exact)/(25 - 3.125)\n",
    "        return narrow_exact, estimate\n",
    "\n",
    "    \n",
    "class particles_generator:\n",
    "    \"\"\"Class for simulation of particles/signal.\"\"\"\n",
    "    def __init__(self, threshold, real_expected_voltage):\n",
    "        global global_counter\n",
    "        self._threshold = threshold\n",
    "        self._real_expected_voltage = real_expected_voltage\n",
    "        #self._guess_expected_voltage = guess_expected_voltage\n",
    "        #self.guess_fraction_above = particles_generator.calculate_guess_fraction_above(threshold, guess_expected_voltage)\n",
    "        \n",
    "        global_counter += 1\n",
    "        self._prng = ROOT.TRandom3(int(time.time()+global_counter))\n",
    "    \n",
    "    @staticmethod\n",
    "    def calculate_guess_fraction_above(threshold, guess_expected_voltage, low=-100, high=1000):\n",
    "        \"\"\"Calculates the size of the fraction of signals expected to be detected above the threshold\"\"\"\n",
    "        signal_func = ROOT.TF1(\"sth\", \"TMath::Landau(x, %f, %f)\" % (guess_expected_voltage, guess_expected_voltage/4), low, high)  # sigma = voltage/4  approximation\n",
    "        guess_fraction_above = signal_func.Integral(threshold, high)/signal_func.Integral(low, high)\n",
    "        return guess_fraction_above\n",
    "    \n",
    "    def generate(self, expected_count):\n",
    "        prng = self._prng\n",
    "        threshold = self._threshold\n",
    "        expected_voltage = self._real_expected_voltage\n",
    "        signal_count = prng.Poisson(expected_count)\n",
    "        sigma = expected_voltage/4  # approximation from document\n",
    "        \n",
    "        above_threshold_count1 = 0\n",
    "        above_threshold_count2 = 0\n",
    "        \n",
    "        # unrolled loop\n",
    "        for i in range(signal_count//2):\n",
    "            above_threshold_count1 += prng.Landau(expected_voltage, sigma) > threshold\n",
    "            above_threshold_count2 += prng.Landau(expected_voltage, sigma) > threshold\n",
    "        if (signal_count&1) == 1:\n",
    "            above_threshold_count1 += prng.Landau(expected_voltage, sigma) > threshold\n",
    "                \n",
    "        return signal_count, above_threshold_count1 + above_threshold_count2  # actuall signal from particles, and the ones above threshold\n",
    "    \n",
    "    \n",
    "class simulation:\n",
    "    \"\"\"Class representing a simulation setup with chosen parameters\"\"\"\n",
    "    def __init__(self, threshold, real_expected_voltage, expected_count, bunch_crossing_count):\n",
    "        self._sim = self\n",
    "        self._threshold = threshold\n",
    "        self._real_expected_voltage = real_expected_voltage\n",
    "        self._expected_count = expected_count\n",
    "        self._bunch_crossing_count = bunch_crossing_count\n",
    "        \n",
    "        channel_count = int(3.6 * 1e6)  # channels of the HGTD: 3.6 million\n",
    "        self._particles_gen = particles_generator(threshold, real_expected_voltage)\n",
    "        self._noise_gen = noise_generator(threshold, channel_count)\n",
    "    \n",
    "    def run(self):\n",
    "        actual_particles_sum = 0\n",
    "        noise_subtracted_activations_sum = 0  # activations = above threshold\n",
    "\n",
    "        noise_gen = self._noise_gen\n",
    "        particles_gen = self._particles_gen\n",
    "        expected_count = self._expected_count\n",
    "        \n",
    "        for i in range(self._bunch_crossing_count):\n",
    "            noise_exact, noise_est = noise_gen.generate()\n",
    "            particles_total, particles_above = particles_gen.generate(expected_count)\n",
    "\n",
    "            actual_particles_sum += particles_total\n",
    "            noise_subtracted_activations_sum += particles_above + noise_exact - noise_est\n",
    "\n",
    "        # note, it is activations, so it is not scaled by fraction_above for the threshold.\n",
    "        return actual_particles_sum, noise_subtracted_activations_sum  # note: in reality we only get to know the value \"noise_subtracted_activations_sum\"\n",
    "    \n",
    "\n",
    "# one could make functions to call the parallel code, however there is no gain in speed and this is cleaner\n",
    "class simulation_parallel:\n",
    "    \"\"\"Class representing a simulation setup with chosen parameters, capable of running simulations in parallel\"\"\"\n",
    "    def __init__(self, pool, threshold, real_expected_voltage, expected_count, bunch_crossing_count):\n",
    "        bunch_crossing_count = bunch_crossing_count // nprocs  # divide bunch crossings equally\n",
    "        self._sim = simulation(threshold, real_expected_voltage, expected_count, bunch_crossing_count)\n",
    "        self._pool = pool\n",
    "        \n",
    "    def run(self):\n",
    "        # parallel simulation\n",
    "        func = self._sim.run\n",
    "        it = [[]]*nprocs # it takes no arguments\n",
    "        results = np.asarray(self._pool.starmap(func, it)).sum(axis=0)\n",
    "        \n",
    "        # actual_particles_sum, noise_subtracted_activations_sum\n",
    "        return results[0], results[1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulations example\n",
    "\n",
    "### Sequential / Single core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_simulate_single_core():\n",
    "    \"\"\"Method for test running the simulation with some example parameters\"\"\"\n",
    "    my_sim = simulation(threshold=25,\n",
    "                        real_expected_voltage=100,\n",
    "                        expected_count=30,\n",
    "                        bunch_crossing_count=1000000)\n",
    "    t = time.time()\n",
    "    print(my_sim.run())  # run the model\n",
    "    print(time.time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallel / Multi core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_simulate_multi_core(pool):\n",
    "    \"\"\"Method for test running the simulation with some example parameters.\n",
    "        Test runing the parallel simulation class. It takes a multiprocessing.Pool as a parameter\"\"\"\n",
    "    my_sim = simulation_parallel(pool=pool,\n",
    "                        threshold=25,\n",
    "                        real_expected_voltage=100,\n",
    "                        expected_count=30,\n",
    "                        bunch_crossing_count=1000000)\n",
    "    t = time.time()\n",
    "    print(my_sim.run())  # run the model\n",
    "    print(time.time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tools / data functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.linspace(start, stop, num=50)  # includes the stop value\n",
    "\n",
    "def create_simulation(threshold, real_expected_voltage, expected_count, bunch_crossing_count):\n",
    "    \"\"\"The current way of creating a simulation. Can swap between single and parallel\"\"\"\n",
    "    #return simulation(threshold, real_expected_voltage, expected_count, bunch_crossing_count)\n",
    "    return simulation_parallel(pool, threshold, real_expected_voltage, expected_count, bunch_crossing_count)\n",
    "\n",
    "\n",
    "def calculate_deviation_from_guess_voltage(threshold, guess_expected_voltage, actual_particles_sum, noise_subtracted_activations_sum):\n",
    "    \"\"\"Formula for calculating the relative deviation between the actual signal and the estimated signal by the provied guess volatge\"\"\"\n",
    "    scale_factor = particles_generator.calculate_guess_fraction_above(threshold, guess_expected_voltage)\n",
    "    estimate = noise_subtracted_activations_sum / scale_factor\n",
    "    return abs( estimate / actual_particles_sum - 1 )\n",
    "\n",
    "def calculate_absolute_difference_from_guess_voltage(threshold, guess_expected_voltage, actual_particles_sum, noise_subtracted_activations_sum):\n",
    "    \"\"\"Formula for calculating the absolute difference between the actual signal and the estimated signal by the provied guess volatge\"\"\"\n",
    "    scale_factor = particles_generator.calculate_guess_fraction_above(threshold, guess_expected_voltage)\n",
    "    estimate = noise_subtracted_activations_sum / scale_factor\n",
    "    return abs( estimate - actual_particles_sum )\n",
    "\n",
    "\n",
    "def compute_deviation_list_from_guess_voltage_list(sim, guess_expected_voltage_list):\n",
    "    \"\"\"Computes the deviation for each provided guess voltage in the list\"\"\"\n",
    "    actual_particles_sum, noise_subtracted_activations_sum = sim.run()\n",
    "    threshold = sim._sim._threshold\n",
    "    deviation_list = [calculate_deviation_from_guess_voltage(threshold, guess_expected_voltage, actual_particles_sum, noise_subtracted_activations_sum) for guess_expected_voltage in guess_expected_voltage_list]\n",
    "    return deviation_list\n",
    "\n",
    "def compute_the_best_threshold_and_deviation_in_terms_of_least_deviation_no_post_processing(threshold_list, guess_expected_voltage_list, number_of_simulation_samples, real_expected_voltage, expected_count, bunch_crossing_count):\n",
    "    \"\"\"Fast way to compute the best threshold from least deviaiton. No post processing, just directly finding the least deviation. \n",
    "        Unreliable if threshold steps are small (<5).\n",
    "        An alternative function is: compute_the_best_threshold_and_deviation_in_terms_of_least_average_deviation_sum()\"\"\"\n",
    "    cur_deviation = 99\n",
    "    cur_threshold_index = 99\n",
    "    \n",
    "    # Create simulation,  run the simulation,  sum over the deviaitons for each guess,  find the threshold with least deviaiton sum.\n",
    "    for index, threshold in enumerate(threshold_list):\n",
    "        sim = create_simulation(threshold, real_expected_voltage, expected_count, bunch_crossing_count)\n",
    "        actual_particles_sum, noise_subtracted_activations_sum = sim.run()\n",
    "        deviation = 0\n",
    "        for guess_expected_voltage in guess_expected_voltage_list:\n",
    "            deviation += calculate_deviation_from_guess_voltage(threshold, guess_expected_voltage, actual_particles_sum, noise_subtracted_activations_sum)\n",
    "        if deviation <= cur_deviation:\n",
    "            cur_deviation = deviation\n",
    "            cur_threshold_index = index\n",
    "    \n",
    "    return threshold_list[cur_threshold_index], cur_deviation  # best threshold for guess voltages and its corresponding deviation\n",
    "    \n",
    "\n",
    "def compute_average_deviation_list_from_guess_voltage_list(sim, guess_expected_voltage_list, number_of_simulation_samples):\n",
    "    \"\"\"Computes the deviation for each provided guess voltage in the list, as an average over simulations\"\"\"\n",
    "    size = len(guess_expected_voltage_list)\n",
    "    deviation_per_guess_list = np.zeros(size)\n",
    "    for i in range(number_of_simulation_samples):\n",
    "        deviation_per_guess_list += np.fromiter(compute_deviation_list_from_guess_voltage_list(sim, guess_expected_voltage_list), dtype=np.float, count=size)\n",
    "    return np.divide(deviation_per_guess_list, number_of_simulation_samples)\n",
    "\n",
    "def compute_the_best_threshold_and_deviation_in_terms_of_least_average_deviation_sum(threshold_list, guess_expected_voltage_list, number_of_simulation_samples, real_expected_voltage, expected_count, bunch_crossing_count):\n",
    "    \"\"\"Computes the best threshold in terms of least 'average deviation over all the different guess voltages'.\n",
    "        There is also the function which does no post processing of the data: compute_the_best_threshold_and_deviation_in_terms_of_least_deviation_no_post_processing()\n",
    "    \"\"\"\n",
    "    deviation_per_threshold_list = np.zeros(len(threshold_list))\n",
    "    \n",
    "    # Performance note: this could be made a little bit faster, if we didn't need to blur the result. You could easily check minimum while looping.\n",
    "    #  see the alternative function without post processing.\n",
    "    for index, threshold in enumerate(threshold_list):\n",
    "        sim = create_simulation(threshold, real_expected_voltage, expected_count, bunch_crossing_count)\n",
    "        # sum over all deviations for each guess voltage,  so this is a messure of how good is this threshold across all guesses\n",
    "        deviation_per_threshold_list[index] = compute_average_deviation_list_from_guess_voltage_list(sim, guess_expected_voltage_list, number_of_simulation_samples).sum()\n",
    "    \n",
    "    # Blur,  assumes deviations follow a smooth curve else this ruins the results\n",
    "    size = 5  # default is 3, not sure what size to use 3 or 5.\n",
    "    deviations_t_blur = uniform_filter(deviation_per_threshold_list, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't matter much\n",
    "    \n",
    "    # Get the minimum\n",
    "    index = np.argmin(deviations_t_blur)\n",
    "    threshold = threshold_list[index]\n",
    "    deviation = deviation_per_threshold_list[index] / len(guess_expected_voltage_list)\n",
    "    \n",
    "    return threshold, deviation  # best threshold for guess voltages and its corresponding deviation\n",
    "\n",
    "\n",
    "def caluclate_SNR(expected_signal_voltage, noise_sigma):\n",
    "    \"\"\"Formula for calculating SNR\"\"\"\n",
    "    return expected_signal_voltage / noise_sigma\n",
    "\n",
    "def calculate_expected_count(collisions_per_BC = 100, particles_per_collision = 75, ratio_of_particles_reaching_HGTD=0.1):\n",
    "    \"\"\"Formula for calculating the expected number of real particles detected\"\"\"\n",
    "    return collisions_per_BC * particles_per_collision * ratio_of_particles_reaching_HGTD\n",
    "\n",
    "\n",
    "def plot_the_data(title, prefix_filename, xlabel, ylabel, x_list, y_list, make_integer_axis = False, make_margin = False):\n",
    "    \"\"\"Help function for plotting list of x,y values as a graph, and saving the graph.\"\"\"\n",
    "    timestr = time.strftime(\"%Y%m%d-%H%M%S\")\n",
    "    FILENAME = prefix_filename + \"_(\" + timestr + \").pdf\"\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    plt.title(title)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    if make_integer_axis:\n",
    "        plt.locator_params(axis=\"both\", integer=True) # integer axis\n",
    "    if make_margin:\n",
    "        plt.margins(y=1)  # add margin to y\n",
    "    \n",
    "    plt.plot(x_list, y_list)\n",
    "    plt.show()\n",
    "    ax = fig.add_subplot(111)\n",
    "    fig.savefig(FILENAME)\n",
    "    #fig.savefig(FILENAME, dpi=200)  # use for png\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoh0lEQVR4nO3deXxcdbnH8c+TpW3aTFvaphPolpYmgVKgpS0U8aosIiKLIopKWcumXEFEuYAL4IaioqBXEC2rRUBANgUuIDsCXSiU7gstXdN9SZesz/1jTiDUJp20OXMmc77v1yuvzpzZHoaT75z8zjO/n7k7IiISH3lRFyAiIpml4BcRiRkFv4hIzCj4RURiRsEvIhIzCn4RkZhR8IukycxmmNmnWrn9BTM7L3MVieyegqgLEMkWZlbd7GpXoAZoCK5f6O4HNLvvtcBQdx+XuQpF2oeCXyTg7sVNl81sEXCeuz8bXUUi4dBQj0iazGyRmR1jZscBVwOnmVm1mb3dwv3PNbNZZrbezJ42s0GZrVhk5xT8Im3k7k8BPwPud/didz94x/uY2cmkPhxOAUqAl4G/ZrRQkRYo+EXCcRFwvbvPcvd6Uh8UI3TUL9lAwS8SjkHATWa2wcw2AOsAA/pFWpUIOrkrsrt2Na3tEuCn7j4xE8WItIWO+EV2TxVQZmYt/Q7dClxlZgcAmFkPM/tSxqoTaYWCX2T3/C34d62ZTd3xRnf/O/AL4D4z2wS8C3w2g/WJtMi0EIuISLzoiF9EJGYU/CIiMaPgFxGJGQW/iEjMdIg+/j59+nhZWVnUZYiIdChTpkxZ4+4lO27vEMFfVlbG5MmToy5DRKRDMbPFO9uuoR4RkZhR8IuIxIyCX0QkZhT8IiIxo+AXEYmZ0ILfzLqY2Ztm9raZzTCz64Ltg83sDTObb2b3m1mnsGoQEZH/FOYRfw1wVLAs3QjgODMbS2rGwt+4+1BgPTA+xBpERGQHoQW/p1QHVwuDHweOAh4Mtt8FfD6sGv41u4o/vDA/rKcXEemQQh3jN7N8M5sGrAKeARYAG4I1SAGW0sJSdGZ2gZlNNrPJq1ev3q3Xf23+Wm56dh4NjZp6WkSkSajB7+4N7j4C6A8cCuzXhsfe5u6j3X10Scl/fOM4LRXJBDX1jSxZt3W3Hi8ikosy0tXj7huA54HDgZ5m1jRVRH9gWVivW1GaAGBO1eawXkJEpMMJs6unxMx6BpeLgE8Ds0h9AJwa3O0s4NGwaijvWwzA3JUKfhGRJmFO0rY3cJeZ5ZP6gHnA3Z8ws5mk1iH9CfAWMCGsArp1LmBAryId8YuINBNa8Lv7O8DInWxfSGq8PyMqkwnmKvhFRD6Q89/cLU8mWLh6C7X1jVGXIiKSFXI++CuTCeobnUVrt0RdiohIVsj54K9IBp09OsErIgLEIPiHlHQjP880zi8iEsj54O9SmE9Z76464hcRCeR88ENquGfequpd31FEJAZiE/yL1m5he11D1KWIiEQuFsFfWZrAHebrqF9EJB7Br84eEZEPxSL4y3p3pVN+njp7RESISfAX5OcxpKSbgl9EhJgEP6TG+edWaYxfRCQ2wV+RTLBswzY2b6+LuhQRkUjFJvgrgxO8OuoXkbiLTfBXfBD8GucXkXiLTfD336uIosJ8Bb+IxF5sgj8vz6hIFiv4RST2YhP8kBrumbNSY/wiEm+xCv7K0gRrqmtYW10TdSkiIpGJVfCXq7NHRCRewd/U0jlvlcb5RSS+YhX8ye6d6d6lQJO1iUisxSr4zSyYukHBLyLxFavgh6bOns24e9SliIhEIpbBv2l7PVWb1NkjIvEUy+AHTd0gIvEVw+AvBhT8IhJfsQv+3sWd6VPcWZ09IhJbsQt+gMpSzdkjIvEVy+Av75tg3qpqGhvV2SMi8RPL4K8sTbC1toFlG7ZFXYqISMaFFvxmNsDMnjezmWY2w8wuDbaPMLPXzWyamU02s0PDqqElTZ09GucXkTgK84i/Hrjc3YcBY4GLzWwYcANwnbuPAH4YXM+ops6eORrnF5EYKgjrid19BbAiuLzZzGYB/QAHugd36wEsD6uGliS6FNKvZ5FO8IpILIUW/M2ZWRkwEngD+BbwtJn9itRfHB9r4TEXABcADBw4sN1rKk8Wa3pmEYml0E/umlkx8BDwLXffBHwduMzdBwCXARN29jh3v83dR7v76JKSknavqzKZYMGqauobGtv9uUVEslmowW9mhaRCf6K7PxxsPgtouvw3IOMndyF1gre2oZFFa7dG8fIiIpEJs6vHSB3Nz3L3G5vdtBz4ZHD5KGBeWDW0prJUc/aISDyFOcZ/BHAGMN3MpgXbrgbOB24yswJgO8E4fqYN7VuMWaql8/gD946iBBGRSITZ1fMKYC3cPCqs101Xl8J8BvXqqmUYRSR2YvnN3SZNi7KIiMRJrIO/sjTBorVb2V7XEHUpIiIZE+vgr0gmaGh0Fq7eEnUpIiIZE+vgV2ePiMRRrIO/rHc3CvJMwS8isRLr4O9UkMeQkm4KfhGJlVgHPwSdPQp+EYmR2Ad/ZTLBknXb2FJTH3UpIiIZEfvgrwhO8M5bpZk6RSQeFPxJdfaISLzEPvgH9upK54I85uobvCISE7EP/vw8ozxZrBO8IhIbLU7SZmbTSS2TuFPuflAoFUWgIpng1flroi5DRCQjWpud84Tg34uDf+8J/j09vHKiUZlM8PDUZWzcWkeProVRlyMiEqoWh3rcfbG7LwY+7e5XuPv04OdK4NjMlRi+D07waopmEYmBdMb4zcyOaHblY2k+rsNoaunUFM0iEgfpLMQyHrjdzHqQWlhlPXBuqFVl2D49ulDcuUAtnSISC7sMfnefAhwcBD/uvjH0qjLMzKhIFuuIX0RiobWunm+3sB2AHRZQ7/AqSxM89e5K3P2D/0YRkVzU2lh9Yhc/OaW8b4L1W+tYU10bdSkiIqFq8Yjf3a/LZCFRa74oS0mic8TViIiEZ5fdOWbW38z+bmargp+HzKx/JorLpKaWTo3zi0iuS6ct8w7gMWCf4OfxYFtO6VPciV7dOqmzR0RyXjrBX+Lud7h7ffBzJ1AScl0Z19TZo+AXkVyXTvCvNbNxZpYf/IwD1oZdWBQqkgnmVlXj3uIURSIiHV46wX8u8GVgJbACOBU4J8yiolKRTFBdU8/yjdujLkVEJDTpfIFrMXBSBmqJ3AedPSs3069nUcTViIiEY5fBb2YlwPlAWfP7u3tOTdsAUNE36Oyp2syR+/WNuBoRkXCkM1fPo8DLwLNAQ7jlRKtH10JKu3fRCV4RyWnpBH9Xd/+f0CvJEuXq7BGRHJfOyd0nzOz4tj6xmQ0ws+fNbKaZzTCzS5vd9k0zmx1sv6Gtzx2mymSCeVXVNDSqs0dEclNrk7RtJrX0ogFXm1kNUBdcd3fvvovnrgcud/epZpYAppjZM0ASOBk42N1rzCyrBtMrShPU1Dfy/rqtDO7TLepyRETaXWsrcCXcvXvwb567FzW7vqvQx91XuPvU4PJmYBbQD/g68HN3rwluW9U+/ynto1JTN4hIFti8vY7fPjuXrbX17f7c6czVc4SZdQsujzOzG81sYFtexMzKgJHAG0AF8F9m9oaZvWhmY1p4zAVmNtnMJq9evbotL7dHypPFAMzTOL+IRKCx0XlwylKO/NWL/PbZebw0d027v0Y6Y/y3AFvN7GDgcmABHy68vktmVgw8BHzL3TeRGl7qBYwFvgs8YDuZAN/db3P30e4+uqQkczNEdO1UwIBeRcxR8ItIhk1bsoEv3PIa3/nb2wzoVcSjFx/BccNL2/110unqqXd3N7OTgd+7+wQzG5/Ok5tZIanQn+juDweblwIPe2pehDfNrBHoA2TusH4XKpMJdfaISMas2rydXz41h79NWUpJojO//tLBfGFkP/LywlkUKp3g32xmVwFnkBqiyQMKd/Wg4Ch+AjBrh9W6HgGOBJ43swqgE9D+f8vsgYpkghfmrKa2vpFOBTm1rryIZJHa+kbuem0RNz03j5r6Bi785BC+eVQ5xZ3Tiebdl86znwZ8DTjX3VcG4/u/TONxR5D6sJhuZtOCbVcDt5NavP1doBY4y7NsVrTK0gT1jc57a7Z8MI2DiEh7emHOKn70xEwWrt7CUfv15fuf258hJcUZee105upZaWYPAeXBpjXA39N43CukWj93ZlzaFUagaVGWuVWbFfwi0q4WrdnCT/4xk2dnrWJwn27cfvZojtovmdEa0pmr53zgAlInZPcl1ZJ5K3B0uKVFZ0hJN/LzTOP8ItJuttTU8/vn5zPh5fcozDeu/Ox+nHNEGZ0L8jNeSzpDPRcDh5JqxcTd52Xbl67aW+eCfMp6d1Uvv4jsMXfnkWnL+PmTs6naVMMph/TjyuP2o2/3LpHVlE7w17h7bVPHpZkVkPpGb06rLE0wc/mmqMsQkQ5s+tKNXPv4DKYsXs9B/Xtwy7hRHDJwr6jLSiv4XzSzq4EiM/s08A1S6+7mtIpkgiffXcm22gaKOmX+TzER6bjWVNfwq6fncP/kJfTu1okbvngQp47qH1p7ZlulE/z/A5wHTAcuBP4J/DnMorJBZTKBOyxYXc3wfj2iLkdEOoC6hkbu+fdifvPsXLbVNjD+iMFcckw53bvssgM+o1oNfjPLB2a4+37AnzJTUnYobzZnj4JfRHbllXlruO7xGcxbVc1/lffhmhOHMbRvdnYFthr87t5gZnPMbKC7v5+porJBWe+udMrPU2ePiLRqybqt/OQfM3l6RhUDe3XlT2eO5pj9+7KTmWiyRjpDPXsBM8zsTWBL00Z3z+l1eAvy89i3b7Hm7BGRndpaW88tLyzgjy8tJN+M736mkvEfH0yXwuw/J5hO8P8g9CqyVGWymDffWxd1GSKSRdydx99ZwfX/nMWKjds5ecQ+XPnZ/di7R1HUpaUtnW/uvpiJQrJReTLBI9OWs3l7HYksOzkjIpk3Y/lGrntsJm8uWscB+3Tn5q+OZExZr6jLarN0vrl7CvALoC+pKRjSXYGrw6v8YOqGakYNir73VkSisW5LLb/+vzn89c336dm1E9efciBfHj2A/Cxpz2yrdIZ6bgBOdPdZYReTbZrm6ZlbtVnBLxJD9Q2NTHzjfW58Zi7VNfWceXgZlx1TQY+uHXsEIJ3gr4pj6AP061lE1075mrpBJIZeW7CG6x6byZyqzRwxtDfXnHjABxM4dnStLbZ+SnBxspndT2oe/Zqm25strJKz8vKMci3KIhIrS9dv5fp/zuYf01fQf68ibh13CJ85oDSr2zPbqrUj/hObXd4KHNvsugM5H/wAFX2LeX5O1iwOJiIh2VbbwK0vLuDWFxdgBpcdU8GFnxzSIdoz26rF4Hf3cyC12Lq7v9r8NjM7IuzCskVlaYK/TVnK2uoaehd3jrocEWln7s6T767kp/+YxbIN2zjhoL256vj96dez47RntlU6Y/y/Aw5JY1tOqmjW2XO4gl8kp8xeuYlrH5vB6wvXsV9pgvsuGMvYIb2jLit0rY3xHw58DCgxs283u6k7kHt/+7SgeWfP4fvm/g4hEgcbttbym2fmcs/ri+leVMiPPz+cr44ZQEF+PNbYbu2IvxNQHNyn+ansTcCpYRaVTfomOtOjqFAneEVyQEOj89c33+fX/zeHjdvqGDd2EN/+dAU9u3aKurSMam2M/0VSc/Hf6e6LM1hTVjEzKpLFCn6RDu6NhWu59vGZzFqxibFDenHNiQew/945/z3UnUpnyobYhn6TimSCx99ejrvnVEuXSBws37CN65+czeNvL2efHl34368dwvEH5lZ7Zlulc3I39ipLE0x8o56qTTWU9ohunUwRSd/2ugb+9NJC/vDCAhrdufToci765L5aUQ8Ff1qaOnvmVG1W8ItkOXfn6RlV/PSfM1mybhufHV7K1cfvz4BeXaMuLWukdQrbzE5o7Xquawr+eRrnF8lq86o2c8aEN7noL1MoKszn3vMO45ZxoxT6O0j3iH8M8EQr13Nar26d6FPcWXP2iGSpjdvq+O2zc7n734vp1imf6046gNMPGxib9sy2Siv43f2a1q7HQWWpOntEsk1Do/PA5CX88uk5rN9ay1cPHch3jq2kV7d4tWe2VTqTtO1UHCZpa64imeC+N5fQ2OjkddA5uEVyyeRF67j28Rm8u2wTY8r24poTD2V4vx5Rl9UhpDtJ245iM0lbk8pkgm11DSxdv42BvTVeKBKVlRu38/MnZ/HItOWUdu/CTV8ZwUkH7xPr9sy22uUkbZJS0WzqBgW/SObV1Dfw55ff43+fn099o/PfRw7lG0fuS9dOak5sq3SWXkwCPwP2cffPmtkw4HB3nxB6dVmkvG8xkGrpPGZYMuJqROLD3Xlu1ip+/I+ZLF67lWOHJfn+54bpAGwPpPNReSdwB/C94Ppc4H4gVsGf6FJIv55FOsErkkHzV1Xzoydm8tLc1QztW8w94w/lv8pLoi6rw0sn+Pu4+wNmdhWAu9ebWcOuHmRmA4C7gSSpcwK3uftNzW6/HPgVUOLua3ar+gyrSBarpVMkAzZtr+N3z83jjlcXUVSYzw9OGMaZhw+iUO2Z7SKd4N9iZr1JhTdmNhbYmMbj6oHL3X2qmSWAKWb2jLvPDD4UjgXe393Co1BRmuDV+Wupb2hUf7BICBobnQenLuWGp2azdkstp40ewHc+U0kfrYXRrtIJ/m8DjwH7mtmrQAlpTMvs7iuAFcHlzWY2C+gHzAR+A1wBPLqbdUeiMpmgtqGRRWu3MjQY8xeR9rG9roFz75zEawvWcsjAntxx9qEc2F/tmWFIZ3bOqWb2SaASMGCOu9e15UXMrAwYCbxhZicDy9z97dbar8zsAuACgIEDB7bl5ULz4WpcmxX8Iu2osdH57oPv8NqCtfzk88M5/bCBas8M0S7HK8ysC3AJ8GPgOuDiYFtazKwYeAj4Fqnhn6uBH+7qce5+m7uPdvfRJSXZcTJnaN9izNA4v0g7+82zc3n87eVccVwl48YOUuiHLJ2B6ruBA0its/v74PI96Ty5mRWSCv2JwTd99wUGA2+b2SKgPzDVzErbXnrmdSnMp6x3N3X2iLSjB6cs5Xf/ms9powfw9U/uG3U5sZDOGP9wdx/W7PrzZjZzVw+y1Ef2BGCWu98I4O7Tgb7N7rMIGN1Runog6OxR8Iu0i38vWMtVD7/Dx/btzU++MFxH+hmSzhH/1KCTBwAzOwyYnMbjjgDOAI4ys2nBz/G7WWfWqEwmWLx2K9vrdtnRKiKtWLC6mov+MoWBvbpyy+mj1KqZQa1N0jadVAtnIfCamb0fXB8EzN7VE7v7K6ROBrd2n7K2FJsNypMJGhqdhau3MGyfeK7XKbKn1m2p5dw7J1GQZ9xx9qH06FoYdUmx0tpQT6wWW0lXZbM5exT8Im1XU9/AhfdMZsXG7fz1/LGaeiECrU3S9pFF1s2sLxD7dQfLenejMN80zi+yG9ydKx58h0mL1vO7r45k1KC9oi4pltJp5zzJzOYB7wEvAouAJ0OuK2t1KshjSJ9i5qqlU6TNfvvsPB6dtpzvHFvBiQfvE3U5sZXO2ZQfA2OBue4+GDgaeD3UqrJcRWmCuasU/CJt8chby7jpuXl88ZD+XHzk0KjLibV0gr/O3dcCeWaW5+7PA6NDriurVfQtZsm6bWypqY+6FJEO4c331nHFg+8wdkgvrj/lQLVtRiyd4N8QfPv2JWCimd0EbAm3rOzWtCjLvFXVEVcikv0WrdnChfdMpv9eRdw6bhSdCtS2GbV0/g+cDGwDLgOeAhbQ+rKMOa+yac4ejfOLtGrD1lTbJsAd54yhZ1ctgp4N0pmkrfnR/V0h1tJhDOjVlS6FeersEWlFbX0jF94zhaXrtzHx/MMY1Ltb1CVJoLUvcG0mmIN/x5sAd/fYNrHn5xnlfROas0ekBe7OlQ+/wxvvreOmr4xgTFmvqEuSZlrr409kspCOpjxZzKvzO8wUQyIZ9ft/zefhqcu47JgKTh7RL+pyZAc6y7KbKpMJqjbVsGFrbdSliGSVR6ct49fPzOWUkf245Gi1bWYjBf9uqvhg6gZ19og0mbJ4Hd998B0OHdyL67+ots1speDfTZXNVuMSEVi8dgvn3z2Ffj2L+OO4UXQuyI+6JGlBOlM2fNPMNKHGDvbu0YVE5wIFvwiwcWsd59w5iUZ3bj97DHt1U9tmNkvniD8JTDKzB8zsONPfbgCYGeXJYi3DKLFXW9/IRX+ZwpJ1W/njuFEM7qO2zWy3y+B39+8D5aRW0zobmGdmPzOz2K+RVlmaaul031nXq0juc3e+9/fp/HvhWn7xxYM4bEjvqEuSNKQ1xu+pZFsZ/NQDewEPmtkNIdaW9SqSCdZvrWN1dU3UpYhE4g8vLOBvU5ZyyVFDOeWQ/lGXI2lKZ4z/UjObAtwAvAoc6O5fB0YBXwy5vqzWdIJ3njp7JIaeeGc5v3x6DicdvA+Xfboi6nKkDdJZbL0XcMqOC7O4e6OZxXqVrqaWzjkrN3PE0D4RVyOSOVPfX8+3H3ib0YP24oZTD1LbZgeTzlw917Ry26z2Ladj6VPcmV7dOqmzR2JlybqtnH/XZEq7d+GPZ4yiS6HaNjsa9fHvoYpksSZrk9jYuC3VtlnX0MjtZ4+hd3HnqEuS3aDg30OVyQTzqqrV2SM5r66hkW9MnMKiNVu49YxRDO1bHHVJspsU/HuoojRBdU09yzduj7oUkdC4Oz945F1enb+W6085kI/tq3NaHZmCfw9pURaJgz++tJD7Ji3h4iP35UujB0RdjuwhBf8eKg+CX+P8kquenL6Cnz85mxMO2pvLP10ZdTnSDhT8e6hHUSGl3bvoiF9y0rQlG/jW/dMYObAnv/rSweTlqW0zFyj420FFaUJH/JJzlq7fynl3TaZv98786czRatvMIQr+dlCZLGb+qmoaGtXZI7lh0/Y6zr1zEjX1Ddxx9hj6qG0zpyj420FFMkFNfSPvr9sadSkie6yuoZGLJ05l4eot3DpuFEP7ahXWXKPgbwcVyQ+nbhDpyNydax6bwcvz1vDTLwzXVCQ5KrTgN7MBZva8mc00sxlmdmmw/ZdmNtvM3jGzv5tZz7BqyJTyZOqLLJq6QTq6P7/8Hve+8T4XfXJfThszMOpyJCRhHvHXA5e7+zBgLHCxmQ0DngGGu/tBwFzgqhBryIiunQoY2KurTvBKh/b0jJX87MlZHH9gKVd8Rm2buSy04Hf3Fe4+Nbi8GZgF9HP3/3P3+uBurwM5MYl3RTLBPAW/dFDTl27k0vve4qD+PbnxyyPUtpnjMjLGb2ZlwEjgjR1uOhd4soXHXGBmk81s8urVq0OucM9VlhazcPUWausboy5FpE2WbdjGuXdNone3zvxZbZuxEHrwm1kx8BDwLXff1Gz790gNB03c2ePc/TZ3H+3uo0tKSsIuc49VJBPUNzrvrdkSdSkiadu8vY7xd05ie20Dd5wzhpKE2jbjINTgN7NCUqE/0d0fbrb9bOAE4HTPkWktKzR1g3Qw9Q2N/Pe9bzFvVTV/GHfIB/uw5L4wu3qM1ALts9z9xmbbjwOuAE5y95xpfB9S0o38PNPUDdIhuDvXPT6TF+eu5scnD+e/yrP/r2ppP+ksvbi7jgDOAKab2bRg29XAzUBn4JlgubbX3f2iEOvIiM4F+Qzu000tndIh3PHqIu55fTEXfGIIXztMbZtxE1rwu/srwM5aA/4Z1mtGrTKZYMbyjVGXIdKqZ2dW8eN/zOQzByS58rj9oi5HIqBv7raj8mQxi9dtZVttQ9SliOzUu8s2csl9b3Fgvx789rSRatuMKQV/O6pMJnCH+auqoy5F5D+s2LiN8XdNomdRIX8+czRFndS2GVcK/nZUUarOHslOW2rqGX/nZLbUNHD7OWPo271L1CVJhMI8uRs7g3p1pVNBnr7BK1mlodG55K9vMadqMxPOGs1+pd2jLkkipiP+dlSQn8e+JcU64pes8uMnZvLc7FVce9IBfKqyb9TlSBZQ8LezymSxevkla9z12iLufG0R4z8+mDPGDoq6HMkSCv52VlGaYPnG7WzaXhd1KRJz/5pdxXWPz+CY/ZNcffz+UZcjWUTB384qg6+9z6tSZ49EZ+byTXzz3rcYtk93bv7qCPLVtinNKPjbWdN8J/oGr0SlatN2xt81ie5FhUw4awxdO6mHQz5Kwd/O+vUsomunfC3DKJHYWlvP+LsmsWlbHRPOGkNSbZuyEzoUaGd5eUZ5MqEjfsm4VNvmNGYu38SEs8YwbB+1bcrO6Yg/BJXJYgW/ZNzP/jmLZ2dVcc2JB3DkfmrblJYp+ENQkUywprqWtdU1UZciMXHP64uZ8Mp7nP2xMs76WFnU5UiWU/CHoLK06QSvOnskfC/MWcW1j83g6P368oMThkVdjnQACv4QqLNHMmX2yk38971vUZlMcPNXR6ptU9Ki4A9B30RnehQVauoGCdWqTds5945JdOucz4SzR9Ots3o1JD3aU0JgZlQmE5q6QUKztbae8+6ezIZtdTxw4eHs3aMo6pKkA9ERf0gqSlOdPTmylrxkkcZG57L7p/Huso3c/JWRDO/XI+qSpINR8IekMplg0/Z6qjaps0fa18+fms3TM6r4/ueGccywZNTlSAek4A9JeVKLskj7u/eN97ntpYWcefggzjmiLOpypINS8Ifkg84ejfNLO3lp7mp+8Oi7fKqyhB+eMAwzdfDI7lHwh6RXt06UJDrriF/axZyVm7l44lTK+xbz+68dQkG+fnVl92nvCVFlMqFlGGWPrd5cw7l3TqKoUz63nz2GYrVtyh5S8IeoIplgblU1jY3q7JHds622gfPunsy6LbVMOGsM+/RU26bsOQV/iCqSxWyra2Dp+m1RlyIdUGOj8+0HpvHO0g3c9JURHNhfbZvSPhT8IaooVWeP7L4bnp7Dk++u5HvH78+xB5RGXY7kEAV/iMr7FgOas0fa7r433+fWFxdw+mEDGf/xwVGXIzlGwR+iRJdC+vUsUvBLm7w6fw3ff+RdPlFRwnUnHaC2TWl3Cv6QVZYmtAyjpG1e1WYu+ssU9i0p5n+/NlJtmxIK7VUhK08Ws3D1FuoaGqMuRbLcmuoazrlzEp0LUrNtJroURl2S5CgFf8gqkwlqGxpZvHZL1KVIFtte18D5d09mTXUNE84aTf+9ukZdkuSw0ILfzAaY2fNmNtPMZpjZpcH2Xmb2jJnNC/7dK6waskHT1A1zVmo1Ltm5xkbn8r+9zbQlG/jtaSM4eEDPqEuSHBfmEX89cLm7DwPGAheb2TDgSuA5dy8Hnguu56yhfYvJM3X2SMt+/cwc/vHOCq48bj+OG7531OVIDIT23W93XwGsCC5vNrNZQD/gZOBTwd3uAl4A/iesOqLWpTCfst7duOPV9/jn9BVRlyNZptGdBau38NVDB3DBJ4ZEXY7EREYm/TCzMmAk8AaQDD4UAFYCO51Q3MwuAC4AGDhwYAaqDM83jx7KMzOroi5DstSnh5Vy+bEVatuUjLGwV4gys2LgReCn7v6wmW1w957Nbl/v7q2O848ePdonT54cap0iIrnGzKa4++gdt4fa1WNmhcBDwER3fzjYXGVmewe37w2sCrMGERH5qDC7egyYAMxy9xub3fQYcFZw+Szg0bBqEBGR/xTmGP8RwBnAdDObFmy7Gvg58ICZjQcWA18OsQYREdlBmF09rwAtna06OqzXFRGR1umbuyIiMaPgFxGJGQW/iEjMKPhFRGIm9C9wtQczW02qA6gj6wOsibqILKL340N6Lz5K78dH7cn7McjdS3bc2CGCPxeY2eSdfYMurvR+fEjvxUfp/fioMN4PDfWIiMSMgl9EJGYU/JlzW9QFZBm9Hx/Se/FRej8+qt3fD43xi4jEjI74RURiRsEvIhIzCv52pkXmd87M8s3sLTN7Irg+2MzeMLP5Zna/mXWKusZMMbOeZvagmc02s1lmdnhc9w8zuyz4PXnXzP5qZl3itG+Y2e1mtsrM3m22baf7gqXcHLwv75jZIbv7ugr+9qdF5nfuUmBWs+u/AH7j7kOB9cD4SKqKxk3AU+6+H3AwqfcldvuHmfUDLgFGu/twIB/4CvHaN+4EjtthW0v7wmeB8uDnAuCW3X1RBX87c/cV7j41uLyZ1C910yLzdwV3uwv4fCQFRsDM+gOfA/4cXDfgKODB4C6xeT/MrAfwCVKLFOHute6+gfjuHwVAkZkVAF2BFcRo33D3l4B1O2xuaV84GbjbU14HejatZthWCv4Q7c4i8znqt8AVQGNwvTewwd3rg+tLSX04xsFgYDVwRzD09Wcz60YM9w93Xwb8CnifVOBvBKYQ332jSUv7Qj9gSbP77fZ7o+APSbDI/EPAt9x9U/PbPNVDG4s+WjM7AVjl7lOiriVLFACHALe4+0hgCzsM68Rl/wjGrk8m9WG4D9CN/xz2iLWw9gUFfwi0yPxHHAGcZGaLgPtI/Rl/E6k/U5tWgOsPLIumvIxbCix19zeC6w+S+iCI4/5xDPCeu6929zrgYVL7S1z3jSYt7QvLgAHN7rfb742Cv51pkfmPcver3L2/u5eROnH3L3c/HXgeODW4W5zej5XAEjOrDDYdDcwknvvH+8BYM+sa/N40vRex3DeaaWlfeAw4M+juGQtsbDYk1Cb65m47M7OPAy8D0/lwTPtqUuP8DwADCRaZd/cdT+rkNDP7FPAddz/BzIaQ+gugF/AWMM7dayIsL2PMbASpE92dgIXAOaQOwmK3f5jZdcBppLrh3gLOIzVuHYt9w8z+CnyK1NTLVcA1wCPsZF8IPhx/T2o4bCtwjrtP3q3XVfCLiMSLhnpERGJGwS8iEjMKfhGRmFHwi4jEjIJfRCRmFPySU8ysehe3lzWfCTHN57zTzE7d9T3b9JzXmtl3gstnm9k+7fn8Iq1R8ItE72xSUxaIZISCXzoEMxsTzEHexcy6BXO4D2/l/sVm9pyZTTWz6WZ2crObC8xsYjAX/oNm1jV4zCgze9HMppjZ0+nOfGhmeWa2yMx6Nts2z8ySwV8Y/wpqf87MBu7w2FOB0cBEM5tmZkVm9kMzmxTMUX9b8MWd5u/BNDP7ZdNfLpZa6+CXwWPeMbML031fJZ4U/NIhuPskUl9Z/wlwA/AXd29tyGY78AV3PwQ4Evh1U4AClcAf3H1/YBPwjWB+pd8Bp7r7KOB24Kc7PqmZ/cjMTtqhtkZSX6v/QnCfw4DF7l4VPOdd7n4QMBG4eYfHPghMBk539xHuvg34vbuPCeaoLwJOCO5+B3Chu48AGpo9zXhSX98fA4wBzjezwa28NxJzBbu+i0jW+BEwiVSoX7KL+xrwMzP7BKmpM/rx4fS2S9z91eDyX4LnegoYDjwTfD7kk5oq+CPc/YctvN79wA9JhfNXgusAhwOnBJfvIfWhtStHmtkVpOan7wXMMLOXgYS7/zu4z718+IFwLHBQs/MQPUgt1vFeGq8lMaTgl46kN1AMFAJdSE1p3JLTgRJglLvXBbODdglu23GeEif1QTHD3Q/fzdr+DQw1sxJSC2f8ZHeexMy6AH8gtSrVEjO7lg/rbvFhwDfd/endeU2JHw31SEfyR+AHpIZMfrGL+/YgtQ5AnZkdCQxqdttAM2sK+K8BrwBzgJKm7WZWaGYHpFtYMG/634EbSc3Muja46TVSfwFA6sPo5Z08fDOQCC43hfyaYE2HU4Pn3wBsDoaRaPacAE8DXw+GqzCzimBxF5Gd0hG/dAhmdiZQ5+73mlk+8JqZHeXu/2rhIROBx81sOqkx9NnNbptDai3k20lNA3yLu9cGQyU3W2p5xAJSK4fN2KGOHwGT3f2xnbzm/aSGos5utu2bpFbb+i6plbfO2cnj7gRuNbNtpIaG/gS8S2r1pUnN7jce+JOZNQIvklqxClIzfZYBU4PzGKvJ4eUKZc9pdk6RDsLMit29Orh8JbC3u18acVnSAemIX6Tj+JyZXUXq93YxH/3LQiRtOuIXEYkZndwVEYkZBb+ISMwo+EVEYkbBLyISMwp+EZGY+X/dHEJdzAKM7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-268-e0c5fec626f9>:76: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  ax = fig.add_subplot(111)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation took 677.5 seconds\n"
     ]
    }
   ],
   "source": [
    "\"\"\"This main function defines what to run, what to simulate, what to process, what to plot...\"\"\"\n",
    "def main():\n",
    " # Put your code here:\n",
    "    \n",
    "    #test_simulate_single_core()\n",
    "    #test_simulate_multi_core(pool)\n",
    "    \n",
    "    # EXAMPLE:\n",
    "    # test plot something:  best threshold vs expected voltage\n",
    "     # parameters:\n",
    "    real_expected_voltage_list = np.linspace(10, 100, num=10+1)\n",
    "    threshold_list = np.linspace(20, 30, num=10+1)\n",
    "    number_of_simulation_samples = 10\n",
    "    expected_count = calculate_expected_count()\n",
    "    bunch_crossing_count = 200\n",
    "    \n",
    "    x_list_voltage = []\n",
    "    y_list_threshold = []\n",
    "    \n",
    "     # loop the expected voltages\n",
    "    for real_expected_voltage in real_expected_voltage_list:\n",
    "         # guess voltages +/- 5 from the real.\n",
    "        guess_expected_voltage_list = np.linspace(real_expected_voltage-5, real_expected_voltage+5, num=10+1)\n",
    "         # get the best threshold from the guess voltage list.\n",
    "        threshold, deviation = compute_the_best_threshold_and_deviation_in_terms_of_least_average_deviation_sum(\n",
    "                                    threshold_list,\n",
    "                                    guess_expected_voltage_list,\n",
    "                                    number_of_simulation_samples, \n",
    "                                    real_expected_voltage,\n",
    "                                    expected_count,\n",
    "                                    bunch_crossing_count\n",
    "                                )\n",
    "         # add data to the plot lists.\n",
    "        x_list_voltage.append(real_expected_voltage)\n",
    "        y_list_threshold.append(threshold)\n",
    "        \n",
    "     # plot the generated lists of data.\n",
    "    plot_the_data(\"Title\", \"filename\", \"x label: voltage\", \"y label: threshold\", x_list_voltage, y_list_threshold)\n",
    "\n",
    "    \n",
    "# === End of main() ===\n",
    "\n",
    "\n",
    "\n",
    "# It is just a startup script, see main() above instead for the simulation code.\n",
    "if __name__ == '__main__':\n",
    "    mp.freeze_support()\n",
    "    \n",
    "    # profiles to determine compute resources\n",
    "    profiles = {\"normal\": 2, \"fast\": 4, \"fastest\": 99}\n",
    "    max_processes = profiles[\"fastest\"]  # edit this\n",
    "    \n",
    "    # setup pool\n",
    "    nprocs = min(mp.cpu_count(), max_processes)\n",
    "    pool = None\n",
    "    with mp.Pool(processes=nprocs) as poool:\n",
    "        pool = poool\n",
    "    \n",
    "        global_counter = 0\n",
    "\n",
    "        # run code\n",
    "        t0 = time.time()\n",
    "        main()\n",
    "        t1 = time.time()\n",
    "        print(\"Simulation took\", round(t1-t0, 1), \"seconds\") # real time not the actuall process_time()\n",
    "\n",
    "        # clean up\n",
    "        pool.close()\n",
    "        pool.join()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test():\n",
    "    my_sim = simulation(threshold=30,\n",
    "                        real_expected_voltage=100,\n",
    "                        expected_count=30,\n",
    "                        bunch_crossing_count=1000)  #500000\n",
    "    my_sim.run()\n",
    "    my_sim.run()\n",
    "    \n",
    "def test2():\n",
    "    my_sim = simulation(threshold=15,\n",
    "                        real_expected_voltage=30,\n",
    "                        expected_count=60,\n",
    "                        bunch_crossing_count=1000)\n",
    "    my_sim.run()\n",
    "    my_sim.run()\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    import timeit\n",
    "    print(timeit.timeit(\"test()\", setup=\"from __main__ import test\", number=1000))\n",
    "    print(timeit.timeit(\"test2()\", setup=\"from __main__ import test2\", number=1000))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.8.2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'cpython'"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Python version test\n",
    "import platform\n",
    "print(platform.python_version())\n",
    "import sys\n",
    "sys.implementation.name"
   ]
  },
  {
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}