Optimizations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code optimizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl; import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import os; import pickle\n",
    "#import sympy\n",
    "#import scipy\n",
    "#from scipy import special\n",
    "fig_dir='/Users/juan/Google Drive/Brown/2017S/Machine Learning Applied to Scientific Data  - ENGN2912Z/NV-Tomography/';"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The line_profiler extension is already loaded. To reload it, use:\n",
      "  %reload_ext line_profiler\n"
     ]
    }
   ],
   "source": [
    "#import NV_generator\n",
    "%load_ext line_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load NV_generator.py\n",
    "# NV-spin-tomography Michael Scheer mgscheer@gmail.com\n",
    "\n",
    "# generating random NV centers\n",
    "\n",
    "import numpy as np\n",
    "import itertools\n",
    "\n",
    "# For explanation of the formulas used here, see chapter 2 of thesis.\n",
    "\n",
    "# Constants\n",
    "# lattice constant of diamond\n",
    "a0 = 3.567 * 10 ** (-10) # in meters\n",
    "# Hyperfine coupling related constants\n",
    "gam_el = -1.760859 * 10 ** 11 # Gyromagnetic ratio for a single electron in rad s-1 T-1\n",
    "gam_n = 67.262 * 10 ** 6 # Gyromagnetic ratio for a single c13 nucleus in rad s-1 T-1\n",
    "hbar = 1.05457173 * 10 ** (-34)\n",
    "pi = np.pi\n",
    "h = hbar * 2 * pi\n",
    "mu0 = 4 * pi * 10 ** (-7)\n",
    "\n",
    "# Carbon Lattice Definition\n",
    "# Basis vectors\n",
    "# i, j, k are basis vectors of Face Centered Cubic Bravais lattice\n",
    "i = a0/2*np.array([0,1,1])\n",
    "j = a0/2*np.array([1,0,1])\n",
    "k = a0/2*np.array([1,1,0])\n",
    "# a and b are the relative positions of the nitrogen and vacancy\n",
    "a = np.array([0,0,0])\n",
    "b = a0/4*np.array([1,1,1])\n",
    "# rotate all the vectors so that the symmetry axis (b) of the NV center is along the z direction\n",
    "# rotation around z axis by pi/4 puts b in yz plane\n",
    "Rz = np.array([[np.cos(pi/4),-np.sin(pi/4),0],[np.sin(pi/4),np.cos(pi/4),0],[0,0,1]])\n",
    "sqrt2angle = np.arctan(np.sqrt(2))\n",
    "# rotation around x axis by sqrt2angle moves b from yz plane to z axis\n",
    "Rx = np.array([[1,0,0],[0,np.cos(sqrt2angle),-np.sin(sqrt2angle)],[0,np.sin(sqrt2angle),np.cos(sqrt2angle)]])\n",
    "b = Rx.dot(Rz).dot(b)\n",
    "i = Rx.dot(Rz).dot(i)\n",
    "j = Rx.dot(Rz).dot(j)\n",
    "k = Rx.dot(Rz).dot(k)\n",
    "\n",
    "hyperfine_prefactor = -1 * mu0 * gam_el * gam_n * hbar/(4 * pi) # Only one hbar instead of hbar ** 2 since we convert to angular frequency units in the hamiltonian.\n",
    "\n",
    "# Takes Carbon Concentration, and gridsize (N) as inputs\n",
    "# Returns a list of 2D numpy arrays with hyperfine couplings of c13 nuclei to the NV center in angular frequency Hz\n",
    "# If sphere is true the c13 nuclei inside a sphere are returned instead of those inside a cube.\t\n",
    "def generate_NV(c13_concentration=0.011,N=25,sphere = True):\n",
    "\t# define position of NV in middle of the grid\n",
    "\tcenter = round(N/2)\n",
    "\tNVPos = center * i + center * j + center * k\n",
    "\n",
    "\tspin_list = []\n",
    "\t#Calculate hyperfine strength for all gridpoints\n",
    "\tfor n, m, l in itertools.product(range(N), repeat = 3):\n",
    "\t\tif (n,m,l) != (center, center, center): # all lattice pairs except the nitrogen and vacancy\n",
    "\t\t\tlattice_position = n * i + m * j + l * k - NVPos # relative position of a lattice point to the NV center\n",
    "\t\t\tfor pos in [lattice_position, lattice_position + b]: # the two lattice positions at this grid location\n",
    "\t\t\t\tx, y, z = pos\n",
    "\t\t\t\tr = np.sqrt(pos.dot(pos))\n",
    "\t\t\t\tcostheta = z/r\n",
    "\t\t\t\tsintheta = np.sqrt(x ** 2 + y ** 2)/r\n",
    "\t\t\t\tA = hyperfine_prefactor * (3 * costheta ** 2 - 1)/(r ** 3)\n",
    "\t\t\t\tB = hyperfine_prefactor * 3 * costheta * sintheta /(r ** 3)\n",
    "\t\t\t\tspin_list.append((r, A, B, costheta, sintheta))\n",
    "\n",
    "\tif sphere:\n",
    "\t\tspin_list.sort() # sort by r value from smallest to largest\n",
    "\t\tspin_list = spin_list[:len(spin_list)/2] # keep the closest half of the spins\n",
    "\n",
    "\t# for each lattice point, let it be a c13 with probability c13_concentration \n",
    "\tc13_inds = np.where(np.random.rand(len(spin_list))<c13_concentration)[0]\n",
    "\tc13_spins = np.array(spin_list)[c13_inds]\n",
    "\tr, A, B, costheta, sintheta = zip(*c13_spins)\n",
    "\treturn np.array(A), np.array(B), np.array(r), np.array(costheta), np.array(sintheta)\n",
    "\n",
    "# Generates C13s near a NV center. The number of C13s is roughly num_spins if N==None.\n",
    "# If enforce_bound = True, the random generation will repeat until there are no couplings greater than coupling_bound.\n",
    "def generate_spins(num_spins, enforce_bound = True, coupling_bound = 150 * 2 * pi * 1e3, c13_concentration = 0.011, N = None, abs_val = False, verbose = False):\n",
    "\tif N == None:\n",
    "\t\tN = int(np.ceil(((num_spins*1.0)/c13_concentration) ** (1.0/3)))\n",
    "\tif verbose:\n",
    "\t\tprint \"N:\", N\n",
    "\tA, B, r, costheta, sintheta = generate_NV(N = N, c13_concentration = c13_concentration)\n",
    "\twhile enforce_bound and (np.any(A ** 2 + B ** 2 > (coupling_bound) ** 2)):\n",
    "\t\tA, B, r, costheta, sintheta = generate_NV(N = N, c13_concentration = c13_concentration)\n",
    "\tif abs_val:\n",
    "\t\tB = np.abs(B)\n",
    "\treturn A, B, r, costheta, sintheta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%lprun -f generate_spins spins = generate_spins(450)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** generate_spins takes around 2 seconds to produce an (A,B) set for 450 spins **\n",
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 16.9578 s\n",
    "File: <ipython-input-2-a9d54bd61d56>\n",
    "Function: generate_spins at line 78\n",
    "\n",
    "Line \\#      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "    78                                           def generate_spins(num_spins, enforce_bound = True, coupling_bound = 150 * 2 * pi * 1e3, c13_concentration = 0.011, N = None, abs_val = False, verbose = False):\n",
    "    79         1            2      2.0      0.0  \tif N == None:\n",
    "    80         1           40     40.0      0.0  \t\tN = int(np.ceil(((num_spins*1.0)/c13_concentration) ** (1.0/3)))\n",
    "    81         1            1      1.0      0.0  \tif verbose:\n",
    "    82                                           \t\tprint \"N:\", N\n",
    "    83         1      2159591 2159591.0     12.7  \tA, B, r, costheta, sintheta = generate_NV(N = N, c13_concentration = c13_concentration)\n",
    "    84         8         1051    131.4      0.0  \twhile enforce_bound and (np.any(A ** 2 + B ** 2 > (coupling_bound) ** 2)):\n",
    "    85         7     14797100 2113871.4     87.3  \t\tA, B, r, costheta, sintheta = generate_NV(N = N, c13_concentration = c13_concentration)\n",
    "    86         1            1      1.0      0.0  \tif abs_val:\n",
    "    87                                           \t\tB = np.abs(B)\n",
    "    88         1            1      1.0      0.0  \treturn A, B, r, costheta, sintheta\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35\n"
     ]
    }
   ],
   "source": [
    "N = int(np.ceil(((450*1.0)/0.011) ** (1.0/3)))\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spins_NV = generate_NV(N=N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%lprun -f generate_NV spins_NV = generate_NV(N=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** generate_NV takes just a little longer **\n",
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 2.60855 s\n",
    "File: <ipython-input-2-a9d54bd61d56>\n",
    "Function: generate_NV at line 47\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "    47                                           def generate_NV(c13_concentration=0.011,N=25,sphere = True):\n",
    "    48                                           \t# define position of NV in middle of the grid\n",
    "    49         1            6      6.0      0.0  \tcenter = round(N/2)\n",
    "    50         1           32     32.0      0.0  \tNVPos = center * i + center * j + center * k\n",
    "    51                                           \n",
    "    52         1            1      1.0      0.0  \tspin_list = []\n",
    "    53                                           \t#Calculate hyperfine strength for all gridpoints\n",
    "    54     42876        43390      1.0      1.7  \tfor n, m, l in itertools.product(range(N), repeat = 3):\n",
    "    55     42875        56643      1.3      2.2  \t\tif (n,m,l) != (center, center, center): # all lattice pairs except the nitrogen and vacancy\n",
    "    56     42874       454674     10.6     17.4  \t\t\tlattice_position = n * i + m * j + l * k - NVPos # relative position of a lattice point to the NV center\n",
    "    57    128622       196505      1.5      7.5  \t\t\tfor pos in [lattice_position, lattice_position + b]: # the two lattice positions at this grid location\n",
    "    58     85748       307066      3.6     11.8  \t\t\t\tx, y, z = pos\n",
    "    59     85748       445021      5.2     17.1  \t\t\t\tr = np.sqrt(pos.dot(pos))\n",
    "    60     85748       101318      1.2      3.9  \t\t\t\tcostheta = z/r\n",
    "    61     85748       309047      3.6     11.8  \t\t\t\tsintheta = np.sqrt(x ** 2 + y ** 2)/r\n",
    "    62     85748       213028      2.5      8.2  \t\t\t\tA = hyperfine_prefactor * (3 * costheta ** 2 - 1)/(r ** 3)\n",
    "    63     85748       139633      1.6      5.4  \t\t\t\tB = hyperfine_prefactor * 3 * costheta * sintheta /(r ** 3)\n",
    "    64     85748       116701      1.4      4.5  \t\t\t\tspin_list.append((r, A, B, costheta, sintheta))\n",
    "    65                                           \n",
    "    66         1            1      1.0      0.0  \tif sphere:\n",
    "    67         1       164342 164342.0      6.3  \t\tspin_list.sort() # sort by r value from smallest to largest\n",
    "    68         1        12740  12740.0      0.5  \t\tspin_list = spin_list[:len(spin_list)/2] # keep the closest half of the spins\n",
    "    69                                           \n",
    "    70                                           \t# for each lattice point, let it be a c13 with probability c13_concentration \n",
    "    71         1          852    852.0      0.0  \tc13_inds = np.where(np.random.rand(len(spin_list))<c13_concentration)[0]\n",
    "    72         1        47060  47060.0      1.8  \tc13_spins = np.array(spin_list)[c13_inds]\n",
    "    73         1          300    300.0      0.0  \tr, A, B, costheta, sintheta = zip(*c13_spins)\n",
    "    74         1          195    195.0      0.0  \treturn np.array(A), np.array(B), np.array(r), np.array(costheta), np.array(sintheta)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of create_diamonds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load learning.py\n",
    "# NV-spin-tomography Michael Scheer mgscheer@gmail.com\n",
    "\n",
    "# Creating datasets and training classifiers, storing objects on disk\n",
    "\n",
    "import numpy as np\n",
    "import sklearn as sk\n",
    "from sklearn import svm, preprocessing, cross_validation, cluster, metrics\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "from sklearn.externals import joblib\n",
    "import analysis\n",
    "import NV_generator\n",
    "import pickle\n",
    "\n",
    "# store an object in a file\n",
    "def store_obj(obj, fname):\n",
    "\tjoblib.dump(obj, fname + \".pkl\")\n",
    "\n",
    "def store_obj_2(obj, fname):\n",
    "\tpickle.dump(obj, open(fname + \".pkl\",'wb'))\n",
    "\n",
    "# load an object from a file\n",
    "def load_obj(fname):\n",
    "\treturn joblib.load(fname + \".pkl\")\n",
    "\n",
    "def load_obj_2(fname):\n",
    "\treturn pickle.load(open(fname + \".pkl\",'rb'))\n",
    "\n",
    "# decides if a fitted A and B are close enough to a spin to be counted as a fit.\n",
    "# the B tends to be harder to fit so the accuracy demanded is less.\n",
    "def is_fit(fitted_A, fitted_B, spin, A_acc = .1, B_acc = .5):\n",
    "\tA_acc = analysis.mag * A_acc\n",
    "\tB_acc= analysis.mag * B_acc\n",
    "\treturn ((spin[0] >= fitted_A - A_acc) and (spin[0] <= fitted_A + A_acc) and\n",
    "\t\t(spin[1] >= fitted_B - B_acc) and (spin[1] <= fitted_B + B_acc))\n",
    "\n",
    "# creates a dataset of diamonds with N measurements already analyzed, saves it to disk\n",
    "def create_diamonds(diamond_num_list, omega_larmor, num_spins = 450):\n",
    "\tverbose, plots = False, False\n",
    "\tnum_subsets = 4\n",
    "\terror_fun = analysis.squared_error\n",
    "\terror_tol = .1/64\n",
    "\tN = 64\n",
    "\tN_vals = np.arange(0, 256, 2)\n",
    "\ttau = analysis.choose_tau_params(N)\n",
    "\tfor diamond_num in diamond_num_list:\n",
    "\t\tprint \"diamond_num: \", diamond_num\n",
    "\t\tA, B, r, costheta, sintheta = NV_generator.generate_spins(num_spins)\n",
    "\t\tdef data_func(N, tau, noise_level = .02):\n",
    "\t\t\tdata = analysis.calc_M(A, B, N, omega_larmor, tau)\n",
    "\t\t\tnoise = np.random.randn(len(data)) * noise_level\n",
    "\t\t\treturn data + noise\n",
    "\t\tdata = data_func(N, tau)\n",
    "\t\tdip_inds, windows = analysis.find_resonances(data, fit_dips_below = None)\n",
    "\t\tsuccessful_fits, good_fits, phis_list, xs_list, scaled_errors = [], [], [], [], []\n",
    "\t\tprint \"number of dips: \", len(dip_inds)\n",
    "\t\tfor dii in range(len(dip_inds)):\n",
    "\t\t\tif dii%50==0:\n",
    "\t\t\t\tprint \"dii: \", dii\n",
    "\t\t\tdip_ind = dip_inds[dii]\n",
    "\t\t\tres_tau = tau[dip_ind]\n",
    "\t\t\tN_data = data_func(N_vals, res_tau)\n",
    "\t\t\tsuccessful_fit, good_fit = True, True # did the fit succeed without throwing an error? was the error in the fit small enough?\n",
    "\t\t\ttry:\n",
    "\t\t\t\tphis, xs, scaled_error = analysis.repeated_spin_fit(N_vals, N_data, error_tol = error_tol, error_fun = error_fun, num_subsets = num_subsets, verbose = verbose, plots = plots)\n",
    "\t\t\t\tgood_fit = scaled_error <= error_tol\n",
    "\t\t\texcept analysis.FitError:\n",
    "\t\t\t\tsuccessful_fit, good_fit = False, False\n",
    "\t\t\t\tphis, xs = None, None\n",
    "\t\t\t\tscaled_error = None\n",
    "\t\t\tsuccessful_fits.append(successful_fit)\n",
    "\t\t\tgood_fits.append(good_fit)\n",
    "\t\t\tphis_list.append(phis)\n",
    "\t\t\txs_list.append(xs)\n",
    "\t\t\tscaled_errors.append(scaled_errors)\n",
    "\t\tprint \"create diamond_dict\"\n",
    "\t\tdiamond_dict = {\"A\" : A, \"B\" : B, \"r\" : r, \"costheta\" : costheta, \"sintheta\" : sintheta, \"N\" : N,\n",
    "\t\t\t\"tau\" : tau, \"data\" : data, \"N_vals\" : N_vals, \"N_data\" : N_data, \"dip_inds\" : dip_inds, \"windows\" : windows,\n",
    "\t\t\t\"successful_fits\" : successful_fits, \"good_fits\" : good_fits, \"error_tol\" : error_tol, \"num_subsets\" : num_subsets,\n",
    "\t\t\t\"phis_list\" : phis_list, \"xs_list\" : xs_list, \"scaled_errors\" : scaled_errors}\n",
    "\t\tprint \"store diamond_dict\"\n",
    "\t\tstore_obj_2(diamond_dict, \"diamonds/gdiamond_\" + str(diamond_num))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "diamond_num:  1\n",
      "number of dips:  1114\n",
      "dii:  0\n",
      "dii:  50\n",
      "dii:  100\n",
      "dii:  150\n",
      "dii:  200\n",
      "dii:  250\n",
      "dii:  300\n",
      "dii:  350\n",
      "dii:  400\n",
      "dii:  450\n",
      "dii:  500\n",
      "dii:  550\n",
      "dii:  600\n",
      "dii:  650\n",
      "dii:  700\n",
      "dii:  750\n",
      "dii:  800\n",
      "dii:  850\n",
      "dii:  900\n",
      "dii:  950\n",
      "dii:  1000\n",
      "dii:  1050\n",
      "dii:  1100\n",
      "create diamond_dict\n",
      "store diamond_dict\n"
     ]
    }
   ],
   "source": [
    "B_field = 0.0403555 # Teslas\n",
    "gam_c = 67.262 * 10 ** 6 # Gyromagnetic ratio for a single c13 nucleus in rad s-1 T-1\n",
    "ms = 1\n",
    "omega_larmor = -1 * ms * gam_c * B_field\n",
    "%lprun -f create_diamonds carbon = create_diamonds([1], omega_larmor, num_spins = 450)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 742.781 s\n",
    "File: <ipython-input-14-bdaeb4fb273b>\n",
    "Function: create_diamonds at line 38\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "    38                                           def create_diamonds(diamond_num_list, omega_larmor, num_spins = 450):\n",
    "    39         1            3      3.0      0.0  \tverbose, plots = False, False\n",
    "    40         1            2      2.0      0.0  \tnum_subsets = 4\n",
    "    41         1            3      3.0      0.0  \terror_fun = analysis.squared_error\n",
    "    42         1            2      2.0      0.0  \terror_tol = .1/64\n",
    "    43         1            1      1.0      0.0  \tN = 64\n",
    "    44         1           20     20.0      0.0  \tN_vals = np.arange(0, 256, 2)\n",
    "    45         1          321    321.0      0.0  \ttau = analysis.choose_tau_params(N)\n",
    "    46         2            9      4.5      0.0  \tfor diamond_num in diamond_num_list:\n",
    "    47         1          259    259.0      0.0  \t\tprint \"diamond_num: \", diamond_num\n",
    "    48         1     12260015 12260015.0      1.7  \t\tA, B, r, costheta, sintheta = NV_generator.generate_spins(num_spins)\n",
    "    49         1            6      6.0      0.0  \t\tdef data_func(N, tau, noise_level = .02):\n",
    "    50                                           \t\t\tdata = analysis.calc_M(A, B, N, omega_larmor, tau)\n",
    "    51                                           \t\t\tnoise = np.random.randn(len(data)) * noise_level\n",
    "    52                                           \t\t\treturn data + noise\n",
    "    53         1       347356 347356.0      0.0  \t\tdata = data_func(N, tau)\n",
    "    54         1         2987   2987.0      0.0  \t\tdip_inds, windows = analysis.find_resonances(data, fit_dips_below = None)\n",
    "    55         1            3      3.0      0.0  \t\tsuccessful_fits, good_fits, phis_list, xs_list, scaled_errors = [], [], [], [], []\n",
    "    56         1          133    133.0      0.0  \t\tprint \"number of dips: \", len(dip_inds)\n",
    "    57      1218         2923      2.4      0.0  \t\tfor dii in range(len(dip_inds)):\n",
    "    58      1217         2656      2.2      0.0  \t\t\tif dii%50==0:\n",
    "    59        25         4049    162.0      0.0  \t\t\t\tprint \"dii: \", dii\n",
    "    60      1217         4196      3.4      0.0  \t\t\tdip_ind = dip_inds[dii]\n",
    "    61      1217         3451      2.8      0.0  \t\t\tres_tau = tau[dip_ind]\n",
    "    62      1217     28244165  23208.0      3.8  \t\t\tN_data = data_func(N_vals, res_tau)\n",
    "    63      1217         5404      4.4      0.0  \t\t\tsuccessful_fit, good_fit = True, True # did the fit succeed without throwing an error? was the error in the fit small enough?\n",
    "    64      1217         2340      1.9      0.0  \t\t\ttry:\n",
    "    65      1217    701259518 576219.8     94.4  \t\t\t\tphis, xs, scaled_error = analysis.repeated_spin_fit(N_vals, N_data, error_tol = error_tol, error_fun = error_fun, num_subsets = num_subsets, verbose = verbose, plots = plots)\n",
    "    66      1185         4564      3.9      0.0  \t\t\t\tgood_fit = scaled_error <= error_tol\n",
    "    67        32          146      4.6      0.0  \t\t\texcept analysis.FitError:\n",
    "    68        32           70      2.2      0.0  \t\t\t\tsuccessful_fit, good_fit = False, False\n",
    "    69        32           58      1.8      0.0  \t\t\t\tphis, xs = None, None\n",
    "    70        32           56      1.8      0.0  \t\t\t\tscaled_error = None\n",
    "    71      1217         4117      3.4      0.0  \t\t\tsuccessful_fits.append(successful_fit)\n",
    "    72      1217         2548      2.1      0.0  \t\t\tgood_fits.append(good_fit)\n",
    "    73      1217         2534      2.1      0.0  \t\t\tphis_list.append(phis)\n",
    "    74      1217         2532      2.1      0.0  \t\t\txs_list.append(xs)\n",
    "    75      1217         2462      2.0      0.0  \t\t\tscaled_errors.append(scaled_errors)\n",
    "    76         1          104    104.0      0.0  \t\tprint \"create diamond_dict\"\n",
    "    77         1            4      4.0      0.0  \t\tdiamond_dict = {\"A\" : A, \"B\" : B, \"r\" : r, \"costheta\" : costheta, \"sintheta\" : sintheta, \"N\" : N,\n",
    "    78         1            3      3.0      0.0  \t\t\t\"tau\" : tau, \"data\" : data, \"N_vals\" : N_vals, \"N_data\" : N_data, \"dip_inds\" : dip_inds, \"windows\" : windows,\n",
    "    79         1            2      2.0      0.0  \t\t\t\"successful_fits\" : successful_fits, \"good_fits\" : good_fits, \"error_tol\" : error_tol, \"num_subsets\" : num_subsets,\n",
    "    80         1            2      2.0      0.0  \t\t\t\"phis_list\" : phis_list, \"xs_list\" : xs_list, \"scaled_errors\" : scaled_errors}\n",
    "    81         1           24     24.0      0.0  \t\tprint \"store diamond_dict\"\n",
    "    82         1       621593 621593.0      0.1  \t\tstore_obj_2(diamond_dict, \"diamonds/gdiamond_\" + str(diamond_num))\n",
    " ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the **create_diamonds** function most of the time (about 95%) is gobbled up by the **analysis.repeated_spin_fit** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "diamond_num:  0\n",
      "number of dips:  1104\n",
      "dii:  0\n",
      "dii:  50\n",
      "create diamond_dict\n",
      "store diamond_dict\n"
     ]
    }
   ],
   "source": [
    "#def create_diamonds(diamond_num_list, omega_larmor, num_spins = 450):\n",
    "diamond_num_list = [0]\n",
    "num_spins = 450\n",
    "verbose, plots = False, False\n",
    "num_subsets = 4\n",
    "error_fun = analysis.squared_error\n",
    "error_tol = .1/64\n",
    "N = 64\n",
    "N_vals = np.arange(0, 256, 2)\n",
    "tau = analysis.choose_tau_params(N)\n",
    "for diamond_num in diamond_num_list:\n",
    "    print \"diamond_num: \", diamond_num\n",
    "    A, B, r, costheta, sintheta = NV_generator.generate_spins(num_spins)\n",
    "    def data_func(N, tau, noise_level = .02):\n",
    "        data = analysis.calc_M(A, B, N, omega_larmor, tau)\n",
    "        noise = np.random.randn(len(data)) * noise_level\n",
    "        return data + noise\n",
    "    data = data_func(N, tau)\n",
    "    dip_inds, windows = analysis.find_resonances(data, fit_dips_below = None)\n",
    "    successful_fits, good_fits, phis_list, xs_list, scaled_errors = [], [], [], [], []\n",
    "    print \"number of dips: \", len(dip_inds)\n",
    "    for dii in range(len(dip_inds)):\n",
    "        if dii%50==0:\n",
    "            print \"dii: \", dii\n",
    "        if dii == 50:\n",
    "            break\n",
    "        dip_ind = dip_inds[dii]\n",
    "        res_tau = tau[dip_ind]\n",
    "        N_data = data_func(N_vals, res_tau)\n",
    "        successful_fit, good_fit = True, True # did the fit succeed without throwing an error? was the error in the fit small enough?\n",
    "        try:\n",
    "            phis, xs, scaled_error = analysis.repeated_spin_fit(N_vals, N_data, error_tol = error_tol, error_fun = error_fun, num_subsets = num_subsets, verbose = verbose, plots = plots)\n",
    "            good_fit = scaled_error <= error_tol\n",
    "        except analysis.FitError:\n",
    "            successful_fit, good_fit = False, False\n",
    "            phis, xs = None, None\n",
    "            scaled_error = None\n",
    "        successful_fits.append(successful_fit)\n",
    "        good_fits.append(good_fit)\n",
    "        phis_list.append(phis)\n",
    "        xs_list.append(xs)\n",
    "        scaled_errors.append(scaled_errors)\n",
    "    print \"create diamond_dict\"\n",
    "    diamond_dict = {\"A\" : A, \"B\" : B, \"r\" : r, \"costheta\" : costheta, \"sintheta\" : sintheta, \"N\" : N,\n",
    "        \"tau\" : tau, \"data\" : data, \"N_vals\" : N_vals, \"N_data\" : N_data, \"dip_inds\" : dip_inds, \"windows\" : windows,\n",
    "        \"successful_fits\" : successful_fits, \"good_fits\" : good_fits, \"error_tol\" : error_tol, \"num_subsets\" : num_subsets,\n",
    "        \"phis_list\" : phis_list, \"xs_list\" : xs_list, \"scaled_errors\" : scaled_errors}\n",
    "    print \"store diamond_dict\"\n",
    "    #store_obj_2(diamond_dict, \"diamonds/gdiamond_\" + str(diamond_num))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%lprun -f analysis.repeated_spin_fit phis, xs, scaled_error = analysis.repeated_spin_fit(N_vals, N_data, error_tol = error_tol, error_fun = error_fun, num_subsets = num_subsets, verbose = verbose, plots = plots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 0.74174 s\n",
    "File: analysis.py\n",
    "Function: repeated_spin_fit at line 215\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "   215                                           def repeated_spin_fit(N_vals, N_data, error_tol, error_fun = squared_error, num_subsets = 4,\n",
    "   216                                           \tfit_fun = spin_fit_fun, extra_params = [], verbose = True, plots = True):\n",
    "   217         1            3      3.0      0.0  \tspin_fits, scaled_errors = [], []\n",
    "   218         9          154     17.1      0.0  \tfor subset in [np.arange(int(len(N_vals) * r)) for r in np.linspace(0, 1, 1 + num_subsets)[1:]]:\n",
    "   219         4            6      1.5      0.0  \t\ttry:\n",
    "   220         4           22      5.5      0.0  \t\t\tspin_ans = spin_fit(N_vals[subset], N_data[subset], error_fun = error_fun, fit_fun = fit_fun,\n",
    "   221         4       689130 172282.5     92.9  \t\t\t\textra_params = extra_params, verbose = verbose, plots = plots)\n",
    "   222         4           16      4.0      0.0  \t\t\tspin_fits.append(spin_ans)\n",
    "   223         4          417    104.2      0.1  \t\t\tscaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *spin_ans[0]), N_data)/len(N_vals))\n",
    "   224                                           \t\texcept FitError:\n",
    "   225                                           \t\t\tpass\n",
    "   226         1            2      2.0      0.0  \tif not spin_fits:\n",
    "   227                                           \t\traise FitError(\"Failed to fit spins with N data.\")\n",
    "   228         1           19     19.0      0.0  \tmin_ind = np.argmin(scaled_errors)\n",
    "   229         1            2      2.0      0.0  \tpopt_spin, num_spins = spin_fits[min_ind]\n",
    "   230         1            1      1.0      0.0  \tif verbose:\n",
    "   231                                           \t\tprint \"best before final pass: \", min_ind\n",
    "   232                                           \t\tphis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "   233                                           \t\tphis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "   234                                           \t\tprint \"answer before final pass: \", xs, phis\n",
    "   235         1            1      1.0      0.0  \ttry:\n",
    "   236         1        51616  51616.0      7.0  \t\tpopts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, popt_spin, error_fun = error_fun, fit_fun = fit_fun)\n",
    "   237         1            3      3.0      0.0  \t\tif valid:\n",
    "   238         1            6      6.0      0.0  \t\t\tspin_fits.append((popts, num_spins))\n",
    "   239         1          157    157.0      0.0  \t\t\tscaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *popts), N_data)/len(N_vals))\n",
    "   240         1          143    143.0      0.0  \t\t\tmin_ind = np.argmin(scaled_errors)\n",
    "   241         1            4      4.0      0.0  \t\t\tpopt_spin, num_spins = spin_fits[min_ind]\n",
    "   242         1            3      3.0      0.0  \t\tif plots:\n",
    "   243                                           \t\t\tplt.figure(figsize = (10,10))\n",
    "   244                                           \t\t\tplt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "   245                                           \t\t\tplt.plot(N_vals, fit_fun(N_vals, *popts), label = 'fit')\n",
    "   246                                           \t\t\tplt.title('Fitted N data')\n",
    "   247                                           \t\t\tplt.xlabel('N')\n",
    "   248                                           \t\t\tplt.ylabel('M')\n",
    "   249                                           \t\t\tplt.legend()\n",
    "   250                                           \t\t\tplt.show()\n",
    "   251                                           \texcept FitError:\n",
    "   252                                           \t\tpass\n",
    "   253         1            6      6.0      0.0  \tif verbose:\n",
    "   254                                           \t\tprint \"best overall: \", min_ind\n",
    "   255                                           \t\tprint scaled_errors\n",
    "   256         1            4      4.0      0.0  \tscaled_error = scaled_errors[min_ind]\n",
    "   257         1           22     22.0      0.0  \tphis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "   258                                           \t#phis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "   259         1            3      3.0      0.0  \treturn phis, xs, scaled_error\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And within **analysis.repeated_spin_fit** most of time (93%) is taken by the **spin_fit** function calls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#def repeated_spin_fit(N_vals, N_data, error_tol, error_fun = squared_error, num_subsets = 4,fit_fun = spin_fit_fun, extra_params = [], verbose = True, plots = True):\n",
    "spin_fits, scaled_errors = [], []\n",
    "for subset in [np.arange(int(len(N_vals) * r)) for r in np.linspace(0, 1, 1 + num_subsets)[1:]]:\n",
    "    one_subset = subset\n",
    "    break\n",
    "#     try:\n",
    "#         spin_ans = spin_fit(N_vals[subset], N_data[subset], error_fun = error_fun, fit_fun = fit_fun,\n",
    "#             extra_params = extra_params, verbose = verbose, plots = plots)\n",
    "#         spin_fits.append(spin_ans)\n",
    "#         scaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *spin_ans[0]), N_data)/len(N_vals))\n",
    "#     except FitError:\n",
    "#         pass\n",
    "# if not spin_fits:\n",
    "#     raise FitError(\"Failed to fit spins with N data.\")\n",
    "# min_ind = np.argmin(scaled_errors)\n",
    "# popt_spin, num_spins = spin_fits[min_ind]\n",
    "# if verbose:\n",
    "#     print \"best before final pass: \", min_ind\n",
    "#     phis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "#     phis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "#     print \"answer before final pass: \", xs, phis\n",
    "# try:\n",
    "#     popts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, popt_spin, error_fun = error_fun, fit_fun = fit_fun)\n",
    "#     if valid:\n",
    "#         spin_fits.append((popts, num_spins))\n",
    "#         scaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *popts), N_data)/len(N_vals))\n",
    "#         min_ind = np.argmin(scaled_errors)\n",
    "#         popt_spin, num_spins = spin_fits[min_ind]\n",
    "#     if plots:\n",
    "#         plt.figure(figsize = (10,10))\n",
    "#         plt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "#         plt.plot(N_vals, fit_fun(N_vals, *popts), label = 'fit')\n",
    "#         plt.title('Fitted N data')\n",
    "#         plt.xlabel('N')\n",
    "#         plt.ylabel('M')\n",
    "#         plt.legend()\n",
    "#         plt.show()\n",
    "# except FitError:\n",
    "#     pass\n",
    "# if verbose:\n",
    "#     print \"best overall: \", min_ind\n",
    "#     print scaled_errors\n",
    "# scaled_error = scaled_errors[min_ind]\n",
    "# phis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "# #phis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "# #return phis, xs, scaled_error  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subset = one_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%lprun -f analysis.spin_fit spin_ans = analysis.spin_fit(N_vals[subset], N_data[subset], error_fun = error_fun, fit_fun = analysis.spin_fit_fun,extra_params = [], verbose = verbose, plots = plots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 0.121627 s\n",
    "File: analysis.py\n",
    "Function: spin_fit at line 166\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "   166                                           def spin_fit(N_vals, N_data, error_fun = squared_error, verbose = True, plots = True, fit_fun = spin_fit_fun, extra_params = []):\n",
    "   167         1          295    295.0      0.2  \tcoefs, freqs = sinusoid_coefs_fft(N_vals, N_data)\n",
    "   168         1         1152   1152.0      0.9  \tlikely_freqs = find_peaks_from_fourier(coefs, freqs)\n",
    "   169         1            2      2.0      0.0  \tparams, errors, spins_found = [], [], []\n",
    "   170                                           \t# the higest number of spins that can be fit with l frequencies\n",
    "   171         1            1      1.0      0.0  \tdef highest_spin(l):\n",
    "   172                                           \t\tif l <= 3:\n",
    "   173                                           \t\t\treturn l\n",
    "   174                                           \t\telif l == 4:\n",
    "   175                                           \t\t\treturn 3\n",
    "   176                                           \t\telse:\n",
    "   177                                           \t\t\treturn 4 # make this 5 if you figure out how to do 5\n",
    "   178         5           14      2.8      0.0  \tfor num_spins in range(1,1 + highest_spin(len(likely_freqs))):\n",
    "   179         4          148     37.0      0.1  \t\tapprox_phis = gen_freqs_from_fourier(likely_freqs, num_spins)\n",
    "   180         4           13      3.2      0.0  \t\tapprox_xs = num_spins * [2.0]\n",
    "   181         4           12      3.0      0.0  \t\tp0s = approx_phis + approx_xs + extra_params[:num_spins]\n",
    "   182         4            7      1.8      0.0  \t\ttry:\n",
    "   183         4       119918  29979.5     98.6  \t\t\tpopts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, p0s, error_fun = error_fun, fit_fun = fit_fun)\n",
    "   184         4           12      3.0      0.0  \t\t\tif valid:\n",
    "   185         3            8      2.7      0.0  \t\t\t\tparams.append(popts)\n",
    "   186         3            6      2.0      0.0  \t\t\t\terrors.append(error)\n",
    "   187         3            6      2.0      0.0  \t\t\t\tspins_found.append(num_spins)\n",
    "   188                                           \t\texcept FitError:\n",
    "   189                                           \t\t\tpass\n",
    "   190         1            1      1.0      0.0  \tif not params:\n",
    "   191                                           \t\traise FitError(\"Failed to fit spins with N data.\")\n",
    "   192                                           \telse:\n",
    "   193         1           28     28.0      0.0  \t\tind = np.argmin(errors)\n",
    "   194         1            1      1.0      0.0  \t\tif plots:\n",
    "   195                                           \t\t\tplt.figure(figsize=(10,10))\n",
    "   196                                           \t\t\tplt.plot(freqs, coefs.real, label = 'real part')\n",
    "   197                                           \t\t\tplt.plot(freqs, coefs.imag, label = 'imag part')\n",
    "   198                                           \t\t\tplt.title('Fourier spectrum')\n",
    "   199                                           \t\t\tplt.xlabel('phi')\n",
    "   200                                           \t\t\tplt.ylabel('Fourier coeficient')\n",
    "   201                                           \t\t\tplt.legend()\n",
    "   202                                           \t\t\tplt.show()\n",
    "   203                                           \t\t\tplt.figure(figsize=(10,10))\n",
    "   204                                           \t\t\tplt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "   205                                           \t\t\tplt.plot(N_vals, fit_fun(N_vals, *params[ind]), label = 'fit')\n",
    "   206                                           \t\t\tplt.title('Fitted N data')\n",
    "   207                                           \t\t\tplt.xlabel('N')\n",
    "   208                                           \t\t\tplt.ylabel('M')\n",
    "   209                                           \t\t\tplt.legend()\n",
    "   210                                           \t\t\tplt.show()\n",
    "   211         1            3      3.0      0.0  \t\treturn params[ind], spins_found[ind]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And within ```analysis.spin_fit()``` 98% of the time is taken by the ```analysis.spin_fit_guess()``` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#def spin_fit(N_vals, N_data, error_fun = squared_error, verbose = True, plots = True, fit_fun = spin_fit_fun, extra_params = []):\n",
    "coefs, freqs = analysis.sinusoid_coefs_fft(N_vals, N_data)\n",
    "likely_freqs = analysis.find_peaks_from_fourier(coefs, freqs)\n",
    "params, errors, spins_found = [], [], []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# the higest number of spins that can be fit with l frequencies\n",
    "extra_params=[]\n",
    "def highest_spin(l):\n",
    "    if l <= 3:\n",
    "        return l\n",
    "    elif l == 4:\n",
    "        return 3\n",
    "    else:\n",
    "        return 4 # make this 5 if you figure out how to do 5\n",
    "for num_spins in range(1,1 + highest_spin(len(likely_freqs))):\n",
    "    approx_phis = analysis.gen_freqs_from_fourier(likely_freqs, num_spins)\n",
    "    approx_xs = num_spins * [2.0]\n",
    "    p0s = approx_phis + approx_xs + extra_params[:num_spins]\n",
    "    break\n",
    "    try:\n",
    "        popts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, p0s, error_fun = error_fun, fit_fun = fit_fun)\n",
    "        \n",
    "        if valid:\n",
    "            params.append(popts)\n",
    "            errors.append(error)\n",
    "            spins_found.append(num_spins)\n",
    "    except FitError:\n",
    "        pass\n",
    "# if not params:\n",
    "#     raise FitError(\"Failed to fit spins with N data.\")\n",
    "# else:\n",
    "#     ind = np.argmin(errors)\n",
    "#     if plots:\n",
    "#         plt.figure(figsize=(10,10))\n",
    "#         plt.plot(freqs, coefs.real, label = 'real part')\n",
    "#         plt.plot(freqs, coefs.imag, label = 'imag part')\n",
    "#         plt.title('Fourier spectrum')\n",
    "#         plt.xlabel('phi')\n",
    "#         plt.ylabel('Fourier coeficient')\n",
    "#         plt.legend()\n",
    "#         plt.show()\n",
    "#         plt.figure(figsize=(10,10))\n",
    "#         plt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "#         plt.plot(N_vals, fit_fun(N_vals, *params[ind]), label = 'fit')\n",
    "#         plt.title('Fitted N data')\n",
    "#         plt.xlabel('N')\n",
    "#         plt.ylabel('M')\n",
    "#         plt.legend()\n",
    "#         plt.show()\n",
    "#     return params[ind], spins_found[ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#def spin_fit_guess(N_vals, N_data, num_spins, param_guess, error_fun = squared_error, fit_fun = spin_fit_fun):\n",
    "def sspin_fit_guess():\n",
    "    try:\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.simplefilter(\"ignore\")\n",
    "            popts, pcovs = opt.curve_fit(analysis.spin_fit_fun, N_vals, N_data, p0 = p0s)\n",
    "        xs = popts[num_spins:2*num_spins]\n",
    "    #     error = error_fun(fit_fun(N_vals, *popts), N_data)\n",
    "    #     #return popts, error, np.all((xs > 0) & (xs <= 2))\n",
    "    except (RuntimeError, TypeError):\n",
    "        raise FitError(\"Failed to fit spins with N data.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%lprun -f sspin_fit_guess sspin_fit_guess()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Timer unit: 1e-06 s\n",
    "\n",
    "Total time: 0.008301 s\n",
    "File: <ipython-input-50-d7129282266a>\n",
    "Function: sspin_fit_guess at line 2\n",
    "\n",
    "Line #      Hits         Time  Per Hit   % Time  Line Contents\n",
    "==============================================================\n",
    "     2                                           def sspin_fit_guess():\n",
    "     3         1            2      2.0      0.0      try:\n",
    "     4         1           17     17.0      0.2          with warnings.catch_warnings():\n",
    "     5         1           12     12.0      0.1              warnings.simplefilter(\"ignore\")\n",
    "     6         1         8266   8266.0     99.6              popts, pcovs = opt.curve_fit(analysis.spin_fit_fun, N_vals, N_data, p0 = p0s)\n",
    "     7         1            4      4.0      0.0          xs = popts[num_spins:2*num_spins]\n",
    "     8                                               #     error = error_fun(fit_fun(N_vals, *popts), N_data)\n",
    "     9                                               #     #return popts, error, np.all((xs > 0) & (xs <= 2))\n",
    "    10                                               except (RuntimeError, TypeError):\n",
    "    11                                                   raise FitError(\"Failed to fit spins with N data.\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, within ```spin_fit_guess``` 99% of the time is taken by the ```opt.curve_fit``` function. It takes about 8000 micro seconds to make a fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting execution times and error statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "B_field = 0.0403555 # Teslas\n",
    "gam_c = 67.262 * 10 ** 6 # Gyromagnetic ratio for a single c13 nucleus in rad s-1 T-1\n",
    "ms = 1\n",
    "omega_larmor = -1 * ms * gam_c * B_field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "global all_them_times\n",
    "all_them_times = []\n",
    "global all_them_errors\n",
    "all_them_errors = []\n",
    "global all_them_num_spins\n",
    "all_them_num_spins = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load analysis.py\n",
    "# NV-spin-tomography Michael Scheer mgscheer@gmail.com\n",
    "\n",
    "# the main code for the analysis\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "import scipy.optimize as opt\n",
    "from scipy import special, fftpack, signal\n",
    "import operator, heapq, itertools, warnings\n",
    "import sklearn as sk\n",
    "import learning\n",
    "\n",
    "# frequency scale -- turns kHz into radial Hz. Used for printing paramaters nicely\n",
    "mag = 2 * np.pi * 1e3\n",
    "\n",
    "# a special error class that is thrown in this code when a fit does not work\n",
    "class FitError(Exception):\n",
    "\tdef __init__(self, value):\n",
    "\t\tself.value = value\n",
    "\tdef __str__(self):\n",
    "\t\treturn repr(self.value)\n",
    "\n",
    "# create a nice looking plot\n",
    "def initialize_data_plot(figsize = (10,5), xlims = None, ylims = [-1.05,1.05]):\n",
    "\tfig,ax = plt.subplots(figsize=figsize)\n",
    "\tax.set_ylim(ylims)\n",
    "\tif xlims is not None:\n",
    "\t\tax.set_xlim(xlims)\n",
    "\treturn fig, ax\n",
    "\n",
    "# calculates the M function for a single spin A, B\n",
    "def calc_M_single(A,B,N,omega_larmor,tau):\n",
    "\ttau = np.array(tau)\n",
    "\tomega_tilde = np.sqrt((A + omega_larmor) ** 2 + B ** 2)\n",
    "\tmx = B/omega_tilde\n",
    "\tmz = (A+omega_larmor)/omega_tilde\n",
    "\talpha = omega_tilde * tau\n",
    "\tbeta = omega_larmor * tau\n",
    "\tcos_phi = np.cos(alpha) * np.cos(beta)- mz * np.sin(alpha) * np.sin(beta)\n",
    "\tvec_num = (mx ** 2) * (1-np.cos(alpha)) * (1-np.cos(beta))\n",
    "\tvec_denom = 1.0 + cos_phi\n",
    "\tfor i in np.where(vec_denom == 0.0)[0]:\n",
    "\t\tvec_denom[i] = .0001\n",
    "\tvec_term = vec_num/vec_denom\n",
    "\tangle_term = np.sin(N * np.arccos(cos_phi) / 2.0) ** 2\n",
    "\treturn 1 - (vec_term * angle_term)\n",
    "\n",
    "# calculates the M function corresponding to many spins given by\n",
    "# A_list, B_list\n",
    "# equivalent to np.array([calc_M_single(A, B, N, omega_larmor, tau) for A, B in zip(A_list, B_list)]).prod()\n",
    "def calc_M(A_list, B_list, N, omega_larmor, tau):\n",
    "\t'''\n",
    "\tcalc_M(A_list, B_list, N, omega_larmor, tau)\n",
    "\tgives back M\n",
    "\t'''\n",
    "\treturn reduce(lambda accum, next: accum * calc_M_single(next[0], next[1], N, omega_larmor, tau), zip(A_list, B_list), 1.0)\n",
    "\n",
    "# calculates what A and B must be given the value of cos(phi) as well as omega_tilde\n",
    "def calc_A_B(cosphi, res_tau, omega_larmor, omega_tilde):\n",
    "\talpha = omega_tilde * res_tau\n",
    "\tbeta = omega_larmor * res_tau\n",
    "\tmz = (np.cos(alpha) * np.cos(beta) - cosphi)/(np.sin(alpha) * np.sin(beta))\n",
    "\tA = mz * omega_tilde - omega_larmor\n",
    "\tB = np.sqrt((1 - mz ** 2)) * omega_tilde\n",
    "\treturn A, np.abs(B)\n",
    "\n",
    "# Firm bounds on omega_tilde arising from an upper bound on hyperfine couplings\n",
    "def omega_bounds(omega_larmor, coupling_bound = 150 * mag):\n",
    "\tupper_omega_bound = omega_larmor + coupling_bound\n",
    "\tlower_omega_bound = omega_larmor - coupling_bound\n",
    "\treturn lower_omega_bound, upper_omega_bound\n",
    "\n",
    "# checks if the A, B pair has a coupling larger than coupling_bound\n",
    "def valid_A_B(A, B, coupling_bound = 150 * mag):\n",
    "\treturn A ** 2 + B ** 2 <= coupling_bound ** 2\n",
    "\n",
    "# for a real sequence y indexed by a sequence x,\n",
    "# writes y = a0 + a1 * cos(omega_0 * x) + b1 * sin(omega_0 * x) + ...\n",
    "# with omega_0 = 2 * np.pi/len(x)\n",
    "# as a \"fourier series\" (though for discrete y not continuous) and finds (very approximately)\n",
    "# the a and b coefficients.\n",
    "# coefs.real = a0, a1, a2, ...\n",
    "# coefs.imag = 0, -b1, -b2, ... (notice the negative signs)\n",
    "# see http://www.feynmanlectures.caltech.edu/I_50.html for fourier series and\n",
    "# http://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.fftpack.fft.html for the formula used for FFT.\n",
    "# notice the 2/N difference between the two formulas for all but the 0 frequency term.\n",
    "def sinusoid_coefs_fft(x, y):\n",
    "\tft = fftpack.fft(y)\n",
    "\tfreqs = fftpack.fftfreq(len(x), d = (x[1]-x[0])/(2*np.pi))\n",
    "\t# get rid of negative frequencies\n",
    "\tft, freqs = map(np.array, zip(*filter(lambda p: p[1] >= 0, zip(ft, freqs))))\n",
    "\tcoefs = 2.0 * np.array(ft)/len(x)\n",
    "\tcoefs[0] = .5 * coefs[0]\n",
    "\treturn coefs, freqs\n",
    "\n",
    "# given the coefs and freqs from sinusoid_coefs_fft, finds the most likely frequencies\n",
    "# present in the signal. It looks for peaks in both the real and imaginary part of the coefficients.\n",
    "# it returns a list of frequencies in decreasing order of likelihood, based on the absolute value of the\n",
    "# fourier spectrum (coefs) at that point.\n",
    "def find_peaks_from_fourier(coefs, freqs):\n",
    "\textremal_inds = signal.argrelextrema(coefs.real, np.greater)[0]\n",
    "\textremal_inds = np.r_[extremal_inds, signal.argrelextrema(coefs.real, np.less)[0]]\n",
    "\textremal_inds = np.r_[extremal_inds, signal.argrelextrema(coefs.imag, np.greater)[0]]\n",
    "\textremal_inds = np.r_[extremal_inds, signal.argrelextrema(coefs.imag, np.less)[0]]\n",
    "\textremal_inds = np.unique(extremal_inds)\n",
    "\tfreq_inds = sorted(extremal_inds, key=lambda ind: -np.abs(coefs[ind]))\n",
    "\treturn freqs[freq_inds]\n",
    "\n",
    "# takes a list of frequencies that are peaks in a fourier spectrum, either of the real\n",
    "# or imaginary part as well as a number of spins to look for. It then approximates the frequencies\n",
    "# of the underlying product cosine functions that were in the signal. See appendix C of thesis.\n",
    "def gen_freqs_from_fourier(freqs, num_spins):\n",
    "\t# sort them in decreasing order\n",
    "\tfreqs = sorted(freqs, key=lambda x:-x)\n",
    "\tif num_spins == 1: # requires 1 freq\n",
    "\t\treturn [freqs[0]]\n",
    "\telif num_spins == 2: # requires 2 freqs\n",
    "\t\tf0, f1 = freqs[:2]\n",
    "\t\treturn [.5 * (f0 + f1), .5 * (f0 - f1)]\n",
    "\telif num_spins == 3: # requires 3 freqs\n",
    "\t\tf0, f1, f2 = freqs[:3]\n",
    "\t\treturn [.5 * (f1 + f2), .5 * (f0 - f2), .5 * (f0 - f1)]\n",
    "\telif num_spins == 4: # requires 5 freqs\n",
    "\t\tphi01sum, phi2, phi3 = gen_freqs_from_fourier(freqs[:4], 3)\n",
    "\t\tf3, f4 = freqs[3:5]        \n",
    "\t\tphi0 = .5 * (f3 + f4)\n",
    "\t\tphi1 = phi01sum - phi0\n",
    "\t\treturn [phi0, phi1, phi2, phi3]\n",
    "\telse:\n",
    "\t\treturn gen_freqs_from_fourier(freqs, 4) # change this if you figure out how to do 5\n",
    "\n",
    "# the form of the signal as a function of N for a fixed tau\n",
    "def spin_fit_fun_single(N, phi, x):\n",
    "\treturn 1 - x * (np.sin(N * phi/2.0) ** 2)\n",
    "\n",
    "# the form of the signal as a function of N for a fixed tau with several spins\n",
    "# there can be any number of spins passed, in the form [phi0, phi1, ..., phin, x0, x1, ..., xn]\n",
    "# where the function for a single spin is 1 - x * (np.sin(N * phi/2.0) ** 2)\n",
    "# equivalent to np.array([single_spin_fit_fun(N, phi, x) for phi, x in zip(phis, xs)]).prod()\n",
    "def spin_fit_fun(N, *args):\n",
    "\tnum_spins = len(args)/2\n",
    "\tphis, xs = args[:num_spins], args[num_spins:2*num_spins]\n",
    "\treturn reduce(lambda accum, next: accum * spin_fit_fun_single(N, next[0], next[1]), zip(phis, xs), 1.0)\n",
    "\n",
    "# for two arrays l1 and l2 returns the squared error between them\n",
    "def squared_error(l1, l2):\n",
    "\treturn sum((l1 - l2) ** 2)\n",
    "\n",
    "# Use scipy.optimize.curve_fit to fit N_data to several spins using fit_fun. Requires a guess of the spin\n",
    "# parameters (phi's and x's) as param_guess\n",
    "def spin_fit_guess(N_vals, N_data, num_spins, param_guess, error_fun = squared_error, fit_fun = spin_fit_fun):\n",
    "\ttry:\n",
    "\t\twith warnings.catch_warnings():\n",
    "\t\t\twarnings.simplefilter(\"ignore\")\n",
    "\t\t\tpopts, pcovs = opt.curve_fit(fit_fun, N_vals, N_data, p0 = param_guess)\n",
    "\t\txs = popts[num_spins:2*num_spins]\n",
    "\t\terror = error_fun(fit_fun(N_vals, *popts), N_data)\n",
    "\t\treturn popts, error, np.all((xs > 0) & (xs <= 2))\n",
    "\texcept (RuntimeError, TypeError):\n",
    "\t\traise FitError(\"Failed to fit spins with N data.\")\n",
    "\n",
    "# fit possibly several spins using a fit_fun to N_data\n",
    "# for reference, see http://stackoverflow.com/questions/13405053/scipy-leastsq-fit-to-a-sine-wave-failing\n",
    "# and appendix C of thesis.\n",
    "def spin_fit(N_vals, N_data, error_fun = squared_error, verbose = True, plots = True, fit_fun = spin_fit_fun, extra_params = []):\n",
    "\tcoefs, freqs = sinusoid_coefs_fft(N_vals, N_data)\n",
    "\tlikely_freqs = find_peaks_from_fourier(coefs, freqs)\n",
    "\tparams, errors, spins_found, times_taken = [], [], [], []\n",
    "\t# the higest number of spins that can be fit with l frequencies\n",
    "\tdef highest_spin(l):\n",
    "\t\tif l <= 3:\n",
    "\t\t\treturn l\n",
    "\t\telif l == 4:\n",
    "\t\t\treturn 3\n",
    "\t\telse:\n",
    "\t\t\treturn 4 # make this 5 if you figure out how to do 5\n",
    "\tfor num_spins in range(1,1 + highest_spin(len(likely_freqs))):\n",
    "\t\tapprox_phis = gen_freqs_from_fourier(likely_freqs, num_spins)\n",
    "\t\tapprox_xs = num_spins * [2.0]\n",
    "\t\tp0s = approx_phis + approx_xs + extra_params[:num_spins]\n",
    "\t\ttry:\n",
    "\t\t\tstart_time = time.time()\n",
    "\t\t\tpopts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, p0s, error_fun = error_fun, fit_fun = fit_fun)\n",
    "\t\t\tend_time = time.time()        \n",
    "\t\t\ttime_taken = end_time - start_time\n",
    "\t\t\tall_them_times.append(time_taken)    \n",
    "\t\t\tall_them_errors.append(error)\n",
    "\t\t\tall_them_num_spins.append(num_spins)\n",
    "\t\t\tif valid:\n",
    "\t\t\t\tparams.append(popts)\n",
    "\t\t\t\terrors.append(error)\n",
    "\t\t\t\tspins_found.append(num_spins)\n",
    "\t\t\t\ttimes_taken.append(time_taken)\n",
    "\t\texcept FitError:\n",
    "\t\t\tpass\n",
    "\tif not params:\n",
    "\t\traise FitError(\"Failed to fit spins with N data.\")\n",
    "\telse:\n",
    "\t\tind = np.argmin(errors)\n",
    "\t\tif plots:\n",
    "\t\t\tplt.figure(figsize=(10,10))\n",
    "\t\t\tplt.plot(freqs, coefs.real, label = 'real part')\n",
    "\t\t\tplt.plot(freqs, coefs.imag, label = 'imag part')\n",
    "\t\t\tplt.title('Fourier spectrum')\n",
    "\t\t\tplt.xlabel('phi')\n",
    "\t\t\tplt.ylabel('Fourier coeficient')\n",
    "\t\t\tplt.legend()\n",
    "\t\t\tplt.show()\n",
    "\t\t\tplt.figure(figsize=(10,10))\n",
    "\t\t\tplt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "\t\t\tplt.plot(N_vals, fit_fun(N_vals, *params[ind]), label = 'fit')\n",
    "\t\t\tplt.title('Fitted N data')\n",
    "\t\t\tplt.xlabel('N')\n",
    "\t\t\tplt.ylabel('M')\n",
    "\t\t\tplt.legend()\n",
    "\t\t\tplt.show()\n",
    "\t\tall_them_times.append('split')\n",
    "\t\tall_them_errors.append('split')\n",
    "\t\tall_them_num_spins.append('split')\n",
    "\t\treturn params[ind], spins_found[ind], times_taken[ind]\n",
    "\n",
    "# repeated attempts to fit N_vals and N_data using several spins. It tries subsets of the data of increasing length\n",
    "# and then uses the best result as a guess for one final fit of the whole data. See appendix C of thesis for more detail.\n",
    "def repeated_spin_fit(N_vals, N_data, error_tol, error_fun = squared_error, num_subsets = 4,\n",
    "\tfit_fun = spin_fit_fun, extra_params = [], verbose = True, plots = True):\n",
    "\tspin_fits, scaled_errors, exec_times = [], [], []\n",
    "\tfor subset in [np.arange(int(len(N_vals) * r)) for r in np.linspace(0, 1, 1 + num_subsets)[1:]]:\n",
    "\t\ttry:\n",
    "\t\t\tspin_ans = spin_fit(N_vals[subset], N_data[subset], error_fun = error_fun, fit_fun = fit_fun,extra_params = extra_params, verbose = verbose, plots = plots)\n",
    "\t\t\tfit_time = spin_ans[-1]\n",
    "\t\t\tspin_ans = spin_ans[0:2]\n",
    "\t\t\tspin_fits.append(spin_ans)\n",
    "\t\t\texec_times.append(fit_time)\n",
    "\t\t\tscaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *spin_ans[0]), N_data)/len(N_vals))\n",
    "\t\texcept FitError:\n",
    "\t\t\tpass\n",
    "\tif not spin_fits:\n",
    "\t\traise FitError(\"Failed to fit spins with N data.\")\n",
    "\tmin_ind = np.argmin(scaled_errors)\n",
    "\tpopt_spin, num_spins = spin_fits[min_ind]\n",
    "\tif verbose:\n",
    "\t\tprint \"best before final pass: \", min_ind\n",
    "\t\tphis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "\t\tphis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "\t\tprint \"answer before final pass: \", xs, phis\n",
    "\ttry:\n",
    "\t\tpopts, error, valid = spin_fit_guess(N_vals, N_data, num_spins, popt_spin, error_fun = error_fun, fit_fun = fit_fun)\n",
    "\t\tif valid:\n",
    "\t\t\tspin_fits.append((popts, num_spins))\n",
    "\t\t\tscaled_errors.append(1.0 * error_fun(fit_fun(N_vals, *popts), N_data)/len(N_vals))\n",
    "\t\t\tmin_ind = np.argmin(scaled_errors)\n",
    "\t\t\tpopt_spin, num_spins = spin_fits[min_ind]\n",
    "\t\tif plots:\n",
    "\t\t\tplt.figure(figsize = (10,10))\n",
    "\t\t\tplt.plot(N_vals, N_data, '.-k', lw=0.4, label = 'data')\n",
    "\t\t\tplt.plot(N_vals, fit_fun(N_vals, *popts), label = 'fit')\n",
    "\t\t\tplt.title('Fitted N data')\n",
    "\t\t\tplt.xlabel('N')\n",
    "\t\t\tplt.ylabel('M')\n",
    "\t\t\tplt.legend()\n",
    "\t\t\tplt.show()\n",
    "\texcept FitError:\n",
    "\t\tpass\n",
    "\tif verbose:\n",
    "\t\tprint \"best overall: \", min_ind\n",
    "\t\tprint scaled_errors\n",
    "\tscaled_error = scaled_errors[min_ind]\n",
    "\tphis, xs = np.array(popt_spin[:num_spins]) % np.pi, np.array(popt_spin[num_spins:2*num_spins])\n",
    "\t#phis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "\treturn phis, xs, scaled_error, scaled_errors, exec_times    \n",
    "\n",
    "# creates a tau vector that matches the tau from real data from Tim and Julia.\n",
    "# should be filled in with realistic tau vectors for other N in the future.\n",
    "def choose_tau_params(N, num_samples = 5100):\n",
    "\tif N != 64:\n",
    "\t\traise FitError(\"N should be 64\")\n",
    "\telse:\n",
    "\t\tmin_tau = 3.0000000000000001e-06\n",
    "\t\tmax_tau = 2.1999999999999999e-05\n",
    "\t\tlength = num_samples\n",
    "\t\ttau = np.linspace(min_tau, max_tau, length)\n",
    "\t\treturn tau\n",
    "\n",
    "# finds local minima in the vector data as well as the intervals in which they are minima.\n",
    "# only finds dips below fit_dips_below and less wide than max_width if the parameters are not None.\n",
    "def find_resonances(data, fit_dips_below = None, max_width = None):\n",
    "\tmin_inds = (signal.argrelextrema(data, np.less)[0]).astype(int)\n",
    "\tmax_inds = (signal.argrelextrema(data, np.greater)[0]).astype(int)\n",
    "\twidths = max_inds[1:] - max_inds[:-1]\n",
    "\tdip_inds = min_inds[1:] if min_inds[0] < max_inds[0] else min_inds[:-1]\n",
    "\twindows = [(max_inds[i], max_inds[i+1]) for i in range(len(dip_inds)-1)] # should this be 3?\n",
    "\t# a dip must go below fit_dips_below, unless fit_dips_below = None and it's width should be at most max_width unless max_width = None\n",
    "\tfilter_fun = lambda p: ((data[p[0]] < fit_dips_below) if (fit_dips_below != None) else True) and (p[1] <= max_width if max_width != None else True)\n",
    "\tdip_inds, widths, windows = map(np.array, zip(*filter(filter_fun, zip(dip_inds, widths, windows))))\n",
    "\t# sort by width\n",
    "\t# sort_inds = np.argsort(-1 * widths)\n",
    "\t# return dip_inds[sort_inds], windows[sort_inds]\n",
    "\treturn dip_inds, windows\n",
    "\n",
    "# Calculates all possible omega_tilde values for a given phi and x at a given res_tau.\n",
    "# See chapter 4 of thesis.\n",
    "def calc_omega_tilde(phi, x, res_tau, omega_larmor):\n",
    "\tbeta = omega_larmor * res_tau\n",
    "\tcosbeta, cosphi_0 = np.cos(beta), np.cos(phi)\n",
    "\tlower_omega_bound, upper_omega_bound = omega_bounds(omega_larmor)\n",
    "\tomega_cosphi = []\n",
    "\tfor cosphi in [cosphi_0, -cosphi_0]:\n",
    "\t\tprod = x * (1+cosbeta) * (1+cosphi)\n",
    "\t\tcoeff = [1, prod - 2 * cosbeta * cosphi, prod + cosbeta ** 2 + cosphi ** 2 - 1]\n",
    "\t\tif coeff[1] ** 2 - 4 * coeff[0] * coeff[2] >= 0: # discriminant is nonnegative\n",
    "\t\t\tcosalphas = np.roots(coeff)\n",
    "\t\t\tfor cosalpha in filter(lambda x: np.abs(x) <= 1, cosalphas):\n",
    "\t\t\t\talpha_0_pi = np.arccos(cosalpha)\n",
    "\t\t\t\tfor alpha_base in [alpha_0_pi, 2 * np.pi - alpha_0_pi]:\n",
    "\t\t\t\t\tlower_int = int(np.ceil((lower_omega_bound * res_tau - alpha_base)/(2 * np.pi)))\n",
    "\t\t\t\t\tupper_int = int(np.floor((upper_omega_bound * res_tau - alpha_base)/(2 * np.pi)))\n",
    "\t\t\t\t\tfor n in range(lower_int, upper_int + 1):\n",
    "\t\t\t\t\t\tomega_cosphi.append(((alpha_base + 2 * np.pi * n)/res_tau, cosphi))\n",
    "\treturn omega_cosphi\n",
    "\n",
    "# for a given dip_ind, measures N data and fits it using repeated_spin_fit. If the fit error is less than error_tol,\n",
    "# it finds all possible A and B values for each fit spin. Each possible A and B value that is within the bounds (is_valid)\n",
    "# is added to spin_dict.\n",
    "def analyze_dip(dip_ind, tau, data_func, omega_larmor, spin_dict = {}, error_tol = .1/64,\n",
    "\t\t\t\tN_vals = np.arange(0,256,2), error_fun = squared_error, fit_fun = spin_fit_fun, extra_params = [],\n",
    "\t\t\t\tnum_subsets = 4, verbose = False, plots = False):\n",
    "\tres_tau = tau[dip_ind]\n",
    "\tN_data = data_func(N_vals, res_tau)\n",
    "\ttry:\n",
    "\t\tphis, xs, scaled_error = repeated_spin_fit(N_vals, N_data, error_fun = error_fun, fit_fun = fit_fun, extra_params = extra_params,\n",
    "\t\t\tnum_subsets = num_subsets, verbose = verbose, plots = plots)\n",
    "\texcept FitError:\n",
    "\t\tif verbose:\n",
    "\t\t\tprint \"Failed to fit spins with N data.\"\n",
    "\t\treturn spin_dict\n",
    "\tif verbose:\n",
    "\t\tprint \"res_tau: \", res_tau\n",
    "\t\tprint \"scaled_error <= error_tol: \", scaled_error <= error_tol\n",
    "\t\tprint \"scaled_error: \", scaled_error\n",
    "\t\tphis, xs = map(np.array, zip(*sorted(zip(phis, xs), key = lambda x:-x[1])))# sort phis in decreasing order by xs\n",
    "\t\tprint \"xs, phis, scaled_error: \", xs, phis, scaled_error\n",
    "\tif scaled_error > error_tol:\n",
    "\t\treturn spin_dict\n",
    "\tfor phi, x in zip(phis, xs):\n",
    "\t\tfor omega_tilde, cosphi in calc_omega_tilde(phi, x, res_tau, omega_larmor): # all the possible omegas with their associated cosphis\n",
    "\t\t\tA, B = calc_A_B(cosphi, res_tau, omega_larmor, omega_tilde)\n",
    "\t\t\tif valid_A_B(A, B):\n",
    "\t\t\t\tif (dip_ind, phi, x) in spin_dict:\n",
    "\t\t\t\t\tspin_dict[(dip_ind, phi, x)] += [(A, B, cosphi)]\n",
    "\t\t\t\telse:\n",
    "\t\t\t\t\tspin_dict[(dip_ind, phi, x)] = [(A, B, cosphi)]\n",
    "\treturn spin_dict\n",
    "\n",
    "# takes in a dict of the type returned by analyze_dip and for each entry which represents a particular spin found at a particular dip\n",
    "# it chooses the A, B value that minimizes the error_fun from the tau data. It returns these A and B guesses along with the calculated errors.\n",
    "def choose_spin_guesses(spin_dict, N, omega_larmor, tau, data, classifier, x_min = 1, error_fun = squared_error):\n",
    "\tguess_As, guess_Bs, dataerrs = [], [], []\n",
    "\tall_guess_As, all_guess_Bs, select_As, select_Bs = [], [], [], []\n",
    "\tfor k in spin_dict.keys():\n",
    "\t\terr = []\n",
    "\t\tfor A, B, _ in spin_dict[k]:\n",
    "\t\t\tall_guess_As.append(A)\n",
    "\t\t\tall_guess_Bs.append(B)\n",
    "\t\t\terr.append(error_fun(calc_M_single(A, B, N, omega_larmor, tau), data))\n",
    "\t\tmin_err_ind = np.argmin(err)\n",
    "\t\tbest_A, best_B, best_cosphi = spin_dict[k][min_err_ind]\n",
    "\t\tselect_As.append(best_A)\n",
    "\t\tselect_Bs.append(best_B)\n",
    "\t\tbest_err = err[min_err_ind]\n",
    "\t\t_, _, x = k\n",
    "\t\tfeatures = [best_cosphi, x, best_A, best_B, best_err]\n",
    "\t\tif classifier(features) and x >= x_min:\n",
    "\t\t\tguess_As.append(best_A)\n",
    "\t\t\tguess_Bs.append(best_B)\n",
    "\t\t\tdataerrs.append(best_err)\n",
    "\treturn (np.array(guess_As), np.array(guess_Bs), np.array(dataerrs),\n",
    "\t\tnp.array(all_guess_As), np.array(all_guess_Bs), np.array(select_As), np.array(select_Bs))\n",
    "\n",
    "# given guesses for As and Bs, clusters these guesses using the DBSCAN algorithm.\n",
    "# for info about DBSCAN see http://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html\n",
    "# for info about clustering with sklearn see http://scikit-learn.org/stable/modules/clustering.html#clustering\n",
    "def cluster_spin_guesses(guess_As, guess_Bs, dataerrs, eps = .075, min_samples = 1):\n",
    "\tX = sk.preprocessing.StandardScaler().fit_transform(zip(guess_As, guess_Bs)) # possibly need to be lists?\n",
    "\tdb = sk.cluster.DBSCAN(eps = eps, min_samples = min_samples).fit(X)\n",
    "\tlabels = db.labels_\n",
    "\tcluster_As, cluster_Bs, cluster_dataerrs = [], [], []\n",
    "\tfor k in np.unique(labels):\n",
    "\t\tif k >= 0:\n",
    "\t\t\tlabel_inds = np.where(labels == k)[0]\n",
    "\t\t\tbest_guess = np.argmin(dataerrs[label_inds])\n",
    "\t\t\tcluster_As.append(guess_As[label_inds][best_guess])\n",
    "\t\t\tcluster_Bs.append(guess_Bs[label_inds][best_guess])\n",
    "\t\t\tcluster_dataerrs.append(dataerrs[label_inds][best_guess])\n",
    "\treturn np.array(cluster_As), np.array(cluster_Bs), np.array(cluster_dataerrs)\n",
    "\n",
    "# an approximation of the background due to the weakly coupled spins\n",
    "#A_background = 4 * mag * (np.random.rand(400) - .5)\n",
    "#B_background = 2 * mag * (np.random.rand(400))\n",
    "#background_dict = {\"A_background\" : A_background, \"B_background\" : B_background}\n",
    "#learning.store_obj(background_dict, \"background_A_B\")\n",
    "background_dict = learning.load_obj(\"background_A_B\")\n",
    "A_background, B_background = background_dict[\"A_background\"], background_dict[\"B_background\"]\n",
    "\n",
    "# given guess_As and guess_Bs, this function considers all ways of removing num_remove spins from the guess list\n",
    "# it compares all of these possibilities along with not taking anything out in terms of the error from this subset to the data\n",
    "# the error function considers the spins along with 400 weakly coupled spins which approximates the background.\n",
    "# it returns the optimal subset as well as a number to remove next time: the same number if a spin was removed, otherwise the number\n",
    "# considered for removal this time plus one (considering taking several out at the same time)\n",
    "def remove_spins(guess_As, guess_Bs, N, omega_larmor, tau, data, num_remove = 1, error_fun = squared_error, verbose = False):\n",
    "\t# creates all subsets of s with size between lower and upper\n",
    "\tdef subset_size_range(s, lower, upper):\n",
    "\t\tans = []\n",
    "\t\tfor i in range(lower, upper+1):\n",
    "\t\t\tfor j in itertools.combinations(s, i):\n",
    "\t\t\t\tans.append(j)\n",
    "\t\treturn ans\n",
    "\tguess_As, guess_Bs = np.array(guess_As), np.array(guess_Bs)\n",
    "\tM_background = calc_M(A_background, B_background, N, omega_larmor, tau)\n",
    "\terr, As, Bs = [], [], []\n",
    "\tfor subset in subset_size_range(range(len(guess_As)), len(guess_As) - num_remove, len(guess_As) - 1):\n",
    "\t\tsubset = np.array(subset)\n",
    "\t\terr.append(error_fun(M_background * calc_M(guess_As[subset], guess_Bs[subset], N, omega_larmor, tau), data))\n",
    "\t\tAs.append(guess_As[subset])\n",
    "\t\tBs.append(guess_Bs[subset])\n",
    "\tbest_ind = np.argmin(err)\n",
    "\torig_error = error_fun(M_background * calc_M(guess_As, guess_Bs, N, omega_larmor, tau), data)\n",
    "\tif verbose:\n",
    "\t\tprint \"new error, old error: \", err[best_ind], orig_error\n",
    "\tif orig_error > err[best_ind]:\n",
    "\t\treturn As[best_ind], Bs[best_ind], num_remove\n",
    "\telse:\n",
    "\t\treturn guess_As, guess_Bs, num_remove + 1\n",
    "\n",
    "# given a data function data_func(N, tau) that represents a diamond, returns an estimation\n",
    "# of the hyperfine couplings of many c13s. The function returns, in order\n",
    "# the final estimation, the total set of estimated parameters, the set remaining after using least squares,\n",
    "# the set remaining after applying the classifier, and then the set after applying clustering. See chapter 4 of thesis.\n",
    "def analyze_diamond(data_func, N, omega_larmor, kernel = \"rbf\", verbose = False, plots = False):\n",
    "\ttau = choose_tau_params(N)\n",
    "\tdata = data_func(N, tau)\n",
    "\tdip_inds, windows = find_resonances(data, fit_dips_below = None)\n",
    "\tlower_ind_cutoff = 3220 # corresponds to 15 microseconds\n",
    "\tif verbose:\n",
    "\t\tprint \"number of dips to measure: \", sum(dip_inds >= lower_ind_cutoff)\n",
    "\tspin_dict = {}\n",
    "\tfor dip_ind in dip_inds:\n",
    "\t\tif dip_ind >= lower_ind_cutoff: # 15 microseconds and on\n",
    "\t\t\tspin_dict = analyze_dip(dip_ind, tau, data_func, omega_larmor, spin_dict, N_vals = np.arange(0,256,2),\n",
    "\t\t\t\terror_tol = .1/64, verbose = verbose, plots = plots)\n",
    "\t# load the classifiers\n",
    "\tguess_scaler = learning.load_obj(\"classifiers/scaler_svm_\" + kernel + \"_di3220_29diamonds_cxABe\")\n",
    "\tguess_clf = learning.load_obj(\"classifiers/clf_svm_\" + kernel + \"_di3220_29diamonds_cxABe\")\n",
    "\tdef guess_classifier(features):\n",
    "\t\treturn guess_clf.predict(guess_scaler.transform([features]))\n",
    "\tguess_As, guess_Bs, dataerrs, all_guess_As, all_guess_Bs, select_As, select_Bs = choose_spin_guesses(spin_dict, N, omega_larmor, tau, data, guess_classifier, x_min = 1, error_fun = squared_error)\n",
    "\tcluster_As, cluster_Bs, cluster_dataerrs = cluster_spin_guesses(guess_As, guess_Bs, dataerrs, eps = .075, min_samples = 1)\n",
    "\tAs, Bs, num_remove = cluster_As, cluster_Bs, 1\n",
    "\twhile num_remove <= 2:\n",
    "\t\tAs, Bs, num_remove = remove_spins(As, Bs, N, omega_larmor, tau, data, num_remove = num_remove, error_fun = squared_error, verbose=verbose)\n",
    "\treturn As, Bs, all_guess_As, all_guess_Bs, select_As, select_Bs, guess_As, guess_Bs, cluster_As, cluster_Bs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# %load learning.py\n",
    "# NV-spin-tomography Michael Scheer mgscheer@gmail.com\n",
    "\n",
    "# Creating datasets and training classifiers, storing objects on disk\n",
    "\n",
    "import numpy as np\n",
    "import sklearn as sk\n",
    "from sklearn import svm, preprocessing, cross_validation, cluster, metrics\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "from sklearn.externals import joblib\n",
    "import analysis\n",
    "import NV_generator\n",
    "import pickle\n",
    "\n",
    "# store an object in a file\n",
    "def store_obj(obj, fname):\n",
    "\tjoblib.dump(obj, fname + \".pkl\")\n",
    "\n",
    "def store_obj_2(obj, fname):\n",
    "\tpickle.dump(obj, open(fname + \".pkl\",'wb'))\n",
    "\n",
    "# load an object from a file\n",
    "def load_obj(fname):\n",
    "\treturn joblib.load(fname + \".pkl\")\n",
    "\n",
    "def load_obj_2(fname):\n",
    "\treturn pickle.load(open(fname + \".pkl\",'rb'))\n",
    "\n",
    "# decides if a fitted A and B are close enough to a spin to be counted as a fit.\n",
    "# the B tends to be harder to fit so the accuracy demanded is less.\n",
    "def is_fit(fitted_A, fitted_B, spin, A_acc = .1, B_acc = .5):\n",
    "\tA_acc = mag * A_acc\n",
    "\tB_acc= mag * B_acc\n",
    "\treturn ((spin[0] >= fitted_A - A_acc) and (spin[0] <= fitted_A + A_acc) and\n",
    "\t\t(spin[1] >= fitted_B - B_acc) and (spin[1] <= fitted_B + B_acc))\n",
    "\n",
    "# creates a dataset of diamonds with N measurements already analyzed, saves it to disk\n",
    "def create_diamonds(diamond_num_list, omega_larmor, num_spins = 450):\n",
    "\tverbose, plots = False, False\n",
    "\tnum_subsets = 4\n",
    "\terror_fun = squared_error\n",
    "\terror_tol = .1/64\n",
    "\tN = 64\n",
    "\tN_vals = np.arange(0, 256, 2)\n",
    "\ttau = choose_tau_params(N)\n",
    "\tall_scaled_errors = []\n",
    "\tall_exec_times = []\n",
    "\tfor diamond_num in diamond_num_list:\n",
    "\t\tprint \"diamond_num: \", diamond_num\n",
    "\t\tA, B, r, costheta, sintheta = NV_generator.generate_spins(num_spins)\n",
    "\t\tdef data_func(N, tau, noise_level = .02):\n",
    "\t\t\tdata = calc_M(A, B, N, omega_larmor, tau)\n",
    "\t\t\tnoise = np.random.randn(len(data)) * noise_level\n",
    "\t\t\treturn data + noise\n",
    "\t\tdata = data_func(N, tau)\n",
    "\t\tdip_inds, windows = find_resonances(data, fit_dips_below = None)\n",
    "\t\tsuccessful_fits, good_fits, phis_list, xs_list, scaled_errors = [], [], [], [], []\n",
    "\t\tprint \"number of dips: \", len(dip_inds)\n",
    "\t\tfor dii in range(len(dip_inds)):\n",
    "\t\t\tif dii%50==0:\n",
    "\t\t\t\tprint \"dii: \", dii\n",
    "\t\t\tdip_ind = dip_inds[dii]\n",
    "\t\t\tres_tau = tau[dip_ind]\n",
    "\t\t\tN_data = data_func(N_vals, res_tau)\n",
    "\t\t\tsuccessful_fit, good_fit = True, True # did the fit succeed without throwing an error? was the error in the fit small enough?\n",
    "\t\t\ttry:\n",
    "\t\t\t\tphis, xs, scaled_error, allerrors, alltimes = repeated_spin_fit(N_vals, N_data, error_tol = error_tol, error_fun = error_fun, num_subsets = num_subsets, verbose = verbose, plots = plots)\n",
    "\t\t\t\tgood_fit = scaled_error <= error_tol\n",
    "\t\t\t\tall_scaled_errors.append(allerrors)\n",
    "\t\t\t\tall_exec_times.append(alltimes)\n",
    "\t\t\texcept FitError:\n",
    "\t\t\t\tsuccessful_fit, good_fit = False, False\n",
    "\t\t\t\tphis, xs = None, None\n",
    "\t\t\t\tscaled_error = None\n",
    "\t\t\tsuccessful_fits.append(successful_fit)\n",
    "\t\t\tgood_fits.append(good_fit)\n",
    "\t\t\tphis_list.append(phis)\n",
    "\t\t\txs_list.append(xs)\n",
    "\t\t\tscaled_errors.append(scaled_errors)\n",
    "\t\tprint \"create diamond_dict\"\n",
    "\t\tdiamond_dict = {\"A\" : A, \"B\" : B, \"r\" : r, \"costheta\" : costheta, \"sintheta\" : sintheta, \"N\" : N,\n",
    "\t\t\t\"tau\" : tau, \"data\" : data, \"N_vals\" : N_vals, \"N_data\" : N_data, \"dip_inds\" : dip_inds, \"windows\" : windows,\n",
    "\t\t\t\"successful_fits\" : successful_fits, \"good_fits\" : good_fits, \"error_tol\" : error_tol, \"num_subsets\" : num_subsets,\n",
    "\t\t\t\"phis_list\" : phis_list, \"xs_list\" : xs_list, \"scaled_errors\" : scaled_errors}\n",
    "\t\t#print \"store diamond_dict\"\n",
    "\t\t#store_obj_2(diamond_dict, \"diamonds/diamond_\" + str(diamond_num))\n",
    "\t\treturn all_scaled_errors, all_exec_times\n",
    "\n",
    "# train and optionally save a support vector machine classifier\n",
    "# scaler preprocesses the data (scales it) and clf is the classifier\n",
    "def train_svm(train_X, train_Y, test_X, test_Y, kernel = 'rbf', C=1.0, verbose = False, pickle = False, suffix = \"\"):\n",
    "\tscaler = sk.preprocessing.StandardScaler().fit(train_X)\n",
    "\tclf = sk.svm.SVC(kernel = kernel, C=C)# svm with rbf kernel (default)\n",
    "\tclf.fit(scaler.transform(train_X), train_Y)\n",
    "\tif verbose:\n",
    "\t\tprint \"training score: \", clf.score(scaler.transform(train_X), train_Y)\n",
    "\t\tprint \"testing score: \", clf.score(scaler.transform(test_X), test_Y)\n",
    "\tif pickle:\n",
    "\t\tstore_obj(scaler, \"classifiers/scaler_svm_\" + kernel + suffix)\n",
    "\t\tstore_obj(clf, \"classifiers/clf_svm_\" + kernel + suffix)\n",
    "\treturn scaler, clf\n",
    "\n",
    "# train and optionally save an extra-trees classifier\n",
    "# scaler preprocesses the data (scales it) and clf is the classifier\n",
    "def train_trees(train_X, train_Y, test_X, test_Y, verbose = False, pickle = False, suffix = \"\"):\n",
    "\tscaler = sk.preprocessing.StandardScaler().fit(train_X)\n",
    "\tclf = sk.ensemble.ExtraTreesClassifier(n_estimators=10, max_depth=None, min_samples_split=1, random_state=0)\n",
    "\tclf.fit(scaler.transform(train_X), train_Y)\n",
    "\tif verbose:\n",
    "\t\tprint \"training score: \", clf.score(scaler.transform(train_X), train_Y)\n",
    "\t\tprint \"testing score: \", clf.score(scaler.transform(test_X), test_Y)\n",
    "\tif pickle:\n",
    "\t\tstore_obj(scaler, \"classifiers/scaler_trees\" + suffix)\n",
    "\t\tstore_obj(clf, \"classifiers/clf_trees\" + suffix)\n",
    "\treturn scaler, clf\n",
    "\n",
    "# creates and optionally saves the dataset used to train the SVM which classifies guesses.\n",
    "def guess_dataset(diamonds, omega_larmor, min_dip_ind = 3220, training_percent = .75, pickle = False, suffix = \"\"):\n",
    "\ttrain_X = [] # feature vectors\n",
    "\ttrain_Y = [] # classifications\n",
    "\ttest_X = []\n",
    "\ttest_Y = []\n",
    "\terror_fun = squared_error\n",
    "\tfor diamond_ind in range(len(diamonds)):\n",
    "\t\tprint \"diamond_ind: \", diamond_ind\n",
    "\t\tdiamond = diamonds[diamond_ind]\n",
    "\t\ttau = diamond[\"tau\"]\n",
    "\t\tdata = diamond[\"data\"]\n",
    "\t\tN = diamond[\"N\"]\n",
    "\t\tdip_inds = diamond[\"dip_inds\"] # the indices for the dips\n",
    "\t\twindows = diamond[\"windows\"] # where was this used?\n",
    "\t\tspin_dict = {} # The keys to the dictionary will be a tuple (dip_ind, phi, x)\n",
    "\t\tfor dii in range(len(dip_inds)):\n",
    "\t\t\tif diamond[\"good_fits\"][dii] and dip_inds[dii] >= min_dip_ind:\n",
    "\t\t\t\tdip_ind = dip_inds[dii]\n",
    "\t\t\t\tres_tau = tau[dip_ind]\n",
    "\t\t\t\tfor phi, x in zip(diamond[\"phis_list\"][dii], diamond[\"xs_list\"][dii]):\n",
    "\t\t\t\t\tfor omega_tilde, cosphi in calc_omega_tilde(phi, x, res_tau, omega_larmor):\n",
    "\t\t\t\t\t\tA, B = calc_A_B(cosphi, res_tau, omega_larmor, omega_tilde)\n",
    "\t\t\t\t\t\tif valid_A_B(A, B): # this checks for the upper bound imposed on (A,B)\n",
    "\t\t\t\t\t\t\tif (dip_ind, phi, x) in spin_dict: # if (dip_ind,phi,x) is already a key in the dictionary then add another value to the list that corresponds to it\n",
    "\t\t\t\t\t\t\t\tspin_dict[(dip_ind, phi, x)] += [(A, B, cosphi)]\n",
    "\t\t\t\t\t\t\telse: # if it is not, then simply initialize the list with the corresponding tuple\n",
    "\t\t\t\t\t\t\t\tspin_dict[(dip_ind, phi, x)] = [(A, B, cosphi)]\n",
    "\t\tfor k in spin_dict.keys(): #now iterate through they dictionary's keys\n",
    "\t\t\terr = []\n",
    "\t\t\tfor A, B, _ in spin_dict[k]:\n",
    "\t\t\t\terr.append(error_fun(calc_M_single(A, B, N, omega_larmor, tau), data)) # calculate the errors from assuming data only comes from a single spin\n",
    "\t\t\tmin_err_ind = np.argmin(err)\n",
    "\t\t\tbest_A, best_B, best_cosphi = spin_dict[k][min_err_ind]\n",
    "\t\t\tbest_err = err[min_err_ind]\n",
    "\t\t\t_, _, x = k \n",
    "\t\t\tfeatures = [best_cosphi, x, best_A, best_B, best_err] # this is the feature vector\n",
    "\t\t\ttag = 0\n",
    "\t\t\tfor spin in set(zip(diamond[\"A\"], np.abs(diamond[\"B\"]))): #iterate through the true spins in the diamond\n",
    "\t\t\t\tif is_fit(best_A, best_B, spin, A_acc = .5, B_acc = 2): # arbitrary thresholds over here\n",
    "\t\t\t\t\ttag = 1 # best_A and best_B have been found to correspond to a one of the true spins\n",
    "\t\t\t\t\tbreak\n",
    "\t\t\tif diamond_ind < training_percent * len(diamonds): # first diamonds for the training set\n",
    "\t\t\t\ttrain_X.append(features)\n",
    "\t\t\t\ttrain_Y.append(tag)\n",
    "\t\t\telse: # the other ones for the test set\n",
    "\t\t\t\ttest_X.append(features)\n",
    "\t\t\t\ttest_Y.append(tag)\n",
    "\tif pickle:\n",
    "\t\tprint \"create dataset\"\n",
    "\t\tdataset = {\"train_X\" : train_X, \"train_Y\" : train_Y, \"test_X\" : test_X, \"test_Y\" : test_Y}\n",
    "\t\tprint \"store dataset\"\n",
    "\t\tstore_obj(dataset, \"datasets/guess_dataset\" + suffix)\n",
    "\treturn train_X, train_Y, test_X, test_Y\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "diamond_num:  0\n",
      "number of dips:  1067\n",
      "dii:  0\n",
      "dii:  50\n",
      "dii:  100\n",
      "dii:  150\n",
      "dii:  200\n",
      "dii:  250\n",
      "dii:  300\n",
      "dii:  350\n",
      "dii:  400\n",
      "dii:  450\n",
      "dii:  500\n",
      "dii:  550\n",
      "dii:  600\n",
      "dii:  650\n",
      "dii:  700\n",
      "dii:  750\n",
      "dii:  800\n",
      "dii:  850\n",
      "dii:  900\n",
      "dii:  950\n",
      "dii:  1000\n",
      "dii:  1050\n",
      "create diamond_dict\n"
     ]
    }
   ],
   "source": [
    "error_sample, times_sample = create_diamonds([0],omega_larmor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1059"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(error_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1059"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(times_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18873\n",
      "18873\n"
     ]
    }
   ],
   "source": [
    "print(len(all_them_times))\n",
    "print(len(all_them_errors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from itertools import groupby"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## split the two lists using the keyword 'split' and then separate out the values with the smallest error in each group\n",
    "all_them_times_split = [list(group) for k, group in groupby(all_them_times, lambda x: x ==\"split\") if not k]\n",
    "all_them_errors_split = [list(group) for k, group in groupby(all_them_errors, lambda x: x ==\"split\") if not k]\n",
    "all_them_num_spins_split= [list(group) for k, group in groupby(all_them_num_spins, lambda x: x ==\"split\") if not k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#pickle.dump(all_them_times_split,open('all_them_times_split.pkl','wb'))\n",
    "#pickle.dump(all_them_errors_split,open('all_them_times_split.pkl','wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_them_best_times_and_errors = []\n",
    "for counter in range(len(all_them_times_split)):\n",
    "    the_times_group = all_them_times_split[counter]\n",
    "    the_errors_group = all_them_errors_split[counter]\n",
    "    the_num_spins_group = all_them_num_spins_split[counter]\n",
    "    the_small_error = min(the_errors_group)\n",
    "    the_small_error_index = np.argmin(the_errors_group)\n",
    "    the_small_time = the_times_group[the_small_error_index]\n",
    "    the_optimal_num_spins = the_num_spins_group[the_small_error_index]\n",
    "    all_them_best_times_and_errors.append([the_small_time, the_small_error,the_optimal_num_spins])\n",
    "all_them_best_times_and_errors = np.array(all_them_best_times_and_errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4185"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(all_them_best_times_and_errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAFgCAYAAACmDI9oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvX98VdWd9/tZ52BsRSFFn9cM/SGEAPmBSKC00tt5jdWK\n0UcR8WlIAkwpONPb8b4erjP2NcWJAw0PecRn+sPL3Mfb27mA+ChJyFwsiqMRW+q85rmNViWKEAkk\n6FjLzLRgoNLWyDnr/rHPWll7nbX3XnuffX4l3/frlddJcs7Z+7vWXmvt7/6u7w/GOQdBEARBEARB\nEA6JYgtAEARBEARBEKUEKcgEQRAEQRAEoUAKMkEQBEEQBEEokIJMEARBEARBEAqkIBMEQRAEQRCE\nAinIBEEQBEEQBKFACjJBEARBEARBKJCCTBAEQRAEQRAKpCATBEEQBEEQhMKkYguQC1dddRWfOXNm\nscUgCIIgCIIgyoBXX33115zz/xD0ubJWkGfOnIlXXnml2GIQBEEQBEEQZQBj7B2bz5GLBUEQBEEQ\nBEEokIJMEARBEARBEAqkIBMEQRAEQRCEAinIBEEQBEEQBKFACjJBEARBEARBKJCCTBAEQRAEQRAK\npCATBEEQBEEQhAIpyARBEARBEAShQAoyQRAEQRAEQSiQgkwQBEEQBEEQCqQgEwRBEARBEIRCWSrI\njLFljLEfnjt3rtiiEARBEIQvHX0daHisAR19HcUWhSAIS8pSQeacP805//rUqVOLLQpBEARB+NIz\n2IMUT6FnsKfYohAEYUlZKsgEQRAEUS40zW1CkiXRNLep2KIQBGEJ45wXW4bILF68mL/yyivFFoMg\nCIIgCIIoAxhjr3LOFwd9jizIBEEAAE63b8HAvGtwun1LsUUhCIIgiKJCCjJBEACAkb17gVTKeSUI\ngiCICQwpyARBAAAqV64EkknnlSAIgiAmMOSDTBAEQRAEQUwIyAeZIAiCIAiCICJACjJBEARBEARB\nKJCCTBAEQRAEQRAKpCATBEEQBEEQhAIpyARBEARBEAShQAoyQRAEURCoGA1BEOUCKcgEQRBEQaBi\nNARBlAukIBMEQRAFgYrREARRLlChEIIgCIIgCGJCQIVCCMKSjr4ONDzWgI6+jmKLQhAEQRBECUAK\nMjHh6RnsQYqn0DPYU2xRCIIgCIIoAUhBJiY8TXObkGRJNM1tKrYoBEEQBEGUAOSDTBAEQRAEQUwI\nyAeZIAiCIBQoDzNBELaQgkwQBEFMCCgPM0EQtpCCTBAEUSDIgllcKA8zQRC2kA8yQRBEgRiYdw2Q\nSgHJJOqOvllscQiCICYc5INMEARRYpAFkyAIojwgCzJBlCmn27dgZO9eVK5ciembNxVbHIIgCIIo\neciCTBDjnFILOKKKhARBEMR4gRRkgihTSm27nioSEgRBEOMFUpAJokyZvnkT6o6+WTLuFVSRkMgV\nyvJBEESpUFI+yIyxOwHcBmAKgB2c8+f9Pk8+yARBEOOHXLN8kF8+QRBBlIwPMmNsJ2Ps3xljb2r/\nv4UxdpwxdpIxthEAOOc/4pz/GYBvAGjOt2wEEZbxYuEif2HCj2KNj1zdhkrNL18wXtYNgphIFMLF\n4lEAt6j/YIwlAfx3ALcCqAfQyhirVz7yQOZ9gigpSvUGHBbyFyb8KNb4yNVtqNT88gXjZd0giIlE\n3hVkzvk/ATir/fvzAE5yzoc556MAugAsZw4PAXiWc/6a6XiMsa8zxl5hjL3yq1/9Kr/ClwlknSgc\npXoDDgv5CxN+lPL48FvvSs0vXzBe1g2CmEgUxAeZMTYTwAHO+TWZv78C4BbO+Z9m/v4TANcBGASw\nFsDPAfRzzn/gd1zyQXag6lwEEQ7yVS1faL0rPDRfiPFEyfggh4Fzvp1z/lnO+TeClGNiDLJOTCxy\n3TGgHYf8bXmXUt+WkixeRJGxlNe7UuhzVYa45CEXEWIiUiwL8hcAfJtz3pj5+34A4Jw/GOa4ZEEm\nJiK5WtDIApc/i1gp9W0pyeKFlBFAZWtr2VsnS6HPVRkAxCJP0HwhCzNRTpS6BfnnAOYwxqoYYxUA\nWgA8VSRZCKKsyNWCFocFrhQsZbmQL1/VUrJulsN1VmUbD9bJUrj+FbNmyde45AmaL2Rhjo9yX1vH\nE3m3IDPGOgF8CcBVAP4NwGbO+Q7G2H8E8DCAJICdnPPQ+YTIgkwQxaEULGVE/inEdSbrY7wUY25O\n9GsYZ/vzff2Kca1KbXyUjAWZc97KOZ/OOb+Ec/5pzvmOzP//kXM+l3NeHVY5ZowtY4z98Ny5c/kR\nmiAIX0rBUkbkn0Jc51LNPFGuFGNuTvRrGKcFPd/XrxjW/nLdYSipSnphIQsyQRATjVKzxhDERKec\n5iRZkO0tyKQgE0SZ0tHXgZ7BHjTNbULbkrZii5MzpbSIlpIsOuTeUhxKeUwQBGFPybhYEAQRTJTS\nvuOtGl4pbcOVkix60E5cW7DjORgo7radbt+Ckc7OkhkTURnP15wg4oYUZIIoAaIou6Vc7SwKpeTX\nbCtLIRQOXVmPy9+zlB4C4ibutqnHsRmf+chFHAfj+Zp7UUr9T5QXZakgU5AeYUu5LI5RlN22JW3o\n/2r/uHCvAEor0MdWlkIoHDbK+ngruJErcbdNHs8yV7M6LkpJKVX7pVzWxlwppf4nygvyQSbGNWH8\nNSe6j+FEb38USqXPCumXXOg2l0ofh0GVGUBB5A/bTxPFl70cxw+RXyhIjyDgvTia/j9RbhheTPT2\nlzOFVAIKPU5oXNoRtiph1DFTKgpnqchBlB8UpEcQ8N4qN227jectZxsmevtLHb8tcTHOAeR921yO\nj3Q61vN4BapWrlyJdIKhdwEPFcQ6nrBxhwhblTCqS5Np7YwSZJwr5DpB5BtSkIkJiUkZjNsHthg3\njVwoJR9gIhsbhUD9TL58TKdv3gQkkwDnsSonXoGq0zdvwuqNFdjRmCjpjC359Om1ufbTN29CZWtr\n3h9yTWtnMTLqlNoD/UTx6Z5IkIJMTEgKoQyOtzRsRHGxUQjUz0ilqrMz9pt2nMqJUCy29s3wDFQt\nh4wt+bRo2va337oWlwJnOkcxrk/UNTxfiixZtMcf5INMEHlivBXyKAbF9DPM9dzF9pGUuXuBkvbf\nVX2MhWJfDn6l+vUt9vUOgny5HfLVD6V+/YkxxrUPcrmkeaMtl4nNeEvDVgyKaZWJcm51ztt8Px9r\nhDgmAOOWey7nzIe8Rqt3gfssCvnKT+1Fru0uNZeEYpGvfiAXtfEHWZDzCD2xE0RulJsF2c8aatpR\nyMca4XfMXK3K+V7TbPo8DhlszhP0mUKPTbqfEEQ8jGsLcrlAT+wEkRu6VaaQgY9RLELqnNe/b/JJ\nz8ca4XdMqRwjuCKcyWKZ7zXNps/jkCFswGNUWeOkEPeTUrHOE0QpQBZkgiBKFt3q2vBYA1I8hSRL\nov+r/cUWLxSl4JM+UFcPcA4whrqBY/6fHccWyzgsyOOR8XzNCUJAFmSCiEC5pWYb7+hW13LIZuCF\njU+6sOANLbvDypIX1uJX2dLiWCFbWoI/O453wGysvxPRp3Q8X3OCCAtZkIkJiZc1jyyUpcV4bJMf\najU0AIGWvIlo8ZuIll0iPmj82DGe+4ksyAThg1eO4nK2UI7HvMvlmAkkFz9OYcGrmDPHypJXaIuf\nzQ5LvndhwmQXIZ9aQofyFdtB/VSmCnK5pHkjShcvRbgcFTJBOSv3haBQylIuCpzY1q9++imr7f1C\nuwHYPITl+0EtzEMB3eTHL1EfxMiNxA7qpzJVkDnnT3POvz516tRii0KUMs/cB7RPc141ylkR9mI8\ntilOCqUsFUKBE4r1gXuWF9Rn3uYhLOqDms0DTNhtX9O1IKtyfBQzZkN/ELO9rvpDJY0HMxPRB1+H\nfJCJkiYnP6j2aQBPASwJbD5bWrIZGG/+tqXWnlL0qVNlAmAtn/A9TjGgdeOksvSZ17Hxp47D53q8\n+G0HjedCzL9ixmzo7Yt6XcfLeCDsIR9kYlyQk9Vv8TpHOV68Ln7BEL9F0m9ruhyza5SaT3Qxcyrb\nyBRmPA1dX40UA16+rlJaa8vdEmZjeY9j2zfqMeLu31yPFzReCjH/iunWpe+YRb2u5EpAeEEKMlHS\n5LR43fZdx3J823fjFwzxL6ymm41Q4rqPd5eUsmlDqftEF0KBCKMEhRlPDyx5B60bJ2H7jb+TSkK+\nXUj0tsStMBYi9VpHXwdurdqHfTtaPa2uXg9Ncfdvru41FbNm+Y6XQsy/UnLrijo2wnyv3B9CiXCQ\niwVBWFIMlwGxhQlA3uxK4WZUDErZpcVLtpy2b5+5D3hll7MDoj3kmWTPtwuJ3pY4t6YLNbeCXAL8\n3o+7f6MeL0q/l6J7UdwUYgyNR3eMUnOFKwTkYkGUBKWwjR0XxXAZEFaglpqWkrHU2BK3tSVuC15Y\n61cU66KtVdjYV6/scnzoX9llJXsu1lWbeaq3JZcdFP18hZpbQVZV8f7WvhlZ10Pt3zjGdtTrFaXf\nJ0I2j0KMofHojlFqrnClBCnIRF4ZT5MvypZlrjfSXLYwi/VwIto80tXleVOO0i/F9h31G8testkq\nQUYFJoQPfa5ttJmneluiKngdfR3oOt7lOl/V1CrXa74Imk/i/eoXh3wVykIqnPo8jtLv41Gx0ymE\nS8l4zOxQ6q5wxYQUZMKTOBSs8TT5oiirtjfSfCizQgnpOt7l+n++/OikYtzZ6WxDcu55U46iYES9\nOcWlzPiN5VxvnEYFJoQPfa5tLNQ8Pd2+BcvXPY71vSl5XgA4de6U67UQ+M2DIIXSVuGMY14HPbzY\nzOdSVOziXvNKyR+6nKB+86YsfZAZY8sALJs9e/afnThxotjijFs8/fF8fCMJN7a+f2HSJdn6jF27\n+1pwcDAwvLH2Dfn/fPnRyeMyBiQSvm0upE9kOftf2speiDbajju/z6np6fbvWiPfP3DPcsw8NIi3\nb5iL2x/Znxf5dQrhT5prGrTT7VtwtrsLBxsYLmxYZez3MPOulPxNi5kijpjYjGsfZCoUUhg8rUo+\nvpGEG1vLTRgLnm2C/OaaZiRZEs01za7/x7ndqlqB5HFbWvKejSAMpWg9s8XWMhyljWEteJO378Hj\nD36Iydv3ZL2njkEbV5QrW1pdClr1i0NIcue1UBTCZSdXy/zI3r1IpDkaX2eeCq1oBzgPHCul5PI2\nnnYXiXEK57xsfz772c9yonD88tvt/Fj9PP7LdTdz/u1PcH7gL4sqz9afbeULdi/gW3+2tahyFBq9\n3cfq5/FjNbX8WP28UMeR1/Pb7VbnMbFg9wJ+zaPX8AW7F4Q6dxR5JiK59onf97/fWs+P1Nby77fW\nW53rzbo6fqymlr9ZV5f1njoGo8xLv3OX2riIOt+iEKbtNp81XRvb/xHEeAHAK9xCxyxLFwsBpXmL\nH7+t2lJLcUNbdA4227Amgq6nTf/GuWVbauOrGMTlKiGOg3Ta5QuuHvtofT0SaY50gmHesWOB/e8n\nWy5yB3231MZFObvsmDDNc1pb42W8jZlyZ1y7WBD5w29Lt9QioYu5RdfR14GHV83D0fr6kkgazzkH\nBw+1dRp0PW361y/AI2wwYKmNr2IQV0ChOI6uHKvHntbcAiSTziuC+9/PjSMXN5agNucyLvIRkOrX\n1nxmjslXcK1pnpP7Q7xMhDR74xGyIBOSjr4OTN6+B0v7OaY1t9CTrg8NjzXg8Qc/RJIjVstWFEuD\nsATqgU/FJhfL30S1uMRtQVaPE+XY+bIa5+M4JuK0PtvImU/La6lZ0gl7CrWeTdR1MyxkQSZC0zPY\ngx2NCazeWEGTK4CmuU14YWEC6QTQu4C7LEZRrEh67uAzXZ3W3z/YwJBiwAsLEy7luBh5kFUrl2r5\nC2v9kjmUu7qM75dTydcwskaxxIrjDy27Q57HdBzxPwDW8vhZvkzvRbku+Sz1G+euhI0VMJ+W11Lf\nYSmnOVloChUoTJbqeCEFmZBEWdz9lLAV+1dg/u75WLF/RZxiFo5n7gPapzmvGm1L2nDvnqNYvfFS\n7GhMuFwbokSKi4UtzR1L8MGFzPr7Fzaswpr7L8WFDatc//eSI5+K89luR7E9293luimEXrjFzpbH\nDle53AhOt2+ReaF1WU+3b8FAXT0GausClYq+Ly/Bsdo69H15SdZ7oi9GT5xwncfLDci27063b3F8\nmAGplBmzligKm2xrZ6fvscOi59i2ve5eiok43oF7lnvOBX2emNqrK4Vx5JT1UjSjPODEhY3yWy5z\nshgU6uGh1B+iyg1SkAlJlMXdTxk8OXLS9RqFuBeWUMqhls7OJEtc/nuVK1cixYDnFzGs2ngJdt9S\nYf19r+vmVZ0sKE2cqZ2210FYsw82sKz2hVq4GXO/agQdL66HAJN1NgyqsqDLOrJ3r3wACFIqprx3\nDizzqlO5ciXEYwRXztMz2IObDqeRSPNsOSyuhZQvmZQKpjp28mEVM1039SEDjMWiAAhlbuahQc/1\nS58npvZ6KYW5rFtBiqbX+/l88LVRfkk586ZQDw/lnNKyFCEFmcgJP2VwduVs12sU4l5YQll3tVK/\nJllMymmUB43pmzdh/6412H1LBZprmmOpbORVnUy/Znq79L87+jpwpsttufNSALys2WEX7sqWFplT\n2UTQ8eLK96pbZ892d6HhsQas2L/CShmRSkNra5asqiJRMWuW73HOf2oqeOZVZ/rmTXj3Kkc5fvcq\nyPOMuQEx17lMfWdSrkwKjxg7W/tmGK9/xZw5rtewmK6bS7m3yLFtg2jb2zfM9Vy/bB50vZTCXNat\nqFX88pnj2Eb5DZqTxXD5KhXo4aFMsckFV6o/lAd5/BN3DtRc8nvGJUuh8rqqbfVr98nbl/FjNbX8\n5O3LjPIt2L2Af6+1jh+prZX/C5sLNte8qmH7LK48ruK8J29fxo/Vz+Pfb63n1zx6jfwx5X8Oc269\nH6OODa8cxbayhM1nna9cwCZ5CzFf4s77m2+Zc8ldXKy80nHmTC8lSi1PNxEMLPMgF13JjfIDYBmA\nH86ePTvOPiOInLFZLHXlohBJ+f1uTkHKThxKS643x3wXZ7Btj+iLO390p+c1C9NW/bx+7YxSTMNW\nlrBjMF9KQZzH9TuW3l69n0pd6cllPokxdqS21vN656P9NmOs1PvdRJi1qdDto4IvZmwV5LJ0seBU\napooUaL46hWi/KvfdnHQ9p/JZSSsy0QUv2yvjBi54LXNa7slLvriyeVPerrB+Lkh6K4pej+aMn+I\nQDIRADnS2ZkVXOZ1PfR+9/I3X/DYy6H6Ll++jnG6VPkdS59zQW5HhSCMC4LXfLLxfRbxDn6BwPlo\nv43rWdB5SzFTRpi1qdDjqpRKi5cjZakgE0QUCuEDF8VXL99J+YOq3RUisCOKX7Z6M4lLRq8bRhQF\n3OtmLdpa/eJQ1s0w6AZpyvwhAsnUwEe/4DKTLKLfvfzNbY7nd7ONS3EJmxrQb077XdOgOSe+WzFr\nlq8McSpsYZQZr/lko4Cp8Q5B7S+0z2zQeaMqmPlUrL3WJlu//nyij/OJ7AceBSoUQhScOMsThyGO\nJP7lmIi9XMvGxt3XUUtyexGlNLPMyMAYKlvcxXjE5ytmzcLo8LB8Hbq+Gg8seQdNc5uwvtdRDtT/\nhWmHLpP42+Z4fvM2H0UsbI4Z19j2Oo6fDPJaAqHbbRobcayLXnOmHNctE1HbUYwiK6W47paiTMWA\nCoUQJUuxtn3isNSWY67Pci0bG9VqfLp9CwZq6zBQV+9Kzzaydy8SaY7G11mgAmK7VW1TmhkYy1s7\nffMmIJkE+FjqNWHVES4UImPG6PAwKleuRPWLQ3j21F1oW9Imj3n7I/vR/9V+rO9NhbKM6X06ffMm\n7NvRir++7hRSPIUFj70sj2eb49eULzksUa1tQWPb6zrq//c6jp8Mfin8gtpmmxEnLF5zphzXLRNR\n1wR5fdLpgrlnFGvd9Vu7yvVeUCzIgjzBKJb1ttRkiCKHn/Uv63MRrTW23y0XS1Ex5JHWIpWMkqNa\naP1kitPiJI4lSoHrVuAUd2S9uzeNxteZSz6h2HjJEYec1+6+FjyTSblz20VZPh2A1bFVGYTMQddb\nn3sNjzVg7XOjWHqY48qW7JR4UfHqn1z7raOvA5O378HSfo5pzd5rAZBttTvdvsWpEMm5Mf2fLWHm\nVqmtC8VgopTqnijtzAWyIBNGSsFpPw5LiR+2flZh+0JaXxIJ401GLxcdxVpja+nx+lypWYrikCdS\neWEAYMzJx8uYVJjrjr6J0eHhQJlEXuKg/MSqbF5y6kFRwgomlGMGhiRL4sKGVag7+iaqn35KWsmi\n5sT1k1UP8FvXexGd2y5ifW8Kb98wVx7P9tiVK1fKPrYZ+x19Heg63pUVJLf0MEeSBxdNCYNXG7z+\nb3M9AWft2NGYwOqNFYEKpzEAUCvAYjPG9c+EmVtxxxmINfbAPctLLmjOizjjDHL9bBRs72uUczk+\nSEGeYORzi6VUIoxtFd+wfWEbQCJuflEWqFBKSYgbfyGJO/uESRHwu1lM37wJdW8NoG7gGKqffgpI\nOMvcSGenlUyn27c4bg4AfnfyhDyHaXyrssnfM+dR5TEFRYnx51cYJkixCaP4eAXk3dwPJLnzevsj\n+1F39E3sbEzi1qp92Lcj2MLpet9i7KvzUvRH25I2XNnSmvNYCcoSosoc5Irgp4CGWTt0g4Bp/Kkl\n2r3Q5alcuRLpBEPvAl7woCuxxs48NJjVR6VyHxAIeQCEfkgI8xCSb+OE7X2NqunFiE0uuFL9oUIh\npUW+c9XaEkfuxyjHKFYOz3zkutQLZIRpU2Au35pafqy2zvqYpn71ygNr+qw8Z+a8tvIfranl32ut\nk+cwtUs9n+s8IeaA3tdR+twkj9d7T//5Ha7xYtu/v1x3Mz9WW8OfXvO5rPEm211TGyhjPnMu57oG\nZV1Py/PafNav3d9vredHamv591vrQ52jWMU3RFue/vM7smQqRp53P3IZE3GPgVwodj+OJ2CZB5l8\nkInYGE9+buUU7ZsPWbP8eEP4s+nZGNTx4DquV2aACD6sWXLrPqe1dfL3urcGPGUWVj09s8P63pSn\n36jpu2HmgNFnGojkQxiX/6GpfwfqagHOpLuI6i+cz7l/tL4eiTSXPtx+rlm5yCHaXDW1CqfOncoa\nW7lm8RC+3gwMb6x9w3hu07Gjvhel7XHEhejXwG99yte4yXVOEuOb2HyQGWNNjLErMr8/wBjbxxhb\nFIeQxPiiHLZ2bP244nZFyee2Yy6y+vnNIpl0fHhDbnuLcWDy9VX9gz0zAyjfCco1rN/MPV1PWjPb\n962tRplNOZdFpoi2JW1Gv1H5XeF329WVlblCzaLhhd7XWX3+zH1A+zTnNQC9/VHznpr6d2jxFKQY\n8PKiS7L8hUU2jFur9kXa7hfX2dRfBxtYVmELr3FhWoPCxiScHDmJFE+h63gXVuxfkfW+vsV9un2L\noxx7jGmBCIQUr2ob/OIy/LbW/b4X5trHGZsSJs97vtwS4sqhXmruIkRhsfFB/hvO+W8YY38E4CYA\nOwD8X/kVixjPnG7fgoG6egzU1hV84bG9EcQdSJjLjcB0o7O9uUaVS9xU1ICxsJiUVdU/2HRM/Tth\n+03IvbMx6VtpTr/x5RQMl9mFS/MxX1Aht0jZFlS4oXLlyrG0bnqfv7IL4Cnn1YehZXfg/c5OvPOJ\nFO5pOAzAPN7V8dTR14Frd1+L+bvnBypSDywdRevGSdjemHD5C4u+nLx9j9XcMlXzG+ns9OyvCxtW\noXXjJOxqnJRd6a6zM1CBCRuTMLtytvzfyZGTWe/ril5Q8K6gpaYFSZZES02Luw0BYzvMA7B6PbuP\nd1srvaZzdPR14OFV83C0vj6nddpvfcpXzERcx82HAl+OSnc5yhwHNgqy2Pu7DcAPOefPAKjIn0hE\nyeFjwYpioZLWuMzvhZx8cViGXfJaWvfCLthqv5pu8HEt3PkM6otiufErvSxYsX8F5u+e77Lu6QQp\nRXr/5RIMV9naihQD3r0KWL7ucVcgoK0F3vd6Ll4HsKTz6sPoiRNgAK7+9ZhiZxrvat/0DPZIi2aQ\nIrW1bwY6t13E1r4ZUqkf2btXWtCX9nPc3ZvGE9tGfeeysZpfhnOfmooUA4aur3Z9RwQz6oFuMnuG\nzzywnfNqGXGhJKvKcthdC6/j+wXr2XzPD/V6cnCrdnu5V/QM9uCmw2kk0hxnujoD1/go94J87TrG\nddx8rI+llmnIhnKUOQ5sFOT3GGP/N4BmAP/IGLvU8nvEeMHHghVla06mhcr8Xm6TbyQTcT7S3WVt\n3Qu7Ba32q+kGH9fCXepuMSb5hPKnWvd0bEsJx3HjE1kqPn0G0u1AtcCLMR656Mht3wU2n3VefaiY\nMwccwL9cNabYre9NofMhJ4WbQO2bprlNYHDm4uRLJuO5P6rDsdo6DC27I+v41S8OIcmdVyA7a8u0\n5hY0vs6QSHPfuay3Vf7d2opvrP09WjdOwgNL3pGfN60xon8rZjvt1NPx+e2w2DyQP7n8SRxZewRP\nLn/S8zO6LDa5iP0ybMRhKFCvZ0tNi5Vi7bWGN81twgsLE1nuLWGPE4Z8GksKpcAHtaEUMg2FpRxl\njoWgKD4AlwG4C8CczN/TAdxsEwGYrx8AywD8cPbs2bFGNhIeHPhLzr/9CedVI47I2jiif/3kUI8f\nR9T3L2+/mh+rreG/vP1q377RCXNu234tVGRzKUVz3/mjO/k1j17D7/zRnQU/v02WCP29N+vq+LGa\nWv5mXXYGjbhk9TqOiOA/UlsbOEauefQaflTJRmHMBqL8zzNjSMj+UTGNZ78x7pWhwC9zgXyvppaf\nvH2ZT4/ES1A2haD38zXXg47r9b7+/zjkyzULiZ8M+cr4oZ+zVLI5Ed7AMouFrUL6CQDXAlgkfmy+\nl+8fSvNWnuRjofdb/NQFK8y5PW/oIZRilaxzH/hLvnX7TL7g0fm+8vgpFgt2L+Dfa63jR2pr5ftx\n9+/Wn23lR2rt05fl8yZng835oyqmUdrml8Ir7PGCFGGZJi6j/P3y2+38SK2Trm7+o/N9+/XOH93J\nn/1irVHGdclRAAAgAElEQVRJ9usr2760bat+/Z/+8zv4kdpa/vSf35H1WfHe91rrXO3yU95F/9im\npYuLXPuxWCndvMiHPPnso3w9YOjnDPvgGAVK+ZYbtgqyTRaL/wLgDQDbAXw38/OdWM3YRNGJGu0O\neFfo8iIf1fz8ttTVqmhh/Pk8XT8st7x1ss79yi70XDEZKXDfvggqWKBnFIh7q7NnsAcHFzqZBGy2\n2KJsxwX5XIfBpkqa3qdDy+7AgIdrQa6FTy5sWIU191+KCxtWWcvqhddYUI8jipyIV1GkZF3vRTz+\n4IeYvH2P8dhPLn8St/zzAOrfGhjL+mHhAmXrImXbVv36zzw0iCR3XnUeWPIOWjdOws7GpNEFQ90a\nl4GTw8OObzggXwtBrkVf8lnkKQr5kCeoD4LGmp9MuQZee7lO6Of0G3txuREWrSJuiIw644HAPMiM\nseMA5nPORwsjkj2UBzk+csmlK3KAphjQunFS4DHizLnphZoHU/pKhswNm/e8zs/ch47hfeiecjk4\nIKPb9b4xyaH24frelOt9m/4N+oya13Xfjta8Xy91/DXNbcrL+dQ2iXEh+swvT7Kap1jPgRwnNuPN\n5jNDy+6QyrE65tV8wm/fMBfVLw5ZjW1xzsTllyN97hwq5sxxKhSGkClMu8XY3No3A9UvDuHcH16O\ny395Dh98ciqm/usHWfOg+3g3ODhaalryniPZNCbLKf97IdbeONHlLWZf55JjPG65i3Yd26c5MTcs\n6RiJyhTbPMg2CvL/C+DPOef/HpdwcUEKcnzkMuHE5FcLKxR78fVThoqB3wKpKocArBRFzweaZ+5z\nAgYXr/O1cD+8ah5uOpzGCwsTuHfP0WxZOzudvK4tLb59prZLWPGiXP9CPzTpbRJKpVD+sgoNdHY6\nH8yxAIefXEinXaWag8ZsWKXtdPsWnOnqlEU+kjy7PX59FFRsJSqexV20/8dVBCUKfgaEYsoVFr92\nlKLyXEoFm8IoueX00GSDbM91n8L0GT8PvL+UOrEVCgHwIIDDjLFexthT4id3EYlSIpftJ1NhhWKj\nbufqW165uJNExZX5QkPPKqAqx17baE1zm3B3bwpPPPh7nF7fOPaGZVaNpf2OgrS0P/sBOSivq9p/\n6tZhLtt+YcafTaT76fYtOFZbh6O1dThwz3IA2VufHX0dmL97Pq7dfS2++adJtP71pehquw5AdqEB\n1eXgdPsWDNTWYaAut/ywAjUbRDrB0LuA46wYLz5bsn79vbMxidZvOa4Hoj9G9u7F2zfMxe5bKvD2\nDXON7g5+W8GqW4IYAwfuWZ51LcLOr6Hrq42p3TyzXeguGgXY9vXbui92hH+Y/vZrR9zb9nGss2pq\nwSDynS40TBGacsvMFIRsz0vvRXIvLFdsFOTdAB4CsA1jPsgTo3eIssXPl60Y/luV1RcAxp1XuBdW\nVTlsW9ImleOqqVW+/nSNh9NIcIaRnzkpsU63b8FA9x/i9KuV7py5BgViWnOLTM2VJWvADV/tP/Wz\n6s03zM0x7I3N5uYzsncvGJwFzuS7KtoBOPli1eppQHYfqONJnpc7iqxoZ9QbtPSRnzMHqzdWYEdj\nAgcbmOv8qo+0OM+Gn3wcgJOeTS/ooI9x0WfVLw6h/6v9uP2R/cb54Xftq59+CpWtrRgdHpZFQWYe\nGsy6FuLck7fvseoP4UespnYD3H3e0deBW6v2Yd+O7DLfA998Bqd/PjnwoTAX/B7gcq0iKPCbM+rY\n0j+nX2v9fduiQnH7FJvW2dPtW3C0vh4Pr5pn1Vd6akE/CqGUBvW9oNgPTXEz3tpji42C/FvO+XbO\n+SHO+YviJ++SEaEp5Wo3qmzFljPKjcBWZs8SuN/4Cupa/h3Tv/EVAP5Kunjv1LlTvlbVyi/McJTu\nLzjWFWebnmNk+Ar3E77Bquz3ABEmWEgtGLG+NyXlVdsXpCwLa+n7mapoerCn/n2bxbpy5UpwAGk4\n/rZe7QAg88YKOvo6/AuDKOWyDzYw2c6oN+jR4WHn9cQJbO2bgSRL4sKGVU6/dnVhoLbOFXgnqs59\n/qURAMD50fOyoIM4tz7GbW9wtkFSS/vTAJz8uOmEu8SyOPfSfm5VNtw4H7WHOq/54hQdYhgZmhxY\nSEUnzp0kvzLUtmud35og5sjZ7q6sz+n95/VwpI9Lvf1xVw81XdeRvXuRSHPcdDhtZaAIo5iF+WzU\ne9Dk7Xtcga5e95JSzy8flvHWHltsfJC/B+BDAE9lXgEAnPPX8itaMOSD7KaUfeFU2QCUrJxe2Pat\n7edM/n7if1VTq3Dq3KksX8AgH0FPvzdLv+SomNqsyipu2F5+hA+vqsdNhzkYd57Y9WBPWz9Er/b7\n+QOqfS4Kj6jnCfIl9AuWtJVT+nwDTgGdRMIdXOrB0I1z8cCSd1A1tQrXdw9iaT+XOwL5CggSbexd\nwMHBsfQwN/qxA2N+3YmpU5E+f15Wz7Sa9+3T0PGJKeiZcjmaalsBuINX9UC+KG2Ny7/1dPsWnO3u\nwsEGhgsbVrnmZpg12W9NWPvcqOzrCxtW+a4DtoFtUdofR0Cm01fAjsZkYGBlvoh6rxSBrukEw7xj\nx/Io4fjzYy4l4gzSO2T4N+ec3xhVuLggBdlNqU2orEAnw++lIKcNtn3rFRxl0/agG1bQ+87NpxPP\nNzD8U3NNVhWwfAXhBCmmIgvB2zfMxe2P7M/6/vzd8wE4Vd9ueT2RFezZ0deBydv3SAXQq/9tg71U\ngrJn5CPY0ySPmnoOwNj5uroAzmX2DJtxKI7PAaQzGStM/a6OhwWPvYyZhwaNnzWNu46+Dixf9ziS\nHJ7Kgpr9Q6WytTU4oPOZ+9Dwq16sfT6NpYc5rmxptVLu9P7Rgy+92p/LfPAbXzYPbaa+EO/3LuDY\n0ZgAA0OCJWKZu0EP4lHbGoZiB99FvVcW8h5bygavcic2BbmUIQW5tBmPEzzKTTXIer5i/wqcHDmJ\nKRVTcOGjC57HPnDPck8lRj0PB9C7iOEv9riVlqg3Jds2658T1hYOgMFRpFZvrAhlGTcpjl5jyevm\nFUZJMlp3DekCo6S+s7F62mQFUS2WsytnZx1LvI80lxb5aways06o4+HxBz9Ekjuf3b9rjdFaa5LD\n9OAn5N7aNwNVPz2BX0zj+MyvnTEgFH2bsagq4YA7xZ4+Fzq6b0XP797FEw+lkFCyc9hm3ogji08Y\npUldE1q/NSmrL8T7UeZMEA2PNUiLtP7gEbatUeXo6OtA1/EuMDA01zSXRGB3qZEvZbzUDGnFgBRk\nIjaiLoKlMhH95AgrYxgl06QQAdkWZGFBBYAja494Hi8oH+/p9i14v7MTDGalKOp1tG2znjru4VXz\ncPNraTBAKu0AQt2YZZszrgcVs2ZhdHjYqMCGtSwHfVa3GKvnurVqX1afeMki+j3N0+Dg1g8oXv2u\n5h0HYEzXpipBLyxkuFd7WNKV7Fk/GQQH8NKSSnz/hg+c42oWY9u82arC97XnPjKmE/TLX2yTYk+/\npg2PXoMUY7i79yIaX09aWZBt+joMpj7ytSBnUikO3TDHtWNik2YxF3ldDx5FtAQX24o8URmPhquw\nxJnmjZjgCB/S7uPdoYJa8uXYH2fWg7CBVV5BGSaZRL89sOQd336YXTnb9ep1XDX4RJVXBNvsbExi\n+Ma5shCETtQgHNugRpE67ubX0jhwz3LsbEyid5ETxPX8ogR2NiZl/t2zhnR3JmTgTUsL6o6+6QS0\nKdfL5vpFDfTRj62OZ68AJF0WYSlL8RQ4OO7uTeOJbaNWadG8+r1y5UqkEwwvLEzIdG0/+9wVmL97\nPr7Y+UU0PNaAqqlVAADGgZtf4zI4VmQQONvdhUSao/F1hupDJ5yMH4xh+42/c51f0H28W64BQf2m\nyq2nExRtXfDYy9iz7SOs7/UIvPNIsWc6HwBs/fFl6Nx2EbOTV7jmWvXTT6HurQGpHIft6zCYguy8\nxuf0zZucHSXOUX3oBDofuoj1vamx7wCeaRZzlbdtSRuubMnu0yjkIkehKgPmIyi82IHmuTBRM1JE\nwdOCzBj7JOf8lwWWJxRkQc4jSmBXx5XT5OIPIPYn/rBW3LBPwHFakG1kEtvvJh+/sEE7up+nSd5c\ng9ii4OlrnbH26YF2pmAjU2CXF17uCbm2Kc4KaSbLn7g2gFMp8a67O7Ouv/6ZoAcYk1zqLgTgzNHO\nhy66XFMASMvzCwsTaHzdyT4hLbQA9j26xtgfNrscNhZUYdlufI3LnY79u9Z4upDoQW9e2M4rfa7E\n6ZcvXKVmV86W/v82a4/uPpTPnbd87+pFuXb5Jh8W00JbYUuxiEs5E4cF+f9hjPUxxrYxxr7EGJsU\no3xEqaOkBhOWx5aaFiRZElVTq6zyddqiWllsUi+FfQLOJaWZLapMfmnahPWvdwEPtMT3DPbg4EKG\nFIOTF9dDXltLTJx5Qk3HUq19b98w1yWTGEO/3bAaa+6/FBc2rAp1Pi9rfJjrJ8bm0LI75OvydY9j\n7XOjxpRTYceGahXU060Jxdc0dtXrpueMNc0lr+u4vjeFzm0XcXdvGk1zm1zp6CpXrnRZni9sWCVz\n90oTCWOeuwxi7oty6CZMFlS9D6umVmHpYS4T613Z4mSoEAVbRIo9kQ/aNle56Neh66t914+glGi5\ncOrcKdcrYLf2mCzkcaKOe5EiMF+5gsOmcSsE+bCYBh0zbgtzMXL3Ez4KMuf8PwL4EoCfAlgBoI8x\nto8x9nXG2NWFEY8oGovXOfXWldyi4uZ56twp42Tt6OvAma7OUMouYFYu/RaCfOdkjJIfNWj7Xf2c\n7c2/aW4TdjYmsWrjJW6FUssRa+s6EeeNQha3yLwCbuuUV1XFfLt5+CEUy9ETJ+Rrkjs+0XFt80bJ\nxbq+NyUVUJPLhlqMxOscLTUt0n1laT/H5O17cLa7E72LEti3azWmb96E6Zs34Uc7V+PRWy4BMHbT\nfX5RQrqx6Ii5ACDwuvldI3GcoZEh+dAn/OjVgi1qbl8vd5Rrd1+L+bvnY/7u+Vixf4V0L6o7+iYe\nWPKO5/phssLFuc0f9Vj6ehZ3wQvXuM+QL2VcfQiL23Uiat5qvX+jHEdXeG3zhcd1DQvljkK4CRWk\nxxirAnArgFsA/CHn/PP5EswGcrEoDl7bPXp0tCmQKeqxC4nYhlVTKwGITa6c29g+zbHusySw+WxR\ntjVNW4ylHvwhFPjE5Zcjfe6ck5/33DljMFSUnMom9C19V+ozxpzcwF7pwTLuGr0LGXY0JrLmkVcG\njt4FjgVPZKVYc/+l8ntBae30doo0Y6Y57Bdop8unu2jpuX5NGQ1M40l1R1ER8pncHLyuRRCn27dg\npKsLnHM8vygReW6JvtTTFwZ93ib3tm1p9rPdXfjFNI5Pn2W+qRLjJCjLSb6ChfNxnDjd+ojik5cg\nPc75Kc75I5zzOwD8UWTpiLLGZAU83b4FT2wbBQPD/l1rAi2pYY6to2+V22xjhdnyEnJzcGmN0i3b\nB+5Zjjfr6nDgnuXWbRPYtFG1cmTJrln3497WtOkrkxVT/1+YPhftFRbBqNXNOvo6ssouq0GMdUff\nRPoDJ0tD+oMPpJI60uUOGvSyAHkF4nmV9tXnwND11dKlgWeU44pZs7L6SQQxphlwYcMq4zzSx6Sw\nas2unA3GncwhuiVPyLO1bwbuursTz566yzgOR7qcym03v5b2nMOTt++RgXZdx7uyrpnwOV++7nFs\n7XOqPTIwYyGMlpoWvLH2DZlWruGxBgxdX210R1ErH86unO2Sz+TmoLfddj1yqvQ57iBR55a6qzbz\n0KDVNnmQdTLsdrvYtfrmn03C6o0VBVPY1Lmi/h7VXSAuK2qU48Tp1leK2K7VcVafLAcozRsRCi/r\nhV/xA6+0V1EtqFkFCEKm77K1bqq5Vl//6uddMr9ZVyctdKYcs2Gq14l+UlOYqdZ3GWzlIXvcFuS4\nLMFhjiOsOn/79xdx9a+Bd68CGv/Z6VdTJTcva1zDYw0yp69MAaZZjNTrWv2TQfldNU+uaezqBUsA\nGK2sXu0+3b4FZ7o6keCZvNAAfvToGmPgXu8f1eEzWj/oMpz60hxjH+jpAE0Wu6A8uzb5g4/W1yGR\ndgdkPnvqLpelUE0nZsr3a7Lm2eZINq0hce5AnW7fgrOZAMbnFzH8dsPqSLl+1/emsPQwx9s3zM1O\n56bNexuFKkobC7kzF5TeMqoFOS65KNAtm6hBruUKpXkj8oLX07/pCdvLChfGgmB6YhXnqpgzJ1L6\nLluqXxxCkjuvutX37Ru806kBAF7ZhdM/n4yBbz7j+1Qut9IV39iRvXtlmq6qqVWuACQ1mEkwffMm\nzDt2DPfuOWq98PtZDML2lX4s8XfFrFnBgSx19RiorcPWvhlIsiSuzhSU+Myvxz6njhcxpryscU1z\nm/DCwoRjpU2njVbcB5a8g9aNk/DAknfGAqRaW13HMVmAegZ7sKMxIa1wQpal/dx1fDk+NcvwyN69\nSGYsuynmKF09gz1Z/d3R14FPnXH64VNn3H3WM9iDmw6nkUhzVP9kEJ0PXURLx0uu8+hz0DTfxDkP\nNjBp5VXHgle/qBxsSDgBpAuZbL8659f3ppDMWHsrV67E1r4Z6Nx2EV97bhQH7lmOgXnX4Ae7P4bO\nbRelhVlcQ7U/TWvA5O178PiDH2Ly9j0umdqWtEmXjqhWLjF+AecBZvVfXxpaOQbGAi53NibxN9tq\nsx5mdJ94W3/VoB0oU3+J76zvTYXefQuLKaC2FCyqFOjmjWnNN42jieYLbaUgM8Y+zhirybcwtjDG\nljHGfnju3LliizLhUBW3ILwULb9JpitbfpHx1U8/Zb3oBi3QJoXRJL9wrQAcy7Gpoh0AYPE6jAxN\nBjjzvfGp76kKv7pVLGR/YMk7ANzBTDq2W2B+QSRhb2b6seSNf3jY8zjywSCzg1X94hD6v9qPS+fM\nAQD5CoyNl6qpVehd4FTmS3DIbA0qbUvacO+eo2BKNgldoVDHX5i2qu4JA/OuQeLyywEAv7lqEpC+\nCLzzz67+Ezmbz3R1YufXvoB0psrhtNZW7N+1BrtvqXDJsLMxiYbHGtB9vFsGsokHMDE+t/bNwAsL\nE0gnMm4GBgVLzx0s+uhrz41K15PpmzehcuVKLO3nuPk1J7jvTFenHDfTN2+SWS68xtKFDauw5n5H\neRTKF9JpKZcok41kEtM3b5IPnDcd5ph5aBBIpTD1vRH5EKoGLarXS+SR7jreJfvi5tfSxvzKwkoY\nVRFSH1hNYyfMFrOawcQU2Bz0oN/R12F8IA7Cr/1RlfIwBClRxVJUoyp3vu5uGco917KXQUC/TlGD\nrMuVQAWZMbYMQD+A5zJ/NzDGvEsSFQDO+dOc869PnTq1mGJMSEw+fupNRb/JqgFcoliBl++jngUD\nKNwTq1faMn3RmHloEEnuvPpy23dR2bIq0BIrrcM3zsV/avoX7NvRiumLz6Hp3DkkM/6aAvE70/6v\nYnvzsbUSq77BD6+qN/pd68eyObZ+Yxaf1Qs7AO7sKTsaE0gzx7ra+DozBpep1muTf2/URV58r/rF\nIcc9IfOAfvmvPkKKMfT87t2sNgkL6+dfGkECjk8xAOMcENduXe9FLD3M8fJ1lXhgyTvo6OtwfJJT\nKVT9dBCP3nIJfrRztetcPHM+3d96+uZNaFvShrt707j5NY5EmrseZBJpR8EUcopx09HXgcu2P2G0\n0nr1o/DZHROKZ40LV4ETxlEx5SOA8azCLKpSInyOxevI3r3SC1m4uqjjPpc1Qx2XpvGrK+t+iP4B\ngDR3HhxMMl22eLHxIc2U3cMGv/ablPK4FbGg+VUsK2TUeW/awbKJTciVfBwzDBPNWmzCxoL8bQCf\nBzACAJzzfgDB5kNiXGKaNGICc7hvsq73DcEaOmreX3FzKtQTq63CGOhakaGjrwO3Vu1zFN5MQQvT\nTUi1DkvF9pVdaDtzFv1v/8LV7rYlbTiy9ogMZjKhXh8/a5et5VTcHE6OnMRNmTRi4uHApIwFHTtL\ngc3cqG0QbRPV40zXSrde65X34kBXMt5efAWSnKPp459xfW765k3Yv2uNkw2FOz7HLyxkUtnVZRIB\naCJd2+dfGpFj4mCDMy+eb2Dyf5WtrUjDOW7vIiZTpon3j1+3BAO1dRiYf62r5LeqsApXj9aNk7Cz\ncWzcdB3vknLcdDjtaz3z6pfK1laXFVp1Bbr9kf2o+85tqL7tfdR95zZp0Tale2yuaUaSJdFc0+w6\nj1puXR33uawZpmOr6Mq6jinYtGewR5YYV2VS10JT39o8EJtY35tyVeZTMe2+FVoRC3t9ih0Ypo4t\nr/tEMXIt55uJZi02ERikxxjr45wvYYwd5pwvzPzvDc75tQWR0AcK0otO1IAFzwpq3V0Yqv8IDyz7\nOJqUVE36dwBkfT9XmUoRz/RehmpZAMaC7JZfjbaBnwHgwOf+FLjtu5H7JY6AirGgsDR+Mc3xiX37\nhrm4/ZH9saRLihIQGKYyok26NsB7TPr1i+01GairBzhHGkDL/ZNwd29aVrETMorrv6MxkRXUVTW1\nCkMjQ+DgmF05G6fOnZIBUL0L0tjRmAQDwxtr3/AMQBRwAPVvDbhSj91/3bArxZq4rkKOgwsZdt9S\ngSe2jSKR5kgnGH60c7VV+0Xqx5tf40gAnopn2P616X/9M7msL+K7puqY6vt6Cjq/dHrqGJQV9QLm\ngU0bxltKsqB1ptTlJ0qPOIP0jjLGVgFIMsbmMMb+DsD/l7OERFGJ6gfm6YrQ/K+4/Zp/z7J66vhZ\nF+N8YjVZZPwsYHFbKXRLu24NEP34fmcn3u/sRCLNHZeBt14CwJ00bpnsF14BSWFl8MLPr+6uuzvR\n+FoaiTRw9ftJl991HOmSvKwkftfD1n/a62FuYN41Lj9T+btScjmon8T80dObGeXOGCEYHGvgjsaE\n3FkQ7RFp+hgcZXT/rjWy2Mr13YPYs+0j3N2bxpPLn3S5eSzth8u6qgaWquhmENGH1YdOYO9/4+h9\n+ytZPtq/3bDa5SstrNgHG5hx/fAK6ll6mMsbjbhmQfMtaC2wWb/0z+h/qzIEySO+a6qOqb7P4AQs\nqunnvNqijlXbCpumduvzd7ylJFPXGdNaFZcFvJA+v2WNVqhqPGOjIP9nAPMAfAigE8B5APfmUygi\n/0T1L/JcfA2V9wTF8KUynTPIvSPOwBH9pqjfhCpXrgQH5LZ3OuFYFE39uLSfewYk2cggIte9Fv4g\nvzowZrx5ty1pw7On7sJdd3d6HluXVe8Hr5uz3/XwUwCGlt2Bgdo6p6yuzxgAY2PHYJntcmbeNjf1\nkzpv1NzYptLVFUrgYYI5S27X8bHqeKpvbnNNc5YypV9/tQ+mNbe4Pq9WOKxsbXWu26IEhm903FI+\nkclKIb7POZfV+gT62E3zNLqOd+HF5rmyTLgesKg+NPQM9rgC7q5saZV9K65ZLsoqELx+dfR1SL9f\nIaPIlKIW/tHznHcd78L83fPx39c0uOZM0PnE++L6iQcZU5yFqW22FTbF+VM8JasJ6i47xjkVUqkp\nJWVRHY+mOR2XK0KxfX5NFNu9xMgru5xCVa/sinyIUhpffgQqyJzz33LO2zjnn+OcL878/vtCCEfk\nj6jWWk9rw23fBTafNeb8LYYvlU0hC2BsAaqaWpVTQELYCT998yYM3zhX+oDKBP6GfhSBSAnuKGEi\nUMhWmQ9a+L2sV0PXVyPFgKEb5njevIOOHfbBQ83W4JXqyzQGxfdEOd3REyf8x0BLy5j1rqXFs9Sy\n3k8i8O+uuzvxF4ecLBYio4sI4BSlq4XMvx866cg0PCwVKYYxX+LpmzdhWnMLGl9nRr/Rac0tUhkW\nfXFr1T4MXV8t/VcFoydPyldV6VJTbqmW9ecXJaRV2ITwnwWAkyMnkeIpdB/vBgCXJVs8NIhrJpS2\ns91dzvUaOIa6twZkpg59vunKqin1nIq+fq3Yv0KWnlblvrs37biaZKzl6gOjOsaE7ILrX/3QNa6D\n1kvb9dRvPtjED7QtaXPJCWSuXdD6GlKpsVUWg9a9OBWh0+1bZJYUNSWiGuuRC8X2+TVRkqnpfIxh\ntpTiw4gJmywWixlj+xhjrzHG3hA/hRCOGB/oCk0hnh711FmqYgUgy+o1NDIUcER/okz42x/Z79rC\nBswWg+mbNzmBbJy7smds7ZuBo/X1eHjVPF/LW9DCrypSk7c/ITNVqPmCTRa0jr4O9C7gY9ZvDVFd\nUU3HJq79gXuWG2WVW/+ZtG9C4TC5majtlZbhDBVz5hgVadMYsN1i1tO3ff6lEQBjGV0++ORU8Mxr\n25I2ObYONoxZ4AFHsWyuaZZKsiq/yQVhZ6NTZGNno6MYieOKVGnqmBMxJeLVM6g2870LG1Zh1cZL\nsKMxYUwnplauE68cHN3HuzF/93z0LkgjnWAYur5a+ug6QYUwKt5ergqqnMItI8kd33x9nOjjvKOv\nAydHnAeDkyMn0dHXId0dVKu7cHXxGmMAcEniEgDAi5+9NC/Kkp8lWlWygxRpNUjwwoZVwePXQqlR\n12VbZTFo3ZPvd3WFrqqpr2kyS0ombSAQrwIZp6tJXPc4m52SgluYfYxhtpTiw4gJGxeLJwA8CuA/\nAVim/BATgHxMwEI+PZoWUH2rXJSVDrJa+aEW8zD1l59FSL1JBxViefuGuTK3avWLQ64S017ftVn4\nRT8szeTErfrJoDEzAABXdL5aOENH+Naq6diCCn14LZwmNwO1vdLafeNcVLa2YnR4OKvUtOh7Pddt\nEK6b3TP3oXLWb5BmTkYK0XcAMPVfPwDLvIr/392bxtJ+jl9Mc66TUPDblrRJdwtRLIQDSKdSOHDP\nctm27uPd6Dre5RqbpmweQsZ3rxJFSMaWduEiIRRgtY9VOUQ6Mf3mnmAJtNS0SD9nBiatyjsak3i+\nwbGer31uVFqZdzQmpTuGiuirNE+7Mj2oBT4AOG4ZyaSjaCtWayD7YUkdQ1MqpsgUbAmWcKzuwr1D\nc0m/7vkAACAASURBVC9Rx5jo7zRP48jaI/jfHu+3UpZMyrrf37Z5lYMU6TfWvoEja4/gyNojdruA\nFkqNui6aDAwmghQd8T4yrjxivvkpkCbXF1NBHcA+a0+hieseF4cvfilS6n7vAhsF+Vec86c456c4\n5++In7xLRpQE+ZiAhXx6NN1o9PN/7bmP0PXQRVkwwTaQCHBX3RLp2kyKdlA/erkWCMSCIgK32pa0\nSaXwhYUMVVOrsPa50ayqZEFyC0su4Fg2hV2KYSxXr+rDrOaBDbJu+Lk4CEVf3NgeXjUPR+vrZT/q\nC6fuZgC4r62wdp8cOYn3M8qv8KsVCtVl259Aw2MNLn9bmzEobnZnujrRMbwP0z87gtXfSmJHo3ur\nW3XBEP63S/ud3MOf+TWyFHxhCUzxFBpn/gPSzFmQZx4adD24AZAWVbV4xe2P7M9K1fXps8ylmKou\nEhwcXce7sLMxKTMnHLhnufTVFenExLHOdne5rnf38W601LTgjbVvoKVm7DoI2ZYe5jI4raVmzDda\nTX3WM9gjFWyhTKtBdGLuPHvqWaz8Kyb7mCthhmMPS2mXu0ZLTQvOj5539a/q3jF98ybXljwwtosQ\nNSYjKBAwaN57ve+nGOVLEfR7cPAKxgxSdMT7auEaAK4x5letLaioj63VXUddc/Kxk1moexzlKs4v\nNmnevgygFcCP4QTqAQA45/vyK1owlOYt/8SZKqkYqGmvTFXvGh5rwOMPfoikmAaZRW365k3u9EJX\nfhmnf/APGBmajMrmFrlIqymVKleudJSvdCZqX0mzFNRvUdKdqfIBGGuH4Rin27fIymaVra3yBpVi\nQOvGSQCAI2uPSP9U6a4gchSnUkgnGFq+5fwt0orphEnBZrwGhjR4Nla8ruNdYGDofOgiEmlHnXp+\nUQL37jmKo/V1SKQh26qmWbun4TBOjpzE7MrZeHL5k8bjn27fgjNdnU66s8ZJ6H/7F+iovQ5dv/8F\n/vbvL+LqXzsBeKK4iXotexc4luP3rnSygJjaLlKDibRqL19XifWP/gyAM36rfjLoZMBgDJUtLcb+\n0PtLTUsm3A/E8d+7Erg6U8pb9ImaQkvMmRcWMnBApnrb2ZhEkiXx7Km75LnuaTiMP+4+LtPS3f7I\nfqzYv8LVp2obBaIyopoyraOvA8vXPY4kh2tcAnBdH9HW3gUcOxoTLtmv3X0tOLjn+FTnjJq2bt6x\nY57jy28sBq2PuaSr85svtukVvY5hK5cprZ16PVtqWiLdB0zX0Csdni1h7k2mNSeuY5cylBLPIc40\nb+sANAC4BWPuFbfnJh5RLti6AJiw9cMK+lzQtqXfcWf9xL/yXdPcJlm6t7K11WWhcG3fDe/DmZOX\nAUolMsBtKRAuBUI5VgNJ/G6CauGMMBYH3dIi2+FVQCPzMHy2u0uWbGbcUZyET6PJ4iPaqPqTVldW\nG/tbTZmmX0+vbUeT7C6rbcB1FmNRbKmnEwzPL0pgduVsDMy7Br+YNlYpjoHhwoZVsoCF6ruqoo8x\nxpzvNtU0A5vPoq35WbTUtODqXzvWdhEcCLjHhCjH/Mp31hgDC5/YNir7f1fjJLRunITzo+flfKh+\ncQgJOOdIM8iiG2HSkgmEpfczvx77nGrJFwhr/I7GpMs6rFuYR/buxeI/WOxKS6f2pXgV41S3LuuZ\nHtqWtMlCPC9+9lJXe9R2iDF6YcOqLNnFuNTHp+ivyZdMBuAEVv5imjMHfjHN30ikj9vT7VswUFeP\ngdq6rLLYqquIcB0R75vWOZPrktd5BWGshl7HsLVsm9LamTK4hMV0DW3uLX73CptASTUo22+91ClH\nV4Z8psSbKNhYkI9zzmsKJE8oyIJceMI8SR+tr7ey0ARZT3WLia0FRRyXAxi+0WxBtm1nw+5rsTZT\nAvjKFnPEdBSLTxTLse15VSvQN/72LXwmo8z1LkpgR2MCndsuSmvd/l1rsL435WtdEJZaAP5tEWgW\ndKfgCMfHqmdjdHjY14rhstreUuF7nb3G5EBtHQAnjV7z/WPWyCNrj2RZNXULsnrNOh+6KC3oqzdW\nuM4ztOwOjJ44gQrFgmyLfu1FO4RlUz54dXYCjOGdKzk+fQbSmitwW0OBlm+5La+L/2AxegZ7sOEn\nH8fnXxrBB5+ciqn/+gEqV67EzsZkVt+p/WkaE+p4a5z5D1kW2/m757v62u8aCYTV+uBCht9uWI3u\n490utwrdUmk6nt/41K+32mdi/AOOIiSKsHgVm9HHuV78xGu++81103fisPblakG2eV/0m/671/0h\nqkzAWB+KaxbWmhu1eFI5WpBN440syA62FmQbBXkXgL/lnHtrOEWCFOTS5uFV83DT4TReWOhsdRt5\n5j6j64JKVDePoMXA7319GxEIXvi98LUgr2/EyM/eQeUXZmD6zt5QxxV43XjVNghlGMkk9u1oNSoC\nNsq6quhO066X6M+KWbOkAgwAI11dSHOO5xc5yq5QOINucrnelEQFO4Cj+f5LMnmnHUUujALQ0vES\nRk+cwL9cBXzzzyaFurnabqEDCPz9TFenfKBZtfESJFgCaZ6W6cwaX2d4bkE6S3kWsn5/Vb10l2ip\nafF0U7BFVUgBf2VYVNTTHzDFZ8V2d4oBa+6/1KXMmtwlTIqOeo4XFjJc2LA6S+lXXQbW96ZcD2GA\nExCozhXTPFDdlZxOdn/O65qHcdUoJ3RXL/G7l8tELkaBMA/OQPm7COYCKcPexKkgDwCoBnAKjg8y\nA8Cp1HTxKJeBr/siGmmf5uTnZEknyjpPGEsQi8pphoU6yFrqR6hFONP+069WYmT4Cke+xeecfKWL\n17mizvXj6grp0PXVMiWbejOomlqV5Sdq00deeN3g1PLFQo677u6U1jbVUm1zk/Pyq7W9uTmKzB5U\nVl/AzsZJ6Jk61dMqqtPR14HF33xcuiOIgi7N90/C7MrZeKR/oZVswicW8H/Q8upTYQUHnAwdqpW1\nbUmbnGNTKqbgwkcXkGAJfJT+CIDjOnNzP3DqS3PwwJJ3pBLK4biMgHNpFReKo1cpZb1vVMUYGLPA\nd/R1SOuvavVV/YvVNprKWv92w2oAcJ2jpaYlS9HRz6P7MK+5/9JQCtSCx17GzEODeO9K4NNnWdYD\noNf8K/Q6HNX6G1XeqNZk4ZagK8s2889WHvV8JtmiWoyLRane28fbg0WcCrIxJL4UMllMVAU5yhN4\n0ADPx8S0Wpyeuc+oDMaN3mfqNqm6RRpm8fVCbbca0GTaohbK8ED3HwJpZy5Wzvktpn92JOuhQe9P\n0YY041i9sQJpAOs83EDi3Fo0uRWYAqxk+zMPIvrnvfrW5GpQd/TNaNvQhvGlKp11bw24jjN0fTXu\nv24YAND94EWpGDMAaQAt909yuV0EySYeBESAGwDX2Djb3YWDDQyzK2dLa77q9nDXuidclsqgIL+b\nX3M+e3BRAo2vM1cg5t2Z96UnOWMYusFRnvUgOgaGBEtg8iWTcX70vOshV1hqxbmeX+S0TyixLnn6\nOQ42JHBhw6osxchLmRZjQL73sU+j53fvIsVYlsVY7XPx0CsCNC9sWBXoj6qOQXX73qRcB82hXN0W\nbD8XJEcUFw+b3bRc1g5VWY5TUY1yTeJQ9vKlyMbhcpcPyu1BI4icg/QYY1Myv/7G44coElFSyAQF\nGeTivB8ll6ckhqTjNuh9Jv/OKMeiDWrOWT0Axxa13X7ljkf27pXtr1TSl40MTTYm9RdpwdY+N4qj\n9fVOYB9zrG4pcHBwVzowL5kA+0BHMW4mb98jAz5Ghx0FUryKzx1cyJBiwMvXVQJwAqFkkRPt821L\nvMtUi3PqFcJM48lv3J5u34KBv+rFadzrHl9KaWkRyDLS1SXzMwv+5SpHOf7gY04p8FM3jgW06ePJ\n1L/CJUINcNPHhshjLSrdAXCVqxZV/gBktVMN8gOc8yTgLOpL+7lMAygC3l5srpGlpgEAiYRROQac\ntGopnpJp00TAXUefU4TDda7DjqJ82fYn8MS2D3G3Kk8aMk+3nqZLpKBLsmRWEJtQjhkY2t56CU3n\nP0CSc9e11/tc9A1LJnHvnqOhqtqJKm0cwAsLE1lrVkffWOlq03p2un2Lscy41/n85l/QWi0qN4pX\nHa91V+YKvz47wFamXevsxNH6OtecVI+nB36pKfz09qhBc0Hrj/j7wD3LQxXYCLrHmAL3RP92He+K\nnCpPXXfiTLtXqgU0Jmo6OU8LMmPsAOf8dsbYKYwZUQSccz6rEAL6MVEtyFHIpwW5EE+XXv6aYWX1\n6gfRBgaGPds+sk7/E8XP2es7QSnphJzqVvkn/persXPhIHqmXIGmmmbr7Usr/3Clv/b81w9dCwC0\nlGN+1iKv9pqsJWrKtmaL9vj1v9fxVR/qs91dMi0cSyZx7g8vx+W/PJflwxrUPybXCi9rvmjjd/7+\nIj7za+Ddq4BXvrMmy4qp+meLdn54aQIVv/0I5z/lBNqpwYM/2P0xTH3vHABnV+TWqn1ZgWn69djZ\nmETX8S7p4vDCQobZlXNQ9dMTeL4B0p/5ksQleO1PXnO5RIxZq518xbqf++TtT+AmxR8YcO/IeFlw\nkUxi5V+NjbYjp/4FAAM+d7fng7SexjBsLIPqCmT6ftAa53XdbNwP9OPlakH2Qv2e7h8sdjsYh8zC\nI2IVVDn0OWUKfrSRySvw2uT/Hbe1Nhf3OYEqk5hn48W6Gjem2JRiu4/E5mJRypCCXBpY+RrniLow\nA4gcyRwULBTWR06Vy3RDCYOqpCdYQvqDqsF0OxuTuGz7E2jMbJUbt5M1twKT7+T7nZ2O20BAhhHZ\nThnw5sakTIQJojzb3YnnGxj+qbnGlTNXdyOJGtCjPlSd7e7Cu9PGMkHsvqVCCepyHhR8z/fMfVjx\nywP440NJ+Z0djc4mnJ9rhZpfunLlStw807EM6llE2pa0udxXutquy/KPTfAxf+jhjE+yyEohFLw0\nACSAU1+aiweWvGPMN6yO9bPdnUB6TDECILO/pBnwiyuBT59xCtL0LPsEzo+ex5SKKTg/eh4MDNWV\n1Tg5clK6VExrdsaEns0iSLHTFfeewR40nTuHtjNnAZZER+N9nrl5bbIb+J3fFJMQxt1Klf3Vf3tV\n+oqLBwz9nLm4W9h8N+h7JgV9xf4V+OPu4/IamhQ/v3zbQb7rfvKJv/XAYWDs2pqyyEQlTp/aKMcq\nVV/jfOCX3ahYxOmD/GPO+ZeD/lcMSEG2J58T0saikev5dWUnTCSzilewUFSiWBK8lHGhFOhb3rpV\nRVhBb34tDQZ3QJIamDb9cxdwGvdm9ZW6YFXMmWP1VK9aAT48cWLMmhyyD3WlQ8+CEDYQyvbmpKb8\nE33WunGSfBhxnc/LEtk+DfNnfBKdD6WyilmIIDnVGqkqa6rfs5NRg7sC0xiY4zOcTgOcy2OLzwjL\nnlitRy+7BJN+9xEOLmTSrYaDI82YywqoK3ptS9qMuyVpAJxlSjzD2aFJp1LynGqfAW6FT7Ui+hXt\nsNnFEj7Z8oFPedhrOPPjLGulUPa29s2QSqnXmhD0IK+PtaiW2jfr6rJ88cMouDa+u1GtzGEU2nwH\nZoXNRuSVcaUclc1S9TXOB+VsQfbzQf4YY2wagKsYY59gjE3L/MwE8Kn4RCUKge6r6eU3ZVvcQ8XG\nPynXBOWq/+L0zZvw9g1zsfQwtyqrrNK2pM1RAmLy85q+eRPq/lsjpuNhNF063bcfRN9O3r5H+hie\n7e6SJVeFv1xLTYursMLbN8x1ybu+N4XGjHIsfCa39s1wfGk7OwHOHD/mxeswsnev9IMVcqn+16PD\nw1bXZfrmTahcuRKjw8N49ypHoUoju1RzkD+e6l+p+ljOrpwt29b54ChWfO1xPLxqHnY2Jn1L2aql\nkP2oXLkS6QRzZGfA2/MuIpmxfurtRMJZFrP6ZPE6zP7oI+lrLRRbADg/el62S/TvlS2tYzd+xe+5\nuaYZSZbE1Iqp8tBL+9OOAp9Rjg8udD4vlF8RJDh841y0PVSLSb/7SF5XMT6eX8Tk9zggS16r4w0Y\nm68cfKwtixLYv2uNnF/iuqbh+GKLkuZiTDbNbZLXumpqlewHdew31zTL3w/cs1yWLvdShlSfbDk2\nlBgFIffsytlIgqHp3Dn0HO9GiqfwwJJ3sH/XGuy+pcJz/oliI2rREWBszMqxtvgc0D4NVYmPAfD2\n9fVCFDs5mOkvLx9Ydfyr80Z932ttDfJTNn1PuBZ4FQDRsSm8EYTf/SSoDQK/wjBAtHtL2PtclPui\nH1F8jeOWoVCI61f99FO+a3kp4ueD/L8DuBfAJwG8hzEf5PMA/p5z/n8WREIfyIJsj+7fGiXSORfC\nPOXb+OKWgpzSAqJsA/ulqjNtFU7evsfKF9h0HADuYhICJduByY3B1F49RZyN/H5laG2iyoHsbWu1\nbcIybsoEIpBWSs7xxn9otAv2bJ+G0z+fjJGhyehdmJTBjVe2tGJnY1L6KIv0aHqf6G1Qi1qo2RhM\nfvOuFHgiQwVjeG4hXGnOxDGdh6FM5onMWJ+/e760LLvchLpvxZ2b33asHow5in7Gav78IoYLX/gd\n2mbdJfvIz4InrzUcy7I6F1UlTmBKYaf6egblFhb99X4mC8WpTGEf3SddnuPcOYBzdE+5HJw5iuj1\n3YPG/Nym62aKP5BjNpN6sWHmZ7IyZ8SBKZe4yT/Yz8qr94uNEhulRHTOuch91um4LNRRLMhh7x+l\nYPEtBRnGCzlbkDnn/wfnvArANznnszjnVZmfBaWgHBPhEOVjH1jiZOfzskyoT7ZxRufqEex+zDzk\nXx5alzNO1FLJQUgLyJQrjFkn1Cd+ESUPANOaW6RlRpQivrBhlbWMss2MyRLXgGM1HLrmoqufp2/e\nhJZvTcLOxmRWOWXxft3RN2U2Ay9rjkn+67sH8fiDH2Ly9j3yc1ZR5Vd+GW2930XbmbNZFqrKlSsB\nxlzZBLIixndfi46/qwKeuc+xxnKO5vMfOFvxCp7jd/E6x8LOGZb2K5k/OjsxefseJx8I5zg5ctLY\nJ6plrW1JGxJszA9ZbYsqt7GfhXGCc1lmWrjZCGVwZ2MSvYvc2TxmV87GzsYk/mZbrSsDC2b80diC\nzrn8PAPQeDjtPMRpffS15z7CneufyLJKCYs7ACQ5UP3ikHyvO2O1Vek63uXKDrBi/wpXv+m7ICam\nb94ElkwigbHziWNwcLnrIOZcT0Y5ZmA4OXISNx1OI6GUgtevv5dFNCvDQu11aJj5GVRNujwvkftt\nS9rQ+DpzyWqaN17jTxxD7ALYZmMQ57BVjgH7zBte+K3TcViogXD3Fhu54vh8PigFGSYangqygHP+\nd4UQhMgv+gLstTipi43tFlgQXil9TAvt6fYtSGR8It++Ya7nMaMsilYoW+FByD6taTamqlMVpJG9\nex2FKJl0yRx0kzClUZq+eRPq3hpA3cAxuR2eYkDvIob7lznbwqpyLtwXxKtvWzyUAZP8S/sdxXJp\n/9gulF975HUf3ucUh9GUNSBzXQeOof6tAZmqSyhrvQscZSAFjp4rJgOv7HJSSH38M+iZcjk6aq9z\nHct0c1+xfwUazvwYQzfUAMkkpjW3SL9b0SahMN90OC3TswlM25xq36ljW7+hdfR1YO1zo+jcdhFb\n+2agsrVVuroItwsAMs3gXxy6HADwT801rrG++A8WI8mSWPwHi7PaOXSjs70/dONcZ2yIc3xhRtZD\nXM9gT5ZSqV6H1Rsr8PwixwWjctZvHH9gwFUCmmFsnqhK88mRky6F7PZH9qPu6JvY2ZjE/N3z5Y+e\nHsyUPk+cx5Umr6YZTbWtUkkEHJeGdILJ79quX/qY7fnwNFKM4VT699KaG4eRQEVvp5ABGHMbYmC+\nc1L9v5+LkXyAAkIrpOrY1vvTZss/7Dpto4Tbuhr4HSusXNM3b8K+HU7gYq5jIaqrRN7ueUWiHFxG\nKIsF4ULd9vrHU/8oI9b/Z+v/jHxMr5Q+pq1L220km+05G1cNnTBp2cIcCzCnprMNthF4bffqAVFe\nhVDCtt2vLWGDLuR1B8Oz/3AOIycu85XPFADJwJAA0HT+N9JdwDSe9CIUukuAX6AP4B0EGjQ+TbKI\ncagG1HllPVmxfwX+y8a3XEFeekU6PUWVX1CXX1DW2FZ/GgcbEq5iJa4MBR/+DqcuuQRNv7mAtv98\nynW+V/7tFePOhFcgnF7YRLiViEIjUVJQegXfBQXlec29qAUubNcJvznv5wZh+p4+703ElYbTKy1f\nlGxCXtjIanuPiDv9aFzHI1cJh2L2Q84uFsTERLUSiCIB4jUqunXSz1ppu41kYx2ycdXQMT2lRw0w\nFMcCnBRjvQu4TPvk2Y5n7kPH31Xh4VX1OFpfj619M2RgEgNDiqdcFgzxFP5fX6pCkiXl1rzaf7Zy\ne7UzSznOuKGMDg/7WjR0C87WvhmO9fSlKowMX+Er34F7luN94e6yd68cM9WV1QBLAJ/7U2mxV8eT\nGgipFqEQn5lS4dQ/0gOv9CBQr4AvdXyKcx24Z7lsp2lsi3EoUsMJK6cpwPDUuVOuIEAArjGu/i7G\ngl9BBjG+TNv1bUvacO+eo1i98VLsaEw480S5/uK4pyouRYoxx5UI7iIv13e75xYDQ0tNi2e6R7Xf\nxcOCKDTiNZeD5p8efCfGnVDc9aA8gdca4tefftiuE35rl58bhOl7YudBDYg0HRMA0jydk/VTt7SL\nnauDC1nOu4yqrEH9bXuPCHPt4pDN1iJKrhIO5dAPgQoyY+zHNv8jxgfqImCzNW9CXyj0hVX8vb43\nlbWg2G4j2Sx+Iprcy1WjUAuaMTo/Q1Y7XtmFnism46bDHIk0R9VPT6D/q/14cvmT0tdVPYa4KVf9\n9AS+9two7lz3P/DwqnrsbEyOba8b5Da1XW+nXmVOuoqon9dYsX8F5u+ejy92ftFVlfB0+xZU/fQE\nkhyo+ukJ6WvsdZyZhwZllo7KlSvHlLVzp4yKnhhfIivI0v600aXowkcXAGQrU/pDx53rn8DXnvsI\n039wAG/W1eHAPcsBuMen6PuZhwaR4il0H+92bZWL46rj8N49RzHvmOMaI9wTVDeFprlN2H1LBf5m\nWy12NTop1aqmVmUp4OI7fsqdON7dvWl0PXQRd/ems+ZLR99YhTiTn3BHX4eTVxlwXIkyiLYvPexY\nO2cnJ6Pl/AeBNxRVWX1hYQJIOlla/OZy0PzTXVzCuCj4va+nYwxaK3JR3GzcIEzfs/HjVf2V41Jk\nAfg+SEbFrz1iPQJgdY9QjxXHdn5QX9s+IMXprlEKRI1VKgeXEb8sFh8DcBmAQwC+hLEsFlMAPMc5\nry2EgH6Qi0VpElQ8Q9x4ntg2ikSaG7dYcnF1sP1crls8YfJ4nunqwsGFDIdub8RTq77nfdBn7kPH\n8D5c9rOPuwpYeJ1PzRF60+G03JoXeZGjtF0cU0/ublPFUC0OAbizFwgZg9oEAN9fVS+33/9izzEp\n088+dwW+f8MHntvmYxUCGR5VbtxemTNMW9oiZ7HI4CD69JqBAde5VNcJsTNwZO0RAO7t2KqpVa6t\nftXt4fruQZnPWnc1EZ9L87S0httmOFDdS2SVOJHxJESuX/19V07Tk4PoXciwo3ESwDOFazJZH7wy\nj6guIn6FhfyynYR1UYiaLUFte+dDF3PeDvbLGJPPaqT6eAD8i56Yvp+PfMhhXddyWa9LLZtGIarP\nFopybEscLhb/K4BXAdQCeC3z+6sA9gOgLBaEJ6olRWwLTt6+Rz7Bi/8dbGBGi4usaqU9jds+oQd9\nTlgTKmbNyskyrEbt60/Q6lP19M2b8JWvfxE7bk7i5G/eQMNjDThwz3Icra/Hw6vmye+dbt+Cgb/q\nxfqz6/DbDaux5v5L8WLzXHkckwVDzRH6wkKGNON470rgiW2jwdYuxoBUKutzUjlmTAaRqe4H4neT\nVUa4MFySuASAkztXuBSYMnZ4bTeL9v82U6JYyPT5l0YAeG+bi3PsaEy68i0LS/b63pSrD9UtfxGM\nlVaOJ9wdTLsQIjOMUI7VnRbV2ie2+sWraPOpc6ewtJ9Ly4M+XsXn1vVelIF9qsuE/tCp5tEVgWtq\nXmaZ8SQzN1bsXyGVyRRP4drd12YFzYn+Ea/i+6PDw6j7zm2YPerItv75NDjg5CfWMo+oqGP35MhJ\nT6uT6n4ismaIMSLHzPFuJyVbJnhQ7XfVRUG9/lHzu8exHazvOPQM9sg5tOEnHwfgHo9+lrmOvg7M\n3z0f1+6+1rMPvdyNgvIw68QVrK0T1nUtl2vg99242hfGIhq3+0cxGU9t0bFN81ZFad4IL3RFSV0o\nxORZ2s+z/Em9XDe8tvFtFsjT7Vukcuf1OddN3mdBC9qWU7e79QV28vY9ePzBD3HZ9ifQ8FgDLpnW\nB8Y4kh/7N6S4c5MUbhfdx7vR8FiDdA8Y2bvX16VAoN4kAeDePcfwo11/gk+dgTEzgYpfQQzZzy0t\nvv1jKj7zlaffR+e2i/jqsx8iyZLY0ZjE6o0VmL55U5aCf7p9C57YNmrc+jf5O9psx+vFVkTmA5nK\nTWurqmgL6+y7VzmuHb/51FTsvqUC+3etwetf/bxLgVBlFwyNDFltNao3lIMNTBZdGbreKVoytOwO\nDNTW4Qe7P4YkS6LxNUf2WT8ZlO4QQLZvslpcggEA52gaZXj1yP9AKp3Cq//2qmsO6cF1HNmp7YZG\nhtyv11c7WTKur8bQD06g+s1J0o+4pbYV/WvfQNuSNtfnVPS+EQ8lurKnup2o1nNX/53/DTo+MQUN\nv+r1TeXmdf2FO9D83fON10w9lp/yI5XL7luzFHYV0ScvX1cp2+P34OenuOnp70xIN5h+j/5T5pCN\nX3TcCpBcn9NpqweXXLbk/b5bDAUvrhR3pcB4aouOXyW9GzO/vscYu0v/iVsQxtgsxtgOxtg/xH1s\nIn5U5dHPEiAmz7TmFnljFkE+1T8ZNOYcViu9qQuazQIpZcgof7qC29HXgd4F3JUOyqt90ord2OdP\n+QAAIABJREFUZU6hJNomAmXUBVamQDvMkeIp3N2bQue2i/iLQ5fL6njpBMMLCxPg4NKiLtKZiRu2\n8CE1WYQnb9+Dzm0Xsa73InoGe7Bi/wp0He+SVk+hnHgp+l4PHKKfdzYmpcJ3un0LBmrrMFBXL4+j\nf79nsGcs8KqfW6WOS6Q5Gl9ngYurkOn2R/bLxdjP6qX6ud91dyc+PnsOkGCudGWif8XDmnBd+PQZ\nx59s6nvn0PnQRazvTWUpEKrsD740S14H1Wotftd9+dUbyoUNq9By/yS03D8J9183jIbHGvDhiRMA\ngCnvnUP/V/tlUjUOyIwF+phwxkkKjz/4IRY89jLeePs9tJz/AD2TPsLMo5NksKo6h4Q8wuqvIiyZ\nwnq9rvciAHc+9dGMnIBjYVfToel51wUmH3yTstdc8/+39+5xVlRnvvdvdXMzYIMQSBhRabl7IU1E\nwZNMoonYGnUwmpYGPSFIkneOZ5IxYd4TTY/xkuwxySe8b2LOOPl4FCRRoO28qFGiLRONmThpIpEW\nkDs0CIYZGbS7gaBg93r/qFrVq9Zeq2rVZe9du/fz5dOfpveuWvWsVauqnnrWc5nnKcmNUxp9D2Bv\n/M69AS01w5wAwgDrX9NsffVM+QUhifXQO9cnDhjTFwJ9Y/LAZ054/ZFf/NQg3KDrR01/p0O0Ledd\nB/y5yLF2SeCLqre9RgGKkr5Tx9i7vwNUVwOc48jqVan749rKU0gFL81aAkTxCXKx+JT7+zrNz7U2\njTPGljHG3maMbVE+v4oxtoMxtpsxdgcAcM73cs4XR+5BhVOqXIKyUmxj2VWVW58yLeUcjptSTSDL\nIpdxFrTsbMEj9VWeVTOwfwLXT9801robrHghEBbPOe2OL+usP3YBAF7/4iU4f+tW3L7yDc8v8JH6\nKjR+qxqP1Fd5D2xdUQGBrITLS/nL6qt9yomp+EnYC4es5HljyDk6V63CoXvvy9u/YXKDl6Vh5LxG\n37joHhRRlkzVsT90732Yu+gxLHz+ZGChBDFPT+zehTX/4wTGXtTpU2DEy9p379iOhc+fREdXh6dM\nCRcUseohK6Wy7BNe3oNqDlzZDp/VWvzfVOJYHF+mh/fgTdeCfeDDzmfrZ49ADwN+O3MwGBg4OK58\nzc1fLJ3TK9yXk/Ev7URu6ixPeVw3A1o3kSfnPonNCzfjlfmveBZ3we7O3aj7eR2udOfYla5rodyv\nQZMmAQAGTZqUV2jGpNzJLhCbF272MkUA+cpeFavSZnPwLM7/tQ61Z0zSHke16o6d2ZU311WXmLh4\nfT3tLG3BoLztpGPJL35qEG6Q4tY0uwmbF27GJtdiryPw+t6w3FPm815U1y4JtIQL1JfGOK4KhciE\nYZIvCkmfq2L+qe5B5UY55CouJIGlpjnnP2GMfZJz/vtYjTP2KQDHAPycc36B+1k1gJ0A5gA4COBV\nAPM551vd73/JOf+CTfsUpFfcXII2eX1tEPlXr3yt11dWFki3P2qgFxCeG1nkTv3GS8NwaVun46rR\n6JSDTSKbOO6/zmB4pL4ai1t7Uf9ar9f+leP9N1C5hKxp3P+yYYNnxRsxfz5uq9uITzXvwJWvOde0\nGNdt087zShpP27bVO0ZYcIr8/YceeBxXuiWPxevMiPnzvcpvQQEuTpCiPq+w+D5sLql9kHOwzr9j\ngDFAxHfs+gFo33fQUWDcFHHP3jYX577oZMxQ87mqcpmCTz/28z8GzikRmGY6p/Xjf+krvjFxxERf\nAJ4cBCMC9lbf/4FXUnrN8puxesfqvLzCE0dM9F6abANo1CDLxa29mNPO8UIdsLx+gK+ssRxsuKdz\njxcEFjaXbCx1QYE/ciCeHLioDXrjHO37DgCsGrn6JZkITgs6XtTy0R5rlzhKrzS3bbc/tGG4//pz\nS22DVTsFkALklccnLO+0bTtpkaTdpM8h3RwtRxeE/pqz2TZIL0hBbuec1zHGXuOcfzyBIOMBPCsp\nyJcCuIdzXu/+fScAcM7vd/8OVJAZY18F8FUAOPvssy/av3+/adOKIKnFNQppXSxBxS+i9CdsW10S\n/bA+yAqCyEhge7wg5HarWXVfBg8AqK7Gt/9pkvdwEX7HcjYAEfkuZ/4A4GWa6K1iWFfHvIwIot1p\nb2wxyi0rIQ2TGwKzC4ixvLW1B/Wuoozqasz/1gDrxP5CAQX8UfQ2EeYr/+l9r19rHr3FLXDB0XGZ\nY700ZXU4dO99eKd5NdbVMRz/+oK8h9SWadNQzR2L7eBJkwILn+gKl1Szat9cVueMbqzbv9ju9bm3\niqHxW/7c2Gob8oNenKNbW3tw1etVWgVbZuKIifh0807MaecYOfssjD3nVce6/P4hbWYZMY7yeApL\nHOC/Vm2L2Oj6H+XlTGdBlhVJWb68zBWDx6Jp+3pg5iLUHflN4FzNteXQsqPZK0Ijb69TwguBqfBH\nqLJnqdRatRVV2XYpdDaDvsI2HCPnNYZeo0meh7ZFk+SCOrrrqZSKcRpjUUz9opikkcViG2NsF4Ap\njLFN0s9mxpi+ZI8dZwI4IP19EMCZjLFRjLGfAZghlGYdnPOHOOczOeczR48enUCM/kExcwnaBsmZ\nlmTEd3LxC9OSo01/TL7PYnlrwogJeUn0TX0Q+wh/TF0AYZKxlnPeej7ZgBdMKJa7n5z7pG8pVriJ\niMh3OfPHiJtuAocT4PVCnSiN7Ch7cpCiSW75OHKxCl1lNOFnvay+Gq2i/PBNN1kFuIgxH9U4XxtF\nHzSvRLDjMaeCNgZNmuS5yayrY5jw8h4813GDL5hRducI83MWOYr3fmYyTu7d680n0UZeiW8l+DQo\nh668/9dfPM3LRCH3eV1dn3uRcD1Q9wX6cuMKV4i/fP1mLyuF8A/W0dHV4bnoHPnDm8idUYOWEwfy\nln3lzBrtX2z3ylgD8PIAy/Ln2nJermHddayijlfQ8neYciHcC4RyrGYiAZwXiFU/+AC3br8Uufol\nqDvyG9QOrw2Us2Vni6+MuVfYRsoeospr8jGN63sqXnTUF57QzBMzF2ndO3RyhLoeXLPUUbIjKMdA\n4YPdWnaaS6ML4hZ1UlHvmeo9RT6GbjyzELiWxliUQ67iQhJYapox9lEArQD+Rv2Oc25lutVYkL8A\n4CrO+Zfdv/87gFmc87+LKjy5WGSLXFsOcxc9hmruWDRvvmOQ7yGXxAK957q/wclduzBo0iRMeOZX\nAMxvt3EsGWlbP8KsDLYWhrC8vi3bV6GHMYBz3PpCb17uZFuE5Qqwz08bOcesm+O5uWYYOPJL6ep4\n47zzUNXrSCas1sK9Qbamyy4P4oGls8B//cXTcMn6Tq07hJzj98TuXb78xqrl83tt53hlmWX3GNna\nKFtX5XzQ8tzPy/e7o9lxo2B+xVk3J1U3k5pBNTh+6riXd1lYWIULyLoZDMuvrMK8087C0Kff9Fni\n1Fy58hjKFm8hf9714p7blprT0WDhHmCTyzjsWqxbMR094E5qOVdZ1q1K2KxyeDJJFuRt/6s1MJ97\nkKxx7yemcdF9bnMM3VI/EC0Pchx5095H7FcsC7KKPNbqPSUoJ3mahLkHqhTb+psFq7ktqZSa5pz/\nh5vWbb/6k0C2twCcJf09zv2MKDFJI25bdrb0lcmty099prMU5tpy+PGC8/HGeecFBgIIX9uTu3Z5\nQSS6wBsgXknQtK0f6tu7alGwtUj15Q6+2dvfl9LrtLNQzTkmDhiG5fUDMP+OAb48w7YIC7FcJvjZ\n2+b6qshF6YeuL6JKIAe8nKwq8rnJteW8zB6D3YAw9PZ6uYw7LpvkZesQwV7CoijOpbCAiCCyS9Z3\nasuPiwwgvKcHJ3ft8oIfVeuoyM8rl2VunNKIajAwzn2p1uRUhnI+aBPCgskB55yGWGblACcAOH7q\nuFd1UQ7gkvM1c8bQNO+5vMDPptlNviAxccwe3oMLV1zoZESRrt2868U9tz3gvpznJuTzpd5v5LZF\nKrbPP/35vDnyvWdOoJpzR6E1pAW0XeUQMrUv3ISmr3UA1yz1tWGSV1ivhw4c6pvrce8nJsujKX2d\nreVeTkWZpnUzTiBc3OC5ptlNuP3mq3D+vP/A2Jld2m2iWjxtg9DksVbvKWLVpdBKoShZr967TBTb\n+pskKDKrhJaaLgCvApjEGKtljA0C0AjgVyWQg1BIOsEbJjtlcvddPhlz2rkvbZDpbdZm2QyAL2Je\njsAG8pPmhz0AdEtPaS+JjbjpJjdlW29entWgtErqOQh6MNYOr0XL+4dQe8YkdPS+h3lT5nmZAXTY\npEWT99XdkFWlWX1IqwUrRLGHZ2+bi23NH8X3nn3fKyahk+/I6r4CMXLGkQnP/MpLCSXOm5pKTHYT\neK7jBtyweJUvDVo1q8bLFw1GDwNevmiwz41hd+duzNnYV7QDAEY1zseTc5/MawtwlVJJeWrfdxDz\nuo85CptrQRGuCre29uCKjRzr6lhegRV5jKpYFarB0Hj0ONpH1+PJuU/mla2Wx+rq2jW+ktQNkxu0\n51h2QRD76l5WhUvBl54/mXduvj2rw/ewFanC8OrDqFsxHbmps9Bw9DgWt/bgytd69Uu7SnYE0/1G\nnotykRXRN+F2VPuG02+MyS/quqzesRwvq6+OfW3LbZjkFZlJuk92e+dR7UPaiHEAzGWp1W3kVJRp\nph6L8yKQyBih3PuTYuuGEPcFJU3kkvVZpNjjUQwCXSwSN87YKjhlqj8M4D8B3M05f4Qx9jkAPwZQ\nDWAZ5zzWlUouFsmRFVebjAQ26FwpTO4VNstmeShBJOoyYljAX9KlJ9v9TXJ5YwF4SooYe9NyXdAS\nqyB0STrisq9uSU8EtelKL6vHkH2bTS4GchDdsvpqLGr9AHM2coxqdDJkyIGDD7bPsHZbkcsry8cT\nLivrZgDLXMVSIAIQAaD7zOH4Te46vysHY+iBkyVief0AfxBV89VoOXEAtQOGoaP3PUemI++g7nAr\nHvtBj+d2dP7Wrb6MHCIDhUDndqI7b6ZgUtP4y5jOv3Bn6WHAgjsGYsKICT5/dE82cf3xHtSNPws9\nbnlpnysG8ktny4FkufolPpcO0/1Gzowg/MwXt/ai/nUnX/gj9VWBmS6SuDiobibGYKwdzejhvV66\nSvHyWqhl9yiuFUFjI8Yx6RJ8Gkvr1vfkmAGEiY9bpHaI4pCKi0VSOOfzOedjOecDOefjOOePuJ//\nmnM+mXM+IY5yzBi7jjH2UFeXfpmlHMlCPuOoVg+TJUJnnTIFYjXNbsLtK9/A+Vu32t9Yrlnq3CA3\nLAfWLvHeWNU8qqJv7zQ7paDbPjsb26ZOw182bEhUQU9t12SJaZjc4AXnyZYbsVwt59lVg6Tk3MGf\nWPUJb2lftl7Jy/hCIVr2pUuxZdo0LPvSpXl9sXnDz7XlMH3FdFy44kK8/sVLcMG2bT5/N2HFWDeD\nedbXZ2+b6xUREQGY8jkR+4nzrwuiu2Kjk75s+ZVVePqe8V7lPcHuzt15S4bqsrc8f03zTeSOvnIj\nB3ONAwwMjVMa8Y2VW9HLHF/nYX/u8sZcBEaC9+WdlgM/AaDl/UPoYQy7e473lR8/3Irak6fwrzPg\n5YYG4OXVBud515pu9aZhspODeeU/vY9tU6fh0L33+YI+5fGUz7HaVpDLRq4thxfq4J1bDo4n5z7p\n5ej2ySaseGBoOHrctyJgKvIDwBdI1rKzxVf+2IQuePX41xd45dWDCmmYvgu6vn3uS9JKzYUrLkTz\njuY8JVCsHjS6qwfClUGtSJgmQq7vtZ1j7EdQ/8V3cnXTJKSxtB4WcO3dY2MGEJowuSFEfR5HCYjr\nj8VD+mu+5KA0b88AhtxBADjneYF7xaY/WZBLlW/Q5s1XDeCJEiTikfKbv01aI9E3YWlqvv8Dbwl9\n2vZ866dAPReqhURt19ZKpY6XHBT23p7deanIVAsxoEljJT2wD917H95dtcrL6XvLnYN9gWymFGom\nq7SsaMm5bkWxCoFnHQZ8xzFZquRx+F7bOV4A2bL6aidv7f4/e+dUWBBFAJqwzIq5pE3FJaX2Uufa\noXvvQ2fzaoyYcBxXf2G4E+QlnT9hNZcD9MSY+66VmV2++SxSjwnEGFVzjvbR9Z4ch+69z6nMyLmn\nRJquL+2cdMfYFJQoz8MoeXXVuSYU45adLaitGoKOD46h4bSz0DTvOe9aDkoXFxbEqdvOl7pvVocX\nLJeWMiSwSSsYdE3kXevSvS03amTgSlDQcdLsRxixVu4C2rINLDSRZsB1GkQd2ygWZF2fTPuXS+Cb\n9XilrQfEJA0L8o8ALA34IVIkSlWxNLFx5BfWHvF/gWqlCHyLTNl37ND+i7GteSwO7b/Y/7kkg+ib\nsDR1nzkcQJ8/swn1XIhUYx964HFfuq8gC5YOdbzG3v0drHlkPk7s3qVNRSa2l1PPBQXHdT7xhJfm\nbd0MlpcWTkXXhmr1Ftvs7tzt5UFe+f1TuLXVURgYGPZdPtlLNyfKW8vyz2nvdS3uq1D38zpfEJ0c\nQMYA1J46henjz/T8toUF8fip41j4/EnMXfQYfvz4r5E7o8acisst+XvoZ7/Mm49j7/4Opm3dirHP\n7EeD65dZO7zWs+hc++DTuOXOwXnKMeD4pN70vxjqx/8Sub1rAN6D3N41nq+nXFLa8/mcOt/3IOh8\n4gmAO24M4hhNs5uwaeEmrf+4sDb94eLT0eue2z2fnuB7aMqWTtkyJdKhBVVbk8+VaX50fHAMDd3H\n0HLigNO2a8Vref+QMV2c+Ez2sZYtZrpgTxGLMP6lnb50a+pYRLW8qfsF3WtNvqZAcFlndX/hP24a\n9yhW16il4tX+6vaX/fuTugSYVh7z+hhQoc8mFWUxifo8jhIQp+uTyQLdsrPFu+9l2TprPV4p6wGF\nxqggc85fDvopppCVQJbzDcpKk/wQVm+MgctMhjydcelc/xbAmfNb/jwgAO+cT30OqK7Gh2YGvziq\n50Iu6azmujRF4gP5Dybdg0TO/KFm9xBK0CvzX8H968/Fd+/YbgyOA/puUmfMn4+/fP1m33K0sFSG\nRdqrypr6YJqzsW8sACdf67UPPo0F3x6Mxjv7Aubk/q6rq0IPA15wM5vIbiRChsYpjdi0cDM+/dIA\nrPz+B7i1NV9xF8e+YiPQUjMMmLnIK/M84eU9ff057SzkRo3Ekd0fClz2FPLJOZTlcRFWVDnoUPS5\nuWYY6safheaaYZ4COPMjMz1/4NU7VqN2eG2e0hC1tK5QMrpOdoEz4IWPV/nKOovsGyK1mxwoZoPI\nErHhPzfkuVN449l9zCtbrXuZ6uE9xjklz5/VO1Z7gbQqDZP7ypTvu3yy47Zx9LjvfhF3KV/dT7yY\nXl27xkrZDn3RiPHQj6L4iaBEr9y7i+mZofZXd08shuKZd4wI42QTiGhD3KX/KM/jqMfQPQdMCqZ8\n30vqClNIrMcrZT2g0IQG6THGJgG4H8B5AIaIzznn5xZWtHD6k4tF1hEKgih1q1vySjtQwbaaEQDt\n/+XAtz2de7Dy+6e0gWI2cpiqsQUtAdosO5mW0NR2w4LjgsZNPXdAeA7ivPy87v9vbXUe1C/UwbX6\nOhW/gpYCVfcDtZSyzBvnTUNVr77s87urVoEDeOHjDL+bNwUdXR3GvMZ1P6/DwudP4srXOKqgCRiz\nOAfyeZDdSlQXE5nGKY2+vm5euNl4PuRleHmMdcvUj3//fVT1Or7MTy272VseF+dBRq3CBkDrBqAG\n8W1euFkv2/7fY/WJA2CMYZ7i4iPmlChJva6O4fjcs30uLur5j7tcbrvcrLqs6MZWd93GXs4u8LJx\nXy50Kb95wDFNLmElDyCLME5puVaIwFMRIFsI4rq62M63zJy/fkTiUtNSQ78HcDeA/xfAdQAWAaji\nnJfsTDHGrgNw3cSJE7+yy82PSwST1Jcpbm15k6IlKz6mi9/2xmPaTvWrvLW1x8uQkMaNJsx3NMmN\nTR23oQ88jis2cm2SePXcbjvvPKCXA1UM07Zu1foyL27tQb1bqljn8xb0IqQ7ZpT+qEqZjFYZgL8o\nxi13Dgbg5OgV2SAAv6+weKGZ0869kt6ykmzrey98eIVCrMtoAMDrk3B1EJkXnpz7pG9+7vn0BC8z\niLAEV7Nqrz+m8T50733oXLXKqZDY2OhZBXurgMZvudk4pLHI8zuWfPbras/2jqXKKpAVFJ1s3pzi\nHNWsCo//4AMvA8b8b1WjsfsYmt7tRq5+iadwi4wYNj7RcZBTCwp846m8JKj+2Wn6u6bpO6q9h3Z1\nOb74ATEY5Uxa42e6n6RJ3Pu8GqMhCo9k3d+4P5CmgvwnzvlFjLHNnPML5c9SkjU2ZEG2J+nNP+4N\nSz6uULoYGFr3fcF7yAPwKbdy8NrJvXtDbzxhAQ5ygJlN9TYAOHRrPd799/3gYOj4zGRc9JGL8qzI\nsuJp224k3Ad63fhxvmAy1conK7PPddyAzlUrAQAjJp3A2Gf2e4FuA6sG4lTvKQB9gXU9DLjr+1M9\na6GqYARV1YuLqNqns3Sa5plQel+oA5bXD8CEERPQ8e4u1J48hd2DBgK6VGPukmXnqlVOIxZpB2XU\nlGlB8z/o+pDn55HVq7xxf3r5LT7FxzpQT/TLbfO2uo3Y/e4uRzmWlOQg5TAskE6kL2voPgqMmaoP\nxtu+Cg2uIrxnwyfx/q5dePPDwP/9lQFecKIaRCm/SKQddCVb+wHkj6US2Ktev0B6Febi3G9t7rFe\nu2Bo33ew5MFOWcFkiJH/nxWlU9wPhEIsAqlFRcxiByRWImkqyP8O4JMAfgngRThV777POZ+ShqBJ\nIAXZHttlN1vLmu1Nx2Q5fOKHToohMAZU+S2ZaWb00JWotpF/27SpAHcetD0MqK6q9lkwhaKaxgPf\nqKC4CggYQ8vw4XnWShmhgHj5fxnHmntqvUAqlfvXn5uXrUFXJtnGhSIqQmEH/O4WQPjDTLhOzNnI\nMerSszH2nFfzsimoc1iXfeLZ3RMxfsPRwLKt2gwaEfquGzNdbmnV8mlSyIMi/cWYvDyjGifhlBy3\nLUnrU+ZGfbZvGVz4jLp5i/POjaR0i7LMaol5XV5hgfZFIIG7QugcVdpOev2GuRV530lZV4L6ZKNU\n21yHOjeZNJTDLGdUCFvxKBRxrMfq801VmLM4vv2NNBXkiwFsAzACwHcBDAfwQ855WxqCJoEU5PiY\nlNColrUoNyC5uMEL+xrybixRLcc2bJs6zfu/SO1mI79sQf7tzME41XsKc9p7sa6uyueHrFN4VMSD\n5Xtt52DCy3usUhl5nylpwuQx1Pnyyjfsq2vXeBa1KlaFoQOHovtkt2cVFnKpn6vH2bxws0/G5zpu\nyDt3UR6eOpcPQEqNpqTBU48zd9FjkX3JvfN0/ge49oK38wpc2Mhpk5ZJt78u/ZqwFk8cMdHneiAs\n6wACXVxk5DHpcfMYRhkf37lrXdpnZRVKspJOz+gCYvliHei+Y5G+MVJ/IrzAR1VKrO+Dln2yTZFn\nK5fAZgUkSrtZtHDqXJ4K5cojE8eYQz7FpSe1QiGc81c558c45wc554s45zdkQTkmkiGiZvd8eoJ1\nCiRB3AjoximNXnYAXdSrcLk4uXev77skidW7zhwO7v7WyW9qe+yyViz49hA03jkAP5vT46ZEGozb\nV77hu+HKWRRMCMvZ+Jd2arMqiFLAPbwHuZ/WesVP5DRhItuAgIFpU0nJ4yramDdlHtq/2I5X5r/i\nFV0AHJ/sVT/4AA3PdALoK5sL+M8V4JQhXn3/B3j8n97Hu6tW5fUjSoYBufSxSGEHwPPzFd/L0fdy\nVPuoxvmRUyJ6ZbPfGAC4mS7kOSCKo8jzQIyfrsBGp5tZoFPKLGDKEqKmX5PTJsqV6iaMmIBqVg3u\n/tNdY7r52jS7yRuTfZdP9rJBqEVZTMhR9bmps1A3/izkps7yFWUIu+ajRP2LfmrbSiHK3ZRaTocp\nRRngH2vdOFrfB0dP9v+2lMWp0NeDlh3Nwe0r6OZt3AwgunaLnXLNBjmjEADvOkrSXxvipGfNcsYq\nwo+NBXkmgCYA5wDwarNyzqcXVrRwyIJsh/d2rSmgEMXPMhUZAto2lb0WN3fhd3do/8XoXP9W7KTs\nsjwma1auLYfV21eDgWNCbxU6BlSF+peaZBFL9d94aRhm/bHL82OWM0IATuq0f50BHL/0PbQMH+4F\nNqlZE2yDJOWl1k8378yzgMuBb/PvcC5tm0IVvQCq3IdCkAU5LEOHCdWCvOfTE3DnrL2+72yQj/+x\nn//R6NogZ/cIvA4enA28vQ0YMw2HnjuEd3afhnUzGI5//RYA0C7XqxkcAOel4OjJo3kWZJtr0BQs\naxpr4e8N6M+t48qz2inKMeCjqBv2XqCVMMhHPKxkfRrln20DQtNwfbJZts+TTeciIizIQJ/yb+E+\nkvtpLVpOH4qGo8fR9LWO0O2N7YQEExebQrlqFOtZRpQ/aZaafhzAcgA3wsliIX6IMkGUzF3tFlCQ\n81DKVoFm12LRHNFiYYONVUd+s5YtHsJi0ct7kTujBp1/2J9nWTS1abJ6iPZNlrqWnS0Acy6QJ9/c\n73N9kI9lYw0QltkHPnMCN98xCI/UV6FlZ4tXZnnORu7L8Sty6wrFSU0pppai1vVdKAmiHacIA3DF\nxl7PqiLn5PX1G/riCqIYSMdnJuf12ZTjWWe1CrNACavu1bVrsOaR+b7cyqZ9hbzC0j59xXSs3rHa\nS7J/0UcuwgXbtuGij1zk5SyV54Cgdnit2dr29jbkRp6Bug8dw7JbzsPN3xqAR+oHoGVni29bWUad\nBav7ZLeXa1qUc5aV3aDzK+azyD8s2jfJLPdN/c6bI8zNK/32tsDVFXlOqW3J1v6g8x5ogTQUkvDu\nX5b5nZtmN/nG1ITNvWNxay8e//7JvhLqg8f6ZMzrqy7Xr7CKg1nlAfZWS8ZMRfv+PzvVBBMgViyq\nWXBZ72KRhjVbhzy/glYGCMIWGwX5MOf8V5zzDs75fvFTcMkCYIxdxxh7qKurq5RiFI29/okoAAAg\nAElEQVQkLgYAPJ9Gxlje8qV8IxEPU1OO1ySEKQ2m7Xt5LwB3yYwxtNQMw4hLz/GWtcJutrobZa4t\n5/nmNk5p1CokDZMbUM3hBMmNcfyYoz6o5b6IG7f8/0HnngsO4OAo4F/dYiHrZjBwt7+iMhsDw8QR\nE73zKC/Vm/oufzZxxET86wygl3Gsm+G+aLQ5FQGfXn4LVlw1yFuOrR1eix8vOB9zFz2Ghc+f9F5o\nxt79HZy3fRv+8QdTceesvfj805/3jadufqoKkUiof2trj6fAiOOK341TGnFra493fGEBB4CBVQOx\nesdq37Hl/qovFQzMe/E4snoVcm05nyIn5BOyAI7Lg1zpz8eYaX0FM94/hIap8/PO68QRE9Gys69w\njHwM4U4iu5UAESqRSdtOGOFULBw6cGhedUIZufiH7Noi2ncHy5vnsizyS23dz+t88149jrzULI5T\nWzXEp0yGKi2GQhLe/QtMt5cWGwXJ5t5R/zpDVS/HhJf3OO1tX++TMU/p17mICFeVixdbuY94cr1/\nyHNxkQl8HmheMrLmGmFyPUpKHKU46bOV6N/YuFh8FsB8AL8B8L74nHO+prCihVMpLhbFStFm2i7O\nkljYkjtgl7EgbMksjmza8Vy7BHWHWwMDt2Qf4DSWKkUSew6Auf7guihmeYkU0Ofi1S2di32+8dIw\np5jG+R/gH687zZcyzkNKKfeYW1BFTjuky0IgchgHzU9ZRi/DRkhAi1osxDuulMJMzZ8su5MIt5R5\nU+bh1tYeHFm9ysvWIXJhq0GVtsvyYfMtaM5GvY6DjqVzU5HnpC57hCoTEHwNqi4oItgz7FrzBZju\nO5AXnKbLsgA4Pre1J99Hx6DBaJCCq5IuySdJAZYXhJu0KIjF/lHmWN71lpJbRjHIQtBfFmQgik+a\nWSweAzAVwBtwVlgBgHPOb00sZUIqRUEudXodq6wPii9u2D5JUxolGZNcWy7/Ydy6FLkzatBSMwwN\nU+dr25TTk2mrlYUdU5FXJLFn3F3KkRTHyOmxNMcS+8g5j+Xcu742XD/J3KiRGPqHD2FOO0fHZZM8\nhV2VYeKIiZj5kZmh6aTCsl/oxudDDzyOORs51s1g+MbKrV7Ku6E9veiursLEMyZ5x7byxXaVhl4A\nXJcj2EWcX1PuZ9sUW+oKg+i7qRpjHNQXAsA/J02KunweTfMoLJuC9cuvEvOg8/kW6KoTpuUvGzsF\nmPTSHPZyYH0/ipnRwva7uhXTnRdgMLRHuD+VgkLd39OSgfDTn8YqTQV5RxZyHuuoFAW5FES19qrp\nbsIupqTfJ37zv3ck6s75qz6LsZz/NcAypKY/06GTPai0rS79m2whnDhiorbQiUmhM1UQDEpFd+jW\nenT+YT9GXDAIYy/IL0Bg2yebsQhrS2elVi1vuiItxjnj7pubOgur3zuYl/7JJmBTZ8E3oY5/45TG\nfOu53B8gklVStcDqxss3FlIeXlFBLqgfYSsCuv1t5oIaYCgr97rP0rLqRb2Xedw70ntp7gl4sZL7\nZpRZnO/Rk4HDO1PJiRzW15IqMQkt7WTZjU6hz31/OidpBun9O2PsvBRkIsoI2T9P+Hbd2trjBTip\nqOluwvzB5LQ8Ov+vMP9AG7864fMq5PX9PXMRGo4eRzWY04aU0kpG9VFT05/pfNh0suvkFWN07YNP\n56W1k5Wejq4OL6hLHk+hSOzu3I3pK6b7/F5lltVX45Y7BwcWjehc/xbAGTq3nHQsXK8+nBcspWJz\nDkQfAWD6igtx4aMXoFny5VZ9dRsmN+SNMQDf+cm15TzfdKDPl9Y4Z9x9m+Y9h80LN2PTwk3aYMKg\ngE01kC8suAuA5+PeNLspPx2U62976Ge/xLZ/WItDrw4NDd5S5d3duduncDZOafR8UJuOvNN3/Um+\nvSKATQT5yantxFw2+mADeXNa7NfrusD08B7fuIigyU+s+oR3zoTfvxyg6H024Fw0dh/ruy41yNec\nyYdU/ly+F0XyU525CE3vdqN9dL0n5/faztHeAz2/a8XP20Ocg8M7tfcZlbh+w3L/1PtfUTH4k9uS\nNb/pcqBQwY+CSjwnNhbkbQAmAOiA44PM4LhYUJq3MiJOqiR1+yQV7nTp0Gx9V+O+DavyBslvOl4U\nVxHhQhBUESksLZy6RC+WnnXLzbLLh0C1pvrcHwIqeh269z50Nq/GiHOPYuzMbuRGnuFzN4lSVjvI\nQgzAXA7ZkqipzsJkBMxWRfV8NE5p9Pn0xi4j7FrYtjV/FOjlAOOY9qNrYlmQfXLLS/ii0IfGaukb\nQ01Z6CgrAl5b7nmV95dXXAAL/26pfLXJBcHGZcLW2hX1PmO6h1hbkItYGjrNiqSRKUF/K53MrB6U\nAWlakK8CMAnAlXDSu10LSvNWdkR9u9RZWmQrmJxWyyYKWFf0IchSZbL0RLGKqMVQ9nx6gjGpuxyx\nr0urZXprlr8XfZzw8h7HZePVh1EnWXZzbTkcWZ1fZENGVcaqWJUnn2hDjLewwKn7y1YzXzGRAKvO\n2Lu/g2nz/gNjZ3Y7xxPZGravAtYu8VnGhNVfd/6FQtnDezD0gZXeuWqY3ODkIOAcjaed5bNi2kSQ\n59r6innoMpAA0pw58o6zPN58dd7cFPNn6AMr81ZIdA8VMe7VrNqz4EexpIhxEikUvevPtWqPmNfo\nzMnGBdaKhHpuAcf3NPfTWkcZlpVjg9XSOR8M4NzJYLFheaQVAXmsGiY3oBoME0+dyrP8ikwsAJxj\nDR5rbLtlZ4sz52qGBWZ6kOU0yWx7jqLeF02FIdTj5Vm2DStUhSROEYtEyBk0StDfSodS26WP0YLM\nGKvhnHczxkbqvuecx6sDmgKMsesAXDdx4sSv7Nq1q1RilBW2AUa2b6BBZXh16EogmzIpBFkd4lhF\nbKyfqm+lKWtGEIfuvc8LxHr58g+we+BAn0Wt7ud1WPj8yUB/YLmwg/DJlK2EJl9mWbGObcWSvs+N\nGumz5tXVnu1aCYHG7qNoqTndK2ssH1OWf/UPelDVy73gQFMAnykAKihQMTCTh2tF1ZWT3jp1mrME\nBmDBtweHnt+0fOWFzFGtO3kBcs1Xo/nEAXDW577hyxqx/89OIZ0Af3KPtUscVxowJwVZgZQZYRmu\nPXkKHYMGGoNgtWOZkiUyzQw9NmTJX7NofU+hTHi5QNba8iYNC/JK9/efAGxwf/9J+rtkcM6f4Zx/\ndfjw4eEbEwDSyQsqkH1AdWV4dSyrr8b8bw3Asvpqn7XFZ2mx8FuLYxWRZfP1TbJ4CN9MNaVZFH+u\nsXd/xysEsnvQIM+NQBy/YXKDMw53DPAVv5CZN2WerxDE7s7dvpuwzlIllGMGFnwuFD9e1bqaGzUS\ndbVnIzdqpDNfRtc7S90zF3nHbew+iqZ33kVD91Ff/mJxTNlXd+S8Ri+3c/OOZmPRCw6uHWv5HHgW\nTxfhBywss778yG4u2obTzgLQl/fZkQ/ebxtrS9h1I2RU+6eukjROafTm1bO3zbVeBVHnYYurHIvv\nvHEEQ8PR48DMRU4hHeFPHmTFE9cZq/K2Mfn0JkHMpY5BAx0LseGa0o51Ql9Wgel61h3Txsc5jGL4\na9rKZup7qj6ra5cAIi4gQZnwcqHQ/r5ENjAqyJzza93ftZzzc93f4ufc4olIFIsoy5IioMm3hB+y\nj25J23ej0SXZdxEPg2X11ZHr2KvKr4fy8JXlivuAE/vJhT1kOcSSsymYR8gwb8o87zP5JiwHvanF\nG+ZNmec7F1ol2LTkD81NX1km7eW9WF1zOnIjz0DTgI/mu3AAnusE4LwUiUIkOiVY9KVxyDhUa5bf\nvXMweCyaWpdi04fnYPPCzdi8cLOXvUFWyD1/bCkoT/hwi+qQHZ+ZjB7m/I6CSRkxKflDH1iJx+5/\nH59u3pk318e/tBPo6fGKlwQhjwHuHYmGDwaCuat+Ym42zW5C+8JNTt7ba5ZixKXnoJdxtM6oCm5f\nXG+jJ3sviroHfypK8zVLfYVVbMlNnYW68Wchd+Y52ip7toS6QEjIYxBXEQp7sUpjTG1lS+qGYsWG\n5QC4M5+y6FZhqNIYl0oMWKtEbIL0fsM5/2zYZ6Wg0oL0srKsExgkZNhepMlSXRys2rIo4hFFbjVQ\nynr5NuL2uly0plRiVvJq2pZR0875lt5H1+elR1NdOABzEQV1P5s0Zzo3EDXFGgDf0uwh3J4fxKhZ\nupXHZsN/bsDuzt2oGVSD46eO+8ZLl5bP+jqSzndeAJuStis3dZZTXc9t841p01DFgV4GnL9tm++4\n32s7B+Nf2ol1MxhWXDXIbj6LMQD6XiID0vCZlvjFdkMHDkX3yW4nRWD7y9745uqXmFPwgaF9X4DL\nRkroXGtMRUfiYhsgDPivg7TuwWmMaamfB6ZUgplUkCvI/YMIJ7GLBWNsiOt//GHG2BmMsZHuz3gA\nZ6YnKmFLVpZ1hHVEWPHC5JEtzrJ188IVF6J5RzMaJjcEt7VhORq6jzlWRs0bu23gnjx+ubYc6o78\nBrn6JdobutbCE3G5V7YyCIutKZWYjiArlGeprh6Kas7ROGSc84VkKWmY3OCMmRuEJe8nlvxFqrBA\nyz78Fu8815iAfst9EZbWPHcLaeVADub00KwsyLKKYMXuk90+mUX7IqBPYH0dSec7r0/iu7e3AbwH\nTdvX+87VyIl/ARh3fksyN0xuwD/O3o+7vj8VK64aZG+BEmMApp2Dap9MFi6xXfdJJxhzd+dunyW5\nqXUp2kd9Nj8Ij1WjoftoKu4Ogaxdgpbtq3yuNdXCXUa3uhTTMhhkAQxKC5fWPdh6TAP6J8tWCNeY\nMHxjkfWgvIDVyYohZSt6JRDkg/x/wfE3ngrgNfT5Hz8N4H8XXjRCJQt+bfL3Qp6wvLC6HKHiAcPB\n/Q9CXd+kfKQ6ZVGrVGmQjxH2oNNmtZBusjYPJPkBJrsCyA+1C1dc6MthbEL1ixRj9uSenWjfdwBN\n29c7OWcPt+LzY0d7OW9lP2JVJrnv6jlUz0dHl1O2VrzkqL63cjYL04Pb6HMsPVzjRt7LfZEtf4Df\nL1r0Vf5tGmf5fOe9rIjvxkzLf/CuXYKxF3U6WUH+9gt5cvbwHuzp3KM/pkmeUSOdMbp4sfZBr54v\n08uV2K5mUA0AN8uEGP/DO7XKmtfWuTcUXslQXoa9Y897Tq+AmV5aQ5QB24h/9dykdQ+2HlPLl/JS\nGE9Ccz+XCO31lHUFvhik5M9fSdi4WHyNc/7TIskTiUpzsSgGpkpxQimTy8ROHDHRV8krSpWpwCV3\nl7CcwVG3kwmL7FYreqlZLaJEqcvZLeQyw6qbhOoi4Vtq3tGMHk153ltbP8AVblnmZfVOXlhwjs2u\nS0WQC4spS4SuMh2gX2pW3TxM+WiFvMK32rQ0HJSVImh5VN1PzRFtqvqW5LwaCZBXHbegXL6qrAXP\nr5uF3LVRM2usXYLc3jVoqTkdDUPO7CttLZSBkCV12ywkxXAvSZLFoxTuFlnK1CGTVblKThau74yQ\nZh7khxlj32SMrWGM/X+MsdsZY0NSkJFQKMUymYoIsGKc+6wv4qGuBkWJNF22uYIFTbObtFXNZGwt\nw2Pv/k5g4J5uXNVgN/Gd6GdHV4cxq0WuLef1O8iSJI77TvNqVPVy1L/O8pau1e1lfFkcuo86hRgk\nODiu2MhRzYE5G/u+m3jGJO8GKNqQXSnUMRBL/wJhNTcFVsr7ijEyZTOR/xarBb6xF/l73bmmtYRZ\nLI+q8qnKsSyHbiVBWMKDcnObyMt6MH6cY/HVyKuOW8PkBi9Dhy4ndKjFMk2rkMbKltY9yTp/+TVL\nnWwq48cht3dN3tdCHm/FYtRItAwfjh5wtJw40DcWlkvqYZbXYrqXaGWxtHzaWsTTJIo1vZjPNgqg\nM0BW9MjYWJCfAHAUwGPuRwsAjOCcl3z29TcLchbefEXOUrWalc6aaBukF0sGN5hpwst7IlmGdQSN\nq8jby8CwaeEmq5yhXuCQZbDa4tZe1L/OtP2QLZ2mgCoRBJPbuwara4b5LLFDH3gMV2yEZ0HevHBz\n7HMlp4yLkws6aOx0qwVq/l7c/U6sALrcqJF5+4hxnThiIp6c+2SozKEBiAbrizpmALTWaSD4PMS+\n9l251CDBtPDJNeqzxjEwnjO1YqAuf7kytnUrpqMH3LHaLtyklUfgG+fBY/ssyJZKgCd72L5FsL6V\nOuiukGTh2UYQAlsLso2CvJVzfl7YZ8WkvxYKycwNssRLMWncTIPcA2RklxLhMmKTlcPmPEVV0mKd\nd3mJ2VU+k4xfksj0KEUxfJlNuo8BY6Z6Ch72/x4tJw6g4bSzHN9THZIbgyhiImfriEquLYeWHc1o\n6D6KpnNv8CndAIwlkNUiNAC0L1IAAgvr2LgcBWF9zm2vbZ3i3brU77bgblM3fpyj0OqO7Z6nQ38a\ngc69p+tfdhWXFO+6OXkST/7VtVplPPWXc8pyUFAy82wjCKSrID8G4H9zztvcv2cB+J+c8y+mImkC\n+psFudKwsdba+sLGtcjlmq/G6hMHvKIegF95sVVyddjKEMeHOgjbcVXJ62tIejXVl1mkDhME9Vu1\n3Pt8lTkHZ8xReL9kqJYYUPUvtoKjU7plq7BGAQ8a06iWfHW+RFEqfMcKesGwVQSV7ZyVpdXOC8SA\njwK3tXnbCDeHIAtyoEKuWpA1KwsFh/wz7cnyWGVZNiIzpOmDfBGAf2eM7WOM7QPwBwAXM8Y2M8Y2\nBe9KZBnVL6zYPtAm/7+g1ErqPqY2VD+0XFsO01dMx4UrLvT7Im9f71gwpRdF2XdNuD/IPq22GH3h\nlAh7o691SCR+nj9mW86pCti6FA2Dx3r+0oIwf8u8vmr8OHVtqKnDgHy/XxXhyy5++3yVXeVYVMLT\nIvnTNc1uQsNpZ6GlZhhyU2eZ9wnD7W9u6iyfj7l3HqfOz3vo2vp+iu2CCuuo82X1jtVeJpWw9Ey+\n7CQnDjgV604cQK75am1KvdzUWcHXulJEpGVHM3oY0FIzzMl2IW3TdO4N5jGw8XtUtmmY7K8MWBTS\n9s/MaEqtVO7xWc6GkGXZiLLDRkG+CkAtgE+7P7XuZ9cCuK5wohGFRg08Ew9knQKV5MZq2tc2mELe\nTt3H1IZOyRbKmJy6LTd1FlpqhmHigGFefmD5QS8q34nfUTApT7m9a1B3zl95QUhyajPfWIXc7LUB\neO4+LScO5J3LsPHO66umNLUuiE20KwL1Gqc0hiqN6rGaZvdVO2yc0oj2L20xu1foxuL9Q45S+P4h\n630Eaio10UYVq9LmwrWWSfdCoihOubacF6jYdOQd33F8AbOaefD5pz+PC1dc2Fdi26XhtLO8/Ne6\neQAALe+9FZwW7JqljnLq5nlu6D6ar7S68yM3amSqL9ZNs/2VAXVoXw6zREYVNW9O7miOr8BnOadw\nlmUjyo5QFwsAYIx9DMBfu3/+G+f89YJKZQm5WOixXZoNTPOl+J8m8WuNs29eHzRLZ8Z+avxyZZ9X\nAF7QUdQKfUEuBrb+dSIIiQGoYtW+ACFfxTZDUJQsi1xqunFKo7e8XjtgGDp634vl86emuhPlmoHw\nKnpGTMvocnsBy6NhY5zEx1ENRFPnTpRjyb7EE0ZMcNwp5AAwJf1YkDuBL2D2nXeBi7/sGxdTSkZ5\nLPOC92zcIgRy9T7l2PI4RKkOqZMxznK4LlivWMcuSPtFcg3w5m5Xl3OPJ59rogJJzcWCMfb3AB4H\nMMb9eYwx9rXkIhKFwjZpvJrmy2dBVSwgNtbepJZiuY08a7bGIqPrZ64th7rDrWg+faiT+sn9rml2\nEzYt3ORLr5Xbuwa9AGCo0Kf2yWRl9+TYvsrKItMwZZ6rdDq+rXJ6Kt9YhSz7ylbXiSMmomVnC5rf\news9jKGj970+a6TFcq/cT9kyDfgLbUQpChBkDdfOiQCrW9icTpLmyifLhuWOJXffwTz/bXV8xMqL\nnKpOyMfB0dHV4ci0fb0x/ZjRncB1lWn/yzDHr1qjoE6sHgpwjhoMyHeVcqtENs17TlvgJNAtQvR5\n6izUjT8LuWn/TTsH5dSPNtd33v0hgZVVXbWInNIraxZeIc+rDxfULcO7Tqo/4nwwenLBjlVyMurm\nQpQPNkF6mwBcyjk/7v49FMAfOOfTiyBfIGRB1pNKxHAMi0Ya2SdEGwwMVawqsgXZsyxxjmpWZbQC\nAghMJ6XrEwC/lV1OKfbuLtT09OJ4dRUaps4PHXdf9oOuo1oLYRRUi5ovqNCicIVsBRQrCnKxFIF3\nbi2svb7UaSHWcACJLMi27YRisLyaiot4mSpcC3Cufkl+Noo41sRXH+77W8wLtR33vNaNPws9jHnX\nTGyLrkLQ9ezLQqK4JVm3l4LVNHZxjUJbbKNmxZDPuW6ftOWthKwdldBHIhZpBukxAHK5rB73MyKj\npJI0PkbQirWlOODNXlgoJ4yY4O+DRh5dPz1L+NT5jh+jq6DkflrrWPqar/aO3TBlniOvm1NYJ2dD\nV5dj4ZP8nz2FwO1HR+cegDF0V1ehhzGs3rE61C9S9HPiyZOOcsyqI+VuVdsXsgm/1T2de7Qlk9X9\ndVZAOaBMXbr3zq1r8crtXZMni2jTV0DmmqXI1S9xLJumcQmYc+q59vognU+PJNZBVwbPR3dHs9dv\nUcwDgH/lRbIAawvgRL2WVLndv1XfdXFeG047y1fC29aiG2ZhC7qehU+/KD2ubffB2b72RXteWXNR\nPjuBwqddWbA5/wUomuC7LqP6wl6z1HkRMu2TtsW7Enx1K6GPREGxsSB/E8BCACLH1fUAHuWc/7jA\nsoVS8Rbkck1pE/BmL6y6AKwtU6HWxXtHou6cv3J8jTlH+74DdlaFMAuE7M95+lDUnjyFjkED0cuq\nQi141iVsNedYzb0r9zmoxLJ6bACef3HQWH9i1SfQfbIbNYNq8Mr8V3xy6XLgJk29Z+q3/J3nO647\nn2lYJn9ai5bTh6Lh6HEnWCyq/EkQ8o+e7GSMsCigAVj4x6t5rRNY2AKvOdl3GchrP81xtLYgF+Fe\naepXQVb0EvaH8hL3H+56agtWrn8TC2adje9ef0GpxSkLUrMgc87/HwCLALzj/izKgnJMIHt+dLbo\n3uxdq1NDV7eXci3Mh1oQ6nM9cxEajh53lkIYM5YBDpJT618t+3OOrseThw6jfXQ95gnLdIAFz7NI\ngiNXv8T8kNOcY7ldtc8dXR3e/00yyNZmoxVQ4vip487v97v6rI2uBa5B7avwnR31WS84UoybsJrb\n+DHnWUtlNixHQ/exvlRw6jmytEwGZWZpOvcGtO//M5rOvcH7zLMWd3UZra6pIKybt7X505+FrHjk\nraisXYKW7avyMpx4cymBhS1wlUq0O2aatv0oMQlhVm6tHDrrcBHulaZ+2caEBCL6BDjj8erDifqT\nikxEJli5/k30cI6V698stSj9DhsXC3DOX+OcP+D+bCy0UIQlbpGEuvHj0k9zFDfAwWa/gAdY0zud\naDx63HNrsCH0gXvNUjR9rQNVrh9xy/DhdlYXSU7tA0VOgzZqJOpqz0Zu1EgrF5em2U2oYs7l5wV7\n6ZQ1jRIjB+fJiinuHYnaqiEAHB9k4SOrtinkUxV5nQy5tpzjLuGmDQt9ICuKiDxuQnmXlXgTLTWn\nO2nbak7P/3LmIjS924320fVOKrigcxR0jKDtTS49+w72WWIjkEb+2ciuU/KLxOCxAO91PlfStKVu\nUTUo+IKmI++gveNNZxwt+pBK8FoRlttN5yfSC0EYYjzAEvUnVZmIkrJg1tmoZgwLZp1dalH6HVYK\nctZgjF3HGHuoq6ur1KKUlmuWomX4cF+2htSIa3Ex7BeqIIgH2MWL0fS1jj7/YSCeFUlD4ENBPYby\nd9gDxadsWb5cyG0alTVZiZHazeuzO+4dHxwD4CihUTM/mIqAAG5O4He78x7IYp/mHc1OIZbxZ/os\n9EE5rAPHRmMtDbMQm9rPteVw4YoLMX3FdK3fdiQlIaaiZTwXAXMlsVItv0hsXw9I2UjySBLxH3Xf\nKPcWeZzjWn9L7IqWOCZEHl/pPpnk5SaVOBUiE3z3+guw5/7PRXOvoAwfVljlQc4qFe+DjML5kuWa\nrzaXqw2SZUezU4723Bt82SZEvt5YfocJo5GD8iV7D04pP22ufom5dLHhYes7RuvSyPJancd7RyJ3\nRg1aaob1ZcpYuwR49REI5Sc37ly0DOLhJYc1x/dlX3CtpHl5dA1yJ8pJa0GSOST7XBdKtrBzZyxZ\nrpvbAf7dsZW9sCwJSa6xsH1t/GdDfM5t+6w9F+WezaBc5S/XGJlKoFznVEqkmcWCyDBpWgJ82Q1i\nVCZr2dniWLMVFwa1mltkolrtlLdjo/VOtmRJx2jZvtpd3h+Wn/lh7xq/9cs9VtORdxy/29alyJ15\njpM/NkLZY6vzOHMRWmqGObLJ+aEly2DTW/u9dpq2r0f7vgOu9TAYn6VYyoPdtH19oFxC7sYpjV4G\njUIs2+bNoQgWECEPi+C2E1W2MNcOo3uJbm67Y9/QfdQ6V3SotVnJkpC3fRIXhLB9VZmj+giLqn4b\nloeeb+25KPdsBoWUv5CWxHKNkakEyv2aKBKkIGeA2EupKd/c5IdLnOVn0z4Ng8eimnM0DhkXT5GP\nkyYroCiFN95TZ/XdJKRjNHQf9Ur1imN6Y1Nzuj9lmhxMJso8Dzhl/3KhO4em83rNUjRMne8f45mL\n4Mu6yHv79hPf8Z5I7h7evkEBigqiEMvmhZtjv6wFHSeskI1xXzdocPOHr/SnXksRm2vFuI1ubgcV\n8zA82MT8DEwxGORTn8Qf2eAGFCazj6hKtgHtOBfK19pE2kpnIeUvpBJLSlh2KfY1UaaQi0UGiJ32\nKOVlkoKl/okrZ5Ll5ID9Qsf7wdnA29ucKPzb2gAEpC6TU299eE6wW4K3zM0cH0JTqq2445WgraBz\nr5ZitnG9iHoMrQtFhGX36Sumg4ODgWGTnALN0P/+luYqatnnVPtvcFVKtdhFOa0587IAACAASURB\nVC3XZ2n5OmzcsjCuWZCBqCjIxaKMiB1RnPIbemx3jTCLSVw541o3At6OvawMQcvtt7UB93Q5v2UX\nCl2EuhxM5h43r7yv3B8AAHcU5XtHOrlu1bGJO166/TSfqdZWUwltr49KKWa45bGjpolSrZY+lx63\nIIc4HoDwZXfpHIty2FwNRAuxuOaVKY+4kpNGdoo0MGUmCdvethR5IAZXpcDtQsgb1whuFiXH5vot\nVpBU2Jir98pSBG+RKwaRUciCTCSnUBYTk2UhgcUhsrU+zb7JFmQAAA+2tCmFIqyPETI26hgEFR4x\nta9akIVFsnZ4LTq6OoxWYtlq6Svf3NWVV5gjcXBWxO3jrOQUrXhIGEmscEnnuO2xI8ioHdckcmbN\nSlksK3PUfhdCrnKwYhMVBVmQieIRZDFJYpEwWYITWBwiW+tFn2Tf3rhcs9SxTN/T6bhYmPyDRf/e\n3hbYT6310mJs1DHI8+8N68Pd7wDnfNL3sbDI7u7cbUx3p65QyHLoCnNEIerqR1CZ8obBYyMFAEbd\npyC8+oibL/iR6PsmXYmy9WeM4PeovU5NctrcY7JmpSyWf25UX9NCyBXVik0QGYEsyBmnrHwldZaA\nFCwS3hicZGg6uBc4bQTw3tHiWRwKZe0xtWtpQdZa2YpkjVGPrbUgx0h356NUvpxBxzWNb6n9Tu8Z\nLv0/Zn74AliCi4LN2GdN5lKTxniU63whKh6yIPcTyqokqM5SkIJFwhuDAaecD050FvdmG9AHbfW5\n5qtR9+gFyDVfHa/dkEpkgqgR+3c9tQUT7vw17npqS7BcAX0zHVtYZJ+c+6SX7g6jJ8Mqi4bJAliq\nKPg4frSljtgXKdwu/nL8NmytrDbbFdOX1WbsC2SlzIoPemTSsKjbtmEaeypWQWQcsiBnnLK3IKdA\nngUZKJ6lLkZGjLpHL0APY6jmHO1f2hJtXExFFeL4IyvtTfj9Z9HDOaoZw577Pxe6q9c3ztE+ut7+\nuLJFD/Bb9wq0ylA0+qk1zFTkR4vNnCyncxoF5fxnxgc9KsW0IJvor3NEpp/eL8odsiD3E8qqJGiB\nrDTeGHx5o6/YQVEIsZLorLgN3cf68iiHtaFaUdRtLf2RbdpbMOtsVDOGBbPOtup6w+SGvn5EKUwh\nW/RU656Q59WH+2QcPdn/G9Gt3UWjP/pLrnUqR+qK/GgRY3B4p3lOltqiXihCcqyXDWnM46Rt9Nc5\nIpM133ciEqQg91MiL/2lvdyVVntqO8VWUEJu4roXmKYJNzoBZxNuDG9DvYGq24q/x0yze5gEtPfd\n6y/Anvs/h+9efwEAaY40X609V02zm9A+ut4pt22RJs1DPkfq+ZLbETIe3un/DWDl+jfRwzlWrn8z\nuL9xWLsEuGeE47f74OzyWuYNuq7iXnMblve91EVR9ILmdX98kQDy+lxWBoys0V/niEwlvAT0Y8jF\nop8SZekv15ZDy/ZVaOg+5ihDaSx3pVCgIko7kShA+rhUZUlSIEUtRBJwLJ8Lxb4DkcY4keuPyY1E\n6u9dT23ByvVvYsGssz2FPjV5xJySycIyr815l68HkRdYbK+7VmzaDLse4rr3EHbQMjxBFJV+7WLB\nGLuOMfZQV1fMaO0KIMrSX8vOFqc0cs2w9JLbizfn0ZMD9wkNQozyBm4rn26Z3/08d0YN6g63phZ0\nE2rJ11lR4i7Lie1Zlf5Bq1sell1BQtL0yX1JZDmT+/zgbOc8jJ7sk1m1dgtkGXyFRqIEs6r9lEuH\n685XVoo6AP7rIWz1wbZNU6noqO49WaAcA79efbjvfkQQRGYoSwWZc/4M5/yrw4cPD984ZcolajmK\nAuMp01PnB1swoihuNn6KUBR53cMtyjKcra+vbpkfAGYuQkvNMOdlIaWsIZGzkKxd4lo3mfmlQDdO\nuv3U7ZSXlqbXnkH7vgNoeqfT8e3WVdOSxjTtjCq5thzqPnQMuZFnAG9vs7q2ZBnk/0fyBb1mqT/r\ngzS/tH00vVAlRYyzcPOQqioax0K+HlSFWHetRF3i1VXFs3XvUftVCiW1LH0+mfI7IuX4UkAQZUBZ\nKsilJEtp19JS1q2V6Tj+VFF8eJM+3Gx9fWUFSd72mqVomDo/1aAbreIW9EALswKrfYE7Dw63Ooqm\nvJ+q2KkvLW9v0x9LpyTNXJR6QJJv5WLMNP21pYyVLIOv0EhUi7bhxctYoEKQRPEyBVAKC+3hnZ5M\nVvcZmzLBSQpFWKYbzKOUSmpWfT6DrvmLF7sva4vjtV2WLwUEkX3IBzkiWUq7ljjFUFq+b1lrp5Bt\np9GObSEKwMpX2piOzfNJht4v1eRbWiSfSPVa0l5bQT63xaQQ5113HgDH1Ucp5W3FPSMAuPdzeUUg\nTZL4NFcyhUxpltXxzqpcRMVj64NMCnIZk1hZT+umnbAdbT8i3FxTe2mJGiQVd8xslOAIxwocP40S\nHDUALlJ/CvliIyxlQOEUwEJic97jVPETyNX0ChV0WAm5awtBJSqLNFeIjNKvg/QIh8QphtJajkzY\nTqDfp8WyYRS3l0C3FPmYttXd4vj/Cf/RDcv7AnTCcskGHEc7D0RfpGV7gTaFWhI/RttzFXesZJ9b\n+ZjlhtwX05gp88s3X8PG+eIvw/NjLZSLQVZdGLJOJaQ0U6n0uUK+4WUPKciVTFo37YTtGP0+LW+u\nkTN2mJTpoAwB4mYH+PuaJOOECKqzySUbdhxTQJ6m3byCIcIVI64fo+25SuorKXzHbcpXJ6EYDzbT\nmCnXkpivq7Y/gbaRc4PH+ZqlwD2dwD1d6Spi8nhUoqJHxKPS5wr5hpc95GJBVBTW7hjqkqhpuTBs\n6TRKzuWg/MBA8HHiyifvCxTWdSGNZWadX3Xay9cZWhrOteWwavsTOPXuJfjg7c9blQhPnTTHoxJd\nDYjKhOZ6ZiEfZCIcuoDtsRir0CCzMOVC3VYOuhozLbhYg0k+m+OHFRjJEjplPm6BDBNZKiSDAviM\nRyXKi1oYGXr5IAiiMiEfZCKcQiwBhZXCFSV+TcvXtsvba5c47dwzojg+XrrlQkVWrfuGrfvB2iUA\n7+3bB4CnHAPxizXojq+O8TVLnW3Ai7IcGCk9ocl1RLZ0xy2QYTqOaWm4REumpqIpRSOKq08Yle6X\nShBE2UAW5EqmEKVkgyxEsvXPZEGytTDZtFVoFFlzbTk073gCJ9+5BA21fx9NofH6w5y8xEKBEJbd\nMVPzz5Euw4NNOrS0La5RWLsEdYdb0cOYXXrCuBbHqP2xtbTHcacpd9K0IBMEQZQYsiAT4VhWu4sU\ntBRkIZI/Gz05+v7atph9Rokw63ZYH0MC4ZpmN+H49vvx3n9e788QYYNoC9xf0OSeLifwSleswVDQ\nI7Tym26MCxlQo1Tma+g+5pS4tik4YpoPYecrSYEME2Ft9tegHLXITiUHXhEOlKGBqADIgkyEW77i\nWPGS+MTGJaxtG+t2kFxyntl7urSbJPYXTRLUJ2+rKxAS9ZhpUYhiH1n0Za0EC7IaPNqf+knY48VH\nMOcFniDKCLIgE/aEWYXi+A2arGnCcmyyICchTM4w63ZoH5n/t8aKEstfNK7/q2lbkQ7N9pwV2hdd\nHtu0LJBZ9GXtr9ZVtV9pzxeyRpYhXPlNEP0PsiAThaEUFuQ0CLKO2aZ+i0oa/q9JKETbWT/PKmQV\ntacfp9UjLKHrhShjKjbN26lTp3Dw4EG89957JZIqnCFDhmDcuHEYOHBgqUUpDEny+ZaaKA/rKA8J\nS8X7rlOLMHnDPbh5wG9QdfGt8cYoiVxRXhDSkiELUJBe6aBxIwoFzS1CQ8UqyB0dHTj99NMxatQo\nMMYMe5YOzjmOHDmCo0ePora2ttTiFAZV2QhSPrJ2A7OVJ2quXMtCHr/4xy/glup1YAzxLWpJci8n\n9dNOg1LMibVLgFcfgbdkrCuYksTHnSBKRdbuscWErklCQ8X6IL/33nuZVY4BgDGGUaNGFc7CXWh/\nPpv2haWY9+b7oKrtvPpIdiL/ozxITH6Yps9NPrPK9jcP+A0Yc9W0uP61Ufxz1W0T+2mngK2Pa5pz\n/ZqlTno9WQaVJD7uOsj3ligG/TW7ig1ZjFUoBXSviUW/syBv27YN06ZNK5FE9hRMzkK/MVtkcrCS\nQ81jHNdimya2YydbG1VLY1Q5o2YIiJLNIi6FqEJXiDZ0uaOT9D/ovBaCYlq3KtmKWOnQuSfIku6j\nYi3IFU/B35iVTA5x5Fi7pE+xufjLdpH/SSqj2RI2drLVG9zZVpdBIkomA3X7OLl207YQJWnPZt+0\nchjrckcn4ZqlTsqqe7qKo0gU07pVyVbESkJ3bfXX7CqEPWRJj0VZKsiMsesYYw91dQVYMEvIrbfe\nijFjxuCCC0pQGrbQN8OLF7ulfhfHl0M8pFlVuoUc1GPoFALdA0R8Btgpp0I5ti1Qkia6cUjj5if6\n8ODsvpeX0ZPTf8kA0lPWxByLktIuS8S5VuPONXpAVgb0IkTooJekWJCLRQH43e9+h2HDhuGLX/wi\ntmzZot0mC3IWlLSyIaR9fN1SUxTXiqynros7trLLC+BaZhHcr7jHenA28PY2YMw0p0JgGm2mtX/W\nKfRc6+/j19+h80cQoZCLRQn51Kc+hZEjR5ZajNISZMlI8W32rqe2YMKdv8ZdTykvIrpjrF3iBA4C\n8SywJrnV4idpl7SOQlwLkhiDMdP6xiJsXOIe6/BO/+802kxrfxPFWCWwDYAtpCVYGT/j9UVkE7IU\nEkRqkIKMfvAQyEKEqiqD6UGubpdQ9pXr30QP51i5/s3wtjYsB8DxAa/CXackuZI+VFSFL0hJs1Hg\nkih5cRUoMQa3tfWNRSEqLIbtl1QBLJQCWYyla5tjyOckreveVPUQyvVFEARRQZCCjH7wEMiC35kq\ng23J5ISyL5h1NqoZw4JZZ/e19erDeqVh5iJ8wKvweM9n0z3XaadKS6LkRVX2kyhZcV8sgvZL+rKS\nlgXN9oUvTdLys0fAS7/say76J7ejjJ/v+iIIgqggyAcZzsNk5fo3sWDW2fju9ekE1u3btw/XXntt\ncXyQs+B3FrfARpqyr13iKMeA0UezEOe6rEnq05qFuVcIsuRXbiJg7Cfc+Wv0cI5qxrDn/s/1faHz\nNZ+5qH+ew2LRX68BguinVGwlvawEvxVVQSb6SFoJrxSUqnLchuWO3/ThnfGPXQ6KZByyND9iYHwR\nTOu8E33012uAIPopFKRXQubPn49LL70UO3bswLhx4/DII4+UWqRsUAxfadsl9iy4pRRKlqBUduIz\ncczDO5O5JERxCyiWr7zFcYQLQttPF+m3LfNgp+9ef4FnOfa5Wuh8zYlkUAo9guiXkIJcAFatWoVD\nhw7h1KlTOHjwIBYvDskZXClkSSlVH2opBw8mkiUGPp9Tm2IiaT3UoyiSKZx/q4Bai+OIuIOZ//VU\nduZkASj7+IpywHQNZCF4miCI2JCCTBSPQlla4jyI1IdaysGDkbhmaZ8faMz0bz5FyKaYSFwLaZKH\nfgrn30rhC3v5QV/w2YYPX9+vrX9pB9mVfcafYpIlgwBBEJEhH+QSUS5ylgVxfQBlP1OgcMGDUHxC\nBy7Pb9umDwHbFC34sMT+lrH6mREf0f4QIGoM/otJfxgTI2Xux04Q/RXyQSYqh7iWyYD0Vmn7oPos\nnzrLUsL0b8LntOBKRon9LWP1MyM+ov3B3SFti3R/GBMjZe7HThCVDinIRPkT90FURMXJp1jojmvT\nhyw8cIPcQbLqc5mFcUP/yCmc9otYfxgTorwgNyHCFnKxKBHlImdqFGC5sV8vz2YZk8tCRlwZCIIg\nTKTtJkSUH+RiQWSLAgSs9Ovl2SxjsrwXwSJP1h8iCJofRBi0akHYQgpyAThw4AAuv/xynHfeeTj/\n/PPxk5/8pNQilZ44ylPIkn3Jb3RZdSlwKZiyYHJZKIIrA70UEUHQ/CDCKFq8BlH2kIJcAAYMGICl\nS5di69ataGtrwz//8z9j69atpRartMRRnkKsziW/0WU8jVNBlIUSvxSU/KWojKhEayrND4Ig0oIU\n5AIwduxYfPzjHwcAnH766Zg2bRreeuutEktVIkqcNzcVTH3IgHxBSlBBlIUSvxSU/KWojKhEayrN\nD4Ig0oIUZKCgVrF9+/Zh48aNmDVrVuptlwVJFKqMZB8w9iED8gUpQQVRFjLwUkDYQdZUgiCI+JCC\nDBTMKnbs2DHceOON+PGPf4yamppU2y4b+oNCleE+FF0JysBLAWEHWVMJgiDiQwoyUBAF6NSpU7jx\nxhtx880344Ybbkit3bKjPyhUGe1D0jR3hfJRDWo3jWOWm29tuclLEARBkILskLICxDnH4sWLMW3a\nNHzzm99Mpc0sQQ/8bJDUx7RQPqpB7aZxzHLzrS03eQmCIAhSkAvCK6+8gl/84hd48cUXUVdXh7q6\nOvz6178utVjWhCnAlfLAz/qLgMm9wlbuQrlnBLWbxjHLzbe23OQlCIIgqJJeySiEnGlVlgurNFQp\nFezKteJSucpdCCplrhIEQRB2UCW9CuOup7bgF237U7Hshlm8TME/Wbe4RiXrlj/TeGddbsFdT21B\n7R1rMf6OtQWbM5Ww2tHfrjuCIIgsQApyP0FWAJIqRnGj3/ubMpL1LACm8c663IKV698El/5fCKK+\nLJSjstnfrjuCIIgsQApyP0EoAv999jmpKkZRFIYgZaQcFY+sy1wulmITC2adDSb9vxDILws257Mc\nlc1ynwcEQRBZhHyQS0S5yJmWP2s5+sUGyVwuvq3lImcxsJmDNF4EQRD9G/JBJlIhLetUOVq5gmQu\nF0tjuchZDGzmYLm4pxAEQRCFhSzIJaJc5OyPpGElLBdLY7nISRAEQRDFoOwsyIyxoYyxFYyx/8MY\nu7nU8iThvffewyWXXIKPfexjOP/883H33XeXWiRCIg2rarlYGoWcADLtT00QBEEQWaKgCjJjbBlj\n7G3G2Bbl86sYYzsYY7sZY3e4H98A4Jec868A+JtCylVoBg8ejBdffBGvv/462tvb8fzzz6Otra3U\nYhEu5ejukRRytSAIgiAIewptQX4UwFXyB4yxagD/DOBqAOcBmM8YOw/AOAAH3M16CixXQWGMYdiw\nYQCAU6dO4dSpU2CMhexFFItysf6mSdKXgqxn9CAIgiCINCmogsw5/x2Ad5SPLwGwm3O+l3N+EsBq\nAHMBHISjJBdcLpVcWw51P69Dri2XWps9PT2oq6vDmDFjMGfOHMyaNSu1tgkiKklfCsgCTRAEQVQS\npfBBPhN9lmLAUYzPBLAGwI2MsX8B8IxpZ8bYVxljGxhjGw4fPpyKQC07W9DDe9CysyWV9gCguroa\n7e3tOHjwIP74xz9iy5b+b3kjK2M+pjEpt7GqRLcUgiAIonLJTJAe5/w453wR5/x/cM4fD9juIc75\nTM75zNGjR6dy7IbJDahm1WiY3JBKezIjRozA5Zdfjueffz71trMGWRnzMY1JuY1VJbqlEARBEJVL\nKRTktwCcJf09zv2sZDTNbkL7F9vRNLsplfYOHz6Mzs5OAMCJEyewbt06TJ06NZW2s0y5WxkLYdU1\njUm5jxVBEARB9GcKngeZMTYewLOc8wvcvwcA2Angs3AU41cBLOCcvxG17azmQd60aRMWLlyInp4e\n9Pb24qabbsJ3vvMd3zZZkDOISsyfW47V/giCIAiCsCcTeZAZY6sA/AHAFMbYQcbYYs75BwD+DkAr\ngG0AnoijHGeZ6dOnY+PGjdi0aRO2bNmSpxyXA1FcAMrNn9ZEoay6/WV8CIIgCKJSKHQWi/mc87Gc\n84Gc83Gc80fcz3/NOZ/MOZ/AOY+cOoIxdh1j7KGurq70hSYARFMWy82f1kQcP1sb5be/jA9BEARB\nVAqZCdKLAuf8Gc75V4cPH15qUfotUZTFSvantVF+5fEhazJBEARBZJ+yVJCJ9EhDYavkDAc2Lwfy\n+JA1mSAIgiCyDynIFQ4pbMmwfTkQLyITxgytWGs7QRAEQZQLA0otAFFaFsw628tWQRQO8SKy5+3j\nlCGDIAiCIDJOWVqQKUgvPSrZPaKYlMJPm/ydCYIgCCIeZakgl0OQXk9PD2bMmIFrr7221KIQGaAU\nLyLkPkMQBEEQ8ShLBbkc+MlPfpLpQiBE/6eSs4sQBEEQRBJIQS4ABw8exNq1a/HlL3+51KIQGaOY\nbg/kPkMQBEEQ8SAFGcChe+/DtvMvwKF770ulvdtvvx0//OEPUVVFw2tLFv1lCyETuT0QBEEQRPYh\nDQ5A5xNPAD09zu+EPPvssxgzZgwuuuiiFCSrHLKoOBZCJnJ7IAiCIIjsQwoygBE33QRUVzu/E/LK\nK6/gV7/6FcaPH4/Gxka8+OKLuOWWW1KQsn+TRcWxEDKR2wNBEARBZB/GOS+1DJFhjF0H4LqJEyd+\nZdeuXb7vtm3blpnguN/+9rf40Y9+hGeffTbvuyzJSRAEQRAEUQkwxv7EOZ8Ztl1ZWpDLIc0bQRAE\nQRAEUZ5QJb0Cctlll+Gyyy4rtRgEQRAEQRBEBMrSgkwQBEEQBEEQhYIUZIIgCIIgCIKQ6JcKctYD\nD7MuH0EQBEEQRCXT7xTkIUOG4MiRI5lVQjnnOHLkCIYMGVJqUQiCIAiCIAgNZRmkJ6V5y/tu3Lhx\nOHjwIA4fPlx8wSwZMmQIxo0bV2oxCIIgCIIgCA1lmQdZMHPmTL5hw4ZSi0EQBEEQBEGUAf06DzJB\nEARBEARBFApSkAmCIAiCIAhCghRkgiAIgiAIgpAoax9kxthhAPtTbnY4gK4SthNlP5ttk25j+u7D\nAP4rpN2skNY5LfQxsjJnbLajOZONY/T3OQOUz7wplzkTtx2aM+lDcyb69knmjPh+BOd8dKg0nHP6\nkX4APFTKdqLsZ7Nt0m1M3wHYUOpzVexzWuhjZGXO2GxHcyYbx+jvc8b9rizmTbnMmbjt0JzJ7vks\nxjEKPWdst08yZ6LKRC4W+TxT4nai7GezbdJt0hqPUlKMPqRxjKzMGZvtaM5k4xg0Z7JDucyZuO3Q\nnEkfmjPRt08yZ2yPAaDMXSyI0sEY28At0qQQhIDmDBEHmjdEVGjOEGlAFmQiLg+VWgCi7KA5Q8SB\n5g0RFZozRGLIgkwQBEEQBEEQEmRBJgiCIAiCIAgJUpAJgiAIgiAIQoIUZIIgCIIgCIKQIAWZIAiC\nIAiCICRIQSZShzF2PWPs/zDGmhljV5ZaHiL7MMbOZYw9whj7ZallIbILY2woY2yFe3+5udTyENmH\n7i1EXEhBJnwwxpYxxt5mjG1RPr+KMbaDMbabMXZHUBuc86c4518B8LcA5hVSXqL0pDRn9nLOFxdW\nUiKLRJw/NwD4pXt/+ZuiC0tkgihzhu4tRFxIQSZUHgVwlfwBY6wawD8DuBrAeQDmM8bOY4xdyBh7\nVvkZI+36j+5+RP/mUaQ3Z4jK41FYzh8A4wAccDfrKaKMRLZ4FPZzhiBiMaDUAhDZgnP+O8bYeOXj\nSwDs5pzvBQDG2GoAcznn9wO4Vm2DMcYAfB/Ac5zz1worMVFq0pgzROUSZf4AOAhHSW4HGXgqlohz\nZmtxpSP6C3SDIWw4E31WG8B5SJ0ZsP3XAFwB4AuMsb8tpGBEZok0ZxhjoxhjPwMwgzF2Z6GFIzKP\naf6sAXAjY+xfADxTCsGIzKKdM3RvIeJCFmQidTjnDwB4oNRyEOUD5/wIHJ91gjDCOT8OYFGp5SDK\nB7q3EHEhCzJhw1sAzpL+Hud+RhAmaM4QSaD5Q0SF5gyRKqQgEza8CmASY6yWMTYIQCOAX5VYJiLb\n0JwhkkDzh4gKzRkiVUhBJnwwxlYB+AOAKYyxg4yxxZzzDwD8HYBWANsAPME5f6OUchLZgeYMkQSa\nP0RUaM4QxYBxzkstA0EQBEEQBEFkBrIgEwRBEARBEIQEKcgEQRAEQRAEIUEKMkEQBEEQBEFIkIJM\nEARBEARBEBKkIBMEQRAEQRCEBCnIBEEQBEEQBCFBCjJBEESRYIzdzhj7kPT3rxljI1Jo9zLG2LOG\n71YxxjYxxr7BGLuPMXaFThaCIAiiD8qDTBAEUSQYY/sAzOSc/1fK7V4G4B8459cqn38UwO855xPT\nloUxNsAtzqD923Y/giCILEIWZIIgiJgwxr7JGNvi/tzufjaeMbadMfY4Y2wbY+yXjLEPMca+DuCv\nALzEGHvJ3XYfY+zD0j6PMsZ2uvtewRh7hTG2izF2ibv9JYyxPzDGNjLG/p0xNiVExBcAnMkYa2eM\n/bXb/hd0sij9uogx9jJj7E+MsVbG2Fj3898yxn7MGNsA4O/d9n7GGFsP4IeMsZGMsadci3UbY2y6\nu989jLFfMMZeAfCLVAafIAiigJCCTBAEEQPG2EUAFgGYBWA2gK8wxma4X08B8CDnfBqAbgC3cc4f\nAPBnAJdzzi/XNDkRwFIAU92fBQA+CeAfAHzb3WY7gL/mnM8A8B0A/xQi5t8A2MM5r+Oc/5v4MEgW\nxthAAD8F8AXO+UUAlgHISZsM4pzP5Jwvdf8eB+C/cc6/CeBeABs559NdmX8u7XcegCs45/NDZCYI\ngig5A0otAEEQRJnySQBPcs6PAwBjbA2AvwbwKwAHOOevuNs9BuDrAH4U0l4H53yz29YbAH7DOeeM\nsc0AxrvbDAewgjE2CQAHMDDF/gimALgAwDrGGABUAzgkfd+sbN/COe9x//9JADcCAOf8RcbYKMZY\njfvdrzjnJwogL0EQROqQgkwQBJE+anCHTbDH+9L/e6W/e9F3r/4ugJc4559njI0H8Nv4IhphAN7g\nnF9q+P54yN8mbLcjCIIoOeRiQRAEEY9/A3C96188FMDn3c8A4GzGmFAwFwD4vfv/owBOT3DM4QDe\ncv//pQTtBMmyA8BoIT9jbCBj7HzLNv8NwM3ufpcB+C/OeXdCOQmCIIoOKcgEQRAx4Jy/BuBRAH8E\nsB7Aw5zzje7XOwD8T8bYNgBnAPgX9/OHADyvC4yz5IcA7meMbUTyFUCt1Yr5twAAAI1JREFULJzz\nkwC+AOAHjLHXAbQD+G+Wbd4D4CLG2CYA3wewMKGMBEEQJYHSvBEEQaSI6/rwLOf8ghKLQhAEQcSE\nLMgEQRAEQRAEIUEWZIIgCIIgCIKQIAsyQRAEQRAEQUiQgkwQBEEQBEEQEqQgEwRBEARBEIQEKcgE\nQRAEQRAEIUEKMkEQBEEQBEFI/P8FRIfX4QCXjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113980b90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "for counter in range(1,5):\n",
    "    plt.loglog(all_them_best_times_and_errors[:,1][all_them_best_times_and_errors[:,2]==counter],10**3*all_them_best_times_and_errors[:,0][all_them_best_times_and_errors[:,2]==counter],'o',ms=2,label=str(counter))\n",
    "plt.ylabel('optimal fit time / ms')\n",
    "plt.xlabel('optimal fit error')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(fig_dir+'errors_and_times.png',dpi=100)\n",
    "plt.savefig(fig_dir+'errors_and_times.jpg',dpi=100)\n",
    "plt.savefig(fig_dir+'errors_and_times.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if False: # In here one would create new diamond datasets.\n",
    "    learning.create_diamonds([20], omega_larmor, num_spins = 450)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:python2]",
   "language": "python",
   "name": "conda-env-python2-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}