Figures_simulations_noise_free.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dipy Simulations\n",
    "\n",
    "Dipy provides tools to generate synthetic signals. These synthetic signals are useful for testing and validating the Dipy procedures during their development. Moreover, for education proposes, these tools are also useful to explore in more detail the properties of different diffusion parameters.\n",
    "\n",
    "To illustrate the sensitivity and specificity of different DKI metrics, DKI is first processed on single voxel synthetic signals for four different sets of ground truth parameters (cases 1 to 4)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's load relevant modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from dipy.core.sphere import disperse_charges, HemiSphere\n",
    "from dipy.sims.voxel import (multi_tensor_dki, single_tensor, multi_tensor, \n",
    "                             _check_directions, all_tensor_evecs)\n",
    "from dipy.data import (get_fnames, get_sphere)\n",
    "from dipy.io.gradients import read_bvals_bvecs\n",
    "from dipy.core.gradients import gradient_table\n",
    "from dipy.reconst.dti import (decompose_tensor, from_lower_triangular)\n",
    "from dipy.sims.voxel import add_noise\n",
    "from dipy.data import get_sphere\n",
    "from dipy.reconst.dti import fractional_anisotropy\n",
    "from dipy.viz import window, actor\n",
    "import dipy.reconst.dki as dki\n",
    "import dipy.reconst.dti as dti\n",
    "import dipy.reconst.msdki as msdki\n",
    "sphere = get_sphere('repulsion724')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The synthetic signals for each set of ground truth parameters are generated according to the gradient directions and b-values of the CFIN dataset. So let's load a object containing the acquisition parameters of the CFIN dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraw, fbval, fbvec, t1_fname = get_fnames('cfin_multib')\n",
    "bvals, bvecs = read_bvals_bvecs(fbval, fbvec)\n",
    "gtab = gradient_table(bvals, bvecs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bval1000 = gtab.bvals == 1000\n",
    "bval2000 = gtab.bvals == 2000\n",
    "bval3000 = gtab.bvals == 3000\n",
    "bval_sel = bval1000 | bval2000 | bval3000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bvals = gtab.bvals[bval_sel]\n",
    "bvecs = gtab.bvecs[bval_sel, :]\n",
    "\n",
    "# add 6 b-values\n",
    "bvals = np.hstack((np.zeros(6), bvals))\n",
    "bvecs = np.vstack((np.zeros((6, 3)), bvecs))\n",
    "gtab = gradient_table(bvals, bvecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voxel Case 1 (Single Tensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define diffusivities of voxel case 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "eval1 = [1.7e-3, 0.3e-3, 0.3e-3];\n",
    "angle = [(130.0, 0.0)]\n",
    "stick = _check_directions(angle)\n",
    "evec1 = all_tensor_evecs(stick[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot tensor component of voxel case 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = np.array([[[eval1]]])\n",
    "evecs = np.array([[[evec1]]])\n",
    "\n",
    "# colors\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = 1;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=1)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "#dir_actor = actor.peak_slicer(np.array([[[stick]]]), np.array([[[1.8]]]), colors=None)\n",
    "#dir_actor.display(y=0)\n",
    "#dir_actor.RotateX(90)\n",
    "#scene.add(dir_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case1.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ground truth total diffusion tensor will be equal to the single component diffusion tensor of this simulation case. In addition to the tensor shape, we also plot the simulations main direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_actor = actor.peak_slicer(np.array([[[stick]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor.display(y=0)\n",
    "dir_actor.RotateX(110)\n",
    "scene.add(dir_actor)\n",
    "\n",
    "window.show(scene)\n",
    "window.record(scene, n_frames=1, out_path='case1_dt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal1 = single_tensor(gtab, S0=1, evals=eval1, evecs=evec1, snr=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "dti_model = dti.TensorModel(gtab, fit_method='NLLS')\n",
    "dti_fit = dti_model.fit(signal1)\n",
    "\n",
    "evals = np.array([[[dti_fit.evals]]])\n",
    "evecs = np.array([[[dti_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = 1;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=1)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case1_dt_dti_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dki_model = dki.DiffusionKurtosisModel(gtab, fit_method='NLS')\n",
    "dki_fit = dki_model.fit(signal1)\n",
    "\n",
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = 1;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=1)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case1_dt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "refmk=0.8\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(110)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case1_kt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voxel Case 2 (Intra - and extra-cellular Tensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define diffusivities of voxel case 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "eval1 = [1.4e-3, 0.1e-3, 0.1e-3]\n",
    "eval2 = [2.0e-3, 0.5e-3, 0.5e-3]\n",
    "angle1 = (130.0, 0.0)\n",
    "angle2 = (130.0, 0.0)\n",
    "stick1 = _check_directions([angle1])\n",
    "stick2 = _check_directions([angle2])\n",
    "evec1 = all_tensor_evecs(stick1[0])\n",
    "evec2 = all_tensor_evecs(stick2[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the two tensor components of voxel case 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = np.array([[[eval1, eval2]]])\n",
    "evecs =  np.array([[[evec1, evec2]]])\n",
    "\n",
    "colors = np.zeros(evals.shape)\n",
    "colors[0, 0, 0, 2] = 1;  # blue\n",
    "colors[0, 0, 1, 0] = 1;  # yellow\n",
    "colors[0, 0, 1, 1] = 1;  \n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=colors, sphere=sphere, scale=1)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ground truth total diffusion tensor and total kurtosis tensors can be computed from the multi_tensor_dki. Since this function ignores high-order-terms, it produced diffusion and kurtosis tensors are analogous to the predicted tensors from equations 22 and 23."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "mevals = [eval1, eval2]\n",
    "S0 = 1\n",
    "angles = [angle1, angle2]\n",
    "fractions = [50, 50]\n",
    "signal2_gt, dt, kt = multi_tensor_dki(gtab, mevals, S0, angles, fractions, snr=None)\n",
    "\n",
    "dki_model = dki.DiffusionKurtosisModel(gtab, fit_method='NLS')\n",
    "dki_fit = dki_model.fit(signal2_gt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the ground truth DT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "reffa = dki_fit.fa\n",
    "refmd = dki_fit.md\n",
    "refmk = dki_fit.mk()\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2_dt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the ground truth KT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(110)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2_kt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's know produce signal considering high order terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal2, st = multi_tensor(gtab, mevals, S0, angles, fractions, snr=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit dti for reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dti_model = dti.TensorModel(gtab, fit_method='NLLS')\n",
    "dti_fit = dti_model.fit(signal2)\n",
    "\n",
    "evals = np.array([[[dti_fit.evals]]])\n",
    "evecs = np.array([[[dti_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dti_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dti_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2_dt_dti_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot dt for dki fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "dki_fit = dki_model.fit(signal2)\n",
    "\n",
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2_dt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot fitted KT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(110)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case2_kt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voxel case 3 (simplistic model of axonal \"degeneration\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "eval1 = [1.4e-3, 0.1e-3, 0.1e-3]\n",
    "eval2 = [2e-3, 0.7e-3, 0.7e-3]\n",
    "angle1 = (130.0, 0.0)\n",
    "angle2 = (130.0, 0.0)\n",
    "stick1 = _check_directions([angle1])\n",
    "stick2 = _check_directions([angle2])\n",
    "evec1 = all_tensor_evecs(stick1[0])\n",
    "evec2 = all_tensor_evecs(stick2[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the two tensor components of voxel case 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = np.array([[[eval1, eval2]]])\n",
    "evecs =  np.array([[[evec1, evec2]]])\n",
    "\n",
    "colors = np.zeros(evals.shape)\n",
    "colors[0, 0, 0, 2] = 1;  # blue\n",
    "colors[0, 0, 1, 0] = 1;  # orange\n",
    "colors[0, 0, 1, 1] = 0.5;  \n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=colors, sphere=sphere, scale=1)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot GT diffusion tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "mevals = [eval1, eval2]\n",
    "S0 = 1\n",
    "angles = [angle1, angle2]\n",
    "fractions = [30, 70]\n",
    "signal3_gt, dt, kt = multi_tensor_dki(gtab, mevals, S0, angles, fractions, snr=None)\n",
    "\n",
    "dki_fit = dki_model.fit(signal3_gt)\n",
    "\n",
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "#reffa = dki_fit.fa\n",
    "#refmd = dki_fit.md\n",
    "#refmk = dki_fit.mk()\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.4]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3_dt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ground truth kurtosis tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(110)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3_kt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot reference DTI diffusion tensor fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal3, dt, kt = multi_tensor_dki(gtab, mevals, S0, angles, fractions, snr=None)\n",
    "\n",
    "dti_model = dti.TensorModel(gtab, fit_method='NLLS')\n",
    "dti_fit = dti_model.fit(signal3)\n",
    "\n",
    "evals = np.array([[[dti_fit.evals]]])\n",
    "evecs = np.array([[[dti_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dti_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dti_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.4]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3_dt_dti_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot DKI diffusion tensor fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "dki_fit = dki_model.fit(signal3)\n",
    "\n",
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(110)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.4]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3_dt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot DKI kurtosis tensor fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(110)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(110)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case3_kt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Voxel Case 4 (Crossing Fibers  )¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "eval1 = [1.4e-3, 0.1e-3, 0.1e-3]\n",
    "eval2 = [2e-3, 0.5e-3, 0.5e-3]\n",
    "eval3 = [1.4e-3, 0.1e-3, 0.1e-3]\n",
    "eval4 = [2e-3, 0.5e-3, 0.5e-3]\n",
    "angle1 = (30.0, 0.0)\n",
    "angle2 = (30.0, 0.0)\n",
    "angle3 = (150.0, 0.0)\n",
    "angle4 = (150.0, 0.0)\n",
    "stick1 = _check_directions([angle1])\n",
    "stick2 = _check_directions([angle2])\n",
    "stick3 = _check_directions([angle3])\n",
    "stick4 = _check_directions([angle4])\n",
    "evec1 = all_tensor_evecs(stick1[0])\n",
    "evec2 = all_tensor_evecs(stick2[0])\n",
    "evec3 = all_tensor_evecs(stick3[0])\n",
    "evec4 = all_tensor_evecs(stick4[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot case 4 components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals = np.array([np.array([[eval1, eval2]]),\n",
    "                   np.array([[eval3, eval4]])])\n",
    "stick = _check_directions(angle)\n",
    "evecs =  np.array([[[evec1, evec2]],\n",
    "                   [[evec3, evec4]]])\n",
    "\n",
    "colors = np.zeros(evals.shape)\n",
    "\n",
    "colors[0, 0, 0, 2] = 1;  #Blue\n",
    "\n",
    "colors[0, 0, 1, 0] = 1;  #Yellow\n",
    "colors[0, 0, 1, 1] = 1;\n",
    "\n",
    "colors[1, 0, 0, 2] = 1;  #Blue\n",
    "\n",
    "colors[1, 0, 1, 0] = 1;  #Yellow\n",
    "colors[1, 0, 1, 1] = 1;\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=colors, sphere=sphere, scale=0.8)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(100)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot gt DT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "mevals = [eval1, eval2, eval3, eval4]\n",
    "S0 = 1\n",
    "angles = [angle1, angle2, angle3, angle4]\n",
    "fractions = [25.0, 25.0, 25.0, 25.0]\n",
    "\n",
    "signal4_gt, dt, kt = multi_tensor_dki(gtab, mevals, S0, angles, fractions, snr=None)\n",
    "\n",
    "dki_fit = dki_model.fit(signal4_gt)\n",
    "\n",
    "evals = np.array([[[dki_fit.evals]]])\n",
    "evecs = np.array([[[dki_fit.evecs]]])\n",
    "\n",
    "#reffa = dki_fit.fa\n",
    "#refmd = dki_fit.md\n",
    "#refmk = dki_fit.mk()\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(100)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(100)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "dir_actor2 = actor.peak_slicer(np.array([[[stick3]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor2.display(y=0)\n",
    "dir_actor2.RotateX(100)\n",
    "scene.add(dir_actor2)\n",
    "\n",
    "dtidir_actor = actor.peak_slicer(np.array([[[dki_fit.directions[0]]]]), np.array([[[1.2]]]), linewidth=3, colors=(1, 0, 0))\n",
    "dtidir_actor.display(y=0)\n",
    "dtidir_actor.RotateX(100)\n",
    "scene.add(dtidir_actor)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4_dt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot gt KT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(100)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(100)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "dir_actor2 = actor.peak_slicer(np.array([[[stick3]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor2.display(y=0)\n",
    "dir_actor2.RotateX(100)\n",
    "scene.add(dir_actor2)\n",
    "\n",
    "\n",
    "dtidir_actor = actor.peak_slicer(np.array([[[dki_fit.directions[0]]]]), np.array([[[1.2]]]), linewidth=3, colors=(1, 0, 0))\n",
    "dtidir_actor.display(y=0)\n",
    "dtidir_actor.RotateX(100)\n",
    "scene.add(dtidir_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4_kt_gt.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal4, st = multi_tensor(gtab, mevals, S0, angles, fractions, snr=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot ref DTI DT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "dti_fit = dti_model.fit(signal4)\n",
    "\n",
    "evals = np.array([[[dti_fit.evals]]])\n",
    "evecs = np.array([[[dti_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dti_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dti_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(100)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(100)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "dir_actor2 = actor.peak_slicer(np.array([[[stick3]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor2.display(y=0)\n",
    "dir_actor2.RotateX(100)\n",
    "scene.add(dir_actor2)\n",
    "\n",
    "dtidir_actor = actor.peak_slicer(np.array([[[dti_fit.directions[0]]]]), np.array([[[1.2]]]), linewidth=3, colors=(1, 0, 0))\n",
    "dtidir_actor.display(y=0)\n",
    "dtidir_actor.RotateX(100)\n",
    "scene.add(dtidir_actor)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4_dt_dti_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot DKI DT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "dki_fit = dki_model.fit(signal4)\n",
    "\n",
    "evals = np.array([[[dti_fit.evals]]])\n",
    "evecs = np.array([[[dti_fit.evecs]]])\n",
    "\n",
    "color1 = np.zeros(evals.shape)\n",
    "color1[0, 0, 0, 0] = dki_fit.md/refmd;  #red\n",
    "\n",
    "scene = window.Scene()\n",
    "\n",
    "tensor_actor = actor.tensor_slicer(evals, evecs, scalar_colors=color1, sphere=sphere, scale=dki_fit.md/refmd)\n",
    "tensor_actor.display(y=0)\n",
    "tensor_actor.RotateX(100)\n",
    "scene.add(tensor_actor)\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(100)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "dir_actor2 = actor.peak_slicer(np.array([[[stick3]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor2.display(y=0)\n",
    "dir_actor2.RotateX(100)\n",
    "scene.add(dir_actor2)\n",
    "\n",
    "dtidir_actor = actor.peak_slicer(np.array([[[dki_fit.directions[0]]]]), np.array([[[1.2]]]), linewidth=3, colors=(1, 0, 0))\n",
    "dtidir_actor.display(y=0)\n",
    "dtidir_actor.RotateX(100)\n",
    "scene.add(dtidir_actor)\n",
    "\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4_dt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot DKI KT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "akc = dki_fit.akc(sphere)\n",
    "\n",
    "scene = window.Scene()\n",
    "dki_actor = actor.odf_slicer(np.array([[[akc]]]), sphere=sphere, colormap='Wistia', scale=dki_fit.mk()/refmk)\n",
    "dki_actor.display(y=0)\n",
    "dki_actor.RotateX(100)\n",
    "scene.add(dki_actor)\n",
    "\n",
    "\n",
    "dir_actor1 = actor.peak_slicer(np.array([[[stick1]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor1.display(y=0)\n",
    "dir_actor1.RotateX(100)\n",
    "scene.add(dir_actor1)\n",
    "\n",
    "dir_actor2 = actor.peak_slicer(np.array([[[stick3]]]), np.array([[[1.2]]]), linewidth=4, colors=(0, 0, 1))\n",
    "dir_actor2.display(y=0)\n",
    "dir_actor2.RotateX(100)\n",
    "scene.add(dir_actor2)\n",
    "\n",
    "dtidir_actor = actor.peak_slicer(np.array([[[dki_fit.directions[0]]]]), np.array([[[1.2]]]), linewidth=3, colors=(1, 0, 0))\n",
    "dtidir_actor.display(y=0)\n",
    "dtidir_actor.RotateX(100)\n",
    "scene.add(dtidir_actor)\n",
    "\n",
    "window.show(scene)\n",
    "\n",
    "window.record(scene, n_frames=1, out_path='case4_kt_dki_fit.png', size=(600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "dwi = np.zeros((4,) + gtab.bvals.shape)\n",
    "dwi[0, :] = signal1\n",
    "dwi[1, :] = signal2\n",
    "dwi[2, :] = signal3\n",
    "dwi[3, :] = signal4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define reconstruction models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "dki_model = dki.DiffusionKurtosisModel(gtab, fit_method='NLS')\n",
    "dki_fit = dki_model.fit(dwi)\n",
    "fa = dki_fit.fa\n",
    "md = dki_fit.md\n",
    "rd = dki_fit.rd\n",
    "ad = dki_fit.ad\n",
    "mk = dki_fit.mk()\n",
    "rk = dki_fit.rk()\n",
    "ak = dki_fit.ak()\n",
    "mkt = dki_fit.mkt()\n",
    "kfa = dki_fit.kfa\n",
    "\n",
    "msdki_model = msdki.MeanDiffusionKurtosisModel(gtab)\n",
    "msdki_fit = msdki_model.fit(dwi)\n",
    "msk = msdki_fit.msk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAANiCAYAAABrea+DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABo0ElEQVR4nO3dfZxkZX3n/c/XQTBREONM1ADjkDhsHBVROmhWIySgARMHE59mEiPkRicmTmLiwwZXFw1mNz5F7+wtrk6UVUkUkUTvVkeJt2hMVMg0EdEZghmRyKAuIyDG+ICjv/uPOo1l091Vw1Sf6urzeb9e9Zo651xV59fV11R/69R1rpOqQpIkSeqau4y7AEmSJGkcDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgPMGS3CfJJ5P8fZJLk9yvhX2ekORTST6e5J1J7rrU+9TKZP/VSmA/liZbvKDG5EqyCqiq+kGSM4Ejq+pPl3if9wO+XlXfTvJnwBVVdfFS7lMrk/1XK4H9WJpsE3tEOMk/JXnQmGv4syR/OK79V9X3q+oHzeKhwM4W9vmVqvp2s3gb8IPF2s9nOfzulotxvhb2X/vvqNiP7cfSxKqqZX0DPgbcAhwyZ/1Tgb9Z5HHX0XuDWD1n/aeBAtbNaftt4N+BrwOfBJ4N3GWR518D3AD82BA/w1ZgBvgu8NZ5tv8E8B7gP4B/A35jP16f44DLgWuA+7f4e7k/8CngrnPW/xSwZ8BjF/3draRbX9/6ZtOPPwAcNexrMefxXwXeCtxjgTb23+H3a//d/9fsN5p+8E3gK8AHgUfvz2vCAu/nzTb78f7/TuzH3gD+HPjdcdcxqbexFzDgl7sO+D5wM/CUOdvu1qy/7wKPva55U/r9vnUPadbNF4RPae7fE9gIfBH434vU9kLgL4f8OX4deCLwvxZ4A34n8C7gHsCjgVuBBzXb7tv88Zh7u++c53gq8MaWfi+HAR8H/tM8254JvHnA4xf93a2k25y+dTfgfOC9w74Wcx5/X+AzwH9fpI391/67FK/Z84Abm75wd+CuwBOAVw/7mrDI+3mz3X68f78T+3EHbgzxARG4H3A9cPACz7GOXu75Zt/tM33bz2y2P23cP+9YXuNxFzCgA5wDfAJ4LfD+ebZ/GDhjkc7zEmBH37rXAC9mkSDct+4Eel83PXiB578UeHrf8qv40YDzauAj/R0T+NO5b8D0/qjcBhzTt+4C4BVDvD79z/3LwGv7ls8CLmne9G8BPg9sAJ4LfAn4GvDrd6LtQcB24OQFavrbvrY/Te8I6NeAbwAfHuZ3t5Juc/sW8Hjg8/vZj/sf/yrgA4u1sf/af0fch+9J7w/nHcLrsP242T7o/dx+bD/2dsff5e3/L1jkA2Lzu3zyAs+xjl7uOWiB7R8FbmLO35au3Jb7GOFnAH/d3H45yX3mbL8aeOgij78MOCzJA5sTGjYBfzXMjqvqn4A9wC8s0GT26PKsVwK/mORhSZ4NnErvjei2Abs6BthXVZ/vW/cZYJixW8c1Zw1/FPhDem/6sx4KTAEXA6uBz9L7KhPgZ4CX0/ugsL9tNwOPAP5bko8ledrsDpszlx9D7z8kwNvpvVnfB/hJ4GV99Q363a04SX4ceBq9ftlvqNciyZHAacDuQW3tv/bfEfp5ekcP3zOg3aDXZND7+R3Yj+3H+qGqurWqpun9HTkjyYP7Nn8M+JX9fc4k9wdOBLbQ+39531HUOkkOGncBC0nyaHrjny6qqq8l+QK9MWqv62v27/S+EljMBfTegP+e3n/6G/ajjC/TGzc2n8Ob/QNQVTcleR3wNnqf2h5dVbcOsY970PuU3u9WeiddLKr5I/GYBTY/FPizqvoIQJJd9Mbl/UWz/Dl++Psfum1VXUDvNZ3PY+h93TL7uvwMsApYVVXfoXc0aNYwv7uV4r1J9tE76rSX3lGjfoNei/cmKXp95VLgpUPu1/6L/XcE7g18rar2DWi34Gsy5Pv5QuzH2I/1Q1X1T0lmPyB+rll9NfCkO/F0zwBmqupvklwN/Ca9McedsZyPCJ8B/F1Vfa1Zfkezrt+h9MbMLOYCem+4Z9L7ZLw/jqA3hmo+t3DHN8lP0ztC8aKqun7IfXyT3livfofR9+Z+Jx0LvL9vecM8y/9yJ9ou5vH0jjzM+k3gdODLSd6SpP+P2TC/u5XiiVV1OL2jaluBv5/zqXvQa/HEqjoUOAn4WXpHi4Zh/71j28XYf+d3E7A6yaADJ4u9JsO8ny/EfnzHtouxH3fD3A+I/07vg+Fivpbk683tBc26Z9D7/0jz7zNGWuUEWJZBOMmP0Tvp4MQkX03yVeCPgIcm6f8a54H0vr5aUFX9G73xNI+nN25q2Bp+jt4b8D8u0OQqel+nzbZ/CL2xXW8D/q9h90NvHNhBSdb3rXsoBzAFT/NVx8HNc886Driyb/lY4Mr9aTvErn/kDbiqLq2qk+m9gT+U3oeRWQN/dytN9aZZ+lt6Jww9um/TUK9FVf09vVkjXjOorf3X/jtCn6I308ITB7Sb9zXZj/fzO7Af24+1oLkfEIf5ULO6qg5vbq9J8ijgaODCZvs7gIckOW7UxS5nyzII03vD/T69/7jHNbcHAv9A82klyd2A4/nhOKjFnAX8UlX9x6CGSQ5L8qv0OsZfVdVnF2i6nd64GpIcAbyP3pmcv0evI53U95wHNfWuAlYludvs0ZWmpr8Fzk1y96Zjns7CX3sN46HAZ6uZ2zLJYfS+lrxqTpvP7GfbBSU5mt7XeFc3y7+eZH2S0PsPei+aN/H9/N2tGOk5nd5rMfs67e9r8X8Dj10oQNh/7b+j1gwtOAc4L8kTk/x4krsmOS3Jq2Dga/JEBryfz2U/th9rYQt8QLwzH2rOAELvw9hX6U0BOLu+O8Z9tt58N+BDwJ/Ps/6p9OZSPQh4CvC3izzHdcw5A7lZfxCLzyN8K70jIM+hN6ZqoedfTe8kjnvS63x/0LftBcAn+pZf1uyz//ayvu0/AbyX3vyVX2I/5q9coLb/Rt8UPvSOPv5L3/JdgG/R+480dNsB+9wKvL5v+bX0vrr5Jr2TWZ7Vt23R391KuvGj8wD/O73xXL857GsxXz+md8Trb+bZh/3X/ruUffk36c3D+x/03oc/APznQa8JQ7yf24/tx94W/J1exw9njTgM+FXgC8Db57T7O+CpCzzHOubMGkFvqN7X6R0ovG/f7TnA/2GBGSZW4m1iL7Gc5HLgrKr63MDGS1fD/wBurKr/e1w1LBdJttN7A94+RNux/+6Wi3G+FvbfH7L/Hhj78fJgP155klxHb8aPffSmEtxFb/arN1bV95s29wN2AD9d88yQkmQdvSGid63mpNckm+idrLq2qr7X1/bH6H24PKOq3j/3uVaiiQ3CWl6S/Bfg/6kfXvZTmhj2X60E9uNuSvLnwBeq6g3jrmUSGYQlSZLUSQNPlktyfpIb05vDcL7tSfI/k+xOclWSh4++TEmSJGm0hpk14q30rs6zkNOA9c1tC70TeSRJkqRlbWAQrqqPs/Bk5tCbYubt1XMZcHgzcFuSJElatkZxieUjgP6r9+xp1n1lsQetXr261q1bN4Lda6W64oorvlZVa8Zdx2Lsx1qMfVgrgf1YK8FC/XgUQXhoSbbQGz7B2rVrmZmZaXP3mjBJ/m3cNQyybt06+7EWZB/WSmA/1kqwUD8exZXlbgCO6ls+sll3B1W1raqmqmpqzZpl/eFSkiRJK9wogvA08Ixm9ohHArdW1aLDIiRJkqRxGzg0Isk7gZOA1Un2AC8F7gpQVW+kd633xwO76V0C8reXqlhJkiRpVAYG4araPGB70bs2tSRJkjQxRjE0QpIkSZo4BmFJktSaJKcmuaa5Iu3Z82xfm+SjST7dXLH28eOoU91gEJYkSa1Isgo4j95VaTcAm5NsmNPsJcBFVfUwYBPwhnarVJcYhCVJUltOAHZX1bVVdRtwIb0r1PYr4LDm/j2BL7dYnzqm1QtqSJKkTpvvarSPmNPmZcDfJfl94O7AKe2Upi4yCC9TyRWt7Kfq+Fb2o+5JfqeV/VS9qZX9qJtyRUvvxcf7XtxnM/DWqvrzJD8PXJDkwVX1g/5Gc69WO69kiUttVC28bTnUoAU5NEKSJLVlmKvRngVcBFBVnwLuBqye+0RerVajYBCWJElt2QGsT3J0koPpnQw3PafNl4CTAZI8kF4Q3ttqleoMg7AkSWpFVe0DtgKXAFfTmx1iZ5Jzk2xsmj0feFaSzwDvBM5sLt4ljZxjhCVJUmuqajuwfc66c/ru7wIe1XZd6iaPCKsThpjA/cwke5Nc2dyeOY46pSTnJ7kxyecW2P7Cvn76uSTfT/ITzbbrkny22TbTbuWSNHkMwlrxhpzAHeBdVXVcc3tzq0VKP/RW4NSFNlbVq2f7KfAi4O+r6ua+Jr/YbJ9a2jIlafIZhNUFw0zgLi0LVfVx4OaBDXs20xtDKUm6EwzC6oL5JnA/Yp52T2qua39xkqPm2Q705q5MMpNkZu9eT2TWeCT5cXpHjv+mb3XRuxDBFc0cqws91j4sSRiEpVnvA9ZV1bHAh4G3LdTQuSu1TDwB+MScYRGPrqqH0xsG9Jwkj5nvgfZhSeoxCKsLBk7gXlU3VdV3m8U3A17mScvdJuYMi6iqG5p/bwTeQ29YkCRpAQZhdcHACdyT3K9vcSO9+S2lZSnJPYETgf+3b93dkxw6ex94HDDvzBOSpB7nEdaKV1X7ksxO4L4KOH92AndgpqqmgT9oJnPfR+9EpTPHVrA6Lck7gZOA1Un2AC8F7gpQVW9smv0a8HdV9R99D70P8J4k0Htvf0dVfaituiVpEhmE1QlDTOD+InpTUUljVVWbh2jzVnrTrPWvuxZ46NJUJUkrk0FYkjSvXHFFK/up4x2SL2k8DMKSpGXJIC5pqXmynCRJkjrJICxJkqROMghLkiSpkwzCkiRJ6iSDsCRJkjrJICxJkqROMghLkiSpk5xHWNKKlfxOK/upelMr+5FWgiSnAn9B75L3b66qV8zZ/jrgF5vFHwd+sqoOb7VIdYZBWJIktSLJKuA84LHAHmBHkumq2jXbpqr+qK/97wMPa71QdYZDIyRJUltOAHZX1bVVdRtwIXD6Iu03A+9spTJ1kkFYkiS15Qjg+r7lPc26O0hyf+Bo4NIW6lJHGYQlSdJytAm4uKq+P9/GJFuSzCSZ2bt3b8ulaaUwCEuSpLbcABzVt3xks24+m1hkWERVbauqqaqaWrNmzQhLVJcYhCVJUlt2AOuTHJ3kYHphd3puoyQ/C9wL+FTL9aljDMKSJKkVVbUP2ApcAlwNXFRVO5Ocm2RjX9NNwIVVVeOoU90x1PRpQ8z5txZ4G3B40+bsqto+2lIlSdKka/LB9jnrzpmz/LI2a1J3DTwi3Dfn32nABmBzkg1zmr2E3qe6h9H7FPeGURcqSZIkjdIwQyOGmfOvgMOa+/cEvjy6EiVJkqTRGyYIDzPn38uApyfZQ+/rjt8fSXWS1DFJzk9yY5LPLbD9pCS3JrmyuZ3Tt+3UJNck2Z3k7PaqlqTJNKqT5TYDb62qI4HHAxckucNzO+efJA30VuDUAW3+oaqOa27nwtDD2CRJfYYJwsPM+XcWcBFAVX0KuBuweu4TOeefJC2uqj4O3HwnHrq/l66VpM4bJggPM+ffl4CTAZI8kF4Q9pCvJC2Nn0/ymSQfTPKgZt3Ql66V1DFJO7cJNHD6tKral2R2zr9VwPmzc/4BM1U1DTwf+Mskf0TvxLkznftPkpbEPwP3r6pvJnk88F5g/f48QZItwBaAtWvXjrxASZoUQ80jPGjOv6raBTxqtKVJkuaqqm/03d+e5A1JVrMfl66tqm3ANoCpqSkPWkjqLK8sJ0kTJMl9k953kElOoPc+fhNDXrpWkvRDQx0RliS1I8k7gZOA1c2UlC8F7gpQVW8Engz8bpJ9wLeBTc1QtHmHsY3hR5CkibHsgnByRSv7qTq+lf1I0v6oqs0Dtr8eeP0C2+4wjE2StLBlF4SXA8O4JEnSyucYYUmSJHWSQViSJEmdZBBWZyQ5Nck1SXYnOXuRdk9KUkmm2qxPkiS1yyCsTkiyCjgPOA3YAGxOsmGedocCzwUub7dCSZLUNk+W07zaOmEQWjtp8ARgd1VdC5DkQuB0YNecdi8HXgm8sI2iJEnS+HhEWF1xBHB93/KeZt3tkjwcOKqqPtBmYZIkaTwMwhKQ5C7Aa4HnD9F2S5KZJDN79+5d+uIkaQUZ5nyNJE9NsivJziTvaLtGdYdBWF1xA3BU3/KRzbpZhwIPBj6W5DrgkcD0fCfMVdW2qpqqqqk1a9YsYcmStLIMc75GkvXAi4BHVdWDgD9su051h0FYXbEDWJ/k6CQHA5uA6dmNVXVrVa2uqnVVtQ64DNhYVTPjKVeSVqTbz9eoqtuA2fM1+j0LOK+qbgGoqhtbrlEdYhBWJ1TVPmArcAlwNXBRVe1Mcm6SjeOtTpI6Y+D5GsAxwDFJPpHksiSntladOsdZI9QZVbUd2D5n3TkLtD2pjZokSXdwELAeOIneMLaPJ3lIVX29v1GSLcAWgLVr17ZcolYKjwhLkqS2DDpfA3pHiaer6ntV9UXg8/SC8Y/wfI0JkrRzuxMMwpIkqS2Lnq/ReC+9o8EkWU1vqMS1LdaoDjEIS5KkVgx5vsYlwE1JdgEfBV5YVTeNp2KtdI4RliRJrRl0vkZVFfC85iYtKY8IS5IkqZMMwpIkSeokg7AkSZI6yTHCkiQtY7niilb2U8cf38p+pOXEI8KSJEnqJIOwJC0jSc5PcmOSzy2w/TeTXJXks0k+meShfduua9ZfmWSmvaolaTIZhCVpeXkrcOoi278InFhVDwFeDmybs/0Xq+q4qppaovokacVwjLAkLSNV9fEk6xbZ/sm+xcvoXaJWknQneERYkibXWcAH+5YL+LskVyTZMqaaJGlieERYkiZQkl+kF4Qf3bf60VV1Q5KfBD6c5F+q6uPzPHYLsAVg7dq1rdQrScuRR4QlacIkORZ4M3B6Vd00u76qbmj+vRF4D3DCfI+vqm1VNVVVU2vWrGmjZElalgzCkjRBkqwF/hb4rar6fN/6uyc5dPY+8Dhg3pknJEk9Do2QpGUkyTuBk4DVSfYALwXuClBVbwTOAe4NvCEJwL5mhoj7AO9p1h0EvKOqPtT6DyBJE8QgLEnLSFVtHrD9mcAz51l/LfDQOz5CkrQQh0ZIkiSpkwzCkiRJ6iSDsCRJkjppqCCc5NQk1yTZneTsBdo8NcmuJDuTvGO0ZUqSJEmjNTAIJ1kFnAecBmwANifZMKfNeuBFwKOq6kHAH46+VEmSNOkGHVxLcmaSvUmubG53ODlUGpVhZo04AdjdnJFMkguB04FdfW2eBZxXVbfA7ZO5S5Ik3a7v4NpjgT3AjiTTVbVrTtN3VdXW1gtU5wwzNOII4Pq+5T3Nun7HAMck+USSy5KcOqoCJUnSinH7wbWqug2YPbgmjcWoTpY7CFhPbxL4zcBfJjl8bqMkW5LMJJnZu3fviHYtSZImxDAH1wCelOSqJBcnOaqd0tRFwwThG4D+Tnhks67fHmC6qr5XVV8EPk8vGP8Ir28vSZIGeB+wrqqOBT4MvG2+Rh5c0ygME4R3AOuTHJ3kYGATMD2nzXvpHQ0myWp6QyWuHV2ZkiRpBRh4cK2qbqqq7zaLbwaOn++JPLimURgYhKtqH7AVuAS4GrioqnYmOTfJxqbZJcBNSXYBHwVeWFU3LVXRkiRpIg08uJbkfn2LG+llD2lJDDNrBFW1Hdg+Z905ffcLeF5zkyRJuoOq2pdk9uDaKuD82YNrwExVTQN/0Bxo2wfcDJw5toK14g0VhCVJkkZhiINrL6J3bQJpyXmJZUmSJHWSQViSJEmdZBBWZwxxWc9nJ/lsc0nPf5x7KXFJkrSyGITVCX2X9TwN2ABsnifovqOqHlJVxwGvAl7bbpWSJKlNBmF1xcDLelbVN/oW7w5Ui/VJkqSWOWuEumK+y3o+Ym6jJM+hNw3gwcAvzfdESbYAWwDWrl07786S3zmwaodU9aZW9iNJ0kpkEJb6VNV5wHlJfgN4CXDGPG22AdsApqamlu1RY8O4JEmLc2iEumLgZT3nuBB44lIWJEmSxssgrK4Y5rKe6/sWfwX41xbrkwBIcn6SG5N8boHtSfI/m9lPrkry8L5tZyT51+Z2h28zJEk/yqER6oQhL+u5NckpwPeAW5hnWITUgrcCrwfevsD204D1ze0RwP8CHpHkJ4CXAlP0TvS8Isl0Vd2y5BVL0oQyCKszhris53NbL0qao6o+nmTdIk1OB95eVQVcluTwJPcDTgI+XFU3AyT5MHAq8M4lLlmSJpZDIyRpssw3A8oRi6y/gyRbkswkmdm7d++SFSpJy51BWJI6pqq2VdVUVU2tWbNm3OVI0tgYhCVpsiw0A8r+zowiSZ1nEJakyTINPKOZPeKRwK1V9RV6J4I+Lsm9ktwLeFyzTpK0AE+Wk6RlJMk76Z34tjrJHnozQdwVoKreSO+Ez8cDu4FvAb/dbLs5ycvpTRUIcO7siXOSpPkZhCVpGamqzQO2F/CcBbadD5y/FHVJ0krk0AhJkiR1kkFYkiS1JsmpSa5pro549iLtnpSkkky1WZ+6xSAsSZJakWQVcB69KyRuADYn2TBPu0OB5wKXt1uhusYgLEmS2nICsLuqrq2q24AL6V0tca6XA68EvtNmceoeg7AkSWrLwCsgJnk4cFRVfaDNwtRNBmFJkrQsJLkL8Frg+UO09VLhOmAGYUmS1JZBV0A8FHgw8LEk1wGPBKbnO2HOS4VrFAzCkiSpLTuA9UmOTnIwsIne1RIBqKpbq2p1Va2rqnXAZcDGqpoZT7la6QzCkiSpFVW1D9hK7/LfVwMXVdXOJOcm2Tje6tRFXllOkiS1pqq207tUeP+6cxZoe1IbNam7PCIsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkThoqCCc5Nck1SXYnOXuRdk9KUkmmRleiJEmSNHoDg3CSVcB5wGnABmBzkg3ztDsUeC5w+aiLlKQuGXTwIcnrklzZ3D6f5Ot9277ft2261cIlacIcNESbE4DdVXUtQJILgdOBXXPavRx4JfDCkVYoSR3Sd/DhscAeYEeS6aq6/T23qv6or/3vAw/re4pvV9VxLZUrSRNtmKERRwDX9y3vadbdLsnDgaOq6gOLPVGSLUlmkszs3bt3v4uVpA64/eBDVd0GzB58WMhm4J2tVCZJK8wBnyyX5C7Aa4HnD2pbVduqaqqqptasWXOgu5aklWjgwYdZSe4PHA1c2rf6bs0Bh8uSPHHJqpSkFWCYoRE3AEf1LR/ZrJt1KPBg4GNJAO4LTCfZWFUzoypUknQHm4CLq+r7fevuX1U3JPlp4NIkn62qL/Q/KMkWYAvA2rVr26tWkpaZYY4I7wDWJzk6ycH03nhvPwGjqm6tqtVVta6q1gGXAYZgSbpzBh186LeJOcMiquqG5t9rgY/xo+OHZ9v47ZwkMUQQrqp9wFbgEuBq4KKq2pnk3CQbl7pAaRSGOAv/eUl2JbkqyUear5ylcVj04MOsJD8L3Av4VN+6eyU5pLm/GngUdzyxWZLUGGZoBFW1Hdg+Z905C7Q96cDLkkZnmLPwgU8DU1X1rSS/C7wKeFr71arrqmpfktmDD6uA82cPPgAzVTUbijcBF1ZV9T38gcCbkvyA3oGOV8zp55KkPkMFYWnCDZwCsKo+2tf+MuDprVYo9Rnm4ENVvWyex30SeMiSFicdoCSnAn9B74Pem6vqFXO2Pxt4DvB94JvAFj/Qaal4iWV1wdBn4TfOAj64pBVJUgcNeZGud1TVQ5r5sF9Fb2YqaUkYhKU+SZ4OTAGvXqSN82FL0p0zcJ7sqvpG3+Ldgf7hP9JIGYTVBUOdhZ/kFODF9GY9+e5CT+YZ95J0pw31DV2S5yT5Ar0jwn/QUm3qIIOwumDgWfhJHga8iV4IvnEMNUqSGlV1XlX9DPDHwEvma+O3cxoFg7BWvCGnAHw1cA/g3UmuTHKH6aokSQdsf+bJht7QiSfOt8Fv5zQKzhqhThh0Fn5VndJ6UZLUPbd/Q0cvAG8CfqO/QZL1VfWvzeKvAP+KtEQMwpIkqRVDzpO9tTln43vALcAZ46tYK51BWJIktWaIb+ie23pR6izHCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKStMwkOTXJNUl2Jzl7nu1nJtmb5Mrm9sy+bWck+dfmdka7lUvSZDlo3AVIkn4oySrgPOCxwB5gR5Lpqto1p+m7qmrrnMf+BPBSYAoo4Irmsbe0ULokTRyPCEvS8nICsLuqrq2q24ALgdOHfOwvAx+uqpub8Pth4NQlqlOSJp5BWJKWlyOA6/uW9zTr5npSkquSXJzkqP18rCQJg7AkTaL3Aeuq6lh6R33ftj8PTrIlyUySmb179y5JgZI0CQzCkrS83AAc1bd8ZLPudlV1U1V9t1l8M3D8sI9tHr+tqqaqamrNmjUjK1waxhAngz4vya7mG4+PJLn/OOpUNxiEJWl52QGsT3J0koOBTcB0f4Mk9+tb3Ahc3dy/BHhcknsluRfwuGadtCz0nQx6GrAB2Jxkw5xmnwammm88LgZe1W6V6hJnjZCkZaSq9iXZSi/ArgLOr6qdSc4FZqpqGviDJBuBfcDNwJnNY29O8nJ6YRrg3Kq6ufUfQlrY7SeDAiSZPRn09llRquqjfe0vA57eaoXqFIOwJC0zVbUd2D5n3Tl9918EvGiBx54PnL+kBUp33nwndD5ikfZnAR9c0orUaQZhSZK07CR5Or05sU9cYPsWYAvA2rVrW6xMK4ljhCVJUluGOqEzySnAi4GNfSeG/ghP+tQoGIQlSVJbhjkZ9GHAm+iF4BvHUKM6xCAsSZJaUVX7gNmTQa8GLpo9GbQ5ARTg1cA9gHcnuTLJ9AJPJx0wxwhLkqTWDHEy6CmtF6XOGuqIsJNfS5IkaaUZGISd/FqSJEkr0TBHhG+f/LqqbgNmJ7++XVV9tKq+1SxeRu8sUEmSJGnZGiYIzzf59RGLtHfya0mSJC17I501om/y61cvsH1LkpkkM3v37h3lrqVFDTHO/TFJ/jnJviRPHkeNkiSpXcMEYSe/1kQbcpz7l4AzgXe0W50kSRqXYYKwk19r0g0zzv26qroK+ME4CpQkSe0bGISd/ForwP6Oc1+UQ3wkSVoZhrqghpNfSz9UVduAbQBTU1M15nIkSdKd5CWW1QVDjXOXJEndYhBWFwwc5y5JkrrHIKwVb5hx7kl+Lske4CnAm5LsHF/FkiSpDUONEZYm3RDj3HfgFRElSeoUjwhLkiSpkwzCkiRJ6iSDsCQtM0NcEvx5SXYluSrJR5Lcv2/b95v53J3TXZIGcIywJC0jfZcEfyy9i7/sSDJdVbv6mn0amKqqbyX5XeBVwNOabd+uquParFmSJpVHhCVpeRnmkuAfrapvNYuX4YmeknSnGIQlaXnZ30uCnwV8sG/5bs0lwC9L8sQlqE86IEMM/XlMkn9Osi/Jk8dRo7rDoRGSNKGSPB2YAk7sW33/qrohyU8Dlyb5bFV9Yc7jtgBbANauXdtavdKQQ3++BJwJvKD9CtU1HhGWpOVlqEuCJzkFeDGwsaq+O7u+qm5o/r0W+BjwsLmPraptVTVVVVNr1qwZbfXS4oYZ+nNdVV0F/GAcBapbDMKStLwMvCR4kocBb6IXgm/sW3+vJIc091cDjwL6j7RJ47a/Q3+kJeXQCElaRqpqX5LZS4KvAs6fvSQ4MFNV08CrgXsA704C8KWq2gg8kN4lwn9A70DHK+Z85SytGA7x0SgYhCVpmRnikuCnLPC4TwIPWdrqpAMy1NCfYVTVNmAbwNTUVB14aeoih0ZIkqS2DBz6I7XJICxJklpRVfuA2aE/VwMXzQ79SbIRIMnPJdkDPIXeUJ+d46tYK51DIyRJUmuGGPqzAy8So5Z4RFiSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJw0VhJOcmuSaJLuTnD3P9kOSvKvZfnmSdSOvVDpA9mNNigPpq0le1Ky/Jskvt1q4NATfi7WcDAzCSVYB5wGnARuAzUk2zGl2FnBLVT0AeB3wylEXKh0I+7EmxYH01abdJuBBwKnAG5rnk5YF34u13AxzRPgEYHdVXVtVtwEXAqfPaXM68Lbm/sXAyUkyujKlA2Y/1qQ4kL56OnBhVX23qr4I7G6eT1oufC/WsjJMED4CuL5veU+zbt42VbUPuBW49ygKlEbEfqxJcSB9dZjHSuPke7GWlYPa3FmSLcCWZvGbSa4Z0VOvBr62f7WMaM/WcED7H1DD/Q+kmKWyvPrxthHtetnUMO79j7qGrvVhuDOv3wh3vgz2vxJr6Fo/3v+/ZaP/g24No9//vP14mCB8A3BU3/KRzbr52uxJchBwT+CmuU9UVduAkf/VSjJTVVOjfl5rmKz9D2A/noAaxr3/ZVLDgfTVYR67ZH0Yxv/6jXv/1jCQ78XWsKz2P8zQiB3A+iRHJzmY3okY03PaTANnNPefDFxaVTW6MqUDZj/WpDiQvjoNbGrOuj8aWA/8U0t1S8PwvVjLysAjwlW1L8lW4BJgFXB+Ve1Mci4wU1XTwFuAC5LsBm6m17GlZcN+rElxIH21aXcRsAvYBzynqr4/lh9EmofvxVpushI+ZCXZ0nxFYg1jrGHc+590y+H1G3cN497/cqlhko379Rv3/q1h8i2H184a2tv/igjCkiRJ0v7yEsuSJEnqpIkLwknum+TCJF9IckWS7UmOSfKhJF9P8v4x1XBCkk8l2ZnkqiRPG0MNJyb55yRXNnU8u+X9H9NsOyzJniSvX6r9Tzr78fj78CI12I+HNO5+PO4+vEgN9uMJMe4+vEgNnerHY+3DVTUxN3rTHH4KeHbfuocCvwCcDDwBeP+YajgRWN8s/xTwFeDwMdRwSLN8D+A64Kfa/D009/8CeAfw+nH3meV4sx+Pvw8P+j009+3Hd/L1a6Mfj7sPD1GD/XiZ38bdh4foQ53ox+Puw61eUGMEfhH4XlW9cXZFVX1m9n6Sk8ZdQ7P85SQ3AmuAr4+jBuAQlu6I/4L7T3I8cB/gQ8BynMNyObAfj78PL1qD/Xgo4+7H4+7DQ9WA/Xg5G3cfHlhDs7zS+/FY+/CkDY14MHDFcq8hyQnAwcAX2q4hyVFJrqJ3ecpXVtWX29p/krsAfw68YAn2uZLYj8ffhxeswX48tHH343H34UVrsB9PhHH34aFq6EA/HmsfnrQgvOwluR9wAfDbVfWDtvdfVddX1bHAA4Azktynxd3/HrC9qva0uE8tgXH24zH3YbAfrwgdfy8G+/GK0PF+3EofnrQgvBM4frnWkOQw4APAi6vqsnHUMKv51PY5emOd2tr/zwNbk1wHvAZ4RpJXLMH+J539ePx9eLEa7MfDGXc/HncfXrSGWfbjZW3cfXjRGjrUj8fahyctCF8KHJJky+yKJMcmWao3mP2p4UTgPcDbq+riMdXwC0l+rFm+F/Bo4Jq29g+8sarWVtU6el9lvL2qzl6C/U86+/H4+/CCNWA/Hta4+/G4+/BiNdiPJ8O4+/BiNXSpH4+1D09UEK7e6YO/BpzSTLGxE/gz4KtJ/gF4N3ByM83GL7dcw2Oa25nNVCNXJjmu5Rp+Frg8yWeAvwdeU1WfbXH/Xx31vlYi+/H4+/CAGuzHQxh3Px53Hx5Qg/14Aoy7Dw+ooTP9eNx92CvLSZIkqZMm6oiwJEmSNCoGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEW5DkCUneNe46pAOR5D5Jrk5yyLhrkQ5Ukj9P8rvjrkPaX0n+KcmDxl3HSmEQHpEk1yX5dpJv9t1eD1BV7wMelOTYAY+/LcnqOes/naSSrGuW35rkT/u2PyjJV5K8YGl+MnXNAn35p6rq/wAfBbYs8ti3Nv34m0luTvLhJD/bt/3MJP/Yt3xYkk8k+ZskBy/tT6YuSbIpyeVJ/iPJjc3930uSpslrgP86qN8luUfTnz84z7brkpwyZ5+3JDlxtD+NumJun2rW/cj7Jr2+e+4iz7Ho++yc9+nZ29OatklybZJdo/7ZliuD8Gg9oaru0Xfb2rftnSwSIBpfBDbPLiR5CPDjCzVO8jB6weRPq+o1B1C3NNfcvvzlZv1fA78z4LGvqqp7AEcANwBvma9RknsBHwH+DXhaVd02otrVcUmeD/wF8GrgvsB9gGcDjwIOBqiqrwD/Amwc8HRPAr4LPDbJfRfZ5xnAecCvVNXfH+jPIC1iGvjFxfrjrEXeZ1815z1+9lvrxwA/Cfx0kp9biuKXG4Nwez4G/MqANhcAz+hbPgN4+3wNk5wAfBj4r1V13igKlIZwOb03yPsPalhV3wYuAo6buy3JGnof4j4HPL2q9o24TnVUknvSO1r2e1V1cVX9e/V8uqp+s6q+29f8Ywx+Xz4DeCNwFfD0Bfb5O8CfA79cVZ884B9CWkRVfQe4AvjlxdrdyffZM4D/F9je3F/xDMLtuRpYl+SwRdpcBhyW5IFJVgGbgL+ap90JwIeAP6qqN4++VGl+zRvpbuChg9omuTu9bzh2z9n0E/QCyKeA/6uqfjDiMtVtPw8cQu+P+SBXs0hfbj7wnUTvm5C/5kcPVMz6XXrB++SqmtnfYqU7adG+y514n03y48CT+WF/39SFIWsG4dF6b5Kv992e1bft35t/Dx/wHLNHhR9Lr6PfME+bRwK3AncYsyaNSH9ffu+cbf/O4v34BUm+3rR7NPBbc7YfBRwDvLWqajTlSrdbDXyt/+hXkk82ffnbSR7T13ZQX/4t4Kqq2gVcSO9cj4fNafNYegcxPjuS6qU5WQJ4wzxtBvXdQe+zL+jbx9eadb9ObxjQ3wEfAO7K4G9MJp5BeLSeWFWH993+sm/boc2/Xx/wHBcAvwGcyQLDIuiNQ5sBPtyM/5FGrb8vP3HOtkNZvB+/pqoOB9YB3wb+05ztnwFeAHxwnlAhHaibgNVJDppdUVX/uemTN/Gjf/cG9eVn0DsyRlXdAPw9d/y6+HfpBY43952IJx2IH8kSwO/N02ZQ3x30Pvuavn3MnqR/BnBRVe1rhl/8DR0YHmEQbs8Dgeuq6huLNaqqf6N30tzjgb9doNn36YXlLwGXDBhuIY1MEy4eQO9NdlFV9SXgucBfJPmxOdv+AngFvQ9zD16KWtVZn6J3VOv0Ido+kAX6cpL/DKwHXpTkq0m+CjwC+I3+kA38H+Bk4BeY/8idtBQW7Luz9ud9NsmRwC8BT+/r708GHj93NquVxiDcnhMZfijDWcAvVdV/LNSgqr4HPAX4GrC9GY8pLbUT6H2g+7dhGlfVh4EvM8+MKVX1Knpn9v9/SeYeNZbulKr6OvAnwBuSPDnJoUnukuQ4YO775GLvy2fQOyF5A70TPo8DHgz8GHDanH1+mV4YPjXJ60byg0gLSHI34Hh6/XNR+/E++1vA5+l9g3dcczsG2EPfbFYrkUF4tN43Z16+9/Rt2wy8aZgnqaovDHPSRTMNyq8D32n2/WMDHiIdqN+kdwb9/ng18F8yz4U4qurlwJuBjyT5mRHUJ83+8X8e8F/oHbH9P/Tef/8Y+CRAkvvRC7nvnfv4Jmg8Ffh/quqrfbcv0hu+doevi5tvQH4JeHKSP1uKn0tqPAH4WN+0losa8n32DOANc/r7V+m936/o4RHxXJWll+QJwG9V1VPHXYt0ZyX5SXpjJB/WjB+TJlaSPwe+UFUOZ9BESXI5cFZVfW7ctawEBmFJkiR1kkMjJEmS1EkGYUmSJHWSQViSJEmdZBCWJElSJx00uMnSWL16da1bt25cu9cEuOKKK75WVWvGXcdi7MdajH1YK4H9WCvBQv14bEF43bp1zMwMnCpXHZZkqIs2jJP9WIuxD2slsB9rJVioHzs0QpIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ00tivLaXnLFWltX3V8tbavLsmftPM7rJf6+9PSuKKl96HjfQ+SllZLf4+4E3+PPCKsTkhyfpIbk3xuge2/meSqJJ9N8skkD227RkmS1C6DsLrircCpi2z/InBiVT0EeDmwrY2iJEnS+Dg0Qp1QVR9Psm6R7Z/sW7wMOHLJi5IkSWPlEWHpjs4CPjjuItRdSU5Nck2S3UnOnmf7IUne1Wy/fPZDXpK7JnlbM8Tn6iQvar14SZogBmGpT5JfpBeE/3iRNluSzCSZ2bt3b3vFqROSrALOA04DNgCbk2yY0+ws4JaqegDwOuCVzfqnAIc0Q3yOB35nsW9CJKnrDMJSI8mxwJuB06vqpoXaVdW2qpqqqqk1a9a0V6C64gRgd1VdW1W3ARcCp89pczrwtub+xcDJSQIUcPckBwE/BtwGfKOdsiVp8hiEJSDJWuBvgd+qqs+Pux512hHA9X3Le5p187apqn3ArcC96YXi/wC+AnwJeE1V3Tx3B36rIUk9niynTkjyTuAkYHWSPcBLgbsCVNUbgXPoBYk39A6ssa+qpsZTrXSnnQB8H/gp4F7APyT5/6rq2v5GVbWNZmaUqakpJ9GV1FkGYXVCVW0esP2ZwDNbKkdazA3AUX3LRzbr5muzpxkGcU/gJuA3gA9V1feAG5N8ApgCrkWSdAcOjZCk5WUHsD7J0UkOBjYB03PaTANnNPefDFxaVUVvOMQvASS5O/BI4F9aqVoa0hCzoqxN8tEkn24udPT4cdSpbhgqCA/qtE2bpybZlWRnkneMtkxJ6oZmzO9W4BLgauCiqtqZ5NwkG5tmbwHunWQ38Dxg9n35POAeSXbSC9T/u6quavcnkBY25KwoL6HX7x9G74PgG9qtUl0ycGhEX6d9LL2TNnYkma6qXX1t1gMvAh5VVbck+cmlKliSVrqq2g5sn7PunL7736E3Vdrcx31zvvXSMnL7rCgASWZnRdnV16aAw5r79wS+3GqF6pRhjggPM5XPs4DzquoWgKq6cbRlSpKkFWCYWVFeBjy9ObF5O/D77ZSmLhomCA/TaY8BjknyiSSXJTl1VAVKkqRO2Qy8taqOBB4PXJDkDnnFaQA1CqM6We4gYD296ak2A3+Z5PC5jey0kiR12jCzopwFXARQVZ8C7gasnvtEXtxIozBMEB6m0+4Bpqvqe1X1ReDz9ILxj7DTSpLUacPMivIl4GSAJA+kF4Q9eqYlMUwQHqbTvpfe0WCSrKY3VMJ5KyVJ0u2GnBXl+cCzknwGeCdwZjM9oDRyA2eNqKp9SWY77Srg/NlOC8xU1XSz7XFJdtG7qtELq+qmpSxckiRNniFmRdkFPKrtutRNQ11ZbohOW/TmsnzeSKuTJEmSlohXlpMkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSVpmkpya5Joku5OcPc/2Q5K8q9l+eZJ1zfrfTHJl3+0HSY5ru35JmhQGYUlaRpKsAs4DTgM2AJuTbJjT7Czglqp6APA64JUAVfXXVXVcVR0H/Bbwxaq6sq3aJWnSGIQlaXk5AdhdVddW1W3AhcDpc9qcDrytuX8xcHKSzGmzuXmsJGkBBmFJWl6OAK7vW97TrJu3TVXtA24F7j2nzdOAd863gyRbkswkmdm7d+9IipakSWQQlqQVJskjgG9V1efm215V26pqqqqm1qxZ03J1krR8GIQlaXm5ATiqb/nIZt28bZIcBNwTuKlv+yYWOBosSfohg7A6Icn5SW5MMu8RsvT8z+Ys/KuSPLztGqXGDmB9kqOTHEwv1E7PaTMNnNHcfzJwaVUVQJK7AE/F8cGSNNBB4y5AaslbgdcDb19g+2nA+ub2COB/Nf9KraqqfUm2ApcAq4Dzq2pnknOBmaqaBt4CXJBkN3AzvbA86zHA9VV1bdu1r0RXXDH3HMSlcfzx1cp+JP0og7A6oao+PjvX6gJOB97eHFW7LMnhSe5XVV9pp0Lph6pqO7B9zrpz+u5/B3jKAo/9GPDIpaxPklYKh0ZIPcOcqS9JklYQg7C0n5x6SpKklcEgLPUMc6Y+4NRTkiStFAZhqWcaeEYze8QjgVsdHyxJ0srmyXLqhCTvBE4CVifZA7wUuCtAVb2R3olJjwd2A98Cfns8lUqSpLYYhNUJVbV5wPYCntNSOZIkaRkYamhEklOTXNNcbODsebafmWRvkiub2zNHX6okSZI0OgOPCCdZBZwHPJbelFI7kkxX1a45Td9VVVuXoEZJkiRp5IY5InwCsLuqrq2q2+hdtvP0pS1LkiRJWlrDBOFhLzTwpCRXJbk4yVHzbHf+VUmSJC0bo5o+7X3Auqo6Fvgw8Lb5Gjn/qiRJkpaLYYLwwAsNVNVNVfXdZvHNwPGjKU+SJElaGsME4R3A+iRHJzkY2ETv4gO3S3K/vsWNwNWjK1GSJEkavYGzRlTVviRbgUuAVcD5VbUzybnATFVNA3+QZCOwD7gZOHMJa5YkSZIO2FAX1Kiq7fSuvNW/7py++y8CXjTa0iRJkqSlM6qT5SRJkqSJ4iWW55Er0sp+6vhqZT+SJEm6I48IS9IyM8Rl7Q9J8q5m++VJ1vVtOzbJp5LsTPLZJHdrtXhJmiAGYUlaRvoua38asAHYnGTDnGZnAbdU1QOA1wGvbB57EPBXwLOr6kHAScD3WipdkiaOQViSlpdhLmt/Oj+8cNHFwMlJAjwOuKqqPgO3z/H+/ZbqlqSJYxCWpOVlmMva396mqvYBtwL3Bo4BKsklSf45yX+Zbwde7l6SegzCkrRyHAQ8GvjN5t9fS3Ly3EZe7l6SegzCkrS8DLysfX+bZlzwPYGb6B09/nhVfa2qvkVv/veHL3nF0n4YdDJo0+apSXY1J32+o+0a1R0GYUlaXgZe1r5ZPqO5/2Tg0qoqelcAfUiSH28C8onArpbqlgYa5mTQJOvpXaTrUc1Jn3/Ydp3qDucRlqRlZMjL2r8FuCDJbnqXtd/UPPaWJK+lF6YL2F5VHxjLDyLN7/aTQQGSzJ4M2v+B7VnAeVV1C0BV3dh6leoMg7AkLTNDXNb+O8BTFnjsX9GbQk1ajuY7GfQRc9ocA5DkE/Q+DL6sqj4094mSbAG2AKxdu3ZJitXK59AISZK0nBwErKc3D/Zm4C+THD63kSd9ahQMwpIkqS3DnAy6B5iuqu9V1ReBz9MLxtLIGYQlSVJbhjkZ9L30jgaTZDW9oRLXtlijOsQxwpKWRP4kreynXlqt7EfSgRvyZNBLgMcl2QV8H3hhVd00vqq1khmEJUlSa4Y4GbSA5zU3aUk5NEKSJEmdZBBWZwy6mlGStUk+muTTSa5K8vhx1ClJktphEFYnDHM1I+AlwEVV9TB6J3C8od0qJUlSmwzC6orbr2ZUVbcBs1cz6lfAYc39ewJfbrE+SZLUMk+WU1cMczWjlwF/l+T3gbsDp7RTmiRJGgePCEs/tBl4a1UdCTweuCDJHf6PJNmSZCbJzN69e1svUpIkjYZBWF0xzNWMzgIuAqiqTwF3A1bPfSIv6ylJ0spgEFZXDHM1oy8BJwMkeSC9IOwhX0mSVqihgvCgaaf62j0pSSWZGl2J0oGrqn3A7NWMrqY3O8TOJOcm2dg0ez7wrCSfAd4JnNlM7C5JklaggSfL9U079Vh6JxjtSDJdVbvmtDsUeC5w+VIUKh2oIa5mtAt4VNt1SZKk8RjmiPAw004BvBx4JfCdEdYnSZIkLYlhgvB8004d0d8gycOBo6rqAyOsTZIkSVoyB3yyXDO91Gvpja8c1NZppyRpgCEuB35Iknc12y9Psq5Zvy7Jt5Nc2dze2HrxkjRBhgnCg6adOhR4MPCxJNcBjwSm5zthzmmnJGlxQ14O/Czglqp6APA6esPSZn2hqo5rbs9upWhJmlDDBOFFp52qqluranVVrauqdcBlwMaqmlmSiiVpZRvmvIzTgbc19y8GTk6SFmuUpBVhYBAectopSdJoDDwvo79N8x59K3DvZtvRST6d5O+T/MJ8O3CYmiT1DJw+DQZPOzVn/UkHXpZyRTsHd+p4p8mVVpCvAGur6qYkxwPvTfKgqvpGf6Oq2gZsA5iamvJNQFJneWU5SVpehrkc+O1tkhwE3BO4qaq+W1U3AVTVFcAXgGOWvGJJmlAGYUlaXoa5HPg0cEZz/8nApVVVSdY0J9uR5KeB9cC1LdUtSRNnqKERkqR2VNW+JLPnZawCzp89LwOYqapp4C3ABUl2AzfTC8sAjwHOTfI94AfAs6vq5vZ/CkmaDAZhSVpmhrgc+HeAp8zzuL8B/mbJC5SkFcKhEZIkSeokg7AkSZI6ySAsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTnIeYUkrVv4kreynXlqt7EeSNFoeEZYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJ2Q5NQk1yTZneTsBdo8NcmuJDuTvKPtGiVJUru8oIZWvCSrgPOAxwJ7gB1JpqtqV1+b9cCLgEdV1S1JfnI81UqSpLZ4RFhdcAKwu6qurarbgAuB0+e0eRZwXlXdAlBVN7Zco3S7Qd9gJDkkybua7ZcnWTdn+9ok30zygtaKlqQJZBBWFxwBXN+3vKdZ1+8Y4Jgkn0hyWZJTW6tO6tP3DcZpwAZgc5INc5qdBdxSVQ8AXge8cs721wIfXOpaJWnSGYSlnoOA9cBJwGbgL5McPl/DJFuSzCSZ2bt3b3sVqiuG+QbjdOBtzf2LgZOTBCDJE4EvAjvbKVeSJtdQQXiIr+meneSzSa5M8o/zHL2QxukG4Ki+5SObdf32ANNV9b2q+iLweXrB+A6qaltVTVXV1Jo1a5akYHXaMN9g3N6mqvYBtwL3TnIP4I+BP1lsB36Y0zgNc/Jy0+5JSSrJVJv1qVsGBuEhv6Z7R1U9pKqOA15F72s5abnYAaxPcnSSg4FNwPScNu+ldzSYJKvpDZW4tsUapVF4GfC6qvrmYo38MKdxGTJTkORQ4LnA5e1WqK4ZZtaI27+mA0gy+zXd7WfcV9U3+trfHahRFikdiKral2QrcAmwCji/qnYmOReYqarpZtvjkuwCvg+8sKpuGl/V6rBhvsGYbbMnyUHAPYGbgEcAT07yKuBw4AdJvlNVr1/yqrVkrrgirezn+ONb+dM9MFM0Xk5v7PsL2yhK3TVMEJ7va7pHzG2U5DnA84CDgV8aSXXSiFTVdmD7nHXn9N0vev33eS2XJs11+zcY9ALvJuA35rSZBs4APgU8Gbi06cO/MNsgycuAbxqCtcwMzBRJHg4cVVUfSGIQ1pIa2clyVXVeVf0MvfFpL5mvjePSJGlxzZjf2W8wrgYumv0GI8nGptlb6I0J3k3vw9uC4yylSZLkLvSGVz5/iLZmCh2wYY4ID/M1Xb8Lgf8134aq2gZsA5iamnL4hCTNY4hvML4DPGXAc7xsSYqTDsygTHEo8GDgY81EKPcFppNsrKqZ/icyU2gUhjkiPPBEo+aqXLN+BfjX0ZUoSZJWiEUzRVXdWlWrq2pdVa0DLgPuEIKlURl4RHjIE422JjkF+B5wC72xa5IkSbcbMlNIrRlmaMQwX9M9d8R1SZKkFWhQppiz/qQ2alJ3eWU5SZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJakZSbJqUmuSbI7ydnzbD8kybua7ZcnWdesPyHJlc3tM0l+rfXiJWmCHDTuAiRJP5RkFXAe8FhgD7AjyXRV7eprdhZwS1U9IMkm4JXA04DPAVNVtS/J/YDPJHlfVe1r+ceQNOtP0s5+Xlrt7GeF8YiwJC0vJwC7q+raqroNuBA4fU6b04G3NfcvBk5Okqr6Vl/ovRvgX0ZJWoRBWJKWlyOA6/uW9zTr5m3TBN9bgXsDJHlEkp3AZ4Fnz3c0OMmWJDNJZvbu3bsEP4IkTQaDsCStIFV1eVU9CPg54EVJ7jZPm21VNVVVU2vWrGm/SElaJgzCkrS83AAc1bd8ZLNu3jZJDgLuCdzU36Cqrga+CTx4ySqVpAlnEJak5WUHsD7J0UkOBjYB03PaTANnNPefDFxaVdU85iCAJPcHfha4rp2yJWnyOGuEJC0jzYwPW4FLgFXA+VW1M8m5wExVTQNvAS5Ishu4mV5YBng0cHaS7wE/AH6vqr7W/k8hSZPBIKzOSHIq8Bf0wsWbq+oVC7R7Er0z8X+uqmZaLFECoKq2A9vnrDun7/53gKfM87gLgAuWvEBJWiEcGqFO6Jub9TRgA7A5yYZ52h0KPBe4vN0KJUlS2wzC6oph5mYFeDm9ixN8p83iJElS+4YKwkNc7vN5SXYluSrJR5qTNKTlZODcrEkeDhxVVR9Y7Imcg1WSpJVhYBAe8ivlT9O7rOex9MZWvmrUhUpLKcldgNcCzx/U1jlYJUlaGYY5IjzwK+Wq+mhVfatZvIzevJfScjJobtZD6c23+rEk1wGPBKaTTLVWoSRJatUwQXiYy332Owv44IEUJS2BRedmrapbq2p1Va2rqnX0PtBtdNYISZJWrpFOn5bk6cAUcOIC27cAWwDWrl07yl1LixpyblZJktQhwwThYS73SZJTgBcDJ1bVd+d7oqraBmwDmJqaqv2uVjoAg+ZmnbP+pDZqkiRJ4zPM0IiBl/tM8jDgTfS+Sr5x9GVKkiRJozUwCFfVPmD2K+WrgYtmv1JOsrFp9mrgHsC7k1yZxK+ZJUmStKwNNUZ4iMt9njLiuiRJkqQl5ZXlJEmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmSJHWSQViSlpkkpya5JsnuJGfPs/2QJO9qtl+eZF2z/rFJrkjy2ebfX2q9eEmaIAZhSVpGkqwCzgNOAzYAm5NsmNPsLOCWqnoA8Drglc36rwFPqKqHAGcAF7RTtTS8IT7oPS/JriRXJflIkvuPo051g0FYkpaXE4DdVXVtVd0GXAicPqfN6cDbmvsXAycnSVV9uqq+3KzfCfxYkkNaqVoawpAf9D4NTFXVsfT696varVJdYhCWpOXlCOD6vuU9zbp521TVPuBW4N5z2jwJ+Oeq+u7cHSTZkmQmyczevXtHVrg0hIEf9Krqo1X1rWbxMuDIlmtUhxiEJWmFSfIgesMlfme+7VW1raqmqmpqzZo17Ranrhvmg16/s4APLmlF6rSDxl2AJOlH3AAc1bd8ZLNuvjZ7khwE3BO4CSDJkcB7gGdU1ReWvlxpaSR5OjAFnLjA9i3AFoC1a9e2WJlWEo8IS9LysgNYn+ToJAcDm4DpOW2m6Z0MB/Bk4NKqqiSHAx8Azq6qT7RVsLQfhvmgR5JTgBcDG+cb3gN+s6HRMAhL0jLSjPndClwCXA1cVFU7k5ybZGPT7C3AvZPsBp4HzJ55vxV4AHBOkiub20+2/CNIixn4QS/Jw4A30QvBN46hRnWIQyMkaZmpqu3A9jnrzum7/x3gKfM87k+BP13yAqU7qar2JZn9oLcKOH/2gx4wU1XTwKuBewDvTgLwparauOCTSgfAICxJklozxAe9U1ovSp3l0AhJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhNUJSU5Nck2S3UnOnmf785LsSnJVko8kuf846pQkSe0ZKggPESIek+Sfk+xL8uTRlyndeUlWAecBpwEbgM1JNsxp9mlgqqqOBS4GXtVulZIkqW0Dg/CQIeJLwJnAO0ZdoDQCJwC7q+raqroNuBA4vb9BVX20qr7VLF4GHNlyjZIkqWXDHBEeJkRcV1VXAT9YghqlA3UEcH3f8p5m3ULOAj640MYkW5LMJJnZu3fviEqUJEltGyYI72+IkCZWkqcDU8CrF2pTVduqaqqqptasWdNecZIkaaRaPVnOI2kakxuAo/qWj2zW/YgkpwAvBjZW1Xdbqk2SJI3JMEF4qBAxDI+kaUx2AOuTHJ3kYGATMN3fIMnDgDfRC8E3jqFGSZLUsmGC8MAQIS1nVbUP2ApcAlwNXFRVO5Ocm2Rj0+zVwD2Adye5Mol9XJKkFe6gQQ2qal+S2RCxCjh/NkQAM1U1neTngPcA9wKekORPqupBS1q5tB+qajuwfc66c/run9J6UZIkaayGGiNcVdur6piq+pmq+u/NunOqarq5v6Oqjqyqu1fVvQ3BknTnDTF3+yFJ3tVsvzzJumb9vZN8NMk3k7y+9cIlacJ4ZTlJWkaGnLv9LOCWqnoA8Drglc367wD/DXhBS+VK0kQzCEvS8jJw7vZm+W3N/YuBk5Okqv6jqv6RXiCWJA1gEJak5WWYudtvb9OcDHorcO9hd+BUlpLUYxCWpI5xKktJ6jEIS9LyMszc7be3SXIQcE/gplaqk6QVxCAsScvLMHO3TwNnNPefDFxaVdVijZK0IgycR1iS1J5h5m4H3gJckGQ3cDO9sAxAkuuAw4CDkzwReFxV7Wr5x5CkiWAQlqRlZogLwHwHeMoCj123pMVJ0gri0AhJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR10lBBOMmpSa5JsjvJ2fNsPyTJu5rtlydZN/JKpQNkP9akOJC+muRFzfprkvxyq4VLQ/C9WMvJwCCcZBVwHnAasAHYnGTDnGZnAbdU1QOA1wGvHHWh0oGwH2tSHEhfbdptAh4EnAq8oXk+aVnwvVjLzTBHhE8AdlfVtVV1G3AhcPqcNqcDb2vuXwycnCSjK1M6YPZjTYoD6aunAxdW1Xer6ovA7ub5pOXC92ItKwcN0eYI4Pq+5T3AIxZqU1X7ktwK3Bv4Wn+jJFuALc3iN5Ncc2eKnsfqufsag/2uIYz8//W4a7hTv4dFarj/AVXzo+zHS1RDXjbePjTi/Y+6hjvThw+krx4BXDbnsUfcod6l68Mw/n58J/Y//j60zGvwvbh9+1/DMngvHPv+F38N5u3HwwThkamqbcC2UT9vkpmqmhr181rDZO2/Lfbjlbv/5VLDUluqPgzjf/3GvX9raI/vxSu7hrb2P8zQiBuAo/qWj2zWzdsmyUHAPYGbRlGgNCL2Y02KA+mrwzxWGiffi7WsDBOEdwDrkxyd5GB6J2JMz2kzDZzR3H8ycGlV1ejKlA6Y/ViT4kD66jSwqTnr/mhgPfBPLdUtDcP3Yi0rA4dGNONztgKXAKuA86tqZ5JzgZmqmgbeAlyQZDdwM72O3aYl+YpvP1nD+Pe/IPvx0MZdw7j3D2Ou4UD6atPuImAXsA94TlV9v+UfYdy/w3HvH6xhQb4XD80aWtp//JAlSZKkLvLKcpIkSeokg7AkSZI6aeKCcJL7JrkwyReSXJFke5JjknwoydeTvH9MNZyQ5FNJdia5KsnTxlDDiUn+OcmVTR3Pbnn/xzTbDkuyJ8nrl2r/k85+PP4+vEgN9uMhjbsfj7sPL1KD/XhCjLsPL1JDp/rxWPtwVU3Mjd5s358Cnt237qHALwAnA08A3j+mGk4E1jfLPwV8BTh8DDUc0izfA7gO+Kk2fw/N/b8A3gG8ftx9Zjne7Mfj78ODfg/NffvxnXz92ujH4+7DQ9RgP17mt3H34SH6UCf68bj7cKsX1BiBXwS+V1VvnF1RVZ+ZvZ/kpHHX0Cx/OcmNwBrg6+OoATiEpTviv+D+kxwP3Af4ELCiJ3M/APbj8ffhRWuwHw9l3P143H14qBqwHy9n4+7DA2tolld6Px5rH560oREPBq5Y7jUkOQE4GPhC2zUkOSrJVfQuT/nKqvpyW/tPchfgz4EXLME+VxL78fj78II12I+HNu5+PO4+vGgN9uOJMO4+PFQNHejHY+3DkxaEl70k9wMuAH67qn7Q9v6r6vqqOhZ4AHBGkvu0uPvfA7ZX1Z4W96klMM5+POY+DPbjFaHj78VgP14ROt6PW+nDkxaEdwLHL9cakhwGfAB4cVVdNo4aZjWf2j5Hb6xTW/v/eWBrkuuA1wDPSPKKJdj/pLMfj78PL1aD/Xg44+7H4+7Di9Ywy368rI27Dy9aQ4f68Vj78KQF4UuBQ5JsmV2R5NgkS/UGsz81nAi8B3h7VV08php+IcmPNcv3Ah4NXNPW/oE3VtXaqlpH76uMt1fV2Uuw/0lnPx5/H16wBuzHwxp3Px53H16sBvvxZBh3H16shi7147H24YkKwtU7ffDXgFOaKTZ2An8GfDXJPwDvBk5uptn45ZZreExzO7OZauTKJMe1XMPPApcn+Qzw98BrquqzLe7/q6Pe10pkPx5/Hx5Qg/14COPux+PuwwNqsB9PgHH34QE1dKYfj7sPe4llSZIkddJEHRGWJEmSRsUgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKS7iDJE5K8a9x1SHdGkmOTfHLcdUha/gzCkgBIcl2SUwCq6n3Ag5IcO6D9bUlWz1n/6SSVZF2z/NYkf9q3/UFJvpLkP5J8s7l9P8l3+pb/69L8lFrJknwsyTOr6irg60meMKBtJXnonPXvadaf1CwfnuT8JF9N8u9JPp/k7L72leQBfcsvaPr3g0b+A0oaOYOwpIW8E9gyoM0Xgc2zC0keAvz4Qo2TPAz4KPCnVXX3qrpHVd0D+Adg6+xyVf2PAy9fHffXwO8MaPN54BmzC0nuDfw8sLevzeuAewAPBO4JbAR2z/dkSV4C/CFwYlXtvLOFS2qPQVjSQj4G/MqANhfQFySAM4C3z9cwyQnAh4H/WlXnjaJAaREfA05Ocsgibf4aeFqSVc3yZuA9wG19bX4OeEdV3VJVP6iqf6mqi+c+UfOtxzOBx1TV50fyE0hacgZhSQu5GliX5LBF2lwGHJbkgU2Y2AT81TztTgA+BPxRVb159KVKP6qqbgC+B/ynRZp9GdgFPK5ZfgZ3/CB3GfDfk/x2kvULPM8rgKfRC8HX3vmqJbXNICxpIf/e/Hv4gHazR4UfSy883zBPm0cCtwIfHFVx0hD+ncH99+3AM5L8LHB4VX1qzvbfp3fkeCuwK8nuJKfNafM44ENV9aUR1CypRQZhSQs5tPn36wPaXQD8BnAmCwyLAM4DZoAPJ7nXKIqThnAog/vv3wK/RC/oXjB3Y1V9u6r+R1UdD9wbuAh4d5Kf6Gu2CXhykj8ZSdWSWmMQlrSQBwLXVdU3FmtUVf9G76S5x9MLFfP5Pr2w/CXgkgHDLaQDluQI4GDgmsXaVdW36H1T8bvME4TntP0G8D+AuwNH9236PHAK8Hv9M0pIWv4MwpIWciLDD2U4C/ilqvqPhRpU1feApwBfA7YnufuBlygt6ETg0qr67hBt/yu9mR6um7shyX9L8nNJDk5yN+C59I4y/0jAbmaJOAV4YZI/PMDaJbXkoHEXIGnZ2gw8fZiGVfWFIdvdluTXgfcD70vyK1X17QOoUZqrmn9/E3jjUA+o+jK9E+cWer7/DawF9gFXAb9SVd+c53k+k+SX6Q0B+k5VDbV/SeOTqhrcSlKnNBci+K2qeuq4a5GGleSfgXOBa4E3VdXPj7kkScucQViSNPGaK7nNAD/bjFuXpIEcIyxJmmhJXgn8HfDHhmBJ+8MjwpIkSeokjwhLkiSpkwzCkiRJ6iSnT1MnJDkV+AtgFfDmqnrFPG2eCryM3nRJn6mq3xj0vKtXr65169aNtlitGFdcccXXqmrNuOtYjH1Yg0xCP5buLIOwVrwkq+hd4vexwB5gR5LpqtrV12Y98CLgUVV1S5KfHOa5161bx8zMzFKUrRUgybI/ccs+rEEmoR9Ld5ZDI9QFJwC7q+raqroNuBA4fU6bZwHnVdUtAFV1Y8s1SpKklhmE1QVHANf3Le9p1vU7BjgmySeSXNYMpZAkSSuYQyOknoOA9cBJwJHAx5M8pKq+Prdhki3AFoC1a9e2WKIkSRoljwirC24AjupbPrJZ128PMF1V36uqLwKfpxeM76CqtlXVVFVNrVnj+SOSJE0qg7C6YAewPsnRSQ4GNgHTc9q8l97RYJKspjdU4toWa5QkSS0zCGvFq6p9wFbgEuBq4KKq2pnk3CQbm2aXADcl2QV8FHhhVd00noolSVIbHCOsTqiq7cD2OevO6btfwPOam7QsDZoPO8la4G3A4U2bs5u+L0mah0eEJWkC9M2HfRqwAdicZMOcZi+h943Hw+gNAXpDu1VK0mQxCEvSZBhmPuwCDmvu3xP4cov1SdLEcWiEpBXr/e9/fyv7+dVf/dU2djPffNiPmNPmZcDfJfl94O7AKfM9kVMATpap/E4r+5mpN7WyH2k58YiwJK0cm4G3VtWRwOOBC5Lc4X3eKQAlqccgLEmTYZj5sM8CLgKoqk8BdwNWt1KdJE0gg7AkTYZh5sP+EnAyQJIH0gvCe1utUpImiGOEJWkCVNW+JLPzYa8Czp+dDxuYqapp4PnAXyb5I3onzp3ZTA04kRwbK2mpGYQlaUIMMR/2LuBRbdclSZPKoRGSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEkTIsmpSa5JsjvJ2fNsf12SK5vb55N8fQxlStLEOGjcBUiSBkuyCjgPeCywB9iRZLqqds22qao/6mv/+8DDWi9UkiaIR4TVCUMcSTszyd6+o2nPHEed0iJOAHZX1bVVdRtwIXD6Iu03A+9spTJJmlAeEdaKN8yRtMa7qmpr6wVKwzkCuL5veQ/wiPkaJrk/cDRw6QLbtwBbANauXTvaKiVpgnhEWF2wv0fSpEm3Cbi4qr4/38aq2lZVU1U1tWbNmpZLk6TlwyCsLpjvSNoR87R7UpKrklyc5Kh2SpOGdgPQ3y+PbNbNZxMOi5CkgQzCUs/7gHVVdSzwYeBtCzVMsiXJTJKZvXv3tlagOm8HsD7J0UkOphd2p+c2SvKzwL2AT7VcnyRNHIOwumDgkbSquqmqvtssvhk4fqEn82tljUNV7QO2ApcAVwMXVdXOJOcm2djXdBNwYVXVOOqUpEniyXLqgtuPpNELwJuA3+hvkOR+VfWVZnEjvaAhLStVtR3YPmfdOXOWX9ZmTZI0yQzCWvGqal+S2SNpq4DzZ4+kATNVNQ38QXNUbR9wM3Dm2AqWJEmtMAirEwYdSauqFwEvarsuSZI0Po4RliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiZEklOTXJNkd5KzF2jz1CS7kuxM8o62a5SkSXLQuAuQJA2WZBVwHvBYYA+wI8l0Ve3qa7MeeBHwqKq6JclPjqdaSZoMHhGWpMlwArC7qq6tqtuAC4HT57R5FnBeVd0CUFU3tlyjJE0Ug7AkTYYjgOv7lvc06/odAxyT5BNJLktyamvVSdIEcmiEJK0cBwHrgZOAI4GPJ3lIVX29v1GSLcAWgLVr17ZcoiQtHx4RlqTJcANwVN/ykc26fnuA6ar6XlV9Efg8vWD8I6pqW1VNVdXUmjVrlqxgSVruDMKSNBl2AOuTHJ3kYGATMD2nzXvpHQ0myWp6QyWubbFGSZooBmFJmgBVtQ/YClwCXA1cVFU7k5ybZGPT7BLgpiS7gI8CL6yqm8ZTsSQtf44RlqQJUVXbge1z1p3Td7+A5zU3SdIABmFJkhYwld9pZT8z9aZW9iPpRzk0QpIkSZ1kEJYkSVInGYQlSZLUSQZhdUaSU5Nck2R3krMXafekJJVkqs36JElSuzxZTp2QZBVwHvBYehcd2JFkuqp2zWl3KPBc4PL2q5SWF08Uk7TSeURYXXECsLuqrq2q24ALgdPnafdy4JXAd9osTpIktc8grK44Ari+b3lPs+52SR4OHFVVH2izMEmSNB4GYQlIchfgtcDzh2i7JclMkpm9e/cufXGSJGlJGITVFTcAR/UtH9msm3Uo8GDgY0muAx4JTM93wlxVbauqqaqaWrNmzRKWLEmSlpJBWF2xA1if5OgkBwObgOnZjVV1a1Wtrqp1VbUOuAzYWFUz4ylXkiQtNYOwOqGq9gFbgUuAq4GLqmpnknOTbBxvdZIkaRycPk2dUVXbge1z1p2zQNuT2qhJkiSNj0eEJUmS1EkGYUmSJHWSQViSJEmdZBCWJElSJxmEJUmS1EkGYUmaEElOTXJNkt1Jzp5n+5lJ9ia5srk9cxx1StKkcPo0SZoASVYB5wGPBfYAO5JMV9WuOU3fVVVbWy9QkiaQR4QlaTKcAOyuqmur6jbgQuD0MdckSRPNICxJk+EI4Pq+5T3NurmelOSqJBcnOaqd0iRpMhmEJWnleB+wrqqOBT4MvG2+Rkm2JJlJMrN3795WC5Sk5cQgLEmT4Qag/wjvkc2621XVTVX13WbxzcDx8z1RVW2rqqmqmlqzZs2SFCtJk8AgLEmTYQewPsnRSQ4GNgHT/Q2S3K9vcSNwdYv1SdLEcdYISZoAVbUvyVbgEmAVcH5V7UxyLjBTVdPAHyTZCOwDbgbOHFvBkjQBDMKSNCGqajuwfc66c/ruvwh4Udt1SdKkcmiEJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMKSJEnqJIOwJEmSOskgLEmSpE4yCEuSJKmTDMLqjCSnJrkmye4kZ8+z/dlJPpvkyiT/mGTDOOqUJEntMAirE5KsAs4DTgM2AJvnCbrvqKqHVNVxwKuA17ZbpSRJapNBWF1xArC7qq6tqtuAC4HT+xtU1Tf6Fu8OVIv1SZKklhmE1RVHANf3Le9p1v2IJM9J8gV6R4T/YL4nSrIlyUySmb179y5JsdJ8Bg3v6Wv3pCSVZKrN+iRp0hiEpT5VdV5V/Qzwx8BLFmizraqmqmpqzZo17RaozhpyeA9JDgWeC1zeboWSNHkMwuqKG4Cj+paPbNYt5ELgiUtZkLSfBg7vabwceCXwnTaLk6RJZBBWV+wA1ic5OsnBwCZgur9BkvV9i78C/GuL9UmDDBzek+ThwFFV9YE2C5OkSXXQuAuQ2lBV+5JsBS4BVgHnV9XOJOcCM1U1DWxNcgrwPeAW4IzxVSztnyR3oTfTyZlDtN0CbAFYu3bt0hYmScuYQVidUVXbge1z1p3Td/+5rRclDW/Q8J5DgQcDH0sCcF9gOsnGqprpf6Kq2gZsA5iamnJ2FEmd5dAISZoMiw7vqapbq2p1Va2rqnXAZcAdQrAk6YcMwpI0AapqHzA7vOdq4KLZ4T1JNo63OkmaTA6NkKQJMWh4z5z1J7VRkyRNMo8IS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTvKCGtIK9f73v7+V/fzqr/5qK/uRJGnUPCIsSZKkTvKIsKQl4RFpSdJy5xFhSZIkdZJBWJIkSZ1kEJYkSVInGYQlSZLUSQZhSZIkdZJBWJIkSZ1kEJYkSVInGYQlaUIkOTXJNUl2Jzl7nu3PTvLZJFcm+cckG8ZRpyRNCoOwJE2AJKuA84DTgA3A5nmC7juq6iFVdRzwKuC17VYpSZPFICxJk+EEYHdVXVtVtwEXAqf3N6iqb/Qt3h2oFuuTpInjJZYlaTIcAVzft7wHeMTcRkmeAzwPOBj4pXZKk6TJ5BFhSVpBquq8qvoZ4I+Bl8zXJsmWJDNJZvbu3dtugZK0jBiEJWky3AAc1bd8ZLNuIRcCT5xvQ1Vtq6qpqppas2bN6CqUpAljEJakybADWJ/k6CQHA5uA6f4GSdb3Lf4K8K8t1idJE8cgrE4YYtqp5yXZleSqJB9Jcv9x1CktpKr2AVuBS4CrgYuqameSc5NsbJptTbIzyZX0xgmfMZ5qJWkyeLKcVry+aaceS+8Eox1JpqtqV1+zTwNTVfWtJL9Lb+qpp7VfrbSwqtoObJ+z7py++89tvShJmmAeEVYXDDPt1Eer6lvN4mX0xl9KkqQVzCCsLphv2qkjFml/FvDBJa1IkiSNnUMjpD5Jng5MAScu0mYLsAVg7dq1LVUmSZJGzSPC6oKhpp1KcgrwYmBjVX13oSdz6ilJklYGg7C6YJhppx4GvIleCL5xDDVKkqSWGYS14g057dSrgXsA705yZZLpBZ5OkiStEI4RVicMMe3UKa0XJUmSxsojwpIkSeokg7AkSZI6ySAsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTjIIS5IkqZMMwpIkSeokg7AkSZI6ySAsSZKkTjIIS9KESHJqkmuS7E5y9jzbn5dkV5Krknwkyf3HUackTQqDsCRNgCSrgPOA04ANwOYkG+Y0+zQwVVXHAhcDr2q3SkmaLAZhSZoMJwC7q+raqroNuBA4vb9BVX20qr7VLF4GHNlyjZI0UQzCkjQZjgCu71ve06xbyFnAB5e0IkmacAeNuwBJ0mgleTowBZy4wPYtwBaAtWvXtliZJC0vHhGWpMlwA3BU3/KRzbofkeQU4MXAxqr67nxPVFXbqmqqqqbWrFmzJMVK0iQwCEvSZNgBrE9ydJKDgU3AdH+DJA8D3kQvBN84hholaaIYhCVpAlTVPmArcAlwNXBRVe1Mcm6SjU2zVwP3AN6d5Mok0ws8nSQJxwhL0sSoqu3A9jnrzum7f0rrRUnSBPOIsCRJkjrJICxJkqROMghLkiSpkwzCkiRJ6iSDsCRJkjrJICxJkqROMghLkiSpkwzCkiRJ6iSDsDohyalJrkmyO8nZ82x/TJJ/TrIvyZPHUaMkSWqXQVgrXpJVwHnAacAGYHOSDXOafQk4E3hHu9VJkqRx8RLL6oITgN1VdS1AkguB04Fdsw2q6rpm2w/GUaAkSWqfR4TVBUcA1/ct72nWSZKkDjMIS/spyZYkM0lm9u7dO+5yJEnSnWQQVhfcABzVt3xks+5OqaptVTVVVVNr1qw54OIkSdJ4GITVBTuA9UmOTnIwsAmYHnNNkiRpzAzCWvGqah+wFbgEuBq4qKp2Jjk3yUaAJD+XZA/wFOBNSXaOr2JJktQGZ41QJ1TVdmD7nHXn9N3fQW/IhCRJ6giPCEuSJKmTDMKSNCG8QqIkjZZBWJImgFdIlKTRc4ywJE0Gr5AoSSPmEWFJmgxeIVGSRswgLEkd49URJanHICxJk2FkV0j06oiS1GMQlqTJ4BUSJWnEDMKSNAG8QqIkjZ6zRkjShPAKiZI0Wh4RliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmFJkiR1kkFYkiRJnWQQliRJUicZhCVJktRJBmF1RpJTk1yTZHeSs+fZfkiSdzXbL0+ybgxlSguyD0vSaBmE1QlJVgHnAacBG4DNSTbMaXYWcEtVPQB4HfDKdquUFmYflqTRMwirK04AdlfVtVV1G3AhcPqcNqcDb2vuXwycnCQt1igtxj4sSSNmEFZXHAFc37e8p1k3b5uq2gfcCty7leqkwezDkjRiB427AGnSJNkCbGkWv5nkmhE99WrgayN6rkmtYdz7H3UN9x/R84zUEvZhuBOvX7JthLsf+/5XYg3Lsh9Lo2AQVlfcABzVt3xks26+NnuSHATcE7hp7hNV1TZg5H+1ksxU1dSon3eSahj3/pdLDQtY9n0Yxv/6jXv/1iBNFodGqCt2AOuTHJ3kYGATMD2nzTRwRnP/ycClVVUt1igtxj4sSSPmEWF1QlXtS7IVuARYBZxfVTuTnAvMVNU08BbggiS7gZvpBQ1pWbAPS9LoGYTVGVW1Hdg+Z905ffe/Azyl7br6LMlX1ftp3DWMe/+wPGqY1wT0YRj/6zfu/YM1SBMjfmsmSZKkLnKMsCRJkjrJICyNQZL7JrkwyReSXJFke5JjknwoydeTvH9MNZyQ5FNJdia5KsnTWt7/iUn+OcmVTQ3PXqr9L1LDMc22w5LsSfL6paxhko27H4+7Dy9Sg/1YmhAOjZBa1lzp65PA26rqjc26hwKHAQcDPw78TlX96hhqOBz4clX9a5KfAq4AHlhVX29x/5dV1XeT3AP4HPCfq+rLo9z/gBoOq6p/SPIXwBrg5qraOur9T7px9+Nx9+EharAfSxPAk+Wk9v0i8L3ZP1oAVfWZ2ftJThp3Dc3yl5PcSO+P6Nfb3j9wCEv7rdWCNSQ5HrgP8CHAuVjnN+5+PO4+PFQN2I+lZc2hEVL7HkzvKNWyriHJCfSO7H2hzf0nOSrJVfQuFfzKpTiKtlgNSe4C/DnwgiXa70ox7n487j68aA32Y2kyGIQl3UGS+wEXAL9dVT9oc99VdX1VHQs8ADgjyX3a3D/we8D2qtrT8n41QuPsw2A/liaFQVhq307g+OVaQ5LDgA8AL66qy9re/6zmCNrngF9ouYafB7YmuQ54DfCMJK9Yohom2bj78bj78KI1zLIfS8ubQVhq36XAIUm2zK5IcmySpfpDuT81nAi8B3h7VV08hv3/QpIfa5bvBTwauKbNGoA3VtXaqlpH72vlt1fV2UtUwyQbdz8edx9erAb7sTQhDMJSy6o3VcuvAac00x3tBP4M+GqSfwDeDZzcTHn0yy3X8JjmdmYz9dOVSY5rcf8/C1ye5DPA3wOvqarPjnr/A2r46lLsb6UZdz8edx8eUIP9WJoQTp8mSZKkTvKIsCRJkjrJICxJkqROMghLkiSpkwzCkiRJ6iSDsCRJkjrJICxJkqROMghLkiSpkwzCkiRJ6qT/H9cFSTevF/LrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x1080 with 12 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cases = [1, 2, 3, 4]\n",
    "\n",
    "fig, axs = plt.subplots(nrows=3, ncols=4, figsize=(12, 15))\n",
    "fig.subplots_adjust(hspace=0.3, wspace=0.3)\n",
    "\n",
    "axs[0, 0].bar(cases, md*1000, linewidth=1.0, color=[0.0, 0.0, 0.8])\n",
    "axs[0, 0].set_title('A) MD ($x10^{-3}mm^2/s$)')\n",
    "axs[0, 1].bar(cases, rd*1000, linewidth=1.0, color=[0.0, 0.0, 0.4])\n",
    "axs[0, 1].set_title('B) RD ($x10^{-3}mm^2/s$)')\n",
    "axs[0, 2].bar(cases, ad*1000, linewidth=1.0, color=[0.0, 0.8, 0.8])\n",
    "axs[0, 2].set_title('C) AD ($x10^{-3}mm^2/s$)')\n",
    "axs[0, 3].bar(cases, fa, linewidth=1.0, color=[1.0, 0.0, 0.0])\n",
    "axs[0, 3].set_title('D) FA')\n",
    "\n",
    "axs[1, 0].bar(cases, mk, linewidth=1.0, color=[0.0, 0.8, 0.0])\n",
    "axs[1, 0].set_title('E) MK')\n",
    "axs[1, 1].bar(cases, rk, linewidth=1.0, color=[0.0, 0.5, 0.0])\n",
    "axs[1, 1].set_title('F) RK')\n",
    "axs[1, 2].bar(cases, ak, linewidth=1.0, color=[0.8, 0.8, 0.0])\n",
    "axs[1, 2].set_title('G) AK')\n",
    "axs[1, 3].bar(cases, kfa, linewidth=1.0, color=[1, 0.5, 0.0])\n",
    "axs[1, 3].set_title('H) KFA')\n",
    "\n",
    "\n",
    "axs[2, 1].bar(cases, mkt, linewidth=1.0, color=[0.7, 0.7, 0.7])\n",
    "axs[2, 1].set_title('I) MKT')\n",
    "axs[2, 2].bar(cases, msk, linewidth=1.0, color=[0.2, 0.0, 0.4])\n",
    "axs[2, 2].set_title('J) MSK')\n",
    "\n",
    "for ii in range(3):\n",
    "    for jj in range(4):\n",
    "        axs[ii, jj].set_xticks(cases)\n",
    "        axs[ii, jj].set_xticklabels(['C1', 'C2', 'C3', 'C4'])\n",
    "\n",
    "axs[2, 0].axis('off')\n",
    "axs[2, 3].axis('off')\n",
    "\n",
    "fig.savefig('Figure2.png', bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}