infMagSim.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.parallel import Client"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rc = Client()"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "import IPython.core.display as IPdisp"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "from IPython.nbformat import current\n",
      "def execute_notebook(nbfile):\n",
      "    with io.open(nbfile) as f:\n",
      "        nb = current.read(f, 'json')\n",
      "    ip = get_ipython()\n",
      "    for cell in nb.worksheets[0].cells:\n",
      "        if cell.cell_type != 'code':\n",
      "            continue\n",
      "        ip.run_cell(cell.input)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#####ScriptStart#####"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "from mpi4py import MPI\n",
      "comm = MPI.COMM_WORLD\n",
      "mpi_id = comm.Get_rank()\n",
      "n_engines = comm.size"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "import io\n",
      "import math\n",
      "import numpy as np\n",
      "if (mpi_id==0):\n",
      "    import matplotlib as mpl\n",
      "    mpl.use('Agg') # non-GUI backend\n",
      "    import matplotlib.pyplot as plt\n",
      "#import ghalton\n",
      "from scipy.stats import norm\n",
      "from infMagSim_cython import *"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "read_ion_dist_from_file = False\n",
      "smooth_ion_velocities = False\n",
      "read_electron_dist_from_file = read_ion_dist_from_file\n",
      "smooth_electron_velocities = smooth_ion_velocities\n",
      "dist_filename = '/home/chaako/analysis/runTmp/dist_func.npz'\n",
      "project_read_particles = False\n",
      "projection_angle = np.pi*5/16.\n",
      "create_electron_dimple = False # Only use with periodic_particles for now, and not with quiet start\n",
      "ion_distribution_gap_center = 0. # Simulate beams by splitting ion distribution\n",
      "ion_distribution_gap_half_width = 0. # Simulate beams by splitting ion distribution\n",
      "periodic_particles = False\n",
      "periodic_potential = periodic_particles\n",
      "prescribe_potential = False\n",
      "def prescribed_potential(grid, t):\n",
      "    z_min = grid[0]\n",
      "    z_max = grid[-1]\n",
      "    return min(1.,t/(z_max-z_min))*np.sin(2.*np.pi*grid/(z_max-z_min))\n",
      "solve_for_electric_field = False\n",
      "simulate_moon = True\n",
      "#v_drift_moon = 25.\n",
      "#moon_radius = 1.\n",
      "#t_moon_center = moon_radius/v_drift_moon\n",
      "boltzmann_electrons = False # Don't move electrons (assume Boltzmann)\n",
      "boltzmann_potential = boltzmann_electrons # Use Boltzmann relation to solve for potential\n",
      "decouple_ions_from_electrons = False # Use full potential to move electrons even if Boltzmann for ions\n",
      "quasineutral = False\n",
      "if quasineutral:\n",
      "    boltzmann_electrons = True\n",
      "quiet_start_and_injection = True\n",
      "use_quasirandom_numbers = False\n",
      "use_quasirandom_dimensions_for_parallelism = False\n",
      "include_object = False\n",
      "use_pure_c_mover = True\n",
      "use_pure_c_solver = use_pure_c_mover"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "def circular_cross_section(grid, t, t_center, v_drift, radius, object_mask):\n",
      "    # object_mask will give fraction of each cell inside object\n",
      "    eps = 1.e-5\n",
      "    radius += eps\n",
      "    n_points = len(grid)\n",
      "    n_cells = len(object_mask)\n",
      "#    if (n_points==n_cells):\n",
      "#        n_cells -= 1\n",
      "    z_min = grid[0]\n",
      "    z_max = grid[n_points-1]\n",
      "    z_grid_center = (z_max+z_min)/2\n",
      "    object_index = n_points/2\n",
      "    z_object = grid[object_index]\n",
      "    dz = grid[1]-grid[0]\n",
      "    if (math.fabs(z_object-z_grid_center)>eps):\n",
      "        print 'object not centered'\n",
      "    t_enter = t_center - radius/v_drift\n",
      "    t_leave = t_center + radius/v_drift\n",
      "    if (t<t_enter):\n",
      "        width = 0.\n",
      "        distance_to = (t_enter-t)*v_drift\n",
      "    elif (t>t_leave):\n",
      "        width = 0.\n",
      "        distance_to = (t-t_leave)*v_drift\n",
      "    else:\n",
      "        x = (t-t_center)*v_drift\n",
      "        width = math.sqrt(max(0,1-x*x)) # Circular cross-section\n",
      "        distance_to = 0.\n",
      "#    if (distance_to>0. and enforce_bdy_cond_outside):\n",
      "#        width = eps\n",
      "    for j in range(n_cells):\n",
      "        object_mask[j] = 0.\n",
      "    if (width>0.):\n",
      "        z_object_left = z_object - width\n",
      "        outside_index_left = int((z_object_left-z_min)/dz)\n",
      "        zeta_left = z_object_left % dz\n",
      "        if (zeta_left==0.):\n",
      "            zeta_left += eps\n",
      "        elif (zeta_left>dz-eps):\n",
      "            zeta_left = dz-eps\n",
      "        object_mask[outside_index_left] = 1.-zeta_left/dz\n",
      "\n",
      "        z_object_right = z_object + width\n",
      "        outside_index_right = int((z_object_right-z_min)/dz)+1\n",
      "        zeta_right = dz - (z_object_right % dz)\n",
      "        if (zeta_right==0.):\n",
      "            zeta_right += eps\n",
      "        elif (zeta_right>dz-eps):\n",
      "            zeta_right = dz-eps\n",
      "        object_mask[outside_index_right-1] = 1.-zeta_right/dz\n",
      "\n",
      "        if (width>dz):\n",
      "            for j in range(outside_index_left+1, outside_index_right-1):\n",
      "                object_mask[j] = 1.\n",
      "    return distance_to"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "z_min = -5.\n",
      "z_max = 5.\n",
      "n_cells = 2000\n",
      "n_points = n_cells+1\n",
      "eps = 1e-6\n",
      "grid, dz = np.linspace(z_min,z_max,num=n_points,endpoint=True,retstep=True)\n",
      "grid = grid.astype(np.float32)\n",
      "if quasineutral:\n",
      "    debye_length = dz\n",
      "else:\n",
      "    debye_length = 0.125\n",
      "number_of_debye_lengths_across = int((z_max-z_min)/debye_length)\n",
      "number_of_uniforming_bins = number_of_debye_lengths_across*4\n",
      "number_per_injection_batch = 100\n",
      "if (mpi_id==0):\n",
      "    print dz, debye_length"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[stdout:0] 0.05\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "seed = 384+mpi_id*mpi_id*mpi_id\n",
      "class uniform_2d_sampler_class:\n",
      "    def __init__(self,seed,low_discrepancy=False,shared_seed=True,rand_dim=1,parallel_dimensions=1,\\\n",
      "\t\tnumbers_per_batch=10000000, mpi_id=0):\n",
      "\tself.low_discrepancy = low_discrepancy\n",
      "\tself.shared_seed = shared_seed\n",
      "\tself.rand_dim = rand_dim\n",
      "\tself.parallel_dimensions = parallel_dimensions\n",
      "\tself.total_dim = rand_dim*parallel_dimensions\n",
      "\tself.numbers_per_batch = numbers_per_batch\n",
      "\tself.mpi_id = mpi_id\n",
      "\tself.total_n = 0\n",
      "\tif low_discrepancy:\n",
      "\t    self.low_discrepancy = True\n",
      "\t    #self.sequencer = ghalton.GeneralizedHalton(self.total_dim,seed) # seed doesn't appear to do anything\n",
      "\t    self.sequencer = sobol_sequencer(rand_dim=self.total_dim)\n",
      "\telse:\n",
      "\t    self.low_discrepancy = False\n",
      "\t    np.random.seed(seed)\n",
      "    def get_n(self,n):\n",
      "\tself.total_n = self.total_n + n\n",
      "\tif self.low_discrepancy:\n",
      "\t    if self.total_n>1e9:\n",
      "\t        print \"WARNING: low-discrepancy sequences may start returning bad values.\"\n",
      "\t    return np.asarray(self.sequencer.get(n))\n",
      "\telse:\n",
      "\t    return np.random.rand(n,self.total_dim)\n",
      "    def get(self,n):\n",
      "\tn_per_batch = int(self.numbers_per_batch/self.total_dim)\n",
      "\tn_batches = max(1,n/n_per_batch)\n",
      "\tnumbers = np.zeros([n,self.rand_dim])\n",
      "\tfor i in range(n_batches):\n",
      "\t    if i==n_batches-1:\n",
      "\t        n_to_get = n-n_per_batch*(n_batches-1)\n",
      "\t    else:\n",
      "\t\tn_to_get = n_per_batch\n",
      "\t    number_batch = self.get_n(n_to_get)\n",
      "\t    for j in range(self.rand_dim):\n",
      "\t\tnumbers[i*n_per_batch:i*n_per_batch+n_to_get,j] = number_batch[:,self.parallel_dimensions*j+self.mpi_id]\n",
      "\treturn numbers\n",
      "    def get_uniform(self,n,n_bins):\n",
      "\tsamples = self.get(n)\n",
      "\tsamples[:,0] = (np.arange(n)+samples[:,0])/float(n)\n",
      "\t# TODO: avoid each engine having extra samples in the same bins (might not matter since just for velocity smoothing)\n",
      "\tsample_bins = np.arange(0.,float(n_bins),float(n_bins)/float(n)).astype(np.int32)\n",
      "\tnumber_in_bins = np.bincount(sample_bins) # TODO: assuming here that there is at least one per bin\n",
      "\tbin_ends = np.cumsum(number_in_bins)\n",
      "\tfor i in range(1,self.rand_dim):\n",
      "\t    for j in range(n_bins):\n",
      "\t\tbin_width = 1./number_in_bins[j]\n",
      "\t\tcurrent_indices = slice(bin_ends[j]-number_in_bins[j],bin_ends[j])\n",
      "\t\tsamples[current_indices,i] = \\\n",
      "\t\t    np.random.permutation((np.arange(number_in_bins[j])+samples[current_indices,i])*bin_width)\n",
      "\treturn samples\n",
      "\n",
      "class injection_sampler_class:\n",
      "    def __init__(self,sampler,n_per_batch):\n",
      "\tself.sampler = sampler\n",
      "\tself.shared_seed = sampler.shared_seed\n",
      "\tself.n_per_batch = n_per_batch\n",
      "\tself.n_left_in_batch = 0\n",
      "\tself.samples = np.zeros([0,sampler.rand_dim])\n",
      "    def get(self,n):\n",
      "\tif self.n_left_in_batch<n:\n",
      "\t    n_extra = n - self.n_left_in_batch\n",
      "\t    n_batches = int(math.ceil(float(n_extra)/float(self.n_per_batch)))\n",
      "\t    n_to_draw = n_batches*self.n_per_batch\n",
      "\t    self.n_left_in_batch = n_to_draw - n_extra\n",
      "\t    bin_width = 1./self.n_per_batch\n",
      "\t    samples = self.sampler.get(n_to_draw)\n",
      "\t    for i in range(1,self.sampler.rand_dim):\n",
      "\t\tfor j in range(n_batches):\n",
      "\t\t    current_indices = slice(j*self.n_per_batch,(j+1)*self.n_per_batch)\n",
      "\t\t    samples[current_indices,i] = \\\n",
      "\t\t\tnp.random.permutation((np.arange(self.n_per_batch)+samples[current_indices,i])*bin_width)\n",
      "\t    samples = np.concatenate((self.samples,samples))\n",
      "\telse:\n",
      "\t    samples = self.samples\n",
      "\t    self.n_left_in_batch = self.n_left_in_batch - n\n",
      "\tif self.n_left_in_batch>0:\n",
      "\t    self.samples = samples[-self.n_left_in_batch:]\n",
      "\telse:\n",
      "\t    self.samples = np.zeros([0,self.sampler.rand_dim])\n",
      "\treturn samples[:n]\n",
      "\n",
      "if use_quasirandom_numbers:\n",
      "    if use_quasirandom_dimensions_for_parallelism:\n",
      "        uniform_2d_sampler = uniform_2d_sampler_class(seed,low_discrepancy=True,rand_dim=2,\\\n",
      "\t\t\t       shared_seed=False,parallel_dimensions=n_engines,mpi_id=mpi_id)\n",
      "    else:\n",
      "        uniform_2d_sampler = uniform_2d_sampler_class(seed,low_discrepancy=True,rand_dim=2,\\\n",
      "\t\t\t       shared_seed=True,parallel_dimensions=1)\n",
      "else:\n",
      "    uniform_2d_sampler = uniform_2d_sampler_class(seed,shared_seed=False,rand_dim=2)\n",
      "\n",
      "if quiet_start_and_injection:\n",
      "    injection_sampler = injection_sampler_class(uniform_2d_sampler,number_per_injection_batch)\n",
      "else:\n",
      "    injection_sampler = uniform_2d_sampler"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "ERROR: Cell magic `%%px` not found.\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "v_th_i = 1\n",
      "v_d_i = 0\n",
      "n_ions = 100000\n",
      "extra_storage_factor = 1.2\n",
      "background_ion_density = n_ions*n_engines*dz/(z_max-z_min)\n",
      "n_bins = n_points\n",
      "v_max_i = 4.*v_th_i\n",
      "v_min_i = -v_max_i\n",
      "inactive_slot_position_flag = 2*z_max"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "largest_ion_index = [n_ions-1]\n",
      "ion_storage_length = int(extra_storage_factor*n_ions)\n",
      "ions = np.zeros([2,ion_storage_length],dtype=np.float32)\n",
      "#ions[0][0:n_ions] = np.random.rand(n_ions)*(z_max-z_min) + z_min # positions\n",
      "if uniform_2d_sampler.shared_seed:\n",
      "    for id in range(n_engines):\n",
      "\tsample = np.asarray(uniform_2d_sampler.get(n_ions)).T\n",
      "\tif id==mpi_id:\n",
      "\t    uniform_2d_sample = sample\n",
      "elif quiet_start_and_injection:\n",
      "    uniform_2d_sample = uniform_2d_sampler.get_uniform(n_ions,number_of_uniforming_bins).T\n",
      "else:\n",
      "    uniform_2d_sample = uniform_2d_sampler.get(n_ions).T\n",
      "ions[0][0:n_ions] = uniform_2d_sample[0]*(z_max-z_min) + z_min # positions\n",
      "if (read_ion_dist_from_file):\n",
      "    dist_file = np.load(dist_filename)\n",
      "    ion_distribution_function = dist_file['ion_distribution_function']\n",
      "    ion_v_edges = dist_file['ion_hist_v_edges']\n",
      "    weights = ion_distribution_function/np.sum(ion_distribution_function)\n",
      "    # TODO: option to sample bins quasi-randomly?\n",
      "    ion_sample_bins = np.random.choice(len(weights),size=n_ions,p=weights)\n",
      "    left_edge_weights = uniform_2d_sample[1]\n",
      "    ions[1][0:n_ions] = left_edge_weights*ion_v_edges[ion_sample_bins] +\\\n",
      "\t(1.-left_edge_weights)*ion_v_edges[ion_sample_bins+1]\n",
      "else:\n",
      "    #ions[1][0:n_ions] = np.random.randn(n_ions)*v_th_i + v_d_i # velocities\n",
      "    #ions[1][0:n_ions] = norm.ppf(uniform_2d_sample[1])*v_th_i + v_d_i # velocities\n",
      "    ions[1,0:n_ions] = uniform_2d_sample[1] # velocities\n",
      "    inverse_gaussian_cdf_in_place(ions[1,0:n_ions]) # velocities\n",
      "    ions[1,0:n_ions] *= v_th_i # velocities\n",
      "    ions[1,0:n_ions] += v_d_i # velocities\n",
      "    #ions[2][0:n_ions] = np.ones(n_ions) # active slots\n",
      "    #ions[2][0:n_ions] = np.ones(n_ions,dtype=np.float32) # relative weights\n",
      "if smooth_ion_velocities:\n",
      "    # TODO: allow quasi-random numbers here?\n",
      "    ion_v_smoothing_scale = v_th_i/10.\n",
      "    if read_ion_dist_from_file:\n",
      "\tion_v_smoothing_scale = ion_v_edges[1]-ion_v_edges[0]\n",
      "    ions[1][0:n_ions] += np.random.randn(n_ions)*ion_v_smoothing_scale\n",
      "if project_read_particles:\n",
      "    if uniform_2d_sampler.shared_seed:\n",
      "\tfor id in range(n_engines):\n",
      "\t    sample = np.asarray(uniform_2d_sampler.get(n_ions)).T\n",
      "\t    if id==mpi_id:\n",
      "\t\tuniform_2d_sample = sample\n",
      "    elif quiet_start_and_injection:\n",
      "\tuniform_2d_sample = uniform_2d_sampler.get_uniform(n_ions,number_of_uniforming_bins).T\n",
      "    else:\n",
      "\tuniform_2d_sample = uniform_2d_sampler.get(n_ions).T\n",
      "    perpendicular_velocities = norm.ppf(uniform_2d_sample[1])*v_th_i\n",
      "    velocity_magnitudes = np.sqrt(ions[1][0:n_ions]*ions[1][0:n_ions]\n",
      "\t\t\t\t  + perpendicular_velocities*perpendicular_velocities)\n",
      "    velocity_angles = np.arctan2(perpendicular_velocities, ions[1][0:n_ions])\n",
      "    ions[1][0:n_ions] = velocity_magnitudes*np.cos(velocity_angles-projection_angle)\n",
      "if ion_distribution_gap_half_width>0.:\n",
      "    ions[1,0:n_ions] += np.sign(ions[1,0:n_ions]-ion_distribution_gap_center)*ion_distribution_gap_half_width\n",
      "ions[0][n_ions:] = inactive_slot_position_flag\n",
      "# List remaining slots in reverse order to prevent memory fragmentation\n",
      "empty_ion_slots = -np.ones(ion_storage_length,dtype=np.int32)\n",
      "current_empty_ion_slot = [(ion_storage_length-n_ions)-1]\n",
      "empty_ion_slots[0:(current_empty_ion_slot[0]+1)] = range(ion_storage_length-1,n_ions-1,-1)\n",
      "ion_hist_n_edges = np.arange(z_min,z_max+eps,(z_max-z_min)/n_bins)\n",
      "ion_hist_v_edges = np.arange(v_min_i,v_max_i+eps,(v_max_i-v_min_i)/n_bins)\n",
      "ion_hist_n, ion_hist_n_edges = np.histogram(ions[0][0:n_ions], bins=ion_hist_n_edges)\n",
      "ion_hist_v, ion_hist_v_edges = np.histogram(ions[1][0:n_ions], bins=ion_hist_v_edges)\n",
      "comm.Allreduce(MPI.IN_PLACE, ion_hist_n, op=MPI.SUM)\n",
      "comm.Allreduce(MPI.IN_PLACE, ion_hist_v, op=MPI.SUM)\n",
      "if (mpi_id==0):\n",
      "    fig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,2))\n",
      "    for ax, data, bins in zip(axes,[ion_hist_n,ion_hist_v], [ion_hist_n_edges, ion_hist_v_edges]):\n",
      "\tax.step(bins,np.append(data,[0]), where='post')\n",
      "    filename='initialIonDistribution.png'\n",
      "    plt.savefig(filename)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "#if not boltzmann_electrons:\n",
      "mass_ratio = 1./1836.\n",
      "n_electrons = n_ions\n",
      "v_th_e = 1./math.sqrt(mass_ratio)\n",
      "v_d_e = 0.\n",
      "background_electron_density = n_electrons*n_engines*dz/(z_max-z_min)\n",
      "n_bins = n_points;\n",
      "v_max_e = 4.*v_th_e\n",
      "v_min_e = -v_max_e\n",
      "dimple_width = v_th_e/10.\n",
      "dimple_velocity = 0.1*v_th_e\n",
      "dimple_height = 0.1\n",
      "def dimple(x, mu=dimple_velocity, sig=dimple_width, height=dimple_height):\n",
      "    return height*np.exp( -np.power(x-mu,2.) / (2*np.power(sig,2.)) )"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "electron_hist_n_edges = np.arange(z_min,z_max+eps,(z_max-z_min)/n_bins)\n",
      "electron_hist_v_edges = np.arange(v_min_e,v_max_e+eps,(v_max_e-v_min_e)/n_bins)\n",
      "if not boltzmann_electrons:\n",
      "    largest_electron_index = [n_electrons-1]\n",
      "    electron_storage_length = int(n_electrons*extra_storage_factor)\n",
      "    electrons = np.zeros([2,electron_storage_length],dtype=np.float32)\n",
      "    n_samples = n_electrons\n",
      "    n_electrons_so_far = 0\n",
      "    while n_electrons_so_far<n_electrons:\n",
      "\tif (n_samples==n_electrons): # First iteration attempt to initialize all electrons\n",
      "\t    remaining_indices = slice(0,n_electrons)\n",
      "\t#electrons[0][remaining_indices] = np.random.rand(n_samples)*(z_max-z_min) + z_min # positions\n",
      "\tif uniform_2d_sampler.shared_seed:\n",
      "\t    for id in range(n_engines):\n",
      "\t\tsample = np.asarray(uniform_2d_sampler.get(n_samples)).T\n",
      "\t\tif id==mpi_id:\n",
      "\t\t    uniform_2d_sample = sample\n",
      "\telif quiet_start_and_injection:\n",
      "\t    uniform_2d_sample = uniform_2d_sampler.get_uniform(n_samples,number_of_uniforming_bins).T\n",
      "\telse:\n",
      "\t    uniform_2d_sample = uniform_2d_sampler.get(n_samples).T\n",
      "\telectrons[0,remaining_indices] = uniform_2d_sample[0]*(z_max-z_min) + z_min # positions\n",
      "\tif (read_electron_dist_from_file):\n",
      "\t    dist_file = np.load(dist_filename)\n",
      "\t    electron_distribution_function = dist_file['electron_distribution_function']\n",
      "\t    electron_v_edges = dist_file['electron_hist_v_edges']\n",
      "\t    weights = electron_distribution_function/np.sum(electron_distribution_function)\n",
      "\t    # TODO: option to sample bins quasi-randomly?\n",
      "\t    electron_sample_bins = np.random.choice(len(weights),size=n_samples,p=weights)\n",
      "\t    left_edge_weights = uniform_2d_sample[1]\n",
      "\t    electrons[1,remaining_indices] = left_edge_weights*electron_v_edges[electron_sample_bins] +\\\n",
      "\t\t(1.-left_edge_weights)*electron_v_edges[electron_sample_bins+1]\n",
      "\telse:\n",
      "\t    #electrons[1,remaining_indices] = np.random.randn(n_samples)*v_th_e + v_d_e # velocities\n",
      "\t    electrons[1,remaining_indices] = norm.ppf(uniform_2d_sample[1]) # velocities\n",
      "\t    #electrons[1,remaining_indices] = uniform_2d_sample[1] # velocities\n",
      "\t    #inverse_gaussian_cdf_in_place(electrons[1,remaining_indices]) # velocities; doesn't work with dimple\n",
      "\t    electrons[1,remaining_indices] *= v_th_e # velocities\n",
      "\t    electrons[1,remaining_indices] += v_d_e # velocities\n",
      "\t    #electrons[2,remaining_indices] = np.ones(n_samples) # active slots\n",
      "\t    #electrons[2,remaining_indices] = np.ones(n_samples,dtype=np.float32) # relative weights\n",
      "\tif smooth_electron_velocities:\n",
      "\t    # TODO: allow quasi-random numbers here?\n",
      "\t    electron_v_smoothing_scale = v_th_e/10.\n",
      "\t    if read_electron_dist_from_file:\n",
      "\t\telectron_v_smoothing_scale = electron_v_edges[1]-electron_v_edges[0]\n",
      "\t    electrons[1,remaining_indices] += np.random.randn(n_samples)*electron_v_smoothing_scale\n",
      "\tif project_read_particles:\n",
      "\t    if uniform_2d_sampler.shared_seed:\n",
      "\t\tfor id in range(n_engines):\n",
      "\t\t    sample = np.asarray(uniform_2d_sampler.get(n_samples)).T\n",
      "\t\t    if id==mpi_id:\n",
      "\t\t\tuniform_2d_sample = sample\n",
      "\t    elif quiet_start_and_injection:\n",
      "\t\tuniform_2d_sample = uniform_2d_sampler.get_uniform(n_samples,number_of_uniforming_bins).T\n",
      "\t    else:\n",
      "\t\tuniform_2d_sample = uniform_2d_sampler.get(n_samples).T\n",
      "\t    perpendicular_velocities = norm.ppf(uniform_2d_sample[1])*v_th_e\n",
      "\t    velocity_magnitudes = np.sqrt(electrons[1,remaining_indices]*electrons[1,remaining_indices]\n",
      "\t\t\t\t\t  + perpendicular_velocities*perpendicular_velocities)\n",
      "\t    velocity_angles = np.arctan2(perpendicular_velocities, electrons[1,remaining_indices])\n",
      "\t    electrons[1,remaining_indices] = velocity_magnitudes*np.cos(velocity_angles-projection_angle)\n",
      "\tif create_electron_dimple:\n",
      "\t    # TODO: unify random number handling?\n",
      "\t    rejection_samples = np.random.rand(n_samples)\n",
      "\t    remaining_indices = rejection_samples < dimple(electrons[1,remaining_indices])\n",
      "\t    n_successful_samples = n_samples - len(electrons[1,remaining_indices])\n",
      "\t    if mpi_id==0 and n_electrons-n_electrons_so_far>0:\n",
      "\t\tprint 'Still need', n_electrons-n_electrons_so_far, 'electrons...doing another iteration.'\n",
      "\telse:\n",
      "\t    n_successful_samples = n_samples\n",
      "\tn_electrons_so_far += n_successful_samples\n",
      "\tn_samples = n_electrons - n_electrons_so_far\n",
      "    electrons[0][n_electrons:] = inactive_slot_position_flag\n",
      "    # List remaining slots in reverse order to prevent memory fragmentation\n",
      "    empty_electron_slots = -np.ones(electron_storage_length,dtype=np.int32)\n",
      "    current_empty_electron_slot = [(electron_storage_length-n_electrons)-1]\n",
      "    empty_electron_slots[0:(current_empty_electron_slot[0]+1)] = range(electron_storage_length-1,n_electrons-1,-1)\n",
      "    electron_hist_n, electron_hist_n_edges = np.histogram(electrons[0][0:n_electrons], bins=electron_hist_n_edges)\n",
      "    electron_hist_v, electron_hist_v_edges = np.histogram(electrons[1][0:n_electrons], bins=electron_hist_v_edges)\n",
      "    comm.Allreduce(MPI.IN_PLACE, electron_hist_n, op=MPI.SUM)\n",
      "    comm.Allreduce(MPI.IN_PLACE, electron_hist_v, op=MPI.SUM)\n",
      "    if (mpi_id==0):\n",
      "\tfig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,2))\n",
      "\tfor ax, data, bins in zip(axes,[electron_hist_n,electron_hist_v], [electron_hist_n_edges, electron_hist_v_edges]):\n",
      "\t    ax.step(bins,np.append(data,[0]), where='post')\n",
      "\t    filename='initialElectronDistribution.png'\n",
      "\t    plt.savefig(filename)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#%%px\n",
      "#for point, j in zip(grid, range(1+0*n_cells)):\n",
      "#    ion_in_cell = np.logical_and(point<ions[0][:],ions[0][:]<=point+dz)\n",
      "#    ions_in_cell = ions.T[ion_in_cell]\n",
      "#    ions_in_cell = ions_in_cell.T\n",
      "#    average_ion_velocity_in_cell = np.average(ions_in_cell[1])\n",
      "#    electron_in_cell = np.logical_and(point<electrons[0][:],electrons[0][:]<=point+dz)\n",
      "#    electrons_in_cell = electrons.T[electron_in_cell]\n",
      "#    electrons_in_cell = electrons_in_cell.T\n",
      "#    average_electron_velocity_in_cell = np.average(electrons_in_cell[1])\n",
      "#    mpi_writeable = np.zeros(1)\n",
      "#    comm.Allreduce(average_ion_velocity_in_cell, mpi_writeable, op=MPI.SUM)\n",
      "#    average_ion_velocity_in_cell = mpi_writeable[0]/n_engines\n",
      "#    mpi_writeable = np.zeros(1)\n",
      "#    comm.Allreduce(average_electron_velocity_in_cell, mpi_writeable, op=MPI.SUM)\n",
      "#    average_electron_velocity_in_cell = mpi_writeable[0]/n_engines\n",
      "#    ions[1][ion_in_cell] -= average_ion_velocity_in_cell\n",
      "#    electrons[1][electron_in_cell] -= average_electron_velocity_in_cell"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "v_drift = 0.125*v_th_i\n",
      "pot_transp_elong = 2.\n",
      "object_radius = 1.\n",
      "t_object_center = (1.+4.*pot_transp_elong*debye_length)/v_drift\n",
      "t = 0.\n",
      "if boltzmann_electrons:\n",
      "    dt = 0.05*debye_length/v_th_i\n",
      "else:\n",
      "    dt = 0.1*debye_length/v_th_e\n",
      "ion_charge_to_mass = 1\n",
      "if (mpi_id==0):\n",
      "    print t_object_center, dt"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[stdout:0] 16.0 0.000291725017506\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "object_mask = np.zeros_like(grid) # Could make 1 shorter\n",
      "object_center_mask = np.zeros_like(grid) # Could make 1 shorter\n",
      "if include_object:\n",
      "    circular_cross_section(grid,1.e-8,1.,1.,1.,object_center_mask)\n",
      "potential = np.zeros_like(grid)\n",
      "if decouple_ions_from_electrons:\n",
      "    potential_electrons = np.zeros_like(grid)\n",
      "else:\n",
      "    potential_electrons = potential\n",
      "previous_potential = np.zeros_like(grid)\n",
      "electric_field = np.zeros_like(grid)\n",
      "\n",
      "if (mpi_id==0):\n",
      "    object_masks = []\n",
      "    potentials = []\n",
      "    potentials_electrons = []\n",
      "    electric_fields = []\n",
      "    ion_densities = []\n",
      "    electron_densities = []\n",
      "    charge_derivatives = []\n",
      "    times = []\n",
      "    ion_distribution_functions = []\n",
      "    electron_distribution_functions = []"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "ion_density = np.zeros_like(grid)\n",
      "accumulate_density(grid, object_mask, background_ion_density, largest_ion_index, ions, ion_density, \\\n",
      "\t\t       empty_ion_slots, current_empty_ion_slot, periodic_particles=periodic_particles, \\\n",
      "\t\t       use_pure_c_version=use_pure_c_mover)\n",
      "electron_density = np.zeros_like(grid)\n",
      "charge_density = np.zeros_like(grid)\n",
      "previous_charge_density = np.zeros_like(grid)\n",
      "charge_derivative = np.zeros_like(grid)\n",
      "if not boltzmann_electrons:\n",
      "    accumulate_density(grid, object_mask, background_electron_density, largest_electron_index, electrons, \\\n",
      "\t\t\t   electron_density, empty_electron_slots, current_empty_electron_slot, \\\n",
      "\t\t\t   periodic_particles=periodic_particles, use_pure_c_version=use_pure_c_mover)\n",
      "\n",
      "initialize_mover(grid, object_mask, potential, dt, ion_charge_to_mass, largest_ion_index, ions, \\\n",
      "\t\t     empty_ion_slots, current_empty_ion_slot, periodic_particles=periodic_particles, \\\n",
      "\t\t     use_pure_c_version=use_pure_c_mover)\n",
      "if not boltzmann_electrons:\n",
      "    electron_charge_to_mass = -1./mass_ratio\n",
      "    initialize_mover(grid, object_mask, potential_electrons, dt, electron_charge_to_mass, largest_electron_index, electrons, \\\n",
      "\t\t\t empty_electron_slots, current_empty_electron_slot, \\\n",
      "\t\t\t periodic_particles=periodic_particles, use_pure_c_version=use_pure_c_mover)\n",
      "injection_numbers = np.zeros(n_engines,dtype=np.int32)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "n_steps = 10\n",
      "storage_step = 1\n",
      "store_all_until_step = 100\n",
      "print_step = 1\n",
      "#mass_correction = math.sqrt(1./1836./mass_ratio)\n",
      "#debye_correction = 0.04/debye_length\n",
      "#n_steps = int(150000*mass_correction*debye_correction/8.)\n",
      "#storage_step = int(max(1,n_steps/100/mass_correction))\n",
      "#store_all_until_step = 0\n",
      "#print_step = storage_step\n",
      "damping_start_step = 1 # don't make zero to avoid large initial derivative\n",
      "damping_end_step = 0 # make <= damping_start_step to disable damping\n",
      "if (mpi_id==0):\n",
      "    print n_steps, n_steps*dt\n",
      "for k in range(n_steps):\n",
      "    if simulate_moon:\n",
      "\tif k==1: # simulate moon by knocking out particles\n",
      "\t    # TODO: fix half-step in velocity?\n",
      "\t    object_mask[0.45*n_cells:0.55*n_cells] = 1.\n",
      "\telif k==2:\n",
      "\t    object_mask = np.zeros_like(grid)\n",
      "    comm.Allreduce(MPI.IN_PLACE, ion_density, op=MPI.SUM)\n",
      "    if periodic_particles:\n",
      "\tion_density[0] += ion_density[-1]\n",
      "\tion_density[-1] = ion_density[0]\n",
      "    else:\n",
      "\tion_density[0] *= 2. # half-cell at ends\n",
      "\tion_density[-1] *= 2. # half-cell at ends\n",
      "    if boltzmann_electrons:\n",
      "\telectron_density = np.exp(potential) # TODO: add electron temp. dep.\n",
      "    else:\n",
      "\tcomm.Allreduce(MPI.IN_PLACE, electron_density, op=MPI.SUM)\n",
      "\tif periodic_particles:\n",
      "\t    electron_density[0] += electron_density[-1]\n",
      "\t    electron_density[-1] = electron_density[0]\n",
      "\telse:\n",
      "\t    electron_density[0] *= 2. # half-cell at ends\n",
      "\t    electron_density[-1] *= 2. # half-cell at ends\n",
      "    if include_object:\n",
      "\tdist_to_obj = circular_cross_section(grid, t, t_object_center, v_drift, object_radius, object_mask)\n",
      "#    elif simulate_moon and t<=2*t_moon_center:\n",
      "#\tdist_to_obj = circular_cross_section(grid, t, t_moon_center, v_drift_moon, moon_radius, object_mask)\n",
      "#\tdist_to_obj = 100. # Don't include in potential\n",
      "    else:\n",
      "\tdist_to_obj = 100.\n",
      "    object_potential = -3.\n",
      "    if prescribe_potential:\n",
      "\tpotential = prescribed_potential(grid,t)\n",
      "    elif quasineutral:\n",
      "\tpotential = np.log(np.maximum(ion_density,np.exp(object_potential)*np.ones_like(ion_density))) # TODO: add electron temp. dep.\n",
      "\tpotential[0] = 0.\n",
      "\tpotential[-1] = 0.\n",
      "    elif solve_for_electric_field:\n",
      "        np.subtract(ion_density,electron_density,out=charge_density)\n",
      "        gauss_solve(grid, charge_density, debye_length, electric_field, periodic_electric_field=periodic_potential)\n",
      "        np.cumsum(-electric_field*dz, out=potential)\n",
      "    else:\n",
      "\tif (dist_to_obj>0.):\n",
      "\t    if include_object:\n",
      "\t\tfraction_of_obj_pot = math.exp(-dist_to_obj/debye_length)\n",
      "\t    else:\n",
      "\t\tfraction_of_obj_pot = 0.\n",
      "\t    object_potential *= math.exp(-dist_to_obj/(pot_transp_elong*debye_length))\n",
      "\t    np.subtract(ion_density,electron_density,out=charge_density)\n",
      "\t    np.subtract(charge_density,previous_charge_density,out=charge_derivative)\n",
      "\t    charge_derivative /= dt\n",
      "\t    if decouple_ions_from_electrons:\n",
      "\t\tpoisson_solve(grid, object_center_mask, charge_density, \\\n",
      "\t\t\t\t  debye_length, potential_electrons, \\\n",
      "\t\t\t\t  object_potential=object_potential, object_transparency=(1.-fraction_of_obj_pot), \\\n",
      "\t\t\t\t  boltzmann_electrons=False, periodic_potential=periodic_potential, \\\n",
      "\t\t\t\t  use_pure_c_version=use_pure_c_solver)\n",
      "\t    np.copyto(previous_charge_density,charge_density)\n",
      "\t    if boltzmann_potential:\n",
      "\t\tnp.copyto(charge_density,ion_density) # TODO: could just make reference\n",
      "\t\tmax_potential_iter = 20\n",
      "\t    else:\n",
      "\t\tmax_potential_iter = 1\n",
      "\t    if k>=damping_start_step and k<damping_end_step:\n",
      "\t\tdamping_factor = dt/2.\n",
      "\t\tcharge_derivative *= damping_factor\n",
      "\t\tcharge_density += charge_derivative\n",
      "\t    else:\n",
      "\t\tdamping_factor = 0.\n",
      "\t    potential_iter = 0\n",
      "\t    potential_converged = False\n",
      "\t    potential_convergence_threshold = 0.0001\n",
      "\t    while (not potential_converged and potential_iter<max_potential_iter):\n",
      "\t\tpoisson_solve(grid, object_center_mask, charge_density, \\\n",
      "\t\t\t\t  debye_length, potential, \\\n",
      "\t\t\t\t  object_potential=object_potential, object_transparency=(1.-fraction_of_obj_pot), \\\n",
      "\t\t\t\t  boltzmann_electrons=boltzmann_potential, periodic_potential=periodic_potential, \\\n",
      "\t\t\t\t  use_pure_c_version=use_pure_c_solver)\n",
      "\t\tnp.subtract(potential,previous_potential,out=previous_potential)\n",
      "\t\tnp.absolute(previous_potential,out=previous_potential)\n",
      "\t\tmax_potential_change = np.amax(previous_potential)\n",
      "\t\tnp.copyto(previous_potential,potential)\n",
      "\t\tif max_potential_change<potential_convergence_threshold:\n",
      "\t\t    potential_converged = True\n",
      "\t\tpotential_iter += 1\n",
      "#\t    if mpi_id==0:\n",
      "#\t\tprint potential_iter, max_potential_change\n",
      "\telse:\n",
      "\t    fraction_of_obj_pot = 1.\n",
      "\t    np.subtract(ion_density,electron_density,out=charge_density)\n",
      "\t    poisson_solve(grid, object_mask, charge_density, debye_length, potential, \\\n",
      "\t\t\t      object_potential=object_potential, object_transparency=0., \\\n",
      "\t\t\t      use_pure_c_version=use_pure_c_solver)\n",
      "    if not decouple_ions_from_electrons:\n",
      "\tpotential_electrons = potential\n",
      "    if (k%storage_step==0 or k<store_all_until_step):\n",
      "\t    occupied_ion_slots = (ions[0]==ions[0])\n",
      "\t    occupied_ion_slots[empty_ion_slots[0:current_empty_ion_slot[0]+1]] = False\n",
      "\t    n_bins = 100\n",
      "\t    ion_hist_n_edges = np.arange(z_min,z_max+eps,(z_max-z_min)/n_bins)\n",
      "\t    n_bins = 100\n",
      "\t    ion_hist_v_edges = np.arange(v_min_i,v_max_i+eps,(v_max_i-v_min_i)/n_bins)\n",
      "\t    ion_hist2d, ion_hist_n_edges, ion_hist_v_edges = \\\n",
      "\t\tnp.histogram2d(ions[0][occupied_ion_slots], ions[1][occupied_ion_slots], \\\n",
      "\t\t\t\t      bins=[ion_hist_n_edges,ion_hist_v_edges]) # could use range=[[z_min,z_max],[v_min,v_max]]\n",
      "\t    ion_hist2d = np.ascontiguousarray(ion_hist2d)\n",
      "\t    comm.Allreduce(MPI.IN_PLACE, ion_hist2d, op=MPI.SUM)\n",
      "\t    if not boltzmann_electrons:\n",
      "\t\toccupied_electron_slots = (electrons[0]==electrons[0])\n",
      "\t\toccupied_electron_slots[empty_electron_slots[0:current_empty_electron_slot[0]+1]] = False\n",
      "\t\tn_bins = 100\n",
      "\t\telectron_hist_n_edges = np.arange(z_min,z_max+eps,(z_max-z_min)/n_bins)\n",
      "\t\tn_bins = 100\n",
      "\t\telectron_hist_v_edges = np.arange(v_min_e,v_max_e+eps,(v_max_e-v_min_e)/n_bins)\n",
      "\t\telectron_hist2d, electron_hist_n_edges, electron_hist_v_edges = \\\n",
      "\t\t    np.histogram2d(electrons[0][occupied_electron_slots], electrons[1][occupied_electron_slots], \\\n",
      "\t\t\t\t\t  bins=[electron_hist_n_edges,electron_hist_v_edges]) # could use range=[[z_min,z_max],[v_min,v_max]]\n",
      "\t\telectron_hist2d = np.ascontiguousarray(electron_hist2d)\n",
      "\t\tcomm.Allreduce(MPI.IN_PLACE, electron_hist2d, op=MPI.SUM)\n",
      "    if ((k%storage_step==0 or k<store_all_until_step) and mpi_id==0):\n",
      "        times.append(t)\n",
      "        object_masks.append(np.copy(object_mask))\n",
      "        ion_densities.append(np.copy(ion_density))\n",
      "\telectron_densities.append(np.copy(electron_density))\n",
      "\tcharge_derivatives.append(np.copy(charge_derivative))\n",
      "        potentials.append(np.copy(potential))\n",
      "\tpotentials_electrons.append(np.copy(potential_electrons))\n",
      "        electric_fields.append(np.copy(electric_field))\n",
      "        ion_distribution_functions.append(np.copy(ion_hist2d))\n",
      "\tif not boltzmann_electrons:\n",
      "\t    electron_distribution_functions.append(np.copy(electron_hist2d))\n",
      "    if (k%print_step==0 and mpi_id==0):\n",
      "\tprint k\n",
      "    move_particles(grid, object_mask, potential, dt, ion_charge_to_mass, \\\n",
      "\t\t       background_ion_density, largest_ion_index, ions, ion_density, \\\n",
      "\t\t       empty_ion_slots, current_empty_ion_slot, periodic_particles=periodic_particles, \\\n",
      "\t\t       use_pure_c_version=use_pure_c_mover)\n",
      "    if not boltzmann_electrons:\n",
      "\tmove_particles(grid, object_mask, potential_electrons, dt, electron_charge_to_mass, \\\n",
      "\t\t\t   background_electron_density, largest_electron_index, \\\n",
      "\t\t\t   electrons, electron_density, empty_electron_slots, current_empty_electron_slot, \\\n",
      "\t\t\t   periodic_particles=periodic_particles, use_pure_c_version=use_pure_c_mover)\n",
      "    expected_ion_injection = 2*dt*v_th_i/math.sqrt(2*math.pi)*n_ions/(z_max-z_min)\n",
      "    n_ions_inject = int(expected_ion_injection)\n",
      "    if not boltzmann_electrons:\n",
      "\texpected_electron_injection = 2*dt*v_th_e/math.sqrt(2*math.pi)*n_electrons/(z_max-z_min)\n",
      "\tn_electrons_inject = int(expected_electron_injection)\n",
      "    # If expected injection number is small, need to add randomness to get right average rate\n",
      "    # TODO: unify random number usage with sampler\n",
      "    if (expected_ion_injection-n_ions_inject)>np.random.rand():\n",
      "\tn_ions_inject += 1\n",
      "    #injection_numbers[:] = 0\n",
      "    #injection_numbers[mpi_id] = n_ions_inject\n",
      "    #comm.Allreduce(MPI.IN_PLACE, injection_numbers, op=MPI.SUM)\n",
      "    injection_number = np.array(n_ions_inject,dtype=np.int32)\n",
      "    comm.Allgather(injection_number, injection_numbers)\n",
      "    if injection_sampler.shared_seed:\n",
      "\tfor injection_number in injection_numbers[:mpi_id]:\n",
      "\t    sample = np.asarray(injection_sampler.get(int(injection_number))).T\n",
      "    if n_ions_inject>0 and not periodic_particles:\n",
      "\tinject_particles(n_ions_inject, grid, dt, v_th_i, background_ion_density, \\\n",
      "\t\t\t     injection_sampler, ions, empty_ion_slots, \\\n",
      "\t\t\t     current_empty_ion_slot, largest_ion_index, ion_density)\n",
      "    if injection_sampler.shared_seed:\n",
      "\tfor injection_number in injection_numbers[mpi_id+1:]:\n",
      "\t    sample = np.asarray(injection_sampler.get(int(injection_number))).T\n",
      "    if not boltzmann_electrons:\n",
      "\tif (expected_electron_injection-n_electrons_inject)>np.random.rand():\n",
      "\t    n_electrons_inject += 1\n",
      "\t#injection_numbers[:] = 0\n",
      "        #injection_numbers[mpi_id] = n_electrons_inject\n",
      "        #comm.Allreduce(MPI.IN_PLACE, injection_numbers, op=MPI.SUM)\n",
      "\tinjection_number = np.array(n_electrons_inject,dtype=np.int32)\n",
      "        comm.Allgather(injection_number, injection_numbers)\n",
      "\tif injection_sampler.shared_seed:\n",
      "\t    for injection_number in injection_numbers[:mpi_id]:\n",
      "\t\tsample = np.asarray(injection_sampler.get(int(injection_number))).T\n",
      "\tif n_electrons_inject>0 and not periodic_particles:\n",
      "\t    inject_particles(n_electrons_inject, grid, dt, v_th_e, background_electron_density, \\\n",
      "\t\t\t\t injection_sampler, electrons, empty_electron_slots, \\\n",
      "\t\t\t\t current_empty_electron_slot, largest_electron_index, electron_density)\n",
      "\tif injection_sampler.shared_seed:\n",
      "\t    for injection_number in injection_numbers[mpi_id+1:]:\n",
      "\t\tsample = np.asarray(injection_sampler.get(int(injection_number))).T\n",
      "    t += dt\n",
      "if (mpi_id==0):\n",
      "    print times[0], dt, times[len(times)-1]"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "ERROR: Cell magic `%%px` not found.\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%px\n",
      "if (mpi_id==0):\n",
      "    times_np = np.array(times, dtype=np.float32)\n",
      "    object_masks_np = np.array(object_masks, dtype=np.float32)\n",
      "    potentials_np = np.array(potentials, dtype=np.float32)\n",
      "    potentials_electrons_np = np.array(potentials_electrons, dtype=np.float32)\n",
      "    electric_fields_np = np.array(electric_fields, dtype=np.float32)\n",
      "    ion_densities_np = np.array(ion_densities, dtype=np.float32)\n",
      "    electron_densities_np = np.array(electron_densities, dtype=np.float32)\n",
      "    charge_derivatives_np = np.array(charge_derivatives, dtype=np.float32)\n",
      "    ion_distribution_functions_np = np.array(ion_distribution_functions, dtype=np.float32) # actuall int\n",
      "    electron_distribution_functions_np = np.array(electron_distribution_functions, dtype=np.float32) # actuall int\n",
      "    filename_base = \\\n",
      "\t'l'+('%.4f' % debye_length)+'_d'+('%.3f' % v_drift)+'_np'+('%.1e' % n_points)+'_ni'+('%.1e' % n_ions)+'_dt'+('%.1e' % dt)\n",
      "    print filename_base\n",
      "    np.savez(filename_base, grid=grid, times=times_np, object_masks=object_masks_np, potentials=potentials_np, \\\n",
      "                    potentials_electrons=potentials_electrons_np, electric_fields=electric_fields, \\\n",
      "\t\t    ion_densities=ion_densities_np, electron_densities=electron_densities_np, charge_derivatives=charge_derivatives_np, \\\n",
      "\t\t    ion_hist_n_edges=ion_hist_n_edges, ion_hist_v_edges=ion_hist_v_edges, \\\n",
      "\t\t    ion_distribution_functions=ion_distribution_functions_np, \\\n",
      "\t\t    electron_hist_n_edges=electron_hist_n_edges, electron_hist_v_edges=electron_hist_v_edges, \\\n",
      "\t\t    electron_distribution_functions=electron_distribution_functions_np, \\\n",
      "\t\t    background_ion_density=background_ion_density, background_electron_density=background_electron_density)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#####ScriptEnd#####"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.parallel import Client"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rc = Client()\n",
      "dview = rc[:]\n",
      "dview.block = True\n",
      "n_engines = len(rc.ids)\n",
      "res = dview.push(dict(n_engines=n_engines))\n",
      "print n_engines"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mpi_ids = np.array(dview.pull('mpi_id'))\n",
      "master_rc_id = np.arange(0,len(mpi_ids))[mpi_ids==0][0]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filename = dview.pull('filename_base', targets=master_rc_id) + '.npz'\n",
      "data_file = np.load(filename)\n",
      "print data_file.files\n",
      "grid = data_file['grid']\n",
      "times = data_file['times']\n",
      "ion_densities = data_file['ion_densities']\n",
      "electron_densities = data_file['electron_densities']\n",
      "potentials = data_file['potentials']\n",
      "electric_fields = data_file['electric_fields']\n",
      "object_masks = data_file['object_masks']"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "k = len(times)-1\n",
      "print times[k]\n",
      "fig, axes = plt.subplots(nrows=3,ncols=2,figsize=(8,6))\n",
      "for ax, data in zip(axes.flatten(),[potentials[k], ion_densities[k]-electron_densities[k], \\\n",
      "                              object_masks[k], object_masks[k], ion_densities[k], electron_densities[k]]):\n",
      "    ax.plot(grid,data)\n",
      "filename='figures/data.png'\n",
      "plt.savefig(filename)\n",
      "IPdisp.Image(filename=filename)"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 8 engines: 14.4s\n",
      "# 16 engines: 18.6s\n",
      "# 32 engines: 29.9s\n",
      "# 1 eng, 10x part: 1m58\n",
      "# 2 eng, 10x part: 1m54\n",
      "# 4 eng, 10x part: 2m31\n",
      "# 8 eng, 10x part: 2m8; 2m6\n",
      "# 15 eng, 5x part: 1m10\n",
      "# 16 eng, 5x part: 1m17; 1m17\n",
      "# 16 eng, 5x part, 10x grid: 1m18\n",
      "# 16 eng, 5x steps: 1m31\n",
      "# 16 eng, 10x grid, 5x steps: 1m37\n",
      "# 16 eng, 5x part, 10x grid, 100x steps: 2h7"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from guppy import hpy\n",
      "h = hpy()\n",
      "print h.heap()"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Partition of a set of 96715809 objects. Total size = 20734669600 bytes.\n",
        " Index  Count   %     Size   % Cumulative  % Kind (class / dict of class)\n",
        "     0 4160324   4 7621713568  37 7621713568  37 IPython.parallel.client.client.Metadata\n",
        "     1 31683369  33 5238958352  25 12860671920  62 unicode\n",
        "     2 13766906  14 3856539632  19 16717211552  81 dict (no owner)\n",
        "     3 15684278  16 1198299824   6 17915511376  86 list\n",
        "     4 3200096   3 896026880   4 18811538256  91 dict of IPython.parallel.client.client.ExecuteReply\n",
        "     5 16641333  17 798783984   4 19610322240  95 datetime.datetime\n",
        "     6      3   0 402654816   2 20012977056  97 collections.defaultdict\n",
        "     7 4344426   4 352809984   2 20365787040  98 str\n",
        "     8 3200096   3 204806144   1 20570593184  99 IPython.parallel.client.client.ExecuteReply\n",
        "     9 3202772   3 76866528   0 20647459712 100 int\n",
        "<1435 more rows. Type e.g. '_.more' to view.>\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import gc"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gc.collect()"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "2082"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,2))\n",
      "for ax, data in zip(axes,electrons):\n",
      "    n_bins = n_points;\n",
      "    occupied_slots = (data==data)\n",
      "    occupied_slots[empty_electron_slots[0:current_empty_electron_slot[0]+1]] = False\n",
      "    ax.hist(data[occupied_slots],bins=n_bins, histtype='step')\n",
      "filename='data.png'\n",
      "plt.savefig(filename)\n",
      "IPdisp.Image(filename=filename)"
     ],
     "language": "python",
     "outputs": [
      {
       "ename": "NameError",
       "evalue": "name 'electrons' is not defined",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-23-37fe8d52c6e1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mfig\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnrows\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mncols\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m8\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[1;32mfor\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxes\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0melectrons\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m     \u001b[0mn_bins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mn_points\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0moccupied_slots\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m==\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0moccupied_slots\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mempty_electron_slots\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mcurrent_empty_electron_slot\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mFalse\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mNameError\u001b[0m: name 'electrons' is not defined"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(nrows=1,ncols=2,figsize=(8,2))\n",
      "for ax, data in zip(axes,ions):\n",
      "    n_bins = n_points;\n",
      "    occupied_slots = (data==data)\n",
      "    occupied_slots[empty_ion_slots[0:current_empty_ion_slot[0]+1]] = False\n",
      "    ax.hist(data[occupied_slots],bins=n_bins, histtype='step')\n",
      "filename='data.png'\n",
      "plt.savefig(filename)\n",
      "IPdisp.Image(filename=filename)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAADICAYAAAAQj4UaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVGX+B/DPDHIZB1QYUAZJLRVDk0W85CVN8YYLdNEB\noowikRZRNFfJ9Le6RmKlJhnlhdzMJDZBUxPT0szLaqIgRWkrXTYFR4UQ5TI4yDy/P45zZoYZ5DYz\nB4bv+/WaF8888z3Pec4AZ+Y557mIGGMMhBBCCCGEEGIFYqErQAghhBBCCOk4qAFCCCGEEEIIsRpq\ngBBCCCGEEEKshhoghBBCCCGEEKuhBgghhBBCCCHEaqgBQgghhBBCCLEaaoAQQgghhBBCrIYaIIQQ\nQgghhBCroQYIIYQQQgghxGqoAUIIIYQQQgixGmqAEEIIIYQQQqyGGiCEEEIIIYQQq7HpBkhubi6C\ngoLQtWtXdOnSBVOnTsX3339vMvbixYsICgqCi4sLZDIZoqKiUFpaajJ269at8PX1hUQigY+PD1JT\nU03GlZeXIzY2Fh4eHnB2dkZgYCDOnz9vtuMjhJCWOHv2LObOnYtBgwbB2dkZvXv3RkREBAoLC41i\nhT43njp1Co899hikUinkcjnmz5+Pqqqqlh88IYQQwYkYY0zoSlhCXl4exowZg969e+Pll19GXV0d\nPvjgA5SVlSEnJwc+Pj58bFFREYYMGQJXV1ckJCSgoqICa9euRa9evZCTkwN7e3s+dvPmzYiLi4NC\nocDUqVNx/PhxfPLJJ3jzzTeRmJjIx2k0GowdOxY//PADEhMTIZPJ8MEHH+DKlSvIzc1Fv379rPp+\nEEKIlkKhwOnTpxEWFgY/Pz8olUqkpqaisrIS3333HQYNGgRA+HNjfn4+Ro0ahUGDBiE2NhZXrlzB\n2rVrMWHCBBw4cMB6bxghhBDzYjbqr3/9K5PJZKysrIzPUyqVzMXFhc2YMcMgNi4ujkmlUnblyhU+\n7/Dhw0wkErEtW7bwedXV1Uwmk7HQ0FCD7WfOnMmcnZ3ZzZs3+bzPPvuMiUQitmvXLj6vpKSEubq6\nsmeffdZsx0kIIc116tQpVltba5BXWFjInJyc2MyZM/k8oc+N06ZNYz179mQVFRV83ocffshEIhH7\n6quvWnj0hBBChGazDRAXFxcWERFhlB8cHMwcHR1ZZWUln9e9e3eTsQMGDGCTJk3in2dnZzORSMS+\n/PJLg7jTp08zkUjEduzYweeFhYUxuVxuVObLL7/MpFIpU6vVLTouQgixlICAADZs2DD+uZDnxlu3\nbjF7e3v26quvGsSp1Wrm4uLCYmJiWnaQhBBCBGezY0DUajUkEolRfufOnaFWq/Hjjz8CAIqLi1FS\nUoJhw4YZxQ4fPtygX7I2XT82ICAAYrEY+fn5BrEBAQEmy6yursalS5dadmCEEGIBjDFcv34d7u7u\nAIQ/NxYUFODu3btGZdrb28Pf35/G0xFCSDtmsw2QAQMG4PTp09BoNHyeWq3GmTNnAABXr14FACiV\nSgCAXC43KkMul6OsrAy1tbV8rJ2dHf8BreXg4ACZTMaXqY1tqEz9/RNCSFuQnp6Oq1evIiIiAoDw\n58b77d/T05POoYQQ0o51EroCljJnzhzExcVh1qxZSExMRF1dHd544w1cu3YNAKBSqQx+Ojo6GpXh\n5OTEx9jb20OlUsHBwcHk/hwdHfmyAKCmpqbRMusrLS3FoUOH0KdPH5N3bwghpLVUKhX+97//YerU\nqXyD4eeff0Z8fDxGjx6NF154gY8DhDs3NrZ/U+dQgM6jhBDLM3UeJc1jsw2Ql19+GVeuXMGaNWvw\n8ccfA+Bu8ScmJmLVqlVwdnYGAP4D6s6dO0Zl1NTUGMRIJBKo1WqT+6upqTH4sJNIJE0qU9+hQ4cw\nc+bMJh8jIYS01I4dO/Dcc8/h2rVrCA4OhqurK7KysiASiQAIf25sbP+dO3c2uT86jxJCrEV7HiXN\nZ7MNEAB44403sGjRIly4cAFdu3bFoEGDsHTpUgDgp+HV3t7X3u7Xp1QqIZPJ+Kkm5XI56urqUFpa\natDiVavVKCsrg5eXF58nl8tNdhHQ7kc/VuvBBx8EwP1B+/r6tuiYhbRgwQKkpKQIXY0WoboLg+pu\nfRcvXsTMmTPx4IMP4tatW5g2bRpu376NEydOwNPTk48T+tzY2P5NnUMBOo8KieouDKq79emfR0nL\n2HQDBAC6deuG0aNH888PHz6MBx54AA8//DAAoGfPnvDw8MDZs2eNts3JyYG/vz//fMiQIQC4Rbym\nTZvG5587dw4ajcYg1t/fHydOnABjjL+iCABnzpyBVCo1WIdES9sFwdfX1+QgzbauW7du7bLeANVd\nKFR34YhEIoSGhuKXX37B4cOH+XOiltDnxkceeQSdOnXC2bNnoVAo+Di1Wo38/Hw888wzJo+LzqPC\noboLg+ouHO35hjSfzQ5CN+Wzzz7DuXPnsGDBAoP8GTNmYP/+/SgqKuLzjhw5gsLCQoSFhfF5gYGB\ncHNzw8aNGw2237hxI6RSKYKDg/k8hUKB69evY/fu3XxeaWkpMjMzERoaarCAFyGEWNtrr72GM2fO\nIDMzE48++qjJGCHPjV27dsWkSZOwY8cOVFZW8rGffPIJqqqqDPZPCCGkfbHZOyDHjx/H66+/jqlT\np8LNzQ3fffcdtm3bhmnTpmH+/PkGsUuXLkVmZiYmTJiA+fPno6KiAmvWrIGfnx+io6P5OCcnJyQl\nJSE+Ph7h4eGYMmUKTpw4gfT0dCQnJ6Nbt258rEKhwMiRIxEdHY0LFy7wq/0yxrBy5UqrvQ+EEGLK\n8ePHERoaitLSUuzYscPgNe0YCqHPjatWrcLo0aPx+OOPY/bs2SgqKsI777yDqVOnYsqUKRZ8dwgh\nhFiUoKuQWNCvv/7Kpk6dyjw8PJiTkxMbOHAge+utt4xW/9X66aef2NSpU5lUKmVubm7s+eefZzdu\n3DAZm5aWxh5++GHm6OjI+vfvz959912TcTdv3mQxMTHM3d2dSaVSNmHCBJabm9tgnXNzcxmA+8a0\nZfVXQW5PqO7CoLpbn/Y8IxaLmUgkMnqIxWKDeKHPjSdPnmRjxoxhEomE9ejRg82bN89gIdmGjo/O\no9ZHdRcG1d362vt5pi2w2TsgDz30EA4ePNjk+IEDBzY5PiYmBjExMY3GdevWDWlpaUhLS2tyPdqz\nyMhIoavQYlR3YVDdhXP27Nkm9b0W+tw4ZswYnDx5skn7twXt+e+K6i4Mqjtpj0SMMSZ0JQgnLy8P\nQ4cORW5ubrselEUIabts/Txj68dHCBEenWdar0MNQieEEEIIIYQIixoghBBCCCGEEKuhBgghhLQh\nx48DNMETIYQQW2bTDZBz587hySefhJeXF6RSKXx9fZGUlASVSmUQd/HiRQQFBcHFxQUymQxRUVEo\nLS01WebWrVvh6+sLiUQCHx8fpKammowrLy9HbGwsPDw84OzsjMDAQJw/f97sx0gIsS0JCcDXXwtd\nC0IIIcRybHYWrIKCAjz22GPw8vLCggUL4ObmhlOnTmHFihXIzc3Fnj17AABFRUUYN24cXF1dsXr1\nalRUVGDt2rUoKChATk6OwYKBmzdvRlxcHBQKBRYtWoTjx48jISEB1dXVSExM5OM0Gg2Cg4Pxww8/\nIDExkZ/nfvz48cjNzUW/fv2s/n4QQgghhBDSFthsA+Tf//431Go1srOz4evrC4CbIlKj0WD79u24\ndesWunbtiuTkZKhUKpw/fx7e3t4AgBEjRmDy5MnYtm0bZs+eDQBQqVRYtmwZQkJCsHPnTgDArFmz\noNFokJSUhNjYWH6xraysLJw+fRpZWVmYPn06ACA8PBw+Pj5YsWIF0tPTrf12EEIIIaSeO3eATp0A\nOzuha0JIx2KzXbAkEgkAoHv37gb5np6esLOzg4ODAwBg165dCAkJ4RsfADBx4kT4+PjwDQ0AOHr0\nKMrKyjBnzhyD8uLj41FVVYXs7Gw+LysrC56ennzjAwDc3d0RHh6OvXv3ora21nwHSgghhJAWcXIC\nJBIgKUnomhDSsdhsA+Sll15Cjx49MGvWLHz//fe4cuUKPvvsM2zatAkJCQmQSCQoLi5GSUkJhg0b\nZrT98OHDDcZsaNP1YwMCAiAWi5Gfn28Qa2pe6OHDh6O6uhqXLl0y12ESQgghpBVqa4Hly4GZM4E+\nfYAbN4SuESG2z2YbIF5eXvjPf/6Dn3/+GUOGDEHv3r0RGRmJhIQErFu3DgCgVCoBAHK53Gh7uVyO\nsrIy/m6FUqmEnZ0d3N3dDeIcHBwgk8lw9epVPk+pVDZYJgCDWEII0ScSCV0DQmxbbS3w7LPG/2vp\n6cAff3CNkD//FKRqhHQYNjsG5Pr165g2bRoAIC0tDTKZDPv378eqVavQo0cPxMfH87NhOTo6Gm3v\n5OQEgBv7YW9vD5VKxXfbqs/R0dFgZq2amppGyySEEEKIdQUGAkeP3j9GpQLc3QEvL6C42Dr1IqSj\nsdkGSFJSEoqLi3Hp0iV4eXkBAJ566iloNBq8+uqriIyM5MeJ3Llzx2j7mpoaALqxJBKJBGq12uS+\nampq+DhtbFPKbMiCBQv4Ae1akZGRiIyMvO92hBCiLyMjAxkZGQZ55eXlAtWGEGHduNF440Pf1avA\n3LlAA7PtE0JawWYbICdPnsSQIUP4xodWaGgotm3bhvz8fAwYMACAriuWPqVSCZlMxk/DK5fLUVdX\nh9LSUoNuWGq1GmVlZQb7kcvlJrtZafdTv071paSkmBxDQgghzWHqwkVeXh6GDh0qUI0IEc4//9n8\nbd5/nxoghFiCzY4Bqa2tRV1dncl8ALh79y569uwJDw8PnD171iguJycH/v7+/PMhQ4YAgFHsuXPn\noNFoDGL9/f2Rl5cHxphB7JkzZyCVSuHj49PyA2shpZK79XztmtV3TQghhAhu48aGX3vyyYZfq6gw\nf10I6ehstgESEBCAvLw8FBYWGuRnZGTAzs4Ofn5+AIAZM2Zg//79KCoq4mOOHDmCwsJChIWF8XmB\ngYFwc3PDxnpnsI0bN0IqlSI4OJjPUygUuH79Onbv3s3nlZaWIjMzE6GhoQaLG1rLF19wt54PHLD6\nrgkhzUCD0AmxHj8/4PffgT17AMYAvY9yXpcu1q8XIbbOZrtgLV68GLt27cLYsWMxd+5cuLm5Yf/+\n/Th48CBmz54NT09PAMDSpUuRmZmJCRMmYP78+aioqMCaNWvg5+eH6OhovjwnJyckJSUhPj4e4eHh\nmDJlCk6cOIH09HQkJycbjNlQKBQYOXIkoqOjceHCBX4ldMYYVq5cafX3ghBCCOnI9K4xon9/YPp0\n4B//AKRSw7jMTGD7duBvf7Nu/QjpaGy2AeLn54dvv/0WK1aswJo1a1BTU4OHHnoIycnJSExM5OO8\nvb1x7NgxLFy4EEuWLIGjoyNCQkKwbt06ozsVcXFxsLe3x7p167Bv3z706tULKSkpSEhIMIgTi8U4\ncOAAFi9ejA0bNkClUmHEiBHYvn07+vfvb5XjJ4QQQgjngQd06fstxSWRAC+/DMTGAmK9PiKXLwO9\nelmufoR0NDbbAAGAESNG4Msvv2w0buDAgTh48GCTyoyJiUFMTEyjcd26dUNaWhrS0tKaVC4hhBBC\n2gaRCKiuBjp35p4PHw5cvy5snQixJTY7BoQQQggh5N7cMwCAFSuavp3+jPm0Ojoh5kUNEEIIaUNo\nEDoh5vX667r0P/7R8nLqTWxJCGkFaoAQQgghxGa98YYubWfXvG0XL9alf/nFPPUhhFADhBBCCCE2\nSv+uxYwZzd9+9WpdWr8rFyGkdWy2AfLiiy9CLBY3+NBf/fzixYsICgqCi4sLZDIZoqKiUFpaarLc\nrVu3wtfXFxKJBD4+PkhtYInU8vJyxMbGwsPDA87OzggMDMT58+ctcqyEEEIIMaa/HvHcuc3fXv+O\nyaBBra8PIYRjs7Ng/e1vf8OUKVMM8jQaDf72t7/hwQcfhFwuBwAUFRVh3LhxcHV1xerVq1FRUYG1\na9eioKAAOTk5BlPxbt68GXFxcVAoFFi0aBGOHz+OhIQEVFdXG0ztq9FoEBwcjB9++AGJiYn8OiDj\nx49Hbm4u+vXrZ503gRBCCOnAQkJ06fHjBasGIaQem22AjBw5EiNHjjTIO3nyJKqrq/Hcc8/xecnJ\nyVCpVDh//jy8vb0BcNP3Tp48Gdu2bcPs2bMBACqVCsuWLUNISAh27twJAJg1axY0Gg2SkpIQGxvL\nL0aYlZWF06dPIysrC9OnTwcAhIeHw8fHBytWrEB6errFj58Q0j7RIHRCzOfQodaXsXUrMGsWl2aM\n/kcJMQeb7YJlyqeffgqRSIRnn32Wz9u1axdCQkL4xgcATJw4ET4+PnxDAwCOHj2KsrIyzJkzx6DM\n+Ph4VFVVITs7m8/LysqCp6cn3/gAAHd3d4SHh2Pv3r2oFbAjKc3iQQghpCNQq3XpzZtbXs6ECbp0\nZmbLyyGE6HSYBkhtbS127tyJMWPGoNe95UyLi4tRUlKCYcOGGcUPHz7cYMyGNl0/NiAgAGKxGPn5\n+QaxAQEBJsusrq7Gpfstw2ohdMWGEKJPpVJhxYoVCAoKgpubG8RiMT7++GOjuIbG0/n6+pos1xLj\n5E6dOoXHHnsMUqkUcrkc8+fPR1VVVcsPnnQIW7fq0rGxLS/nwQd16dY0ZAghOjbbBau+Q4cOoays\nzKD7lXYgunY8iD65XI6ysjLU1tbC3t4eSqUSdnZ2cHd3N4hzcHCATCbD1atXDcodb6KzqXY/V69e\nxSAazUYIEdDNmzeRlJSE3r17w9/fH99++y1EDVypcHR0xFb9b3MAunbtahRniXFy+fn5mDhxIgYN\nGoT169fjypUrWLt2LQoLC3HgwAEzvRvEFplz5fJJk4DDh4FvvjFfmYR0ZB2mAfLpp5/CwcEB4eHh\nfJ5KpQLAfbjW5+TkxMfY29tDpVLBwcHBZNmOjo58WQBQU1PTaJmEECIkDw8PXLt2Dd27d0dubi6G\nDx/eYKy9vb1B11VTLDVObunSpZDJZPj222/h7OwMAOjTpw9mz56Nr7/+GpMnT27V+0Bs18qV5isr\nKwu49+dLCDGDDtEAqaysxN69ezF16lS4urry+RKJBABw584do21qamoMYiQSCdT6HUrrxWrjtLFN\nKbMhCxYs4D+otSIjIxEZGXnf7Qgh7Z85u0tmZGQgIyPDIK+8vBwA16jo3r07AIA1MjiMMQaNRoPK\nykp06dLFZMz9xsmlp6cjOzubvwN9v3FyO3bs4O883759G4cPH8bChQv5xgcAREVF4ZVXXsHOnTup\nAUJM0uuUYBb6N/xOnADGjjVv+YR0NB2iAbJnzx6oVCqD7leArkuU/pogWkqlEjKZjJ+GVy6Xo66u\nDqWlpQbdsNRqNcrKyuDl5WVQ7lUTZz/tfvRjTUlJSTE5hoQQQprD1IWLvLw8DB06tFnlVFdXo0uX\nLqiuroarqysiIyPx1ltvQSqV8jFNGSenPQffb5zcli1bcOnSJQwaNAgFBQW4e/euUZn29vbw9/en\ntZVIg65c0aWLisxb9owZwI0b5i2TkI6mQwxCT09Ph4uLC5544gmD/J49e8LDwwNnz5412iYnJwf+\n/v788yFDhgCAUey5c+eg0WgMYv39/ZGXl2d0VfHMmTOQSqXw8fFp9TERQog1eHl54dVXX8W2bdvw\n73//G0888QQ++OADBAUFoU5vlbfmjpNraOwdAD72fuP0PD09TV7oIQQANm3SpXv2NG/ZJSXmLY+Q\njsjmGyAlJSU4fPgwnn76aX4Mhr4ZM2Zg//79KNK7RHLkyBEUFhYiLCyMzwsMDISbmxs2btxosP3G\njRshlUoRHBzM5ykUCly/fh27d+/m80pLS5GZmYnQ0FCDxQ0JIaQtS05ORnJyMhQKBcLDw/HRRx9h\n1apV+M9//oOsrCw+zhLj5Bobp0fj6UhDtm0zf5kLFpi/TEI6KpvvgvXZZ5+hrq7OqPuV1tKlS5GZ\nmYkJEyZg/vz5qKiowJo1a+Dn54fo6Gg+zsnJCUlJSYiPj0d4eDimTJmCEydOID09HcnJyQZjNhQK\nBUaOHIno6GhcuHCBn+GFMYaV5hwVRwixWW15wbNXXnkF//jHP3DkyBFEREQAsMw4ucbG6XXu3LnB\nOtJYOgIAejPkt9rbbwMpKVz66lWgkd7UxEbcbywdaTmbb4B8+umn6NGjByZNmmTydW9vbxw7dgwL\nFy7EkiVL4OjoiJCQEKxbt87oTkVcXBzs7e2xbt067Nu3D7169UJKSgoSEhIM4sRiMQ4cOIDFixdj\nw4YNUKlUGDFiBLZv347+/ftb7FgJIe1fW2106HNycoKbmxvKysr4PEuMk2tsnN79xtPRWLqO69df\ndemBA81Xrv5Xgp49aWHfjsJcY+mIIZtvgJw6darRmIEDB+LgwYNNKi8mJgYxMTGNxnXr1g1paWlI\nS0trUrmEENJeVFRUoLS0FB4eHnye/ji5adOm8fkNjZM7ceIEGGMGa4/UHyf3yCOPoFOnTjh79iwU\nCgUfp1arkZ+fj2eeecZix0jar1df1aXt7ISrByGkYTY/BoQQQkjL3LlzBxUVFUb5SUlJAICgoCA+\nzxLj5Lp27YpJkyZhx44dqKys5GM/+eQTVFVVGYzTI0Rr1y5dWmzmbznx8eYtj5COyubvgBBDdMuY\nEKKVmpqK8vJyvjvUvn37cPnyZQBAQkICysrKMGTIEDz77LMYMGAAAODQoUP48ssvMW3aNDz55JN8\nWZYaJ7dq1SqMHj0ajz/+OGbPno2ioiK88847mDp1KqZMmWLpt4i0Y5Zon6akAO+/z6U1GvM3cAjp\nKKgB0kG0h37lhBAdawxCX7duHf744w8AgEgkwueff47du3dDJBIhKioKrq6uCA0Nxddff42PP/4Y\ndXV16N+/P1avXo1FixYZlWeJcXJDhgzB4cOH8eqrr2LhwoXo0qULYmJisHr1asu9MaTdqq3Vpf/+\nd/OX30nvW9NzzwH1xiYTQpqIGiCEENKGWPNiwe+//95ozPbt25tVpiXGyY0ZMwYnT55sVj1Ix7Rj\nhy49YoRl9/Xvf1MDhJCWsvmbh3l5eXjiiScgk8kglUoxePBgvPfeewYxFy9eRFBQEFxcXCCTyRAV\nFYXS0lKT5W3duhW+vr6QSCTw8fFBamqqybjy8nLExsbCw8MDzs7OCAwMpFV7CSGEEAt66SVd2lKN\n+XqzOxNCWsCm74B89dVXCA0NxdChQ7F8+XI4Ozvjl19+QXFxMR9TVFSEcePGwdXVFatXr0ZFRQXW\nrl2LgoIC5OTkGEzFu3nzZsTFxUGhUGDRokU4fvw4EhISUF1djcTERD5Oo9EgODgYP/zwAxITE/n+\nzePHj0dubi769etn1feBEEII6Uj0B6KbW04OcG+iNpSUAHqTwRFCmshmGyC3b99GVFQUQkNDDVbr\nrS85ORkqlQrnz5+Ht7c3AGDEiBGYPHkytm3bhtmzZwPgVuRdtmwZQkJCsHPnTgDArFmzoNFokJSU\nhNjYWH6QZVZWFk6fPo2srCxMnz4dABAeHg4fHx+sWLEC6enpljx0QgghpEN74gnLla0/TOm994DX\nX7fcvgixVTbbBevTTz/FjRs3sGrVKgBAVVUVNBqNUdyuXbsQEhLCNz4AYOLEifDx8eEbGgBw9OhR\nlJWVYc6cOQbbx8fHo6qqCtnZ2XxeVlYWPD09+cYHALi7uyM8PBx79+5Frf4oOUIIMYFmrCOkefbt\n06U7Weny6r0ZqQkhzWSzDZDDhw+jS5cuuHLlCgYMGAAXFxd07doVc+bMwZ07dwAAxcXFKCkpwbBh\nw4y2Hz58uMGYDW26fmxAQADEYjHy8/MNYk2twDt8+HBUV1fj0qVLZjlGQojtoRnrCGmZqCiha0AI\naSqbbYAUFhbi7t27eOqppzBt2jTs3r0bL730EjZt2oTo6GgAgFKpBADI5XKj7eVyOcrKyvi7FUql\nEnZ2dnB3dzeIc3BwgEwm4+fR18Y2VCYAg1hCCCGEtN6tW9bbV9++1tsXIbbIZseAVFZWorq6GnFx\ncUhJSQEAPPXUU1Cr1di8eTNef/11qFQqAICjo6PR9k5OTgC4sR/29vZQqVRwcHAwuS9HR0e+LACo\nqalptMz7WbBggcGiXQAQGRmJyMjI+25HCCH6MjIykFFvntDy8nKBakOI5ej3sL73kW9RBw4A99bm\npAUJCWkBm22ASCQSADD60h4ZGYnNmzfju+++w8MPPwwAfJcsfTU1NQblSCQSqNVqk/uqqanh47Sx\nTSmzISkpKSa7cBFCSHOYunCRl5eHoUOHClQjQixDbxgm4uMtvz/tLFgAtxbIc89Zfp+E2BKbbbN7\neXkBAHr06GGQ3717dwDAzZs3+RhtVyx9SqUSMpmMn4ZXLpejrq7OaH0QtVqNsrIyvixtrKluVtr9\n6McSQogpNAidkKYrK9OlrTUA/dlnuZ/Hj1tnf4TYEpttgGgHixcVFRnkaxsGHh4e8PLygoeHB86e\nPWu0fU5ODvz9/fnnQ4YMAQCj2HPnzkGj0RjE+vv7Iy8vD6zeN4gzZ85AKpXCR//SiZXRlxpC2jYa\nhE5I8734ovX3qZ0Uc8sW6++bkPbOZhsg4eHhALiVy/V9+OGHsLe3x/jx4wEAM2bMwP79+w0aKkeO\nHEFhYSHCwsL4vMDAQLi5uWHjxo0G5W3cuBFSqRTBwcF8nkKhwPXr17F7924+r7S0FJmZmQgNDTVY\n3NBa6EsNIYQQW+fra719PfqoLr1/v/X2S4gtsNkxIP7+/njppZfwr3/9C3fv3sW4cePw7bffIisr\nC0uXLoWnpycAYOnSpcjMzMSECRMwf/58VFRUYM2aNfDz8+NnywK4AeRJSUmIj49HeHg4pkyZghMn\nTiA9PR3JyckGg8YVCgVGjhyJ6OhoXLhwgV8JnTGGlStXWv29IIQQQjqCd9+13r70u3rt2gWEhFhv\n34S0dzbbAAGATZs2oVevXvjoo4/w+eefo0+fPkhJSUFCQgIf4+3tjWPHjmHhwoVYsmQJHB0dERIS\ngnXr1hnLYC0LAAAgAElEQVTdqYiLi4O9vT3WrVuHffv2oVevXkblAYBYLMaBAwewePFibNiwASqV\nCiNGjMD27dvRX38JVUIIIYS0iv70u5MmCVOHbduAjz4SZt+EtEc23QDp1KkTli9fjuXLl983buDA\ngTh48GCTyoyJiUFMTEyjcd26dUNaWhrS0tKaVC4hhOij8VqENM3Onbq0tbsb29kBdXXW3SchtsBm\nx4AQQkh7ROO1CGme2Fju55NPWn/f+ndf6KIBIU1HDRBCCCGEtHuJidbfp1SqS//+u/X3T0h7RQ0Q\nQgghhLRLxcW69OjRwtUDMFyNnRByf9QAIYQQQki7pDcDvuBojhlCms5mGyDffvstxGKxyUdOTo5B\n7MWLFxEUFAQXFxfIZDJERUUZrXiutXXrVvj6+kIikcDHxwepqakm48rLyxEbGwsPDw84OzsjMDAQ\n58+fN/txEkJsE/UnJ6Rx338vdA2AI0eErgEh7Y9Nz4IFAPPnz8fw4cMN8vr27cuni4qKMG7cOLi6\numL16tWoqKjA2rVrUVBQgJycHIOpeDdv3oy4uDgoFAosWrQIx48fR0JCAqqrq5Go1/lUo9EgODgY\nP/zwAxITE/l1QMaPH4/c3Fz069fP8gdOCGmXaBA6Ic0n5J2Qe+saE0KaweYbIGPHjsX06dMbfD05\nORkqlQrnz5+Ht7c3AGDEiBGYPHkytm3bhtmzZwMAVCoVli1bhpCQEOy8N+ffrFmzoNFokJSUhNjY\nWH4xwqysLJw+fRpZWVn8vsPDw+Hj44MVK1YgPT3dkodMCCGEdCj/+pdw+xbr9SU5ehSYMEG4uhDS\nXthsFywtxhgqKipw9+5dk6/v2rULISEhfOMDACZOnAgfHx++oQEAR48eRVlZGebMmWOwfXx8PKqq\nqpCdnc3nZWVlwdPT06Dh4+7ujvDwcOzduxe1tbXmOrxmo24dhBBCbMHt27p09+7C1UNfYKDQNSCk\nfbD5Bkh0dDS6du0KiUSCwMBA5Obm8q8VFxejpKQEw4YNM9pu+PDhBmM2tOn6sQEBARCLxcjPzzeI\nDQgIMFlmdXU1Ll261Orjai7q1kEIIcSWvPSS0DXQobsehDSPzTZAHB0doVAosGHDBuzbtw9vvPEG\nCgoKMHbsWL6xoFQqAQByudxoe7lcjrKyMv5uhVKphJ2dHdzd3Q3iHBwcIJPJcPXqVT5PqVQ2WCYA\ng1hCCCGENN+uXULXQEevEwQhpAlsdgzIqFGjMGrUKP55SEgIFAoF/Pz88Nprr+HLL7+ESqUCwDVW\n6nNycgLAjf2wt7eHSqWCg4ODyX05OjryZQFATU1No2Xez4IFC/jxJFqRkZGIjIy873ak/fvjD25V\n3717gXt/LqSDMkd3yYyMDGRkZBjklZeXt75gQgSm///RFqbilUh0acao1wEhjbHZBogpffv2xRNP\nPIHPP/8cjDFI7p0x7ty5YxRbU1MDAHyMRCKBWq02WW5NTQ0fp41tSpkNSUlJMdmFi9i2H34A/vIX\nLp2fD4wcaf59XL4M9O7N7Uev1yBpQ8z5xcXUhYu8vDwMHTrUfDshRAAlJbp0ZqZw9TDl1ClgzBih\na0FI22azXbAa8sADD0CtVqOqqorvEqXtiqVPqVRCJpPx0/DK5XLU1dUZrQ+iVqtRVlYGLy8vPk8u\nl5vsZqXdj34sIVraxgcA6N28Mxu1mmt8AG1j7nwiLJVKhRUrViAoKAhubm4Qi8X4+OOPTcYKvVbS\nqVOn8Nhjj0EqlUIul2P+/Pmoqqpq2YETm3D9ui7dyDU9q9F+tD/2mLD1IKQ96HANkN9++w0SiQTO\nzs7o2bMnPDw8cPbsWaO4nJwc+Pv788+HDBkCAEax586dg0ajMYj19/dHXl4eWL0+FGfOnIFUKoWP\nj485D4nYqBkzzFfW//4HeHoa5hUVma980v7cvHkTSUlJ+O9//8ufv0Qmbr9o10r67bffsHr1aixa\ntAjZ2dmYPHmy0Yx+mzdvxuzZszF48GCkpqZi1KhRSEhIwNtvv20Qp10rKSMjg3/9xo0bGD9+PH75\n5ReD2Pz8fEycOBE1NTVYv349YmJisGXLFoSFhZn5HSHtiZ+f0DUw9u67QteAkHaE2agbN24Y5eXn\n5zN7e3v21FNP8XlxcXGsc+fO7MqVK3ze4cOHmUgkYps3b+bzVCoVk8lkLDQ01KDMmTNnMmdnZ3bz\n5k0+77PPPmMikYhlZWXxeSUlJaxbt24sMjKywTrn5uYyACw3N7d5B9sEW7cyBjCmd0ikjeF6Dhs+\namvNU3b//sZlDx1qnrKJeY0Zw/1+amosU772PPPdd9+x69evM8YYO3fuHBOJROzjjz82io+Li2NS\nqdTkOXLLli18XnV1dbPPkbt27eLzSkpKmKurK3v22WcNtp82bRrr2bMnq6io4PM+/PBDJhKJ2Fdf\nfdXg8VniPEraDu15bOxYoWuis2+frl56f67EBtF5pvVs9g5IREQEQkJCsGrVKqSlpeGVV17B6NGj\n4ezsjDfffJOPW7p0KTp37owJEyYgNTUVq1evRlhYGPz8/BAdHc3HOTk5ISkpCfv370d4eDg+/PBD\nvPDCC0hPT8eyZcsMBo0rFAqMHDkS0dHRSEpK4ldBZ4xh5cqVVn0fSPvQ0N0Ie3vzDEYuLDTO05uR\nmnRA9vb26H5v8QR2nz8yIddKun37Ng4fPoyZM2fC2dmZj42KioKzs7PB/knHod/D+dtvBauGkZAQ\nXfrXX4WrByHtgc02QJ5++mmUlpZi/fr1iI+PR2ZmJhQKBc6dO4cBAwbwcd7e3jh27Bj69u2LJUuW\nYO3atQgJCcHXX3/Nj//QiouLw5YtW1BQUIC5c+fi9OnTSElJwZIlSwzixGIxDhw4gIiICGzYsAGJ\niYno3r07vvnmG/Tv398qx0/al1OnGn5NLAZ++qnlAy3/+KNl2xEi9FpJBQUFuHv3rlGZ9vb28Pf3\nb3DMCLFtPXvq0uI29C1GvwejXq9sQogJNjsL1rx58zBv3rwmxQ4cOBAHDx5sUmxMTAxiYmIajevW\nrRvS0tKQlpbWpHKt7cIFwM4OUCq5BZR27wZqa7k8c449IE3TwDhd3iOPcD+1F6orKgAXl6aV3aeP\n4fNRo4DTp5tVPWImVVXAuXPA+PFAp07AwYNAejqwZQv3vK1p6lpJ9vb2zV4rafz48SbLBLi1kgYN\nGnTf/Xt6euLkyZMtPjZCCCHCaUPXDoi1bNwIDBoEPPywbvXW6dOBiAhAoQDefhu4fZvLv3OH+3Jc\nb6wpMbMTJ5oWJxIB06YBXboAo0e3bF9ff92y7Ujz3LrFTXv8/fdAZSUwbhzQrRvX+ACAu3eBSZOA\njz7i1n1pi5q6VpL2p7nXSmps/42tqURsj/4M9w38uQlKv34mZuMnhNzTBq+5EUtiDKjXRdvIq69y\nD33z5gHXrnHjBvr3B+RyQK9LNjETmYyb3/7RRwETk7MB4K6aA9xdDJEIePFFbparzz4DunfnrrJ3\n7sy9Vv/G3h9/cK+R1mMMOH8eCAjg0j/+CLz/PlBWBqxaBTRnsjuFgutWkpdnufq2hNBrJTW2/873\n+WOmBV1t0//+p0s3MBO0oPQbRVeuAP36CVcXYh60oKtlUAOkg9D2TW3NgOb607gC3J2UrVu5L7UD\nB3JduEjLrV/P/a5ycoABA4B7XeHva9s27uemTcDy5dwX2eBgrqFYf86DXr3MXuUOJzMTeO89wNsb\nyMgAVqwwfp9bMl6nuBjo0aNtLWDW0rWS9LthtWatpMb2f781lWhBV9v08MO6dFO7oQqlf3/zTCJC\nhEULulpGh+qCtWrVKojFYgwePNjoNaEX2rKWigrzlvfTT9yK3X5+XB/2hx/mxpA8/TQwdiz3ZfqZ\nZ7irwkolrT3RmOee06UvXmzetitWcO/3rVvAp58afykmpqlUwLBhXKPvhx+APXuANWu4v9k//wQK\nCrgrrSIR9wgP57rMaS+Imft91mjMW15rCL1W0iOPPIJOnToZlalWq5Gfn29QJiFtxf79QteAkHZA\n4GmArebKlSusc+fOzNnZmQ0ePNjoNXd3d9a/f3/23nvvseTkZObm5sb8/f2ZWq02iN20aRMTiUQs\nLCyMffjhhywqKoqJRCL21ltvGcTV1dWx0aNHM2dnZ/b666+z999/nw0aNIh16dKFFRYWmqyjJeeV\n/te/TK8zIcRj8mTG/v1vLn3mjK6OV69yeRcvmv3w2zzte2NKVZV53vcTJ3RlRkUxNmCAdY5NX2Pr\nmty9y9hPPzW/3Js3Gfv1V8O8nTsZ01vyh3f5Mve3lpDArYvz66/C/0+Yelh6HRD988zZs2fvuw6I\nkGslTZs2jXl5eZlcB+TQoUNNOj5iO7T/H0uWCF2Tht29e/9zOmn/6DzTeh3m3yMiIoJNmjSJjR8/\nnj3yyCMGrwm90JZWR2mANOWxdStj2dmMOTgwNns2l1dZyX3RnD+fW/Dps88Y++tfGZs+nXtdo2Hs\nnXcY+/RT7phv3WIsI4MxpZKx3FzGYmMZ++YbxpKTGTt/nrH//IexujrDBaOUSm6bu3e556dOMfbd\nd1zZlnS/D6u6OsY8PVv3fl67Znp/dXWWOR7t+6fv/Hlunzk5jFVXc3nvv8/YmjXc+3vnDmNz53Ix\n2n9FtZprgF2+zNjt21zeyZOMjRun+53U1jIWGMht99e/cn8XI0bojnHAAO7nRx8J/3fd1hog7733\nHktKSmJxcXFMJBKxGTNmsKSkJJaUlMRu3brFGNNdoOnXrx9/gcbV1ZX95S9/MbpA88EHH/AXaNLS\n0vgLNKtXrzaIq6urY6NGjWIuLi4GF2i6du3KLl26ZBCbl5fHnJycWEBAANu4cSNbtmwZk0gkLCgo\nqNHjI7YlJUX3/3HnjtC1uT9tPfUvshHbQeeZ1usQDZBjx46xTp06sR9//JE9/vjjRndAunfvziIi\nIoy2GzBgAJs0aRL/PDs7m4lEIvbll18axJ0+fZqJRCK2Y8cOPi8sLIzJ5XKjMl9++WUmlUqNPrgZ\nE6YBEh7OWHAwY8eOcV8K588X/ouXOR5iceu2Dw1lbMgQ3fPRoxkrKOAaKJag3U9TfPKJ6ZXNG3qY\nWvFc+9o335j3OPLzGdu0iSv74kXuS3RL3v++fZsW5+Ul/N9aY4+nn9alV67kfn7+Ofdz5EjG0tMZ\nc3Iy3s4aDZA+ffowkUjERCIRE4vFTCwW8+k//viD3+ann35iU6dOZVKplLm5ubHnn3+e3bhxw2T5\naWlp7OGHH2aOjo6sf//+7N133zUZd/PmTRYTE8Pc3d2ZVCplEyZMaPDcd/LkSTZmzBgmkUhYjx49\n2Lx581hlZWWjx0fapz//NH3RR///o61rT3UlzUfnmdaz+X+Nu3fvMj8/PxYXF8cYY0YNkKKiIiYS\nidiaNWuMtp05cyaTyWT88zfeeIOJRCJWUlJiEHfnzh1mZ2fHFi1axOf169ePBQcHG5Wp7Trw448/\nGr1m7QZIQ2pqGJs3j7F//pN7ro3v00f4L3Nt4XH7NmNPPslYWhpjH3xgnrsjzf2gqqtj7O23GVu3\nrvH63m9/u3e3rt63b3N3iU6eZOy334T/3bS1R2Uldzdo+XLGGuh5afQ70T6s2QXLltj68dmiu3cZ\nKy3l0kFB3N+/t7dhjEbT+HmtLVm4sOG6/vknd06w1B1oYnl0nmk9mx+EvmnTJly+fBlJSUkmX2/q\nQlva2OYstNVQmQBMzgBjTfdbH9HREdiwgRvUDABLl3KDzX//nTudfvwxlz50iFvjID+fi6+/4J2t\n6tKFW7dh9mxuSmOxGPjb36xbB7EYWLwYWLiQ+518+KHpuN9+M53fo4eunKb6/XduPQutZct065E8\n9hjw0ENNL8vWMMatlfPdd8DkydyUyMXFgFTKzQy3cmXj03G++KJVqkpIm9OpE+Dubjh1eFER9zw9\nnXu+a5cuvoHZntuURYt0af31Mp98kptu/fXXuXPDokXc+YOQjsamGyB//vknli9fjuXLl0Mmk5mM\nEXqhLaHMmtX02FWruCl2taKiuMbGlCnc7Fd/+Qu3Tsjvv3MzB129Crzyii6+oAD49ltz1bxt2ryZ\nO8Zff23+tub48Jk1iyvn3Xd1eZWVwIMPmo7XTpfcnA/yhx7ipr1kDLC3B5KTW17f9ig8HPjkEyAm\nhlvAMzWVa8hrF+3s1Ilbv+Wrr7jY+8wQa9LmzeavMyFt3d2793995kzuZ1iYLu/ezM9tmv71x4UL\nuZ/V1cC+fYZx69Zxn5uEdDQ2vQ7I//3f/8Hd3R3z5s1rMEbohbZMscQCWvpfcquqdOuCmJurK/dz\n3TrdehTaeds//5ybnlfLxcX80wILSbuqPMDdJYqKatp25lxlfuJE7uevv3JX3xuiXVYhPLz5DaDm\n3DVpT5YvB55/nvsblsmAo0e56Xft7IDQUN2XHu1Uyeb+H3JwAD74oPGFQpuDFtAibd2bbzYe094/\nJ7SzSDd0TnZ35xYhvTeTNSEdgs02QAoLC5GWloaUlBQU6S0+UVNTA7VajT/++ANdunQRfKEtUyyx\ngFZ1tS5tjZWwRSLg7383zHvqKaCmhrtSrL9gYW4ucOYMEB/PfQl8/XXda6+/zuW1Ny+80PQGiPbq\nmDkMGtS0BoVC0fTF8u7csc4XgK5duTVMWmPLFiA2Vvf8nXeAt97iGk1KJbdS/I0bXHeOY8e4916j\nASQS44U29RuU+izVeAeAuDjzNkBoAS3S1v3jH43H6F+Pe+cdy9XF3C5f1i3+2tidHu1HvkZj2XMM\nIW2FjV7LBIqLi6HRaJCQkICHHnqIf+Tk5ODSpUt48MEHkZSUJPhCW9ZSUmLV3TXI0dF4tfShQ7kv\nXcXFXF957VumUHAfTtqhfKWlXL9gxriTdH4+171r3z7uJJ+RAWRnc9t+/DF3zDNncl1iHnuMy586\nlds+MRE4cgS4dg2oqzP+gh0Swv185RVgxIiWHWtqKreQXWPef79l5bfGBx80PfappwAPj+bvQyQC\n9P/M7ex0XcS03bd+/JFrgDIGlJdzv4dFi7hV3fXvtMyezf2cOZNrLDg5cd2h9P+9Dh/m4vSHqr7y\nCvc7Li7m7jRdvcr9nD6dq0vv3lw3tfqND0KIdWm/gG/dyt2l11+QUz+dkGDderXGAw/o0r6+hq9N\nmmR6m5EjLVcfQtoUoUfBW0ppaSnbs2cP27t3L//Ys2cPe+SRR1ifPn3Y3r17+ZmohF5oS8uSsyoc\nOsR9JUtONnvRFnH7NrcGhLloNIwtXcqtJ9GQmze5qYhv3TKe2erHHxl75hlufYr167npeJs6E9L6\n9fevl1Azu2j3WVzccMywYS2b+UmjMXyY0sAsqry6OsNtd+82Pff/zZuM/fe/jR9ve0GzYLWOrR+f\nLSkouP+577//NT6/tDemzpHbt3Ov1T//t9dj7IjoPNN6He5P/fHHHzdaiFDohba0LPkHfeYMd2L7\n/nuzF91h7dzJvaeTJjX+pbyh9/3UKV1MdrZ166/db0xM4zENPV5/nbHFixlLTWVsxgzGioqsV39b\n1bkzNUBaw9aPz5Y05Ut3e/9y3pQGhouLdabgJuZD55nWs9kuWA0RiUQQ1etg6e3tjWPHjqFv375Y\nsmQJ1q5di5CQEHz99df8+A+tuLg4bNmyBQUFBZg7dy5Onz6NlJQULFmyxCBOLBbjwIEDiIiIwIYN\nG5CYmIju3bvjm2++Qf/+/S1+nMTywsK4j4uvv9ZNT9yQv/yFG0tRn/6g/L/+1fx1bIoPP+S6Julr\nrB+ydjDlP/4BvP02N34nK4ubeIC0Ds2GRToC/XPO3LlN2+aXXyxTF0vas8fw+eXLxjG3bxt+Pqxf\nb9k6EdIW2Owg9IYcPXrUZP7AgQNxUDsBeSNiYmIQExPTaFy3bt2QlpaGtPstukFsRlQU4O3Nzfx1\nb7IzA05OwKlTwKhR3PNLl4Dr161bx4b07Mmt91JYyI350M70ZMr169xEBj/+aL36EUJsi/5kbO+9\n17Rt+va1TF0s6cknDZ/rjwvRpz/D/2uvAfWuaRJiczrcHRBCLCkwEFCpjD90tEaP5u4sPP44MGCA\ndetW38WLhs+Tk7mZse7X+Nizh5tJytmZBksSQlpu0KCmxQkxSYdQ9O+W0OKExNZRA4QQC9iz5/4f\nIMePGz4vKLBsfUxpzkJ5X37JHU9DDStCCLGEhx7ifn70kbD1aI365/uG6J9fbXW9JUK06E+cEAtq\n6HZ7fY88Ytl6mOLi0rS4kBAgKMiydSGEdBz64+V2775/7NSpwM6dTV9XqS0aO5YbV0d3NQjRoQYI\nIRZkasBhfVOnWr4epohEjfe9fvJJ4IsvrFMfQkjH8OKLurT+RBymiETchB/t/Y5AUxcX3LdPl24r\n63cRYgnt/F+6YT/99BPCwsLQt29fSKVSyGQyjB49Gunp6UaxFy9eRFBQEFxcXCCTyRAVFYXS0lKT\n5W7duhW+vr6QSCTw8fFBamqqybjy8nLExsbCw8MDzs7OCAwMxPnz5816jKR9aGwl9717rVMPU+bO\nBcaNM/3amDHGM7gQQkhrXLumS0+ZIlw92qrQUF26e3fh6kGIpdlsA+Ty5cuorKzEiy++iA0bNuCf\n//wn7O3t8fzzz2PVqlV8XFFREcaNG4fffvsNq1evxqJFi5CdnY3JkyejtrbWoMzNmzdj9uzZGDx4\nMFJTUzFq1CgkJCTg7bffNojTaDQIDg5GRkYG//qNGzcwfvx4/NIe5xEkrbJyJfD999wq3h9+aPy6\no6P166Tv2DFu4PzOnVxjpKqK6ypw8qSw9SKE2B65XJfev1+4erRlQt0VJ8SqhF6IxJrq6uqYv78/\n69WrF58XFxfHpFKpyVXQt2zZwudVV1c3exX0Xbt28XklJSXM1dWVPfvssw3WjxYi7Bj27+d+F6NH\nm3e1d2IbPvmEFiJsDVs/vvZOu9jeX/4idE3arspK3ft0+rTQtSGm0Hmm9Wz2DogpYrEY3t7eBosL\n7tq1CyEhIfD29ubzJk6cCB8fH+zcuZPPO3r0KMrKyjBnzhyDMuPj41FVVYXs7Gw+LysrC56enpg+\nfTqf5+7ujvDwcOzdu9fozgrpWKZMAd55Bzh8GKi3ziUhhNgs/WW48vOFq0dbp13oFdCtG0WIrbH5\nBkh1dTVKS0vx66+/Yv369Th06BASExMBAMXFxSgpKcGwYcOMths+fLjBmA1tun5sQEAAxGIx8vXO\npufPn0dAQIDJMqurq3Hp0iWzHBtpn+ztgVdeASQSoWtCCCHWExgodA3aD/2vGhqNcPUgxFJsvgGy\ncOFCdO/eHf3798err76KDRs2IDY2FgCgVCoBAHL9Tqn3yOVylJWV8XcrlEol7Ozs4O7ubhDn4OAA\nmUyGq1ev8nlKpbLBMgEYxBJCCCG2Tn8K2uasQdRRHTqkS4eFCVcPQiylk9AVsLRXXnkF4eHhuHr1\nKtLT0zF37lxIJBK88MILUKlUAABHE6OAnZycAAAqlQr29vZQqVRwcHAwuQ9HR0e+LACoqalptExC\nCCGko9CfRre4WLh6tBeurrp0Y2ulENIe2XwDZMCAARgwYAAAYObMmZg6dSoWLFiAiIgISO71gblz\n547RdjU1NQDAx0gkEqjVapP7qKmp4eO0sU0psyELFixAt27dDPIiIyMRGRl53+0IIURfRkYGMjIy\nDPLKy8sFqg0hpKnqrxty+zbQpYswdSHEEmy+AVLfjBkz8PXXX+Pnn3/mu0Rpu2LpUyqVkMlk/IB1\nuVyOuro6lJaWGnTDUqvVKCsrg5fePWW5XG6ym5V2P16N3H9OSUkxOYaEEEKaw9SFi7y8PAwdOlSg\nGpGO6MIFXfrZZ4WrR3uj0ejuHHXtSiupE9ti82NA6tN2fxKLxejZsyc8PDxw9uxZo7icnBz4+/vz\nz4cMGQIARrHnzp2DRqMxiPX390deXh5YvbPFmTNnIJVK4ePjY7bjIYQQQtqyQYN0aRNrAZMG1L8L\nQhNoEltisw2QkpISo7za2lps374dMpkMg+6dEWfMmIH9+/ejqKiIjzty5AgKCwsRpjfyKzAwEG5u\nbti4caNBmRs3boRUKkVwcDCfp1AocP36dezW67hZWlqKzMxMhIaGGkwDTAghhNgq/QVNo6OFq0d7\ndeuWLn3kiHD1IMTcbLYLVmxsLCoqKjBu3Dh4eXnh2rVrSE9Px6VLl/DRRx/Bzs4OALB06VJkZmZi\nwoQJmD9/PioqKrBmzRr4+fkhWu9s6eTkhKSkJMTHxyM8PBxTpkzBiRMnkJ6ejuTkZIMxGwqFAiNH\njkR0dDQuXLgAmUyGDz74AIwxrFy50urvBSGEECKEsWN16a1bhatHe6U/7mPaNOqGRWyHzTZAnnnm\nGWzduhUbN27En3/+iS5duuDRRx9FamoqJk6cyMd5e3vj2LFjWLhwIZYsWQJHR0eEhIRg3bp1Rncq\n4uLiYG9vj3Xr1mHfvn3o1asXUlJSkJCQYBAnFotx4MABLF68GBs2bIBKpcKIESOwfft29O/f3yrH\nTwghhAhJf76DPn2MuxSRpvm//wPeeINLq9VAAxNyEtKu2GwXrIiICHz11VdQKpVQq9UoLS1Fdna2\nQeNDa+DAgTh48CAqKyvx559/Yvv27fDw8DBZbkxMDC5evIiamhpcunTJqPGh1a1bN6SlpaGkpASV\nlZX45ptvaGA5IaTd+fbbbyEWi00+cnJyDGIvXryIoKAguLi4QCaTISoqCqWlpSbL3bp1K3x9fSGR\nSODj44PU1FSTceXl5YiNjYWHhwecnZ0RGBhosEgsabumTdOlf/9duHq0d0lJurSJGf4JaZds9g4I\nIYQQ85k/fz6GDx9ukNe3b18+XVRUhHHjxsHV1RWrV69GRUUF1q5di4KCAuTk5BjcUd68eTPi4uKg\nUCiwaNEiHD9+HAkJCaiurkZiYiIfp9FoEBwcjB9++AGJiYl8d9bx48cjNzcX/fr1s/yBkxbRaIDv\nvj+lCnkAACAASURBVBO6FoSQtooaIIQQQho1duxYTJ8+vcHXk5OToVKpcP78eXh7ewMARowYgcmT\nJ2Pbtm2YPXs2AG4mwmXLliEkJAQ7d+4EAMyaNQsajQZJSUmIjY3lx9RlZWXh9OnTyMrK4vcdHh4O\nHx8frFixAuk0pVKb9dVXuvS1a8LVw1aUlwPaoaavvQasXi1sfQhpLZvtgkUIIcR8GGOoqKjA3bt3\nTb6+a9cuhISE8I0PAJg4cSJ8fHz4hgYAHD16FGVlZZgzZ47B9vHx8aiqqkJ2djafl5WVBU9PT4OG\nj7u7O8LDw7F3717U0rykbZZ+96sePYSrh63o2lWXfvNN4epBiLlQA4QQQkijoqOj0bVrV0gkEgQG\nBiI3N5d/rbi4GCUlJRg2bJjRdsOHDzcYs6FN148NCAiAWCxGfn6+QaypsXPDhw9HdXU1Ll261Orj\nIuZ3+bIuTVPvms/33+vSmzYJVw9CzMFmGyBnz57F3LlzMWjQIDg7O6N3796IiIhAYWGhUSwNnCSE\nENMcHR2hUCiwYcMG7Nu3D2+88QYKCgowduxYvrGgVCoBAHK53Gh7uVyOsrIy/m6FUqmEnZ0d3N3d\nDeIcHBwgk8lw9epVPk+pVDZYJgCDWNJ29O6tS9PUu+bj56dLx8UJVw9CzMFmx4C89dZbOH36NMLC\nwuDn5welUonU1FQEBATgu+++4xcipIGThBDSsFGjRmHUqFH885CQECgUCvj5+eG1117Dl19+CZVK\nBYBrrNTn5OQEgBv7YW9vD5VKBYcG5hF1dHTkywKAmpqaRsskbUtdneFzmnrXvB57TLe445UrwAMP\nCFsfQlrKZhsgf//73zF8+HB06qQ7xIiICAwePBhvvvkmPvnkEwA0cJIQQpqrb9++eOKJJ/D555+D\nMQaJRAIAuHPnjlFsTU0NAPAxEokEarXaZLk1NTV8nDa2KWWasmDBAoMFYgEgMjISkZGR9zs00kpf\nfqlLV1YKVw9bdfw4IL7Xd6VXL1qY0BoyMjKQkZFhkFeuv8gNaRGbbYDoX7HT6tevHwYOHIiff/6Z\nz2ts4KS2AXK/gZPp6enIzs7Gc889B+D+Ayd37NiB2tpao0UOCSGkPXnggQegVqtRVVXFd4nSdsXS\np1QqIZPJ+HOeXC5HXV0dSktLDbphqdVqlJWVwcvLi8+Ty+Umu1lp96MfW19KSgqtvSSA0FBdWioV\nrh62SiTiHtqGR2Ul4OwsbJ1snakLF3l5eRg6dKhANbINNjsGxBTGGK5fv85/6NHASUIIaZnffvsN\nEokEzs7O6NmzJzw8PHD27FmjuJycHPj7+/PPhwwZAgBGsefOnYNGozGI9ff3R15eHli9y7xnzpyB\nVCqFj4+POQ+JtNL167o0jf2wHP1ubi4uwtWDkNboUA2Q9PR0XL16FREREQBo4CQhhDSmpKTEKO/7\n77/Hvn37MGXKFD5vxowZ2L9/P4qKivi8I0eOoLCwEGFhYXxeYGAg3NzcsHHjRoMyN27cCKlUiuDg\nYD5PoVDg+vXr2L17N59XWlqKzMxMhIaG0p3kNsbTU5d+6SXh6mHraFwNsQU22wWrvp9//hnx8fEY\nPXo0XnjhBQBoswMnqe8yIcQczNF3OSIiAp07d8aoUaPQvXt3XLhwAVu2bIGzszPe1FuQYOnSpcjM\nzMSECRMwf/58VFRUYM2aNfDz80O03lysTk5OSEpKQnx8PMLDwzFlyhScOHEC6enpSE5ONjj3KRQK\njBw5EtHR0bhw4QI/oQdjDCtXrmzhu0Is4ccfdem5c4WrR0eRmwtoewDt2QM89ZSw9SGkuTpEA+Ta\ntWsIDg6Gq6srsrKyILp3+aAtDpwEqO8yIcQ8zNF3+emnn0Z6ejrWr1+P27dvo3v37lAoFFixYgUe\neughPs7b2xvHjh3DwoULsWTJEjg6OiIkJATr1q0zulMRFxcHe3t7rFu3Dvv27UOvXr2QkpKChIQE\ngzixWIwDBw5g8eLF2LBhA1QqFUaMGIHt27ejf//+LXhHiKUMHqxLp6QIV4+OQv8rwtNP02B00v7Y\nfAPk1q1bmDZtGm7fvo0TJ07AU+8ecVscOEkIIW3JvHnzMG/evCbFDhw4EAcPHmxSbExMDGJiYhqN\n69atG9LS0pCWltakcomwEhMBOzuha9ExaDS6GbEKCwFqk5P2xKbHgNTU1CA0NBS//PIL9u/fj4cf\nftjgdRo4SQghhLSO/piE5cuFq0dHo/+++/gAJjpeENJm2WwDpK6uDhEREThz5gwyMzPx6KOPmoyj\ngZOEEEJI63l60tS71va//+nSJ04IVg1Cms1mu2D9/e9/xxdffIHQ0FCUlpZix44dBq/PnDkTAA2c\nJIQQQlpq3Dhd+o8/hKtHR9W7ty49eTKNBSHth802QL7//nuIRCJ88cUX+OKLLwxeE4lEfAOEBk4S\nQgghLaN/1b2BiSKJhV24AAwcyKVrawHqZEHaA5ttgBw9erTJsTRwkhBCCGmeTz/VpWn8gXB8fXVp\nBwe6C0LaB5sdA0IIIYQQy6irA557Tvec7n4Ia9EiXZrWOibtATVACCGEENIs+jMz370rXD0IZ80a\nXbpnT+HqQUhTUQOEEEIIIc2inRDS35/W/Wgr3ntPly4pEa4ehDQFNUAIIYQQ0mRvvaVLm1hGiwhk\n7lxdunt34epBSFPYdAOkqqoKK1asQFBQENzc3CAWi/Hxxx+bjL148SKCgoLg4uICmUyGqKgolJaW\nmozdunUrfH19IZFI4OPjg9TUVJNx5eXliI2NhYeHB5ydnREYGIjz58+b7fgIIYQQa1uyRJfuZLNT\n2bRPnp669K+/ClcPQhpj0w2QkpISJCUl4b///S+/UrlIf+nQe4qKijBu3Dj89ttvWL16NRYtWoTs\n7GxMnjwZtbW1BrGbN2/G7NmzMXjwYKSmpmLUqFFISEjA22+/bRCn0WgQHByMjIwM/vUb/9/evUc1\nceVxAP8mEiBEUAF5RIt4UOxWRUSkitZSqW9xffDQrW/Bc6gV0T2Ki117rF2trVakbLXqLqUrtSq+\nWFdXFx+o6wsVdtXKgh4VgaAiokESELj7R5pASICASWaCv885OczcuQzfhMmd3MzMncePERQUhDt3\n7pjuSRNCCCEm8rvf1U832j0SHigqqp/u1Yu7HIS0pF1/dyGVSlFSUgIXFxdcu3YNgwcP1ltv3bp1\nUCgUyM7ORvfu3QEAAQEBGDVqFH744QdERUUBABQKBVatWoWJEydi7969AIAFCxagrq4Oa9euxcKF\nCzU3JExLS8PFixeRlpaGqVOnAgDCw8Ph7e2Nzz77DKmpqaZ++oQQQojRMAbs3l0/T0c/+EcoBHbu\nBNR3C/juO+Djj7nNRIg+7foIiLW1NVx+PRGSNTMw9v79+zFx4kRN5wMAgoOD4e3treloAKp7i5SV\nleHjRu/mRYsW4eXLl/jHP/6hKUtLS4Obm5um8wEAzs7OCA8Px+HDh3WOrBBCCCF8JmzwieH5c+5y\nkOYtWFA/vWgRdzkIaU677oAYoqioCE+ePIG/v7/OssGDB2tds6GeblzXz88PQqEQOTk5WnX9/Pz0\nrrOyshJ5eXnGegqEEEKISa1YUT89dizg4MBdFtKy27frpwMCuMtBSFPe+A6ITCYDALi7u+ssc3d3\nR1lZmeZohUwmQ4cOHeDs7KxVz9raGk5OTihucPcfmUzW5DoBaNUlhBBC+CopSfs+E8eOcZeFGObt\nt+uns7KA997jLgsh+rzxZ3AqFAoAgI2Njc4yW1tbTR2RSASFQgHrJm73amNjo1kXACiVyhbX2ZTY\n2FjNtSRqM2bMwIwZM1p4NoQQUm/37t3Y3fCkfahG5yOkNRredJAGcrQcZWWAo6Nq+vx5brMQ0tgb\n3wERi8UAgKqqKp1lSqVSq45YLEZ1dbXe9SiVSk09dV1D1qlPQkKC3tO3CCGkNfR9cXH9+nUMGjSI\no0TEkk2apLrxILEMXboAly8D776rmq+sBOzsuM1EiNobfwqW+pQo9alYDclkMjg5OUEkEmnq1tbW\n6twfpLq6GmVlZZBKpVrr1XealfrvNKxLCCGE8FFBQf304cPc5SBtExAA9O6tmp4/n9sshDT0xndA\nunXrhq5duyJLz+1cr1y5orl/CAAMHDgQAHTqXr16FXV1dVp1fX19cf36dZ3Rty5fvgyJRAJvb29j\nPg1CCCHE6Hr0UP3s2JHbHKTttm5V/dyzRzWUMiF88MZ3QABg2rRpOHLkCAoLCzVlJ0+eRH5+PsLC\nwjRlI0eOhKOjI7aq382/2rp1KyQSCSZMmKApCw0NxaNHj3DgwAFNWWlpKfbt24eQkBDNURVCCCGE\nj/7+9/rp//2Puxzk9YwcWT/9623NCOFcu78GJCkpCeXl5ZrTodLT01Hw6zHlmJgYODg4ID4+Hvv2\n7cMHH3yAJUuWQC6X4+uvv4aPjw/mzZunWZetrS3Wrl2LRYsWITw8HKNHj8a5c+eQmpqKdevWaV04\nHhoaiiFDhmDevHn45Zdf4OTkhO+++w6MMaxZs8a8LwIhhBDSCnV1qms+1OisYcslENRP/+UvQEIC\nHdEi3Gv3HZBNmzbhwYMHAACBQICDBw/iwIEDEAgEmD17NhwcHNC9e3dkZmZi2bJlWLlyJWxsbDBx\n4kRs2rRJ50hFdHQ0RCIRNm3ahPT0dHh4eCAhIQExMTFa9YRCIY4ePYrly5cjMTERCoUCAQEB+PHH\nH9FbfUImIYQQwkMdOtRPV1Zyl4MYB2P1HRF7e1UHs2HHhBBza/cdkHv37hlU75133sE///lPg+pG\nRkYiMjKyxXqdO3fGjh07sGPHDoPWSwghhHDNw6N+OjgYaGbQRmJBLlwAAgNV00IhdUIIt+gaEEII\nIYQAANasAR4+rJ/PyOAuCzGuoUMBH5/6eaEQqK3lLg95s7X7IyCEEEIIaVlsLLBlS/08fThtf/7z\nH+2jHlZWdCSEcIM6IIQQQsgbrvEH0Npa1TfkpP1peD0IoPo/l5WpblxIiLlQ82JCVVVViIuLg1Qq\nhZ2dHYYMGYIMOp5NCCEGo3bUtDIzdTsfdXXU+WjvGANsbevnHR2BqVO5y0PePNTEmNDcuXOxefNm\nzJo1C4mJiejQoQPGjx+Pf//731xHM4ndu3dzHaHNKDs3KDtpCbWjpqE+7SYoqL6sb1+gpqbtp+NY\n8nviTcyuUABHjtTPHzyo+t8/fWqkYAaw5NedvB7qgJjIlStXsGfPHnz55ZfYsGEDIiMjcerUKfTo\n0QMrVqzgOp5JWHJDQtm5QdlJc6gdNT51B6PhMLsAcO0acPOmbnlrWPJ74k3NPmEC8OyZdpmzs2ob\naTgYgalY8utOXg91QEwkLS0NVlZWWLhwoabMxsYGCxYswMWLF1FUVMRhOkII4T9qR42juhr48UfV\nh8pGt7ZCSorqdBw/P26yEe517qzaBjZt0i738FBtMwIBIJNxk420X9QBMZHs7Gx4e3ujY6PbjQ4e\nPBgAkJOTw0UsQgixGNSOtk1tLbBjR/2HRxsbYM4c7Tp2dqrTsGbP5iYj4Z9ly1QdkTVrdJdJpfXb\n04YNQEWF+fOR9oU6ICYik8ng7u6uU64uKy4uNnckQgixKNSO1mMMqKoCysuB//4XSEpS3VRO/aGw\n4cPKCmhw0EhLdrZqXS9f0tCrRL/Vq1XbiFKpf/nKlaq7qevb9vz9gY0bgbNngdu3gefPVesipDEa\nhtdEFAoFbGxsdMptfx12QqFQ6CxT/vpuv337ttHznDyp+jlmDNC7t9FXDwC4dascI0ZcN83KTYyy\nc4Oy6zp3TvUzOxuwtjb66jXti742iG+M3Y7OmaO6zoHfyiEQGG+76tgRGDiwfj4mxmir1kHvZ26Y\nMvt776l+1tYCeXlAaWnz9a9dUz0MZ9zt3RjCw4G4uObrWFI7ylfUATERsViMqqoqnXL1zlEsFuss\nu3fvHgBg5syZJstVUqJ6mMq5c4NMt3ITo+zcoOz6DR1qslUDAO7fv49hw4aZ9o+8Jr62o6ZnvO2q\noqK+U2sO9H7mhiVnN+b2bgx796oehrCEdpSvqANiIu7u7npPD5D9eiWXVCrVWTZmzBjs2rULnp6e\neneshBDyupRKJe7du4cxY8ZwHaVF1I4SQvjIktpRvqIOiIkMHDgQZ86cgVwuh729vab88uXLAABf\nX1+d33F2dsZHH31ktoyEkDdTYGAg1xEMQu0oIYSvLKUd5Su6CN1EQkNDUVtbi+3bt2vKqqqqkJyc\njCFDhqBbt24cpiOEEP6jdpQQQtonOgJiIgEBAQgLC8Mf/vAHPH78GF5eXkhJSUFBQQGSk5O5jkcI\nIbxH7SghhLRPAsZogDRTqaqqwh//+Efs2rULz549w4ABA7B27VqMGjWK62iEEGIRqB0lhJD2h07B\nMiEbGxt89dVXKC4uhkKhwKVLl/TuNM+ePYtJkybBw8MDYrEYrq6uCA4OxrFjx/Su98KFCxg+fDgk\nEgnc3d2xZMkSvHz50tRPR6+TJ09i/vz58Pb2hkQigZeXF6KiolDSxFBbfMpeUlKClStX4oMPPoC9\nvT2EQiEyMzObrM+n7FVVVYiLi4NUKoWdnR2GDBmCjIwMTrI05eXLl/jss88wduxYODo6QigUIiUl\nRW/d27dvY+zYsbC3t4eTkxNmz56N0pbGezSRrKwsfPLJJ+jbty86duyIHj16ICIiAvn5+Tp1+ZQb\nAG7duoWwsDB4eXlBIpHAyckJgYGBSE1N1anLt+xNMbQdbU5GRgZGjhyJzp07w8HBAf7+/thr6DA3\nPBEVFQWhUIiQkBCuoxiktfsGLlhCO6pPa9ooS/CnP/0JQqEQ/fv35zqKQa5fv45JkybByckJEokE\n/fv3x7fffst1LItDp2DxQH5+PqysrBAdHQ03NzeUlZVh165dmDBhAv72t79pXVCZk5OD4OBg9O3b\nF5s3b8bDhw+xceNG5Ofn4+jRo2bPHhcXh/LycoSFhaF37964e/cukpKScOTIEeTk5MDV1ZW32XNz\nc/HVV1/B29sbPj4+uHjxIgRN3JmLb9nnzp2L/fv3Y+nSpejduzeSk5Mxfvx4nD59mjdDAj558gRr\n165Fjx494OvrizNnzuh9fQsLCzFixAh06dIF69evh1wux8aNG3Hjxg1cuXIFIpHIrLk3bNiAixcv\nIiwsDD4+PpDJZEhKSoKfnx8uXbqEvn378jI3ABQUFKCiogJz586FVCpFZWUl0tLSMGvWLNy/fx+r\nVq3ibXZTSU5ORmRkJEaPHo3169ejQ4cOyM3NRWFhIdfRDHb16lWkpKTA1ta2yTaKb1qzb+CKJbSj\n+hjaRlmCwsJCrFu3DhKJxCK27RMnTiAkJASDBg3C6tWr0bFjR9y5cwdFRUVcR7M8jPBSZWUlc3Nz\nYyNGjNAqHzduHOvWrRuTy+Wasp07dzKBQMBOnDhh7pjs3LlzOmVnz55lAoGAffrpp1rlfMsul8vZ\ns2fPGGOM7du3jwkEApaZmam3Lp+yX758mQkEArZp0yZNmVKpZL169WKBgYFmzdKcqqoq9ujRI8YY\nY1evXmUCgYClpKTo1IuOjmYSiYQ9fPhQU5aRkcEEAgHbvn272fKqXbhwgb169UqrLD8/n9na2rKZ\nM2dqyviWuym1tbXM19eXeXh4aMosJfvrunfvHhOLxSw2NpbrKG1WV1fHhg4dyiIjI5mnpycLCQnh\nOpJBWrNv4IKltKP6GNpGWYKIiAj24YcfsqCgINavXz+u4zTr+fPnzNXVlU2bNo3rKO0CdUB4rF+/\nfiw4OFgz//z5cyYSiVhcXJxWverqamZvb88iIyPNHbFJTk5OLDQ0VDPP9+zNdUD4ln358uVMJBJp\ndYYYY2z9+vVMIBCwwsJCs+YxRFZWVpMdEBcXFxYREaFT3qdPH/bhhx+aI55B/Pz8mL+/v2beUnIz\nxtjEiROZl5eXZt6Ssr+OuLg4Zmtry168eMEYU33pUFdXx3Gq1klJSWGdOnViJSUlrEePHhbTAWlK\n430DVyyxHW1J4zaK7zIzM5mVlRW7efMme//991n//v25jtSsrVu3MoFAwHJzcxljjFVUVLDa2lqO\nU1kuugaER168eIHS0lLk5uYiPj4eeXl5+P3vf69ZfuPGDdTU1MDf31/r90QiEXx9fZGdnW3uyHpV\nVFRALpfD2dlZU2Yp2fXhW/bs7Gx4e3ujY8eOWuWDBw8GoDpdzFIUFRXhyZMnOq8toHo+fNkuGGN4\n9OiRZpvme+7KykqUlpbi7t272Lx5M44fP44VK1YA4H92Y8rIyMDbb7+NI0eOoHv37nBwcICzszNW\nr14NZgHjr8jlcsTFxSE+Pp4Xpyy9Ln37Bq60p3YU0G2j+K62thaLFy9GVFSUxZwylpGRAQcHBzx8\n+BB9+vSBvb09OnXqhI8//hhVVVVcx7M4dA0Ij4SHh+PEiRMAAIlEgr1792LcuHGa5eq7/7q7u+v8\nrpubG86fP2+eoC1ISEjAq1evEBERoSmzlOz68C27TCbTm0Vdpu/O0XzV3Gvr7u6OsrIyvHr1ivNr\nElJTU1FcXIwvvvgCAP9zL1u2THPvDCsrKyQmJmLhwoUA+J/dmNTX182fPx9xcXEYMGAA9u/fjy++\n+AI1NTVYt24d1xGb9fnnn0MikWDp0qVcRzEKffsGrrSndhTQbaP4btu2bSgoKMCpU6e4jmKw/Px8\n1NTUYPLkyYiMjMSGDRtw+vRpfPvttygvL8dPP/3EdUSLQh0QI2OMGdwTtrW11ZrfsGEDli9fjoKC\nAmzbtg3Tp09Henq6ZsQXhUIBQDUqjL51qZdzkV3t7NmzWLNmDSIiIhAUFKQpt4TsTTF19tZSKBRN\nZlEvtxQtvbbqOlx+GM7NzcWiRYsQGBiIOXPmaDIB/M29dOlShIeHo7i4GKmpqfjkk08gFosxZ84c\n3mdvSlve4xUVFWCMadpWAJgyZQrKysqwZcsWxMfH63wDbgptyZ6Xl4fExET8/PPPnP8vTLlv4Ep7\nakf1tVF89vTpU6xevRqrV6+Gk5MT13EMVlFRgcrKSkRHRyMhIQEAMHnyZFRXV+P777/H559/jl69\nenGc0nLQKVhGlpmZCTs7O4MeeXl5Wr87YMAABAcHY968eTh//jx69eqFRYsWaZaLxWIA0LsjUCqV\nsLOz4yw7oGoEp0yZAh8fH+zcuVNrGd+zN8fU2VtLLBY3mUW93FK09No2rMOFkpISTJgwAV26dEFa\nWppmlBa+5+7Tpw9GjhyJmTNn4tixYwgODkZsbCyUSiXvszelLe9xsVgMgUCAGTNmaK1r+vTpUCgU\nZjvNpi3ZlyxZgmHDhmHKlClmydgcU+4buNJe2tGm2ig++/TTT+Hs7IzFixdzHaVV1NtE4/ZEPX/p\n0iWzZ7JkdATEyH7zm9/ghx9+MKium5tbk8tEIhFCQkLw5Zdfory8HJ07d9YcGlafQtGQTCaDVCpt\nU2a118n+8OFDjB49Gl26dMHRo0chkUi0lvM5e0tMnb213N3d9Z4eoM5n7jyvo6XX1snJibNvf58/\nf45x48bhxYsXOHfunNZ2w+fc+kybNg3/+te/kJuba3HZ1dryHpdKpbh7967O9RMuLi4AgGfPnhk1\nY1Nam/3UqVM4fvw4Dhw4gPv372uW1dTUoLKyEg8ePICjoyPs7e1NE7gRU+4buNIe2tHm2ii+ys/P\nx44dO5CQkKA1FLZSqUR1dTUePHgABwcHdOnShcOU+kmlUvzyyy+ctyftBXVAjMzV1RWzZ882yrrU\nh4CFQtWBqn79+sHKygpZWVkIDQ3V1KuurkZOTg6mT5/+Wn+vrdmfPn2K0aNH49WrVzh9+rTeiyX5\nmt0Qps7eWgMHDsSZM2cgl8u1PoBcvnwZAODr62vWPK+jW7du6Nq1K7KysnSWXblyhbPnolQqERIS\ngjt37mguZG6Ir7mb0rAtsbTsam15j/v7++POnTsoLCxEz549NeXqD55du3Y1asamtDZ7QUEBAGDq\n1Kk6y4qLi9GzZ08kJCQgJibGaBmbY8p9A1csvR1tqY3iq6KiItTV1SEmJkbv9tuzZ0/Exsbim2++\n4SBd8/z9/ZGRkYHCwkL07t1bU27u9qTd4Gr4LVJPfa+Ehp49e8beeust5uPjo1U+btw4JpVK9d6P\n4vjx4ybP2lhFRQULCAhgnTp1YtevX2+2Lt+yN2TIfUD4kl09fv3GjRs1Zerx64cOHWrWLIZqbhje\n6OhoZmdnp/eeFN9//705YzLGGKupqWGTJk1i1tbW7NixY03W41tuxhh7/PixTll1dTXz8/Njzs7O\nrKamhjHGz+ymcOjQISYQCNiqVas0ZbW1tWz48OHM2dmZVVdXc5iuaQUFBezw4cNaj0OHDjEXFxcW\nEBDADh8+zO7evct1zGa1Zt/ABUtsR9UMbaP4qLS0lB06dEhn2+7Xrx/z9PRkhw8fZjdv3uQ6pl7Z\n2dlMIBCwjz76SKt8xowZzNramslkMo6SWSYBYxYwFmE7N2jQILz11lsICAiAi4sLCgoKkJycjCdP\nnuDo0aMYOXKkpm52djYCAwPxzjvvICoqCoWFhfjmm2/w/vvv49ixY2bPPnnyZKSnp2P+/Pk6Fxba\n29vjt7/9rWaeb9kBaEYMuXXrFvbs2YP58+fD09MTgOo8VTW+ZY+IiMDBgwexdOlSeHl5ISUlBVev\nXsXJkycxfPhws+dpSlJSEsrLy1FcXIxt27Zh6tSpmm8WY2Ji4ODggMLCQgwcOBCdO3fGkiVLIJfL\n8fXXX8PDwwNZWVlmPx0oNjYWiYmJCAkJQVhYmM7ymTNnAgDvcgOqC6zlcjlGjBgBqVSKkpISpKam\nIi8vD8nJyZpvsfmY3VRGjRqFU6dOISoqCj4+Pjh06BAyMjKwfft2REZGch2vVTw9PeHj44P09HSu\no7SoNfsGrlhKO9qYoW2UJQkKCsLTp09x48YNrqM0KzIyEn/9618RHh6OESNG4MyZM0hLS0N8aup+\n/gAAATJJREFUfLzFjEDGG1z3gAhjf/7zn9l7773HunbtykQiEXNzc2OhoaHs2rVreuufP3+eDRs2\njInFYubq6soWL17MKioqzJxaxdPTkwmFQiYQCHQePXv21KnPp+yMMSYQCDT5Gz4PoVCoU5dP2ZVK\nJVu+fDlzd3dntra27N133+XkbvIt8fT01HpNG77WDx480NS7desWGzNmDJNIJMzR0ZHNmjVL77f5\n5hAUFNTkNt14u+BTbsYY+/nnn9moUaOYm5sbE4lEzMnJiY0fP55lZGTo1OVbdlOpqKhgsbGxzN3d\nndnY2LABAwawn376ietYbWJJd0Jv7b6BC5bSjjbWmjbKUgQFBfH+RoSMMfbq1Su2Zs0a5unpyayt\nrZm3tzfbsmUL17EsEh0BIYQQQgghhJjN/wGG0EYMzfD9ywAAAABJRU5ErkJggg==\n",
       "prompt_number": 13,
       "text": [
        "<IPython.core.display.Image at 0x60e83350>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# python: 10.3s\n",
      "# cython: 9.8s\n",
      "# cimport and cdef: 13.6s\n",
      "# dtype and ndim: 452ms\n",
      "# no inj or plot: 13.5ms\n",
      "# 100x part: 1.25s (might have been fluke; more like below)\n",
      "# no bounds check: 973ms\n"
     ],
     "language": "python",
     "outputs": [],
     "prompt_number": 122
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#dir(MPI.COMM_WORLD)\n",
      "help(MPI.COMM_WORLD.Allreduce)"
     ],
     "language": "python",
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Help on built-in function Allreduce:\n",
        "\n",
        "Allreduce(...)\n",
        "    Comm.Allreduce(self, sendbuf, recvbuf, Op op=SUM)\n",
        "    \n",
        "    All Reduce\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "outputs": []
    }
   ]
  }
 ]
}