analysis_ellipsoids_kappa_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal components analysis for order parameter construction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from analysis_tools.read_gsd import ReadGSD\n",
    "from analysis_tools.feature_creation import FrameToFeatures, TrajectoryToFeatures, TrajectoryToFeaturesBatch, TrajectoryToFeaturesBatchR\n",
    "from analysis_tools.feature_creation import FrameToFeaturesComposition, TrajectoryToFeaturesComposition\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array, arange, floor, floor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hard disks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Kappa=4\n",
    "particle_inc = 104 #skip this many particles in generating features; #wp: because it selects index 0 be default, rounds up to integer if not a mulptiple of Ntot\n",
    "N_nn = 3120 #number of nearest nighbors for pca analysis #wp:Should not be bigger than Ntotal\n",
    "nn_inc = 10 #reduces the number of nearest neighbors to include as features\n",
    "split = 1 #chunks the data up so it can be processed by the pca tool if really large\n",
    "N_batch=1\n",
    "frame_skip=4\n",
    "\n",
    "#method = ['distance']\n",
    "method = ['orientation']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify what data to read in and process\n",
    "#wp:there are 51 eta points given range below\n",
    "folder='ellipsoids/kappa_4.0_3N'; \n",
    "file_data = [(arange(0.6100, 0.6800, 0.01), folder),\n",
    "             (arange(0.6850, 0.7200, 0.005), folder),\n",
    "             (arange(0.7300, 0.90001, 0.01), folder)] \n",
    "            \n",
    "             "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#wp: the `_' variable is a dummy return variable like `~' in matlab\n",
    "#wp: Here, the for loop returns two arguments but only keeps the first that contains all the numbers\n",
    "etas = []\n",
    "[etas.extend(etas_) for etas_, _ in file_data]\n",
    "etas = array(etas)\n",
    "#print(etas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perfom PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array_split\n",
    "from sklearn.decomposition import IncrementalPCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train on ideal gas structure first to identify trivial variance contributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#n_components_ig = N_nn/nn_inc-1  #for this the more the better, however, 50 is probably enough\n",
    "n_components_ig = N_batch*N_nn/nn_inc  #for this the more the better, however, 50 is probably enough\n",
    "#n_components_ig = N_nn/nn_inc*(N_batch+2) #for this the more the better, however, 50 is probably enough\n",
    "#print n_components_ig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "from analysis_tools.feature_creation import FrameToFeatures, TrajectoryToFeatures\n",
    "\n",
    "###wp: Ideal gas pca structure\n",
    "incpca_ig = IncrementalPCA(n_components=n_components_ig, whiten=True)  #NOTICE: whitening is on here\n",
    "\n",
    "frames = None\n",
    "features = None\n",
    "split_features = None\n",
    "#split=1\n",
    "\n",
    "#features_partial_ig=[];\n",
    "\n",
    "for etas_, file_base in file_data:\n",
    "    for eta in etas_:\n",
    "        print 'eta={}'.format(eta)\n",
    "        filename = \"{}/trajectory_{:.4f}.gsd\".format(file_base, eta)\n",
    "        frames = ReadGSD(filename, shuffle_data=True, randomize=True)        \n",
    "        #print frames[0]['orientation'][0:5]\n",
    "        \n",
    "        frames_split = array_split(frames, split)\n",
    "        \n",
    "        #frames = None\n",
    "        count = 0\n",
    "        \n",
    "        for frames_partial in frames_split:\n",
    "        #for frames_partial in frames_split[0:-1:2]:\n",
    "            #features_partial = TrajectoryToFeatures(frames_partial, N_nn=N_nn, method=method, \n",
    "            #                                        particle_inc=particle_inc, nn_inc=nn_inc)\n",
    "            #wp: Rather than using all frames by default, allow setting to skip some\n",
    "            \n",
    "            #if eta > 0.8200:\n",
    "            #    N_nni=N_nn2; particle_inci=particle_inc2; nn_inci=nn_inc2;\n",
    "            #else:\n",
    "            N_nni=N_nn; particle_inci=particle_inc; nn_inci=nn_inc;\n",
    "                \n",
    "            frames_partial=frames_partial[::frame_skip]\n",
    "            #frames_partial=frames_partial[0:200]\n",
    "        \n",
    "            features_partial = TrajectoryToFeaturesBatch(frames_partial, N_nn=N_nni, method=method, \n",
    "                                                         particle_inc=particle_inci, nn_inc=nn_inci,N_batch=N_batch)\n",
    "            \n",
    "            #features_partial = TrajectoryToFeaturesBatchR(frames_partial, N_nn=N_nni, method=method, \n",
    "             #                                            particle_inc=particle_inci, nn_inc=nn_inci,N_batch=N_batch)\n",
    "\n",
    "\n",
    "            \n",
    "            print features_partial.shape\n",
    "            incpca_ig.partial_fit(features_partial)\n",
    "            count = count + 1\n",
    "            print count\n",
    "            \n",
    "        clear_output()\n",
    "        frames_split = None\n",
    "        #features_partial_ig.append(features_partial)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train on actual data with trivial data variance removed by the ideal gas model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#n_components = N_nn/nn_inc-1 #this number will not affect the quality of the result just the computation time and amount of info returned\n",
    "n_components = N_batch*N_nn/nn_inc  #for this the more the better, however, 50 is probably enough\n",
    "#n_components = 2*N_nn/nn_inc  #for this the more the better, however, 50 is probably enough\n",
    "#n_components = N_nn/nn_inc*(N_batch+2) #for this the more the better, however, 50 is probably enough"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "incpca = IncrementalPCA(n_components=n_components , whiten=False)  #NOTE: not whitening here\n",
    "\n",
    "frames = None\n",
    "features = None\n",
    "split_features = None\n",
    "#split=1\n",
    "features_partial_ellipse=[]\n",
    "\n",
    "for etas_, file_base in file_data:\n",
    "    for eta in etas_:\n",
    "        print 'eta={}'.format(eta)\n",
    "        filename = \"{}/trajectory_{:.4f}.gsd\".format(file_base, eta)\n",
    "        frames = ReadGSD(filename, shuffle_data=True, randomize=False)\n",
    "        \n",
    "        frames_split = array_split(frames, split)\n",
    "        frames = None\n",
    "        count = 0\n",
    "        for frames_partial in frames_split:\n",
    "            #features_partial = TrajectoryToFeatures(frames_partial, \n",
    "            #                                        N_nn=N_nn, \n",
    "            #                                        method=method, \n",
    "            #                                        particle_inc=particle_inc, \n",
    "            #                                        nn_inc=nn_inc)\n",
    "            \n",
    "\n",
    "            #if eta > 0.8200:\n",
    "            #    N_nni=N_nn2; particle_inci=particle_inc2; nn_inci=nn_inc2;\n",
    "            #else:\n",
    "            N_nni=N_nn; particle_inci=particle_inc; nn_inci=nn_inc;\n",
    "                \n",
    "            frames_partial=frames_partial[::frame_skip]\n",
    "            #frames_partial=frames_partial[0:200]\n",
    "               \n",
    "            features_partial = TrajectoryToFeaturesBatch(frames_partial, N_nn=N_nni, method=method, \n",
    "                                                         particle_inc=particle_inci, nn_inc=nn_inci,N_batch=N_batch)\n",
    "            \n",
    "            #wp:project *hard disk* features onto ig eigenvectors\n",
    "            incpca.partial_fit(incpca_ig.transform(features_partial))\n",
    "            #count = count + 1\n",
    "            #print count\n",
    "        clear_output()\n",
    "        frames_split = None\n",
    "        #features_partial_ellipse.append(features_partial)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### System train without ideal gas correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "########wp: System fit *without ideal gas*\n",
    "incpca = IncrementalPCA(n_components=n_components , whiten=False)\n",
    "\n",
    "frames = None\n",
    "features = None\n",
    "split_features = None\n",
    "\n",
    "split=1\n",
    "for etas_, file_base in file_data:\n",
    "    for eta in etas_:\n",
    "        print 'eta={}'.format(eta)\n",
    "        filename = \"{}/trajectory_{:.4f}.gsd\".format(file_base, eta)\n",
    "        frames = ReadGSD(filename, shuffle_data=True, randomize=False)\n",
    "        \n",
    "        frames_split = array_split(frames, split)\n",
    "        frames = None\n",
    "        count = 0\n",
    "        for frames_partial in frames_split:\n",
    "            #features_partial = TrajectoryToFeatures(frames_partial, \n",
    "            #                                        N_nn=N_nn, \n",
    "            #                                        method=method, \n",
    "            #                                        particle_inc=particle_inc, \n",
    "            #                                        nn_inc=nn_inc)\n",
    "            \n",
    "            frames_partial=frames_partial[::frame_skip]\n",
    "            features_partial = TrajectoryToFeaturesBatchR(frames_partial, \n",
    "                                                    N_nn=N_nn, \n",
    "                                                   method=method, \n",
    "                                                    particle_inc=particle_inc, \n",
    "                                                    nn_inc=nn_inc,N_batch=N_batch)\n",
    "            \n",
    "            #wp: fits pure system features (no idea gas reduction)\n",
    "            incpca.partial_fit(features_partial)\n",
    "            count = count + 1\n",
    "            print \"Frame split: \", count\n",
    "        #wp:resets frame split\n",
    "        clear_output()\n",
    "        frames_split = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving trained models to files and reloading "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saves trained model to avoid retraining every time as well as memory issues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['trained_models/kappa_4.0_pos.pca']"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.externals import joblib\n",
    "folder_model='trained_models'\n",
    "model_ig=folder_model+'/kappa_4.0_pos_ig.pca'; \n",
    "model=folder_model+'/kappa_4.0_pos.pca';\n",
    "#wp: Saves model as file\n",
    "#ideal gas trained model\n",
    "joblib.dump(incpca_ig,model_ig)\n",
    "#full trained model\n",
    "joblib.dump(incpca,model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loads saved models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.externals import joblib\n",
    "#paths to where models are located\n",
    "folder_model='trained_models'\n",
    "model_ig=folder_model+'/kappa_4.0_ig.pca'; \n",
    "model=folder_model+'/kappa_4.0.pca';\n",
    "\n",
    "incpca_ig=joblib.load(model_ig);\n",
    "incpca=joblib.load(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting of models and analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Kappa=4\n",
    "particle_inc = 3 #skip this many particles in generating features; #wp: because it selects index 0 be default, rounds up to integer if not a mulptiple of Ntot\n",
    "N_nn = 3120 #number of nearest nighbors for pca analysis #wp:Should not be bigger than Ntotal\n",
    "nn_inc = 10 #reduces the number of nearest neighbors to include as features\n",
    "split = 1 #chunks the data up so it can be processed by the pca tool if really large\n",
    "N_batch=1\n",
    "frame_skip=400\n",
    "\n",
    "#method = ['distance']\n",
    "method = ['orientation']\n",
    "\n",
    "###########\n",
    "file_data = [(arange(0.6100, 0.9000, 0.01), folder)]\n",
    "etas = []\n",
    "[etas.extend(etas_) for etas_, _ in file_data]\n",
    "etas = array(etas)\n",
    "#print(etas)\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the order parameters based on two nested models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = None\n",
    "features = None\n",
    "raw_OPs = []\n",
    "#split=1\n",
    "\n",
    "for etas_, file_base in file_data:\n",
    "    for eta in etas_:\n",
    "        print 'eta={}'.format(eta)\n",
    "        filename = \"{}/trajectory_{:.4f}.gsd\".format(file_base, eta)\n",
    "        #wp: randomize undefined here?\n",
    "        frames = ReadGSD(filename, shuffle_data=True, randomize=False)\n",
    "        \n",
    "        frames_split = array_split(frames, split)\n",
    "        frames = None\n",
    "        count = 0\n",
    "        means = []\n",
    "        raw_OPs_single = []\n",
    "        for frames_partial in frames_split:\n",
    "            #features_partial = TrajectoryToFeatures(frames_partial, \n",
    "            #                                        N_nn=N_nn, \n",
    "            #                                        method=method, \n",
    "            #                                        particle_inc=particle_inc,            \n",
    "            #                                       nn_inc=nn_inc)\n",
    "            \n",
    "            #if eta > 0.8200:\n",
    "            #    N_nni=N_nn2; particle_inci=particle_inc2; nn_inci=nn_inc2;\n",
    "            #else:\n",
    "            N_nni=N_nn; particle_inci=particle_inc; nn_inci=nn_inc;\n",
    "                \n",
    "            frames_partial=frames_partial[::frame_skip]\n",
    "            #frames_partial=frames_partial[0:50]\n",
    "        \n",
    "            features_partial = TrajectoryToFeaturesBatch(frames_partial, N_nn=N_nni, method=method, \n",
    "                                                         particle_inc=particle_inci, nn_inc=nn_inci,N_batch=N_batch)\n",
    "            \n",
    "            #wp: This returns the 'Z=X*U' matrix with mxn dimensions, where m=samples and n=feature vectors\n",
    "            transformed_features_partial = incpca.transform(incpca_ig.transform(features_partial))\n",
    "            #transformed_features_partial = incpca_ig.transform(features_partial)\n",
    "            #wp:No ideal gas correction\n",
    "            #transformed_features_partial = incpca.transform(features_partial)\n",
    "            raw_OPs_single.append(transformed_features_partial)\n",
    "            #count = count + 1\n",
    "            #print count\n",
    "        clear_output()\n",
    "        raw_OPs.append(raw_OPs_single)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute order parameter stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array, mean, std, abs, asarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "OPs = []\n",
    "OPs_std = []\n",
    "OPs_abs = []\n",
    "\n",
    "for single_eta in raw_OPs:\n",
    "    merged_data = []\n",
    "    \n",
    "    #wp: here single_eta is just the feature vector of samples; the for loop below extracts it as array 'split'\n",
    "    for split_eta in single_eta:\n",
    "        merged_data.extend(split_eta)\n",
    "\n",
    "    #wp: averages over all samples m\n",
    "    OPs.append(mean((split_eta), axis=0))\n",
    "    OPs_std.append(std(split_eta, axis=0))\n",
    "    OPs_abs.append(mean(abs(split_eta), axis=0))\n",
    "\n",
    "OPs = array(OPs)\n",
    "OPs_std = array(OPs_std)\n",
    "OPs_abs = array(OPs_abs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The feature matrix input for each eta is of dimensions 30000 312\n",
      "Number of rows:  33\n",
      "Array size for each eta:  (30000, 312)\n",
      "[ 1.11947629e+01  9.65937508e-01 -1.37175728e+00  1.73796322e+00\n",
      "  1.27214079e+00 -1.57220045e+00 -1.85895768e-01 -2.62107327e-01\n",
      " -2.21219680e-01 -1.97656154e+00  2.47747421e+00 -6.45705281e-01\n",
      " -6.78707255e-01  2.23544623e+00  1.02465918e+00  3.25675861e-01\n",
      "  1.59229822e-01  5.02755892e-01 -2.14336439e+00 -2.23007973e-01\n",
      " -2.06817325e+00 -1.30418100e+00 -6.03491298e-01 -7.85823688e-01\n",
      " -8.37460238e-01 -1.57256774e+00 -4.97127330e-01 -5.88099560e-01\n",
      "  4.68604527e-02  6.38406926e-01 -2.39266226e+00  6.15138687e-01\n",
      " -1.93793227e+00 -8.73091048e-01 -7.51114798e-01 -9.28567745e-01\n",
      "  1.70516763e+00  5.10937699e-01 -6.85135928e-01  1.31490111e+00\n",
      " -4.78726308e-01 -9.93300792e-02  7.08630586e-01 -1.67750055e+00\n",
      "  2.00044644e-01  2.50105819e-03 -2.06217842e+00  2.46159166e+00\n",
      "  1.32675495e+00 -1.74389826e+00 -4.05346150e-01  7.27944600e-01\n",
      "  1.96425236e-01 -9.51987704e-02  5.48833525e-01  8.36747295e-01\n",
      "  1.21889172e+00 -1.90029888e+00 -8.28054663e-01  7.74033315e-01\n",
      "  8.69392373e-01 -2.44219471e+00  5.00224251e-02 -7.52930843e-01\n",
      " -1.50635660e+00  1.25632466e+00  2.34628785e-01 -1.78776962e+00\n",
      " -7.57990019e-01  1.62489510e-01  1.03204875e+00 -2.81829521e-02\n",
      "  9.82100794e-01 -8.82238729e-01  5.92202466e-01 -2.75467352e-02\n",
      "  1.81607088e-01 -1.52320093e-02  7.58797413e-01  7.46719575e-01\n",
      " -4.72209793e-01 -1.30827218e+00 -7.90130633e-01  1.34718629e-01\n",
      "  1.37336281e+00 -3.33040934e-01  1.62485376e-01 -6.14500806e-01\n",
      " -1.09681735e+00 -2.40739997e-01  1.95145404e-01  1.09462949e+00\n",
      " -7.98075655e-01 -1.99753034e-01 -2.56182119e-01  1.76237461e+00\n",
      "  1.02405814e+00 -1.86002229e+00  2.80264200e-01  3.82665486e-01\n",
      " -9.44921723e-01  1.12961251e+00  1.52981969e+00  1.15898030e+00\n",
      "  6.35529412e-02 -1.31769558e+00 -8.00067547e-01 -1.17429113e+00\n",
      "  2.49074951e-01  7.60437846e-02 -7.92428044e-01 -4.23051059e-01\n",
      " -5.65934202e-02 -3.79244414e-01  1.90520807e+00 -1.55802091e-01\n",
      "  4.09494078e-01  2.23314071e+00 -1.03283896e+00 -5.22892932e-01\n",
      " -1.53119478e+00  7.83510848e-01  1.53762316e+00 -6.56064382e-02\n",
      " -3.35182923e-01  1.49910846e-01  1.34492664e+00 -1.90719581e-01\n",
      "  2.65175184e-01  1.57183467e+00  1.42978580e+00 -1.83817849e-01\n",
      " -1.65156618e+00  3.38527737e-01 -7.65455819e-01 -1.54887263e+00\n",
      "  1.26121818e+00 -7.60723007e-01  7.82477784e-01 -8.29476870e-01\n",
      " -1.75769329e+00 -8.17469438e-01 -3.40081624e-01  8.43621427e-01\n",
      " -5.39233953e-01 -1.04141021e+00  1.77874469e-01 -6.21967464e-01\n",
      " -1.62266350e+00  1.07847207e+00  3.64331910e-01  1.33434550e+00\n",
      "  4.26801008e-01 -4.31320501e-01  9.73454015e-01  3.17721408e-01\n",
      "  1.27907119e+00 -1.40423687e-01  2.59532068e-01 -9.76147205e-01\n",
      " -4.41216135e-01 -2.10358486e+00  5.04524773e-01 -3.01516180e-01\n",
      "  8.39324432e-01 -7.44234729e-01 -4.94395849e-01  6.79122410e-01\n",
      " -7.34736607e-01 -5.23415853e-01  1.14428554e+00  6.03139469e-01\n",
      "  2.37956239e-01 -6.95945599e-01 -6.00317372e-01 -4.98920869e-01\n",
      " -3.44308060e-03 -1.75021484e-01  1.03881191e+00  8.32472117e-01\n",
      " -2.26906346e+00  7.33830140e-02  1.29075412e+00  3.95220535e-01\n",
      "  8.79675039e-02  5.58455878e-01  7.80965456e-01 -9.48286098e-01\n",
      " -4.63399552e-01 -8.64649343e-01 -1.09644830e+00  1.42339549e+00\n",
      "  1.21839688e+00 -4.35664257e-01  1.15611893e+00 -1.99032220e+00\n",
      "  1.10764736e+00  5.23987119e-01 -1.87513640e-01  1.45608090e+00\n",
      "  1.30012584e+00  8.31744393e-03  1.07335411e+00  3.36642553e-01\n",
      "  1.89090916e-01 -6.00177697e-01  2.21312242e-01 -1.17472657e-01\n",
      " -1.17958713e+00 -1.00508721e-02 -4.48467363e-02  4.81386316e-01\n",
      "  2.37569231e-01 -5.18082889e-01  6.21375954e-01  2.51181001e+00\n",
      "  4.37682221e-01 -6.58160941e-01  2.18242499e-01 -4.90493222e-01\n",
      "  1.39494950e-02  5.51190510e-01 -7.92028559e-01  5.90418396e-01\n",
      " -5.62029038e-01  1.02844752e+00  1.37642691e+00  2.08596689e+00\n",
      " -7.25128754e-01 -7.19250211e-01  1.27113813e+00  5.08617106e-01\n",
      "  1.35201789e-02 -3.72393657e-01  3.52023728e-01  1.70198116e+00\n",
      " -6.88790509e-01 -1.15996560e+00 -1.08358184e+00 -2.41112671e-01\n",
      "  2.90521882e-01  2.33562781e-01  5.95263345e-01 -2.29742765e+00\n",
      " -2.06745306e-01  8.72355973e-01  9.06504129e-01  4.03343529e-01\n",
      " -1.01060736e-01 -1.00317665e+00 -4.65646198e-01 -4.47616462e-01\n",
      "  5.03606793e-01 -8.50569646e-01 -7.31458571e-02 -6.30880388e-01\n",
      " -7.44227654e-01  1.79906388e+00 -8.58267436e-01 -1.20217145e+00\n",
      " -7.39114311e-02 -7.80217321e-01 -2.56115122e-01 -8.85526510e-01\n",
      "  1.07979042e+00 -8.85859436e-01 -2.35623661e-01  2.93992737e-01\n",
      "  7.35846382e-01  5.44997100e-02  5.38977209e-01  1.64414571e-01\n",
      "  6.61156596e-01 -9.76434622e-01 -6.93711609e-02  3.16788183e-01\n",
      "  1.42896477e+00  4.82764577e-01 -9.36093048e-01 -5.98346163e-01\n",
      " -1.92339966e-01 -5.79673654e-01  1.07991031e+00 -2.36215856e+00\n",
      "  1.64474115e-02  8.51729281e-02 -1.68421880e+00 -1.27272006e+00\n",
      "  1.57314880e+00  6.62927921e-01 -3.17037302e-01  2.42195248e+00\n",
      " -1.17046765e+00  1.29874755e-01 -2.28510123e-01  1.31685620e-01\n",
      "  1.65359660e-01 -2.56683562e-01 -7.81547601e-01 -5.91360284e-01\n",
      " -2.84874860e-01 -1.03182215e+00 -1.15960245e+00 -2.41835941e+00\n",
      "  5.80827097e-01  6.87125909e-01 -4.56057661e-01  1.46423989e-01\n",
      "  1.05422765e+00  6.01580072e-01 -1.02228581e+00  7.04625792e-02]\n",
      "Shape of average order parameter (33, 312)\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'plt' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-c7c494c0fd0f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mraw_OP_samples\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 22\u001b[0;31m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0metas\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mraw_OP_samples\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'o-'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"
     ]
    }
   ],
   "source": [
    "#wp:figuring things out...\n",
    "print \"The feature matrix input for each eta is of dimensions\",len(features_partial),len(features_partial[0])\n",
    "\n",
    "eta_index=0\n",
    "print \"Number of rows: \", len(raw_OPs)\n",
    "print \"Array size for each eta: \", raw_OPs[eta_index][0].shape\n",
    "print raw_OPs[eta_index][0][0]\n",
    "\n",
    "raw_OP_sample=[]\n",
    "raw_OP_samples=[]\n",
    "samplenum=20\n",
    "for mi in range(samplenum):\n",
    "    for eta_index in range(len(raw_OPs)):\n",
    "        raw_OP_sample.append(raw_OPs[eta_index][0][mi][0])\n",
    "    raw_OP_samples.append(asarray(raw_OP_sample))\n",
    "    raw_OP_sample=[]\n",
    "\n",
    "#wp: for the averaged value\n",
    "print \"Shape of average order parameter\", OPs.shape\n",
    "\n",
    "for i in range(len(raw_OP_samples)):\n",
    "    plt.plot(etas,raw_OP_samples[i], 'o-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Check out the results of PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from random import randint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa1762b4250>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAGRlJREFUeJzt3X+MH/Wd3/HnaxfbWAkNJF6l8S9s50x6Tokwt3FScUfbhB8mOWH3CsU5nUokJIsrVlKloWeUq0kdnY6ASnWt3IJ7sZSLjjoEKF2pTl0SSKrcFdglNnbs1MfacPFuuOBgIHeNz3jtd/+Yz9rD19/v7Kz3u/vd/ezrIa2+M5+Z+fo9jHnt+DMzn1FEYGZms0NXpwswM7Op49A3M5tFHPpmZrOIQ9/MbBZx6JuZzSIOfTOzWcShb2Y2izj0zcxmEYe+mdksclGnC2i0YMGCWLZsWafLMDObUV544YWfR0TPWOvVCn1Ja4E/ArqBP46I+1qs90+Bx4CPRsRAarsHuAM4DXwuInZX/VnLli1jYGCgTllmZpZI+ss6640Z+pK6gW3A9cAQ0C+pLyIONqx3CfB54LlS2ypgA/BhYCHwHUlXRMTpujtiZmbtU6dPfw0wGBFHIuJtYCewrsl6XwG+CvxtqW0dsDMiTkbEy8Bg+j4zM+uAOqG/CDhamh9KbWdJuhpYEhH/Y7zbpu03ShqQNHDs2LFahZuZ2fhN+O4dSV3Ag8C/utDviIjtEdEbEb09PWNehzAzswtU50LuMLCkNL84tY26BPj7wPckAfxdoE/SzTW2NTOzKVQn9PuBlZKWUwT2BuC3RxdGxFvAgtF5Sd8DvhgRA5JOAI9IepDiQu5K4Pn2lX/Ok3uGeWD3IX765gkWXjqfu2/8EOtXn9eTZGY2q40Z+hExImkTsJvils0dEXFA0lZgICL6KrY9IOlR4CAwAtw1GXfuPLlnmHue2M+JU8VXD795gnue2A/g4DczK9F0e11ib29vjPc+/Wvue5rhN0+c177o0vn82eZPtKs0M7NpS9ILEdE71npZDMPw0yaBX9VuZjZbZRH6Cy+dP652M7PZKovQv/vGD3HxnHfuyvw53dx944c6VJGZ2fSUReivX72ILb+56uz8okvn84e/daUv4pqZNcgi9AE+/ZGFAPyb31zFn23+hAPfzKyJbEK/eC4MptvdSGZm00k+oZ8+nflmZq1lE/pd6VQ/cOqbmbWSTeiPdu+cceabmbWUT+inDh5375iZtZZP6I9eyHX3jplZS/mFvjPfzKylfEL/bPeOU9/MrJV8Qt9n+mZmY8om9M/dsmlmZq1kE/qjD2ed8am+mVlLtUJf0lpJhyQNStrcZPmdkvZL2ivpB5JWpfZlkk6k9r2SHmr3Dpyrofh05puZtTbm6xIldQPbgOuBIaBfUl9EHCyt9khEPJTWvxl4EFiblh2OiKvaW3bTOgF375iZValzpr8GGIyIIxHxNrATWFdeISJ+UZp9Fx3KXsl375iZVakT+ouAo6X5odT2DpLuknQYuB/4XGnRckl7JH1f0m9MqNoxCHfvmJlVaduF3IjYFhEfBH4P+P3U/CqwNCJWA18AHpH0dxq3lbRR0oCkgWPHjl1wDV2Sn8g1M6tQJ/SHgSWl+cWprZWdwHqAiDgZEa+n6ReAw8AVjRtExPaI6I2I3p6enrq1n6fo3rngzc3Mslcn9PuBlZKWS5oLbAD6yitIWlma/TTwUmrvSReCkbQCWAkcaUfhzQh5lE0zswpj3r0TESOSNgG7gW5gR0QckLQVGIiIPmCTpOuAU8AbwO1p82uBrZJOAWeAOyPi+GTsCADygGtmZlXGDH2AiNgF7Gpo21Ka/nyL7R4HHp9IgePRJXzPpplZhWyeyIXR7h2nvplZK3mFvi/kmplVyir0i1s2zcyslaxCX3jANTOzKlmFPu7eMTOrlFXoj46pb2ZmzWUV+h5wzcysWl6hD34i18ysQl6h7wHXzMwqZRX6Xb6Qa2ZWKavQxwOumZlVyir0i5t3nPpmZq1kFfru3jEzq5ZV6HvANTOzanmFvs/0zcwqZRX6HnDNzKxaVqEPPtM3M6tSK/QlrZV0SNKgpM1Nlt8pab+kvZJ+IGlVadk9abtDkm5sZ/Hn1+FhGMzMqowZ+unF5tuAm4BVwGfKoZ48EhFXRsRVwP3Ag2nbVRQvUv8wsBb4T6MvSp8Mkm/YNDOrUudMfw0wGBFHIuJtYCewrrxCRPyiNPsuzmXvOmBnRJyMiJeBwfR9k6JL8pm+mVmFOi9GXwQcLc0PAR9rXEnSXcAXgLnAJ0rbPtuw7aIm224ENgIsXbq0Tt1NecA1M7NqbbuQGxHbIuKDwO8Bvz/ObbdHRG9E9Pb09FxwDfLdO2ZmleqE/jCwpDS/OLW1shNYf4HbTogv5JqZVasT+v3ASknLJc2luDDbV15B0srS7KeBl9J0H7BB0jxJy4GVwPMTL7s54Vs2zcyqjNmnHxEjkjYBu4FuYEdEHJC0FRiIiD5gk6TrgFPAG8DtadsDkh4FDgIjwF0RcXqS9sXj6ZuZjaHOhVwiYhewq6FtS2n68xXb/gHwBxda4Hh4wDUzs2pZPZEr5NA3M6uQV+gLj7JpZlYhs9D3LZtmZlXyCn3cp29mViWv0Pd9+mZmlfIL/U4XYWY2jWUV+h5wzcysWlah7wHXzMyqZRX6+O4dM7NKWYV+ly/kmplVyir01ekCzMymubxCX/ITuWZmFbIKfQ+4ZmZWLavQ94BrZmbVsgp9POCamVmlrEJf+IlcM7MqtUJf0lpJhyQNStrcZPkXJB2UtE/SdyVdXlp2WtLe9NPXuG07dXkcBjOzSmO+OUtSN7ANuB4YAvol9UXEwdJqe4DeiPilpN8F7gduS8tORMRVba67Ra3u3jEzq1LnTH8NMBgRRyLibWAnsK68QkQ8ExG/TLPPAovbW2Y9PtE3M6tWJ/QXAUdL80OprZU7gG+X5i+WNCDpWUnrL6DG2jzgmplZtVovRq9L0u8AvcA/LDVfHhHDklYAT0vaHxGHG7bbCGwEWLp06YRqcOSbmbVW50x/GFhSml+c2t5B0nXAl4CbI+LkaHtEDKfPI8D3gNWN20bE9ojojYjenp6ece1AQw0eZdPMrEKd0O8HVkpaLmkusAF4x104klYDD1ME/mul9sskzUvTC4BrgPIF4Lbq8vsSzcwqjdm9ExEjkjYBu4FuYEdEHJC0FRiIiD7gAeDdwLckAfwkIm4GfhV4WNIZil8w9zXc9dNWvk/fzKxarT79iNgF7Gpo21Kavq7Fdn8OXDmRAsfDA66ZmVXL6olcD7hmZlYtq9DHA66ZmVXKKvT9RK6ZWbW8Qr/TBZiZTXNZhX7xRG6nqzAzm76yCv1i7B2nvplZK9mFvp/INTNrLbPQ94BrZmZV8gp9/ESumVmVvELfF3LNzCplFfrFE7lOfTOzVrIKfXfvmJlVyyv03b1jZlYpr9DHwzCYmVXJK/R9pm9mVimz0O90BWZm01teoY+7d8zMqtQKfUlrJR2SNChpc5PlX5B0UNI+Sd+VdHlp2e2SXko/t7ez+EYecM3MrNqYoS+pG9gG3ASsAj4jaVXDanuA3oj4CPAYcH/a9r3AvcDHgDXAvZIua1/5jbV6wDUzsyp1zvTXAIMRcSQi3gZ2AuvKK0TEMxHxyzT7LLA4Td8IPBURxyPiDeApYG17Sj+fB1wzM6tWJ/QXAUdL80OprZU7gG+PZ1tJGyUNSBo4duxYjZKa8907ZmbV2nohV9LvAL3AA+PZLiK2R0RvRPT29PRc+J9ffNsFb29mlrs6oT8MLCnNL05t7yDpOuBLwM0RcXI827aLhM/0zcwq1An9fmClpOWS5gIbgL7yCpJWAw9TBP5rpUW7gRskXZYu4N6Q2iaFkG/ZNDOrcNFYK0TEiKRNFGHdDeyIiAOStgIDEdFH0Z3zbuBbKp6Q+klE3BwRxyV9heIXB8DWiDg+KXtCGmVzsr7czCwDY4Y+QETsAnY1tG0pTV9Xse0OYMeFFjgevpBrZlYtqydywU/kmplVySr0u+T+HTOzKlmFvjPfzKxaXqGPu3fMzKpkFfpdXb6Qa2ZWJavQL96R69Q3M2slq9DHT+SamVXKKvQ9nr6ZWbWsQt/dO2Zm1fIKfXfvmJlVyiv0PeCamVmlrELfA66ZmVXLKvTxhVwzs0pZhb7SZzj5zcyayir0u5RemOjMNzNrKqvQT5nvfn0zsxZqhb6ktZIOSRqUtLnJ8msl/VDSiKRbGpadlrQ3/fQ1bttO7t4xM6s25puzJHUD24DrgSGgX1JfRBwsrfYT4LPAF5t8xYmIuKoNtY6pq6uI/TPOfDOzpuq8LnENMBgRRwAk7QTWAWdDPyJeScvOTEKN4+ancs3MmqvTvbMIOFqaH0ptdV0saUDSs5LWj6u6cTrbp+/MNzNrqtaL0Sfo8ogYlrQCeFrS/og4XF5B0kZgI8DSpUsv+A8SvnvHzKxKnTP9YWBJaX5xaqslIobT5xHge8DqJutsj4jeiOjt6emp+9Xn6Tp7945T38ysmTqh3w+slLRc0lxgA1DrLhxJl0mal6YXANdQuhbQbu7eMTOrNmboR8QIsAnYDfwYeDQiDkjaKulmAEkflTQE3Ao8LOlA2vxXgQFJLwLPAPc13PXTVqPdOx50zcysuVp9+hGxC9jV0LalNN1P0e3TuN2fA1dOsMba/HCWmVm1zJ7I9YVcM7MqeYX+6IRD38ysqaxCf/TuHffpm5k1l1Xon+3e6XAdZmbTVWahX3x6wDUzs+byCv306QHXzMyayyv0z3bvOPXNzJrJLPTThDPfzKypvEIfj6dvZlYlq9D3gGtmZtWyCn0PuGZmVi2v0Mf36ZuZVckr9EefyHWnvplZU5mFvsZeycxsFssr9NOn+/TNzJrLKvS70t54wDUzs+Zqhb6ktZIOSRqUtLnJ8msl/VDSiKRbGpbdLuml9HN7uwpvWqcv5JqZVRoz9CV1A9uAm4BVwGckrWpY7SfAZ4FHGrZ9L3Av8DFgDXCvpMsmXnarWotPD7hmZtZcnTP9NcBgRByJiLeBncC68goR8UpE7APONGx7I/BURByPiDeAp4C1bai7km/eMTNrrk7oLwKOluaHUlsdE9l23Lo8+I6ZWaVpcSFX0kZJA5IGjh07NoHvKT7du2Nm1lyd0B8GlpTmF6e2OmptGxHbI6I3Inp7enpqfvX5fCHXzKxandDvB1ZKWi5pLrAB6Kv5/buBGyRdli7g3pDaJoXfkWtmVm3M0I+IEWATRVj/GHg0Ig5I2irpZgBJH5U0BNwKPCzpQNr2OPAVil8c/cDW1DYp3L1jZlbtojorRcQuYFdD25bSdD9F102zbXcAOyZQ4zik7h2HvplZU9PiQm67uHvHzKxaVqHvAdfMzKrlFfrp0yf6ZmbN5RX67t4xM6uUVeiPPpHryDczay6r0McDrpmZVcoq9D3yjplZtaxC/2z3js/0zcyayir0/USumVm1vELfA66ZmVXKKvTPPpHrt6iYmTWVVeifvXuns1WYmU1bWYW+POCamVmlvELf9+mbmVXKKvT9RK6ZWbWsQt+3bJqZVcsr9NNn+FzfzKypWqEvaa2kQ5IGJW1usnyepG+m5c9JWpbal0k6IWlv+nmoveWfVwcAvmPTzKy5MV+XKKkb2AZcDwwB/ZL6IuJgabU7gDci4lckbQC+CtyWlh2OiKvaXHeLWotPX8g1M2uuzpn+GmAwIo5ExNvATmBdwzrrgK+n6ceAT6oDr7HygGtmZtXqhP4i4Ghpfii1NV0nIkaAt4D3pWXLJe2R9H1JvzHBeit5wDUzs2pjdu9M0KvA0oh4XdKvAU9K+nBE/KK8kqSNwEaApUuXXvAf5rt3zMyq1TnTHwaWlOYXp7am60i6CHgP8HpEnIyI1wEi4gXgMHBF4x8QEdsjojcient6esa/F4mfyDUzq1Yn9PuBlZKWS5oLbAD6GtbpA25P07cAT0dESOpJF4KRtAJYCRxpT+nn8ztyzcyqjdm9ExEjkjYBu4FuYEdEHJC0FRiIiD7ga8A3JA0Cxyl+MQBcC2yVdAo4A9wZEccnY0eg1L0zWX+AmdkMV6tPPyJ2Absa2raUpv8WuLXJdo8Dj0+wxtrcvWNmVi2vJ3LP3iTq1Dczayar0O/yE7lmZpWyCn3fsmlmVi2v0E+fHnDNzKy5vELf3TtmZpUyC/3i08MwmJk1l1fod7oAM7NpLqvQP3f3js/0zcyaySr0ffeOmVm1vELfT+SamVXKK/Q99o6ZWaUsQ999+mZmzWUV+k8d/BkA//qxfVxz39M8uadx2H8zs9ktm9B/cs8w9337/56dH37zBPc8sd/Bb2ZWkk3oP7D7ECdHzryj7cSp03y570CHKjIzm34m+x25U+anb56ghzd4eO6/5yJGzi04Dfu23M28rtOs6Po5cxasgO4ZtNunT8Ebr8Bly133VHDdU8t1v1P3XLjtT+GS97fvOxvUqlbSWuCPKN6c9ccRcV/D8nnAnwC/BrwO3BYRr6Rl9wB3AKeBz0XE7rZVX7Lw0vnc+Tc7WK3B5itEUUH8bIae+bvuqeW6p5brLghefnwLKz77cHu/t2TM0E/vuN0GXA8MAf2S+iLiYGm1O4A3IuJXJG0AvgrcJmkVxasTPwwsBL4j6YqION3uHbn3mvn84+88XXqRipnZzLPw5ce571ufYfOt/2hSvr9On/4aYDAijkTE28BOYF3DOuuAr6fpx4BPqhjych2wMyJORsTLwGD6vra74fVvMEdnxl7RzGwam8MIC1/8j5N2E0qd0F8EHC3ND6W2putExAjwFvC+mttO3F//Fex/rO1fa2Y21boV/LPu7/O1//nspHz/tLh7R9JGSQOSBo4dOzb+L/j+/cWFFTOzDMxhhFv/3yOT8t11Qn8YWFKaX5zamq4j6SLgPRQXdOtsS0Rsj4jeiOjt6empX/2ooecBd+2YWR66FXx8zuFJ+e46d+/0AyslLacI7A3Abzes0wfcDvwf4Bbg6YgISX3AI5IepLiQuxJ4vl3Fn3XnD1ouenLPMPc8sY8Tp/xLwcxmhvlzuvnD37qSKybhu8cM/YgYkbQJ2E1xy+aOiDggaSswEBF9wNeAb0gaBI5T/GIgrfcocBAYAe6ajDt3qqxfvYj1q9t/GcHMrF2e3DPMA7sP8dM3T7Dw0vncfeOHJi23NN1eLdjb2xsDAwOdLsPMbEaR9EJE9I613rS4kGtmZlPDoW9mNos49M3MZhGHvpnZLOLQNzObRabd3TuSjgF/OYGvWAD8vE3ldIr3ofNmev3gfZgupmofLo+IMZ9unXahP1GSBurctjSdeR86b6bXD96H6WK67YO7d8zMZhGHvpnZLJJj6G/vdAFt4H3ovJleP3gfpotptQ/Z9embmVlrOZ7pm5lZC9mEvqS1kg5JGpS0udP11CXpFUn7Je2VNJDa3ivpKUkvpc/LOl1nmaQdkl6T9KNSW9OaVfgP6bjsk3R15yo/p8U+fFnScDoWeyV9qrTsnrQPhyTd2Jmq30nSEknPSDoo6YCkz6f2GXMsKvZhxhwLSRdLel7Si2kf/m1qXy7puVTrNyXNTe3z0vxgWr5sSguOiBn/QzHk82FgBTAXeBFY1em6atb+CrCgoe1+YHOa3gx8tdN1NtR3LXA18KOxagY+BXwbEPBx4LlO11+xD18Gvthk3VXp79Q8YHn6u9Y9DfbhA8DVafoS4C9SrTPmWFTsw4w5Fum/57vT9BzgufTf91FgQ2p/CPjdNP0vgIfS9Abgm1NZby5n+nVe3j6TlF80/3VgfQdrOU9E/G+K9yaUtap5HfAnUXgWuFTSB6am0tZa7EMr64CdEXEyIl4GBin+znVURLwaET9M038N/JjiHdQz5lhU7EMr0+5YpP+ef5Nm56SfAD4BjL68u/E4jB6fx4BPStIUlZtN6E/NC9gnRwD/S9ILkjamtvdHxKtp+q+A93emtHFpVfNMOzabUtfHjlK32rTfh9RFsJriLHNGHouGfYAZdCwkdUvaC7wGPEXxL5A3I2IkrVKu8+w+pOVvAe+bqlpzCf2Z7Ncj4mrgJuAuSdeWF0bxb8AZdYvVTKw5+c/AB4GrgFeBf9fZcuqR9G7gceBfRsQvystmyrFosg8z6lhExOmIuIriPeBrgL/X4ZJayiX0a72AfTqKiOH0+Rrw3yj+wvxs9J/d6fO1zlVYW6uaZ8yxiYifpf95zwD/hXPdBtN2HyTNoQjLP42IJ1LzjDoWzfZhJh4LgIh4E3gG+AcU3Wejr6Qt13l2H9Ly9wCvT1WNuYT+2Ze3pyvkGyhe1j6tSXqXpEtGp4EbgB9x7kXzpM//3pkKx6VVzX3AP093jnwceKvU9TCtNPRv/xOKYwHFPmxId10sB1YCz091fY1SP/DXgB9HxIOlRTPmWLTah5l0LCT1SLo0Tc8Hrqe4NvEMcEtarfE4jB6fW4Cn07/IpkYnr3q384fizoS/oOhL+1Kn66lZ8wqKOxFeBA6M1k3Rv/dd4CXgO8B7O11rQ93/leKf3Kco+irvaFUzxZ0N29Jx2Q/0drr+in34RqpxH8X/mB8orf+ltA+HgJs6XX+q6dcpum72AXvTz6dm0rGo2IcZcyyAjwB7Uq0/Arak9hUUv5AGgW8B81L7xWl+MC1fMZX1+olcM7NZJJfuHTMzq8Ghb2Y2izj0zcxmEYe+mdks4tA3M5tFHPpmZrOIQ9/MbBZx6JuZzSL/H0JVnuEE8SNJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "####wp: Comparing the component variance for the hard disk and ideal gas system\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.plot(incpca.explained_variance_ratio_,marker='o')\n",
    "plt.plot(incpca_ig.explained_variance_ratio_,marker='^')\n",
    "\n",
    "#plt.ylim(-0.1,0.1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explained variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI4CAYAAAB3HEhGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3XuYXGWZ7v/7qa7uqnSnKyEh3ZAESAAJhyAmHJTBaGUrA4yowROizDgwwOwfoqKjjnuD0mxw0HEOig7+BnDAPQwi6oBs1OAWUYMMAiacAsFDQjiEpMmJdDpJH5/9R1V3KqtP1Z1atWqt/n6uq6+utWqtp56EJOTO+673NXcXAAAAAGD/paJuAAAAAACSgoAFAAAAABVCwMKEmdkJZna9mT1jZq+Z2U4z+52Z3Wxmb466PwAAAKDajGewMBFm9kVJfyupbpTLbpZ0qbv3VKcrAAAAIFoELIybmX1Z0meLh72Sbpf0C0l9kt4o6S8lNRbf/w93P7/KLQIAAACRIGBhXMzsNEkPFg+7JJ3l7g8ErjlehcA1o3jqXHe/s2pNAgAAABHhGSyM1xdLXv+vYLiSJHd/StLHS079XehdAQAAADWAEawEM7MjJZ0oaa6kKZI6JK2V9LC7vzqBeodIWi/JJO2WdLC7vzbCtXWSni9+tiSd6u4Pj/czAQAAgDhJR91A0pnZXEknSTq5+P0k7Z06N2Cpu/+iQp+XknSBpE9KOm6Ey/rN7H5J1w03AjWKs1QIV5K0YqRwJUnu3mdmP5F0cfHUOyURsAAAAJBoBKwQmNkXtDdUHVTFz50j6XuSTh3j0pSk0yWdbmY3SvqYu3eX8REnlLx+pIzrH9HegPX6Mq4HAAAAYo2AFY6rq/2BxXD1oKR5w7ztKkwPzA3z3iWS5pjZOWUsp76g5PXaMtpaV/L66DKuBwAAAGKNRS4SwMzqJf1AQ8PVckl5SVPcfZqkaZL+XNKawHXvkPSlMj7qgJLXm8u4vvSa6WVcDwAAAMQaASt83ZJWSrpRhely7wzhMy5RYf+pUl9x97Pc/Zfu3iVJ7r7D3W+TdIqk4LNXl5vZSM9sDWgueb27jL52lbwebvQMAAAASBSmCIbjFkmPSXpU0hOlzzeZ2bxKfpCZNUq6MnD6Z+7+2eGulyR37zCzD0haLamleDqlwhLsyyrZHwAAADCZELBC4O4XVvHjLtC+C2m4pMvGusndN5vZ5yX9a8npd5vZse7+zAi3dZS8nlJGb40lr3eUcT0AAAAQa0wRjL/3BI4fcPfnyrz3dkk7x6hXalvJ6wPLqF96zfYyewIAAABii4AVY2Z2gKS3BE5/t9z73X2npHsDp0ebIlga3A4v4yPml7wOLqwBAAAAJA4BK97eqqHTPH89zhrB6080s5FW/Hui5PUpZdQuvebJcXUFAAAAxBABK96Cq/7tkDTS81MjeXiYc8eMcO1PVHjGS5LebGbTRipqZnWSzio5FRwpAwAAABKHgBVvwSD0vLv7sFeObLgNg4cNWO7+oqQVxcNGSR8dpe4HJc0t+YzhghwAAACQKASseFsQOH5hvAXcfaukzjHqlipdEv4LZpYPXmBmx0v6WsmpKyYQ/AAAAIDYYZn2eJsROH5lgnVekXTkKHUHufsKM/tHSX8jKSPpp2Z2u6RfSOpTYcPjC7R3ifY73P2OCfYFAAAAxAoBK96mBo53TbDO7jHqBn1GUq+kT0uql/SR4lfQLZL+eoI9AQAAALFDwIq3psDxngnWGVfAKk73+5yZ3SHpIklvkzRHUp0Ko2G/knSLu68YuUqBmTF1EAAAAFXl7hZWbQJWvE0JHHdPsE7XGHWH5e6PS7psgp9ZWmd/SwAAAABlMQstW0likYu4C45YNUywTmaMugAAAADKQMCKt52B4+wE6wRHrIJ1AQAAAJSBgBVvwSDUOOxVYwsGrI4J1gEAAAAmNQJWvG0LHB80wTrB+4J1AQAAAJSBRS7i7TlJJ5YcHzbeAmY2Q0NXDXxuf5oar7a2tsHX+Xxe+Xy+mh8PAAAAVAwBK96eDRzPMzPz8S3LN7+MuqEqDVgAAABAnDFFMN5WB45zko4ZZ403DXOuqgELAAAASAoCVrz9UlJf4Nybx1kjeP0qd+cZLAAAAGACCFgx5u5bJa0InD633PvNbKqkswOn79rfvgAAAIDJioAVf/8ZOF5qZgvKvPdDGrrARbAeAAAAgDIRsOLvFkntJccm6Rtj3WRmMyVdEzh9r7sHn+sCAAAAUCZWEYw5d99pZtdKur7k9NvN7Evu/rnh7jGzZkl3SmopLSXpivA6HRnLtAMAACApbHwreqMcxSl6S0d4e6akawPn/lnS70a4/jF3f2yMz6uX9JCkkwJvLZd0naSH3b3bzHIqPHN1pYauNvg1d798tM8Jw/hXlQcAAAAmzszk7hZaff5yW3lm9pcqTN2rhKvdva2MzzxE0oOSDh3m7X5JHZKmjXD7cknvcveeiTY5UQQsAAAAVFPYAYtnsBLC3V+UdJqk3wzzdkojh6tvSVoWRbgCAAAAkoaAlSDu/pIKIesSSWtGu1TSzyW9zd0vcveuavQHAAAAJB1TBBPMzI6SdKKkuZKyknZKWivpv9y9fbR7q4UpggAAAKgmnsFCohGwAAAAUE08gwUAAAAAMcE+WIgc+2ABAAAgKZgiiEgxRRAAAADVxBRBAAAAAIgJAhYAAAAAVAgBCwAAAAAqhIAFAAAAABVCwAIAAACACiFgAQAAAECFELAAAAAAoELYaBiRY6NhAAAAJAUbDSNSbDQMAACAamKjYQAAAACICQIWAAAAAFQIAQsAAAAAKoSABQAAAAAVQsACAAAAgAohYAEAAABAhRCwAAAAAKBCCFgAAAAAUCEELAAAAACokHTUDQBtbW2Dr/P5vPL5fGS9AAAAAPvD3D3qHjCJmZnzaxAAAADVYmZydwurPlMEAQAAAKBCCFgAAAAAUCEELAAAAACoEAIWAAAAAFQIAQsAAAAAKoSABQAAAAAVQsACAAAAgAohYAEAAABAhaSjbgAYzfotnbppxVrdvWqDOrt61ZRJa9mi2bp4yeE6bGZT1O0BAAAA+zB3j7oHTGJm5iP9GnzguXZdettK9fT1q7d/7zXplKm+LqUbzl+spQtaqtUqAAAAEsDM5O4WWn0CFqI0UsBav6VTZ351hXb39I1475T6Oi2/fAkjWQAAAChb2AGLKYKIXFtb2+DrfD6vfD6vm1asVU9f/6j39fT16+YV63TNsoUhdwgAAACUhxEsRGqkEayFV92nnV29Y94/NZPW01efEUZrAAAASKCwR7BYRRA1qbOMcCVJnd3lXQcAAABUAwELNakpU97s1aYGZrkCAACgdhCwUJOWLZqtdGr0kdt0ynTOojlV6ggAAAAYGwELNeniJYervm70X571dSldtGR+lToCAAAAxkbAQk06bGaTbjh/sabU1w0ZyUqnTFPq63TD+YtZoh0AAAA1hVUEEanRNhqWCvth3bxinb772Ivq7u3X1Exa5yyao4uWzCdcAQAAYNzYaBiJNlbAGnD3qpd1+Xcf1/1/81YdMWtqFToDAABAErFMOyCpJZeRJG3asSfiTgAAAICREbAQC625rCSpfUdXxJ0AAAAAIyNgIRYGAhYjWAAAAKhlBCzEwtRMWk0NddrECBYAAABqGAELsdGay2pTByNYAAAAqF0ELMRGSy6jdqYIAgAAoIYRsBAbLc1ZtXcwRRAAAAC1i4CF2GjNZbRpxx6xdxsAAABqVTrqBoC2trbB1/l8Xvl8ftjrWnNZ7enp1449vZo2pb46zQEAAADjQMBC5EoD1mhaBvfC2kPAAgAAQE1iiiBio7U5I0ks1Q4AAICaRcBCbLDZMAAAAGodAQux0ZIrjmCxFxYAAABqFAELsdHYkFZzNq12pggCAACgRhGwECutuSxTBAEAAFCzCFiIlZbmDJsNAwAAoGYRsBArjGABAACglhGwECstuYzad3TJ3aNuBQAAABiCgIVYaW3OqruvX9t39UTdCgAAADAEAQuxMrgXFku1AwAAoAYRsBArrQN7YbFUOwAAAGoQAQuxMjiCxUIXAAAAqEEELMTKrObCCFY7AQsAAAA1iICFWMnW12l6Yz1TBAEAAFCTCFiIncJmw4xgAQAAoPYQsBA7hc2GGcECAABA7SFgIXZamrM8gwUAAICaRMBC7LTmMmrv6FJ/v0fdCgAAALAPAhZipzWXVW+/a+uu7qhbAQAAAPZBwELs7N1smGmCAAAAqC0ELMROS3Gz4XYWugAAAECNIWAhdlqLAYsRLAAAANSadNQNAG1tbYOv8/m88vn8qNfPmjowRZARLAAAANQWAhYiVxqwytGQTmlGUwObDQMAAKDmMEUQsdTSnGEECwAAADWHgIVYas1lGcECAABAzSFgIZZacxkWuQAAAEDNIWAhllpzWb3a0aW+fo+6FQAAAGAQAQux1JLLqt+lLTt5DgsAAAC1g4CFWGptZql2AAAA1B4CFmKJzYYBAABQiwhYiKXBgMVKggAAAKghBCzE0oFTG2QmtTNFEAAAADWEgIVYStelNLMpw15YAAAAqCkELMRWYS8sRrAAAABQOwhYiK3WXJZFLgAAAFBTCFiILUawAAAAUGsIWIitluastnR2qaevP+pWAAAAAEkELMRYay4rd2nzTkaxAAAAUBsIWIit1lxGkpgmCAAAgJpBwEJsDW42zEIXAAAAqBEELMRWS3NhBKudgAUAAIAakY66AcSPmZmkIyUtLn6dKGmRpBkD17i7hd3HzKkZpUxq72CKIAAAAGoDAQsT8Q+SPhV1E3Up06zmDFMEAQAAUDOYIoiJqAsc75L0RBSNFDYbZgQLAAAAtYERLEzEM5K+Kmll8WuNpEMkrat2Iy3NWb20bVe1PxYAAAAYFgErAmZ2pArPLc2VNEVSh6S1kh5291ej7K0c7n5j8Fzhsazqa81ltPKFbZF8NgAAABA0aQOWmc2VdJKkk4vfT1LJIg1FS939FxX6vJSkCyR9UtJxI1zWb2b3S7rO3R+oxOcmXWsuq62d3erq7VMmHZy5CAAAAFTXpApYZvYF7Q1VB1Xxc+dI+p6kU8e4NCXpdEmnm9mNkj7m7t1h9xdnA5sNv9rRpbkHNEbcDQAAACa7SRWwJF1d7Q8shqsHJc0b5m1XYXpgbpj3LpE0x8zOcfee8DqMt5bBzYYJWAAAAIgeqwiGyMzqJf1AQ8PVckl5SVPcfZqkaZL+XIXFIkq9Q9KXwu0y3thsGAAAALVkMgesbhVWwLtR0sWS3hnCZ1wi6Y2Bc19x97Pc/Zfu3iVJ7r7D3W+TdIqk4LNXl5vZSM9sTXqtxREsNhsGAABALZhsUwRvkfSYpEclPVH6fJOZzavkB5lZo6QrA6d/5u6fHeked+8wsw9IWi2ppXg6JemLkpZVsr+kmNHYoHTK2GwYAAAANWFSjWC5+4XufoO7P1qFxSMu0L4Labiky8a6yd03S/p84PS7zezY0e4zs1vNzCv4VbVFQPZHKmVqac6w2TAAAABqwqQKWFX2nsDxA+7+XJn33i5p5xj1UNSSy6q9gxEsAAAARG+yTRGsCjM7QNJbAqe/W+797r7TzO6V9MGS08skXTvKbT+VtL3sJsfWWcFaoWrNZbRuc2zaBQAAQIIRsMLxVg39uf31OGv8WvsGrBPNbLq7Dxui3P12FUa+Jp3WXFYPr90adRsAAAAAUwRDElz1b4ekZ8ZZ4+Fhzh0zsXaSrTWX1Wu7e7Snpy/qVgAAADDJEbDCEQxCz7u7j7PG2jLqQqV7YbHQBQAAAKJFwArHgsDxC+Mt4O5bNfQ5qGBdqLDIhSRtYqELAAAARIyAFY4ZgeNXJlgneF+wLlRY5EJiBAsAAADRY5GLcEwNHO+aYJ3dY9SNhJlNl/TpwOlpgWuGW/HwH0ZapGN/tDYXR7DYbBgAAAARI2CFoylwPNG/+ddkwJI0XdIVY1wz3Ps3a5il5Nva2gZf5/N55fP58TXTWK+GuhRTBAEAABA5AlY4pgSOuydYJzjnLVg3EUoD1kSYmVpyGaYIAgAAIHIErHDskdRYctwwwTqZYepGzt2fl2RR91GqNZdliiAAAAAixyIX4dgZOM5OsE5wxCpYF0WtuQwBCwAAAJEjYIUjGIQah71qbMGA1THBOonX0pxliiAAAAAiR8AKx7bA8UETrBO8L1gXRa25rDq6etXZ1Rt1KwAAAJjECFjheC5wfNh4C5jZDA1dNTBYF0UtzcW9sDoYxQIAAEB0WOQiHM8GjueZmbm7j6PG/DLqJsL+LtMuFUawJKl9xx7NPzC4Sj4AAABQHQSscKwOHOckHSPpmXHUeNMw5xIfsCaqNVcYwdrECBYAAAAixBTBcPxSUl/g3JvHWSN4/Sp35xmsEbSUjGABAAAAUSFghcDdt0paETh9brn3m9lUSWcHTt+1v30lWS6bVrY+xVLtAAAAiBQBKzz/GTheamYLyrz3Qxq6wEWwHkqYWXGzYaYIAgAAIDoErPDcIqm95NgkfWOsm8xspqRrAqfvdffgc10IaG3OMoIFAACASLHIRUjcfaeZXSvp+pLTbzezL7n754a7x8yaJd0pqaW0lKQrwus0epVYRVCSWnIZrd6wozJNAQAAABMwqQJWcYre0hHenjnMuXeZ2dEjXP+Yuz82xkf+/5L+QtJJJef+1sxOkHSdpIfdvdvMcio8c3WlCqsNlrre3Z8c43NirRKrCEqFpdp/vqZd7i4zq0hNAAAAYDwmVcCSdKqkb47j+k+O8t7VkkYNWO7eY2bvkfSgpENL3jqz+NVvZh2Spo1QYrmkz5Tf7uTW0pzRru4+7ezqVXO2Pup2AAAAMAnxDFbI3P1FSadJ+s0wb6c0crj6lqRl7t4TVm9JM7jZMHthAQAAICIErCpw95dUCFmXSFoz2qWSfi7pbe5+kbuTFMahZWCzYRa6AAAAQEQm1RRBd79V0q0RfXafpJsk3WRmR0k6UdJcSVlJOyWtlfRf7t4+chWMZnAEi6XaAQAAEJFJFbBqhbv/TtLvou4jaQYCFiNYAAAAiApTBJEYUzNpNTXUsdkwAAAAIsMIFiJXqX2wpMIo1qYORrAAAAAQDQIWIlepfbCkwkIX7UwRBAAAQESYIohEac1lmSIIAACAyBCwkCgtzRlt2rFH7h51KwAAAJiECFhIlNZcVl29/dqxpzfqVgAAADAJEbCQKC2De2HxHBYAAACqj4CFRGltzkgSz2EBAAAgEgQsJAqbDQMAACBKBCwkSkuuOILFXlgAAACIAPtgIXKV3Gi4sSGt5mxa7UwRBAAAQAQIWIhcJTcalgb2wmIECwAAANXHFEEkTmsuQ8ACAABAJAhYSJyW5iyrCAIAACASBCwkTksuo1c7uuTuUbcCAACASYaAhcRpbc6qu69f23f1RN0KAAAAJhkCFhJncC8slmoHAABAlRGwkDitA3th8RwWAAAAqoyAhcQZHMFiJUEAAABUGQELiTOruTCC1U7AAgAAQJURsJA42fo6TW+sZ4ogAAAAqi4ddQNAW1vb4Ot8Pq98Pr/fNVubs0wRBAAAQNURsBC50oBVKS25jDZ1MIIFAACA6mKKIBKppTnLM1gAAACoOgIWEqk1l9GrHV3q7/eoWwEAAMAkQsBCIrXmsurtd23d1R11KwAAAJhECFhIpL2bDTNNEAAAANVDwEIitRQ3G25nqXYAAABUEQELidRaDFiMYAEAAKCaCFhIpFlTB6YIMoIFAACA6iFgIZEa0inNbGrQpg5GsAAAAFA9BCwk1qzmDHthAQAAoKoIWEis1lyWKYIAAACoKgIWEqs1l1E7UwQBAABQRemoGwDa2toGX+fzeeXz+YrUbc1l9WpHl/r6XXUpq0hNAAAAYDQELESuNGBVUksuq36XtuzsGtwXCwAAAAgTUwSRWK3NLNUOAACA6iJgIbHYbBgAAADVRsBCYg0GLBa6AAAAQJUQsJBYB05tkBlTBAEAAFA9BCwkVroupZlNbDYMAACA6iFgIdFacxmewQIAAEDVELCQaK25rNo7mCIIAACA6iBgIdEKI1gELAAAAFQHAQuJ1tKc1ZbOLvX09UfdCgAAACYBAhYSrTWXlbu0eSejWAAAAAgfAQuJ1prLSGKpdgAAAFQHAQuJNrjZMCsJAgAAoAoIWEi0luIIFnthAQAAoBqqHrDM7PPV/kxMXjObMkoZUwQBAABQHVGMYF1lZmdF8LmYhOpSplnNbDYMAACA6khH8JkpSbeZ2Snu/sfx3GhmR7v7mpD6QkTa2toGX+fzeeXz+YrWZ7NhAAAAVEsUAWu5pDMl3WVmb3L3XeXcZGZ/Kum7kg4IszlUX2nACkNLc1YvbSvrlxkAAACwX6KYIvghSc9LOk7SLeXcYGaXSvqRpFx4bSGpWnMZRrAAAABQFVUPWO6+XdJ7JXVJep+ZfXaka63geklfl1QnaVxTCgGpMEVwa2e3unr7om4FAAAACRfJMu3uvkrSpZJM0hfN7O3Ba8ysWYVRq48Wr/ulpDdVs08kw8Bmw68yigUAAICQRbYPlrvfKulGFUam7jCzwwbeK75+SNIZKoSrmyWd7u5bI2gVMdcyuNkwAQsAAADhinqj4Y9LekzSDBUWvcia2WmSHpF0rCSX9Cl3v8TdeyPsEzHW2lwIWGw2DAAAgLCFtoqgmX1Z0uOSVkl6zt09eI27d5vZeyWtlHSCpPslLZaUkdQh6Tx3/3FYPWJyaClOEWQvLAAAAIQtzGXaP6PCCJQk7Tazp1UIW48Xv550993u/qKZfVCF5dvfpMKUwHWS3unuz4TYHyaJGY0NSqdMm3gGCwAAACELM2DtljSl+LpR0imSTi55v9/Mfq+9oev/qrA/1gpJ73H3LSH2hkkklTK1NGfUzjNYAAAACFmYAWuqpAWS3hD4aim+Xyfp6OI1Hyyec0lHSLrVzB6X9ISkx939DyH2iUmgJZdVewdTBAEAABCu0AJW8ZmrNcWvOwbOm9nBGhq6jtDeBTdmSzpY0p+V3NMp6SlJq9z9srB6RnK15jJat7kz6jYAAACQcGGOYA3L3V+R9IqknwycM7MmFRa5KA1dCyVli5dMlXSqCs9oEbAwbq25rB5eyyr/AAAACFfVA9Zw3L1ThX2vHho4Z2YDUwhLQ9cJkTSI2GvNZfXa7h7t6elTtr4u6nYAAACQUDURsIbj7n2SVhe//iPidhBzLc2Fpdrbd3Tp0JmNEXcDAACApIp6o2GgKlpzhdmmm1joAgAAACEiYGFSYLNhAAAAVAMBC5NCa3NxBIu9sAAAABAiAhYmhemN9WqoS6mdESwAAACEqKYClpnNMrP/ZmafiLoXJIuZqSWXUXsHI1gAAAAITySrCJrZFEnHSTo+8DWr5LKvRdAaEqw1l+UZLAAAAIQq1IBlZilJr9PeALWw+P1wSRa8vOS1h9kXJqfWXEbPbeyIug0AAAAkWMUClpkdrKEjUsdIygQvDRx7yffge0DFtDRnteJ3m6NuAwAAAAk27oBlZs3aOxJVOip1QPDSwLEHvlvg+6uSnix+PTHevoCxtOay6ujqVWdXr5oyNbvHNgAAAGJszL9lmlla0gckvV/SGyQdOtxlw5wbKUj1SFqjQogaCFRPuvvG8ttGkrS1tQ2+zufzyufzoXxOa3EvrPaOLs0nYAEAACAEo/4t08wyku6W9KcDpwKXBEelVHKdSdqofUelnpT0rLv37kfPSJjSgBWmlsG9sPZo/oFNVflMAAAATC5j/TP+X0k6o+R4pFGp4HsvSnqPu/92vzsEKmRgBIuVBAEAABCWsfbBOr/4vTQ8PSvpNkmfk/ROSYsl3aO9YcslzZX0czO7tKLdAvuhJVcYwWrfwV5YAAAACMdYI1iHa+/qfmskvd/dVw9z3TIze6+k6yUdXLynWdLXzeyDki5y999Vrm1g/HLZtLL1KbV3MIIFAACAcIw1gjWw8a9L+psRwlXhAvcfSDpW0s0l95ikN0t63Mz+h5nV7We/wISZWXGzYUawAAAAEI6xAlbpM1bPjlXM3V9z90skLZX0+4HTkrKSrpX0qJktmkijQCW0Nmd5BgsAAAChGStg3aPCghUmqb/cou7+K0mvl/R3knq19xmuN0j6jZl9qbhCIVBVLbmM2jsYwQIAAEA4Rg1Y7r7M3edJmilpw3gKu3u3u18p6URJj6oQ0lyF574+I+lJM3vLRJoGJqowRXCP3IM7CwAAAAD7b6wRLEmSu2+b6N5V7v60pFMlXS6ps+St16mw0uA3zax5IrWB8WrNZbSru087u9iKDQAAAJVXVsDaX15wvQqLYPxYe0ezUpIukfSMmZ1djV4wue3dbJhpggAAAKi8qgSsAe7+krufLelDkl4teWuOpB+a2e1mdmA1e8Lk0lLcbLidhS4AAAAQgqoGrAHufoekYyR9W3tHs0zSuZKeNbMPR9EXkq+1uNnwJvbCAgAAQAgiCVjS4HNdF0g6XdK6krdmSvrf0XSFpBsIWO1MEQQAAEAIIgtYA9z9fkkLJX1FUl/E7SDhpmbSamqo4xksAAAAhCLygCVJ7r7H3f9W0imSfht1P0i21lyWKYIAAAAIRTrqBkq5++Nm9kZJn4i6FyRXSy7DIhcAAAAIRU2MYJVy9353/+eo+0ByFTYbZoogAAAAKq/mAhYQtkLA2iN3j7oVAAAAJAwBC5NOS3NGXb392rG7N+pWAAAAkDA19QwW4sHMmiS9TdJbJZ0o6SgVltfvk7RF0pOS7pP07+6+Lao+R9JSshfWtMb6iLsBAABAkjCChXExs09JelXSDyV9SoWQdbCkBklTJM2V9GeSviZprZn9RUStjqi1OSNJ2sRCFwAAAKgwRrAwXkepEKSkQtD6uaRHJG2U5JKOkXS+pPmSpkv6tpk1u/u/RNDrsNhsGAAAAGEhYGG8XIVQ9c+Slrv7kAeZzOxaSd+UdGHx1D+Y2Y/dfV312hxZS644gsVeWAAAAKgwpghivK5w97e5+73DhStJcvduSRdLeqJ4Kivpz6vV4FgaG9JqzqYZwQIAAEDFMYIVATM7UoXFIeaqMN2uQ9JaSQ+7+6tR9jYWd99a5nX9ZvZdSScUT50w2vXVNrBUOwC6ymmkAAAgAElEQVQAAFBJkzZgmdlcSSdJOrn4/SRJMwKXLXX3X1To81KSLpD0SUnHjXBZv5ndL+k6d3+gEp8bsddKXjdG1sUwWnMZAhYAAAAqblIFLDP7gvaGqoOq+LlzJH1P0qljXJqSdLqk083sRkkfK063i6vXl7x+PqomhtPanNVv1pU1GAcAAACUbVIFLElXV/sDi+HqQUnzhnnbVZgemBvmvUskzTGzc9y9J7wOw2FmB0r6YMmp/xNVL8NpyWXV3rFH7i4zi7odAAAAJASLXITIzOol/UBDw9VySXlJU9x9mqRpKiwCsSZw3TskfSncLkPzLyr8uCTpN5J+EmEvQ7Q0Z9TT59q2K3bZFQAAADVsMgesbkkrJd2owop37wzhMy6R9MbAua+4+1nu/kt375Ikd9/h7rdJOkVS8Nmry81spGe2apKZXSHpA8XDTkkXubtH2NIQA3th8RwWAAAAKmmyTRG8RdJjkh6V9ETp801mNq+SH2RmjZKuDJz+mbt/dqR73L3DzD4gabWkluLplKQvSlo2xufdKukjE254qIPdfeN4bzKzCyRdUzx0FcLV0xXsqyJai3thtXd06ZiDI24GAAAAiTGpApa7Xzj2VRVzgfZdSMMlXTbWTe6+2cw+L+lfS06/28yOdfdnKtxjRZnZRyTdLMlU+PFe4u53RNvV8BjBAgAAQBgmVcCqsvcEjh9w9+fKvPd2Sf8oaWqg3mgB66eStpff3pg6x3OxmV0o6SYVRtz6VRi5uqWC/VTUrObiCBYBCwAAABVEwAqBmR0g6S2B098t935332lm92rfVfiWSbp2lHtuVyGYVZ2Z/X8qLGphkvok/WXxmbKala2v0/TGem3a0RV1KwAAAEiQybzIRZjeqqHh9dfjrBG8/kQzmz7xlsJhZp+UdIMK4apX0odqPVyt39KpK+9+Sjt29+jfH16vhVfdpyvvfkrrt4xr0A4AAAAYgoAVjuCqfzs0+vS+4Tw8zLljJtZOOMzsf0r6p+Jht6T3u/udEbY0pgeea9eZX12hOx55Uf3FdQ13dvXqjkde1JlfXaEHnmuPtkEAAADEGgErHMEg9PwElilfW0bdyJjZ1SqsbihJeySd4+53R9jSmNZv6dSlt63U7p4+9fbv+5+jt9+1u6dPl962kpEsAAAATBgBKxwLAscvjLeAu2/V0IUmgnUjYWZflvSF4uEuSWe7+48jbKksN61Yq56+/lGv6enr180r1lWpIwAAACQNi1yEY0bg+JUJ1nlF0pGj1K06M/tfkkr38vpXSc1mNuo+XZJ2uftPw+tsbHev2jBk5Cqot99116qXdc2yhVXqCgAAAElCwArH1MDxrgnW2T1G3SgEV0f8ZPFrLOslzRvujba2tsHX+Xxe+Xx+Yp2NobOrt7zrusu7DgAAAAgiYIWjKXA80c2WajFgVVxpwApTUyatnWWErKYGflsAAABgYngGKxxTAsfdE6wT3KQpWLfq3D3v7jaBr3lR975s0WylUzbqNemU6ZxFc6rUEQAAAJKGgBWO4IhVwwTrZMaoi3G4eMnhqq8b/Zd8us500ZL5VeoIAAAASUPACsfOwHF2gnWCI1bBuhiHw2Y26YbzF2tKfd2Qkay64uHxc6bp0BmNEXQHAACAJCBghSMYhCb6N/ZgwOqYYB0ULV3QouWXL9F5pxyqqZm0zKSpmbQ+9MbDdPFb5uvR57fpO4+8GHWbAAAAiCme5g/HtsDxQROsE7wvWBcTcNjMJl2zbOGQpdj7+11rXunQ1f9ntRYfNl1HH5SLqEMAAADEFQErHM9JOrHk+LDxFjCzGRq6auBz+9NUrarWMu1jSaVM//SBN+jPrl+hy25fpXsuO02NrCgIAACAcTD30TdenSzMbJ6kdYHTS939FxOodaWka0pO7ZA03cfxk21mJ0p6LHD6NHd/aLz91DIzG89PS1X8+g+bdf63fqP3nzhXf/++E6JuBwAAABVkZnL30ZeW3g88gxWO1YHjnKRjxlnjTcOce3Zi7WA8TjvyQH00f6TufOwl3b3q5ajbAQAAQIwQsMLxS0l9gXNvHmeN4PWr3J1nsKrk8re/TifPO0BX3PWU1m3ujLodAAAAxAQBKwTuvlXSisDpc8u938ymSjo7cPqu/e0L5UvXpfS1Dy5SfTqly25fqa7eYF4GAAAAhiJghec/A8dLzWxBmfd+SEMXuAjWQ8hmT5+ir7zvBK3esEPX/XhN1O0AAAAgBghY4blFUnvJsUn6xlg3mdlM7btAhiTd6+7B57pQBacf26oLTpunWx96Xvet3hh1OwAAAKhxrEEdEnffaWbXSrq+5PTbzexL7v654e4xs2ZJd0pqKS0l6YrwOo1erSzTPpLPnXW0Hn1+qz77/Se1cM40zZke3P8ZAAAAKJhUy7QXp+gtHeHtmZKuDZz7Z0m/G+H6x9w9uIx68PPqJT0k6aTAW8slXSfpYXfvNrOcCs9cXamhqw1+zd0vH+1z4qwWl2kfzvObO3X21x/UgoOadcclb1J9HYO/AAAAcRT2Mu2TLWD9pQpT9yrhandvK+MzD5H0oKRDh3m7X1KHpGkj3L5c0rvcvWeiTda6uAQsSbrniQ36+HdW6dL8EfrsmUdH3Q4AAAAmgH2wYs7dX5R0mqTfDPN2SiOHq29JWpbkcBU37zphtj548iH65i//qBW/fzXqdgAAAFCDCFhV4O4vqRCyLpE02nJ0Lunnkt7m7he5e1c1+kP5rnrncXpdy1R98ruPq71jT9TtAAAAoMZMqimCtcLMjpJ0oqS5krKSdkpaK+m/3L19tHuTJk5TBAf8blOH3vWNB3XSYTP0vy88RalUaCPMAAAAqDCewUKixTFgSdIdj7ygz/3nU/rMGQv00aVHRt0OAAAAyhR2wGKZdkSu1pdpH865Jx+iX/9xi/7xp8/plPkzdPK8GVG3BAAAgBrACBYiFdcRLEnq2NOjs7/+oLp7+/Xjjy/RAU0NUbcEAACAMbCKIFCjmrP1+vp5i7R5Z5c+8/0nFdegCAAAgMphBAuRivMI1oBvPbhO19z7jC5beqS27+7W3as2qLOrV02ZtJYtmq2Llxyuw2Y2Rd0mAAAAxCIXSLgkBCx31zk3/FqPv/ia6lJSX//e99IpU31dSjecv1hLF7RE1yQAAAAkMUUQqHkvbN2lNRs7JO0briSpt9+1u6dPl962Uuu3dEbQHQAAAKqJgAXsp5tWrFVv3+ijcD19/bp5xboqdQQAAICoELCA/XT3qg3q7R89YPX2u+5a9XKVOgIAAEBUCFjAfurs6i3vuu7yrgMAAEB8sdEwIhfHjYZLNWXS2llGyGpq4LcbAABA0rGKICKVhFUEr7z7Kd3xyIujThNMp0znnXKorlm2sIqdAQAAIIhVBIEad/GSw1VfN/pvpT53vfsNs6vUEQAAAKJCwAL202Ezm3TD+Ys1pb5O6dS+/xiSTpka0ik11KV02e2r9MyGHRF1CQAAgGpgiiAilYQpggPWb+nUzSvW6a5VL6uzu1dNDWmds2iOLloyX51dfbrw1kfVsadH3/gwmw4DAABEJewpggQsRCpJAWssG1/bo7/69qNas7FDbe86Tn/+psOibgkAAGDSIWAh0SZTwJIKS7p/7Dur9PM17bp4yXz9j7OOUSoV2u9vAAAABBCwkGiTLWBJUm9fv6659xl9+7/W64zjWvXVcxdpSkNd1G0BAABMCgQsJNpkDFgD/u3BdbrmR8/o9XOm6aaPnKSW5mzULQEAACQeAQuJNpkDliT9dPVGfeKOxzWjqUG3XHCyjmptjrolAACARCNgIdEme8CSpKdeek0XfvtR7enu0zfPP1Fvft2BUbcEAACQWAQsJJqZ+VVXXTV4nM/nlc/no2soIi9v360Lb3lUf3x1p/7unOP1gZMPibolAACARCJgIdEYwdqrY0+PLv2PlVrx+826NH+EPv2nC1hhEAAAoMIIWEg0Ata+evr69YUfrtZ3HnlBZ7/+YP3D+0/Qph17dNOKtbp71QZ1dvWqKZPWskWzdfGSw3XYzKaoWwYAAIgVAhYSjYA1lLvrxl+t1XU/WaMjW6bqpW271Nvn6u3f+/OUTpnq61K64fzFWrqgJcJuAQAA4oWAhUQjYI3s2w89r6vuWT3qNVPq67T88iWMZAEAAJQp7ICVCqswgP3z+/YO1Y3xW7+nr183r1hXnYYAAAAwJgIWUKPuXrVBfWMM7vX2u+5a9XJ1GgIAAMCYCFhAjers6i3vuu7yrgMAAED4CFhAjWrKpMu7rqG86wAAABA+AhZQo5Ytmq10GftgLT50ulgoBAAAoDYQsIAadfGSw1VfN/pvUTPpV7/frAtufVQvbt1Vpc4AAAAwEpZpR6RYpn10DzzXrktvW6mevv5h98H6xnmL9PzWXfqnnz6nPnd9/G2vKyuYAQAATFbsg4VEI2CNbf2WTt28Yp3uWvWyOrt71dSQ1jmL5uiiJfMH97965bXdartnte5bvUlHtU7V351zvE6aNyPizgEAAGoPAQuJZmZ+1VVXDR7n83nl8/noGoq5nz2zSVfds1ovb9+t8045RH975tGa3tgQdVsAAAA1g4CFRGMEq/I6u3r1tft/r289uE7Tp9TryrOP0bI3zJFZaH+OAAAAxAYBC4lGwArPMxt26H/e9ZQef3G7Tjtypq5590IdPmuqpMK0w5tWrNXdqzaos6tXTZm0li2arYuXHD447RAAACCJCFhINAJWuPr7Xbc/8oK+vHyNunr6denSI3Ts7Jw+8Z3HR1w444bzF2vpgpYIuwYAAAgPAQuJRsCqjvaOPbr23md1zxMbZJJG+xmfUl+n5ZcvYSQLAAAkUtgBi7WcgUmgpTmr689bpP92dMuo4UqSevr6dfOKdVXpCwAAIGkIWMAk8si6rWNe09vvumvVy1XoBgAAIHkIWMAk0tnVW9513eVdBwAAgH0RsIBJpCmTLuu6lJn+/eH1at+xJ+SOAAAAkoVFLhApFrmorivvfkp3PPLiPqsHBpmk5kxaO7p6ZSYtOmS6zlx4kM447qAxF75g+XcAAFDrWEUQiUbAqq71Wzp15ldXaHdP34jXTKmv008+8Wb19LnuW71Ry1dv1NMv75AkHX1Qs85ceJDOXHiQFrQ277N58QPPtevS21ay/DsAAKhpBCwkGgGr+iYShF7cuks/fWaT7nt6ox5dv1Xu0mEzG3XmcQfpjIUH6YAp9fqz6x8cM7ix/DsAAIgaAQuJRsCKxvotnbp5xTrdtepldXb3qqkhrXMWzdFFS+aPGYBe7ejSz57dpOVPb9RDf9ysnj5Xtj6lrp7+UZeAT6dM551yqK5ZtrCyPxgAAIBxIGAh0QhY8bZjT48eWNOuT935hPpGea5rwNRMWk9ffUYVOgMAABgeGw0DqFm5bL3e/YY56i8jXEks/w4AAJKPgAVgv41n+fc7HnlBm3d2hdwRAABANJgiiEgxRTAZyl3+vSlTp51dfUqZdNK8GYOLZMyZPmXU+iz/DgAAKoVnsJBoZuZXXXXV4HE+n1c+n4+uIUzIeJZ/39Xdr/tWb9R9qzdqzcYOSdLxc6YV99pq1ZEtzfvcx/LvAACgkghYSDRGsJJjIkHo+c2dg3ttrXphuyTpiFlNOuO4wl5buWxaZ32N5d8BAEDlELCQaASsZNmf5d83vrZH//eZQth6eO1W9fW7GhvqtLu7j+XfAQBAxRCwkGgELAxnW2e37l/Trr/9wZMs/w4AACoq7IBV3tJfAFBFBzQ16H0nztVnvvdEWdd3dvWqq7dPmXRd2Z/BwhkAACAMjGAhUoxgYTQLr7pPO7vK2zuroS6lhXNyWnzoAVp82AE68bAD1JrLDnstC2cAADB5MUUQiUbAwmjKWf49nTKdduSBWnBQs1au36YnX35N3b39kqQ506do0aHTtfjQQuA6dnZOG7bvLmvFQxbOAAAgmQhYSDQCFkZT7vLvpWGou7dfqze8ppUvbNfKF7Zp5fpteuW1PZKkTDql3JR6be7oCm3hDKYeAgBQ2whYSDQCFsZSiel8r7y2WyvXb9dv12/TrQ+tUxnrZkxo4QymHgIAUPsIWEg0AhbKsT/LvwfN/9yPRh29KvXWo2bpiFlTdfisJh0xa6qOmNWkWc0ZmQ39M3kio20AAKD6CFhINAIWqq3chTPSKdNRrc1at7lzn9DUnEnr8FlNOnzWVB1+YJOOaCkEsG8/9Ly+99hLYz4vxtRDAACiRcBCohGwUG3lLpwxEIT6+12v7Nijta/u1B/bd2rt5k6tfbVTf3x15+CzXePB1EMAAKJFwEKiEbBQbZWcytfZ1at1mwth6xN3PF52D5e85XDNmT6l8HVA4SuXrQ+9XwAAQMBCwhGwEIUwRoTKnXpoJtWnUuru69/nfHMmXQhbA6Fr+hTNnj5F9z65Qfc/2x7K1EOmHQIAJiMCFhKNgIWoVHLhDGl8Uw+vftdx2tzZpZe37dbL23drw/bdg69f2lY43rGnvA2WB4x36mGY0w4JbgCAWkbAQqIRsJAUlZ7K17GnRy8XN0UuV3M2rdnTpmj29KwOnj5Fs6dldfC0KTp4elZzpk/RQdOyyqTrQp12SHADANQ6AhYSjYCFJIly6mFDOqUPnnyINmzfo1de261XXtujrZ3dQ647cGqD+l3a1tld8c2W4xjcCG0AMPkQsJBoBCwkTZRTD4NhaHd332DY2rB97/fv/fYl9ZWx27KZ9Po50zStsUEHNNZr+pR6TW9s0PTGeh1Q/D598L0G/f19a/TdRyfW62jCCm5xHG0Lo26ceg2rbpx6DatunHqNW13+Iaf2ELCQaAQsYHRhhIvxbLb8lqNm6bVd3dq2q0fbd3WP+9mw4WTSKX3l/SeoqaFOTZm0mhrSasrUaWomrcZMWo31dUql9v5/b39C5kjiONoWRt049RpW3Tj1GlbdOPUat7r8Q05t9krAQqIRsICxVfp/0OVOOxxu4Yzevn7t2NOrbbu6tb0Yurbv6tG2Xd269kfPlv+DGkPjYPiq0wtbd6mMATdl61O6+S9OVuNAWGsY+J5WQzq1z7VhhDYpvOAWRt049RpW3Tj1GlbdOPUat7r8Q07t9hp2wEqHVRgAUBlLF7Ro+eVLKjb1cNmi2WWFi3MWzRl6vi6lGU0NmtHUMOS9r/7s92UFt6aGOv3wstPU2dWnzq5edXYPfO8tfO8qPe7T81t2lfXj2tPTr/O/9Zth32uoS6kxUzc4WvaH9p1jhrbeftedj72oNxwyXQ3p1OBXpq74PV23z/mGupT+5YE/qCewBH9QT1+/bl6xblzB7aYVayteN4yacasbp17DqhunXuNWN6xe12/p1KW3rRw2uPX2u3r7+3TpbSsnFDIrXTdOvVYSI1iIFCNYQPWF9a+qYY0KlTvi1thQp1svOGUwqO3q6tPOrl7t6u7Vzq6+4vfC+eWrN5b9+WEwSTOaGmS294xZ4byk4msbfL3xtT1lTes0k+bNbCrcaXs/y8yG1P59e0dZI4MpkxbOmTZYa7DIYG3t897jL24vq26dmU6ef0DxXttbZ7Ce7XP80B+2qK+M/1/UpUz5o2bt83M7UKf050CSfvZse1nPI6ZTprOOP3jI/XuP9/58yKR7Ht8w6u+DAfUp03tPnDvk18G+tQff0R2PvqCevjLq1pn+4tR5w9cp+bUgk/7twXVl1WyoS+mv33r4Pr8Ghutx4PgbP//DkP3+hq2bTulTpx81fK8lPx8D/v6+59TdO3bdTDqlK99xzD6/Vkfq1SRddc9qdZVRN5tO6YvnHL/Pr9fgr9UBn/3+k+XVrE/pq+e+YbDL0X59mUn//vB6rfjdqxrtP1udSfkFLfrIn8wb9ucz8NtDJtO/PbhWP1/TPmbd049t1cVvOWLEXkt/rr/5iz/qp89s1Gi/FOpS0lkLD9ZHlx457M9n8L+dZPra/b/Tj596ZdS66ZTpna+frU/96VH7nDeTDpnRxBRBJBcBC4hGGFM24hTcyg1tTQ11+vEnlqi7t19dvf3q7utXV0/he3dv8auvb/D153+4uuwf14ffeKgkySUV/hgs/PjcB44ll8td+t5vXyq77jtPmK2BP1cHyvowtccTMpcumDUY8Pb2JgX//HaXHvzD5rLrnjJ/xpD+Bvsu1h94veqF7WXXPW52rqxen9vUUXbN+Qc2Dfa0b48Dn+GDr1/atrvsui3NmZIeB7sb0r+kYVcGHcnUTHrYXoM/1+UEgAEDf8Hlf9uIu/VfPpspggCAyqr0tENJOmxmk244f/GYwW28tS9ecrh+8NuX1ds/cnCrr0vpoiXzy65Z7jTJ9yyeO65+v7z8ubKfb/viOceXXfcnT28su+7Xz1tUVs3xPIt3ywWnlFVzvHXv/OtTQ6n7o48vqXjNBz6dL6vmeOs+csXbQ6lb7sbj+1tzMMQFwuDrr75PnV0j/54trfvoFW8fNmQPFxD/5Lr71dk9dt2mTJ1+8emlxfv3bW64wPm2f/yldpVTt/iPLnvvH6g1tNd3fePBsmo2NtTp+//9T/bpNRhiS3t99zd+Xd6ItqTv/fdTB/sa/Lkd4R8Jzrvp4TKqFureeuEpQ+oM948lf/Xtx8quecOHFw/7Dxf7Hhf6/8Qdj5dd9yvvP2HIPzqd++Wybp8wAhYi19bWNvg6n88rn89H1gswmRw2s0nXLFs4rul6Y4lLcAsjtEn793xbtevGqdew6sap17Dq7m/NwelrgbGAcxbNKbvulIa6svs9Z3F5dd+zaK5mFUcHy/GecuuO4x9dyq353sVzdWxx1LUcTZl0eSPwmbROmjej7LpTx1H3rUfNqnjNs44/uKyaknTFXU+XXfd9J84dcv7csj9pYlJjXwKEq62tbfCLcAXE30Bwe/rqM7Tuunfo6avP0DXLFu7Xg8YDwe28Uw7V1ExaZoX/cZ93yqFafvmScU9nHAhtU+rrlE7t+zfDdMo0pb5uwqNt9XWj/691IsEtjLpx6jWsunHqNay6ceo1bnXD6nXZotlD/twKmmiAr3TdOPVaSQQsAEAsVDq4VTq0DfQYRnALo26ceg2rbpx6DatunHqNW13+ISdevVYSi1wgUixyASCJ1m/prOg0yTDrxqnXsOrGqdew6sap17jVDaNmLe4tFade2WgYiUbAAgAAGL/JHjL3py4BC4lGwAIAAEA1hR2weAYLAAAAACqEgAUAAAAAFULAAgAAAIAKIWABAAAAQIUQsAAAAACgQghYAAAAAFAhBCwAAAAAqBACFgAAAABUCAELAAAAACqEgAUAAAAAFULAAgAAAIAKIWABAAAAQIUQsAAAAACgQghYAAAAAFAh6agbQPyY2UmSTpK0WNLxkmZJOlBSk6QOSeslPSrpu+5+f1R9AgAAANVm7h51D4gZM9suaVqZl/9K0ofd/aURajm/BgEAAFAtZiZ3t9Dq85dbjFcxYPWpMEr1hKS1krZLckmtkv5E0jmSMsVbXpC0yN23DlOLgAUAAICqIWCh5pjZ8ZJWu3v/KNccKemnkuYXT13v7p8Y5joCFgAAAKqGgIXYMrP3Svp+8XC9u88b5hoCFgAAAKom7IDFIhcRKI7unChprqQpKiwMsVbSw+7+apS9VdjTJa8PjqwLAAAAoEombcAys7kqrIR3cvH7SZJmBC5b6u6/qNDnpSRdIOmTko4b4bJ+M7tf0nXu/kAlPjdiryt5/UpkXQAAAABVMqkClpl9QXtD1UFV/Nw5kr4n6dQxLk1JOl3S6WZ2o6SPuXt32P2Fofhj/krJqTuj6gUAAAColkn1DJaZjfcHu98jWMWg8aCkecO87SpMD8yNcPuPJJ3j7j3700OYintizS0eplXYD+sUSR9QYV8sSXpI0hnuvnOY+3kGCwAAAFXDM1gxZmb1kn6goeFquaQvqfDMVZeZ5SS9S9IVko4uue4dxev+JvxuJ+zTks4d4b0XJH1LhSmPNRsSAQAAgEpJRd1AhLolrZR0o6SLJb0zhM+4RNIbA+e+4u5nufsv3b1Lktx9h7vfpsLIT/DZq8vNbKRntmpZn6T7JT1IuAIAAMBkMdmmCP6bpMdU3CC39PkmM5snaV3glglPETSzRkl/1L7Pev3M3U8f474DJa2W1FJy+ofuvmyM+26V9JGJ9DqCg91943huMLOMpNmS8ios5nF88a1vSfroQKAM3MMUQQAAAFRN2FMEJ9UIlrtf6O43uPujVVg84gLtG65c0mVj3eTumyV9PnD63WZ2bAV7C4W7d7n7One/RdJi7d0D669UGCkEAAAAEu3/tXfvwZJV1QHGvwUDMzC8MQjIU14iKZUCEUVlRtFEQUG0fIJiioCPELTKSBKVghICMVERIRoVFMsQQAFRUFRUVBweToLGAAI6KppCgeIxM8DwmpU/To/07Nv39uOevqf73u9XNTWzd5+9z6qpu+r26t5nb5/BGp7Di/b3M/PWHseeD3wU2KiY7+YpxnwbuL/38Lp6cDqDM/PxiDgaeDnVJh5HRsQ/9fF/IEmSJI0dC6whiIjNgRcX3Rf2Oj4zV0bE5cAb27oPA06ZYsz5VIXZyMjMByLiGuCVQACLAQssSZIkzVpzaongDDqQicXrj/uco7x+n4jYbPCQGrOi7d+bNxaFJEmSNAMssIaj3PVvOVMv7+vkug59ew4WTqN2a/v33Y1FIUmSJM0AC6zhKAuh3wywVd6yHuYdaRHxAqrNLtb4YVOxSJIkSTPBAms49ijad/Q7QWbey8SNJsp5Z1xEvDMiDoqIKbe2jIgXA5e0dV2ZmbcNNzpJkiSpWW5yMRxbFO07B5znTmDXKeZtwvOAfwPujIjvAD8H/gA8DCykivcg4PltY5ZRHeYsSZIkzWoWWMOxUdF+aMB5Hu4yb5O2Ad7aw3VfBo7LzD8OOR5JkiSpcRZYw7GwaK8acJ5RLLCOo1r69yJgb2AXYCtgAVUheR9wC3AtcGFm3tJtwpNOOulP/+iJ4Y4AABKmSURBVF60aBGLFi2qO2ZJkiRpRlhgDccGRfvRAed5pMu8My4zVwBfa/2pRXuBJUmSJI0zN7kYjvIbq/UHnGd+l3klSZIkjRALrOFYWbQXDDhP+Y1VOa8kSZKkEWKBNRxlIbThgPOUBdaKAeeRJEmSNAMssIbjvqK99YDzlOPKeSVJkiSNEDe5GI5bgX3a2jv2O0FEbMHEXQNvnU5Qo8pdBCVJkjRbWGANR7k1+U4REZmZfcyxcw/zzgruIihJkqTZwiWCw3FT0d4E2LPPOfbv0DcrCyxJkiRptrDAGo4fAE8UfS/sc47y+hsz02ewJEmSpBFmgTUEmXkv8KOi+w29jo+IjYBDiu5LpxuXJEmSpOGywBqeS4r24ojYo8exb2biBhflfJIkSZJGjAXW8HweuKutHcBZ3QZFxJbAh4vuyzOzfK5LkiRJ0ohxF8EhycyVEXEKcGZb90ERcXpm/n2nMRGxMXARsFX7VMAHhhdp89ymXZIkSbNF9Ldz+HhrLdFbPMnLWwKnFH0fB26b5Pqlmbm0y/3WA5YA+xYvXQmcBlyXmY9GxCZUz1x9kIm7DX4iM98z1X3GWf+710uSJEmDiwgyM4Y2/1x6cxsRR1Et3avDyZl5Ug/33B64Btihw8urgRXAppMMvxJ4dWY+NmiQo84CS5IkSTNp2AWWz2ANWWb+DjgAuL7Dy+sweXF1DnDYbC6uJEmSpNnGAmsGZObvqYqsY4BfTHUp8D3gpZl5dGY+MhPxSZIkSarHnFoiOCoiYndgH2A7YAGwElgGXJuZd001drZxiaAkSZJmks9gaVazwJIkSdJM8hksSZIkSRoTnoOlxnkOliRJkmYLlwiqUS4RlCRJ0kxyiaAkSZIkjQkLLEmSJEmqiQWWJEmSJNXEAkuSJEmSamKBJUmSJEk1scCSJEmSpJpYYEmSJElSTTxoWI3zoGFJkiTNFh40rEZ50LAkSZJmkgcNS5IkSdKYsMCSJEmSpJpYYEmSJElSTSywJEmSJKkmFlgaG1dffXXTIfTMWIdnnOI11uEZp3iNdXjGKV5jHZ5xinecYoXxineUYrXA0tgYpcTpxliHZ5ziNdbhGad4jXV4xileYx2ecYp3nGKF8Yp3lGK1wJIkSZKkmlhgSZIkSVJNPGhYjYoIfwAlSZI0o4Z50LAFliRJkiTVxCWCkiRJklQTCyxJkiRJqsm8pgOQJhMRuwL7ANsBGwArgGXAdZl5d5OxSXNZRKwH7Ac8E3hKq/se4Gbghsx8rKnYpLq0fs53p/o5fyqwCbAKuBf4LdXP+oNDuq/5pVkvIjalyrEdgG2AhVS1yXLgPqqf+Z/X/TM/EznmM1gaKRGxDvB24L3AXpNcthr4LnBaZn5/pmKTRklEbAfsCzy39fe+wBbFZYsz8+oa77k18I/A26jebHbyAPBF4NTM/GNd95ZmQkTsDbwKeAmwPzB/isufAJYAZwMXZ+bj07y3+aVZLSL2pMqvA4C9ge17GPYIcDnw2cz81jTvP2M5ZoGlkRERTwO+DDy/j2GfAY7LzEeHE5U0OiLiRJ4sqrbuYUhtBVZEHAZ8Adi0xyH3A0dl5mV13F8apog4HPhnYNcBp/gZ8NbM/J8B729+adaLiLOAd09jim8Cx2bm7wa494zmmM9gaSS0iqtr6FxcJdXXxZ0cA1zS+rpXmu1Opvr0r5fiqjYRcQRwCZ1/Ma2i+oSxtBlVbr5lmLFJNdmP7sXV/cDDk7z2bOCGiDio3xubXxIAj1MtC3yAaqVSJ68AlkTEbv1M3ESOWWCpca3i6GJgp+KlK4FFwAaZuSlVYhwJ/KK47mDg9OFGKc1NEbEvcA7Qfl7Iw8BpwNOBDamekdyFKg9XtV23DnBuaw5pnDwInA8cTfWMyPzM3DwzNwS2Bf4KuKUYMx/4akQ8q9ebmF+ao+6lKnjeCxwIbJOZ62XmFpm5GbAAeA7wIeDOYux2wBURsaCXGzWVYy4RVOMi4t3AWUX3v2Tm+ye5fmPgMmBxW/dq4FmZedNwopSaVxzM/Sjwv8BS4CfAH4CvF0OmvUQwIm6gWpK4xv3AwZm5ZJLrDwCuYO1PCq/PzP2nE4c0TBFxOnACcBvwceD8zJxs5cSaMetTPX91dPHSNZn5oh7va35pzoiIV1IVNz/MzCd6HLMJ1YcdBxcvfTAzT+1hfCM5ZoGlRkXEhsCvWHvJ01WZ+bIu454C3ARs1dZ9WWYeVn+U0miIiHN5sqD6WfuzhxGxE/DrYsi0CqzWcykXF91HZuaXuox7K3Be0f2azPzqoLFIwxQR76JaJvSFXt/4tcYF8DXgkOKlF2XmNV3Gml9SDyJiPtUHiu3LeG/LzD26jGssxyyw1KgO314lsGdm3trD2GOAfy+698rMm2sMURoLQyqwfkK1qcYaSzLzgB7HLmHtZypvyMznDRqLNKoiYhfgdtZegvSvmfl3XcaZX1KPIuJvgE8W3dtmZrmEsH1MYznmM1hq2uFF+/u9FFct5wMru8wnaQARsQNr/2IC+HQfU5QffuzX2lpemlUy81dUuwi223uqMeaX1LdrO/RtO9nFTeeYBZYaExGbAy8uui/sdXxmrqQ6G6GdSwSlehxatB+neii5Vxe3xrQzPzVb/bJod9vp0/yS+rOiQ9+6U1zfaI5ZYKlJB1Kd2N3ux33OUV6/T0RsNnhIklpeWrR/mpkP9jq49QFI+al+31tYS2Ni/aLd7fkL80vqz9M69N0xxfWN5pgFlpq0V9FeDvT7/NR1Hfr2HCwcSW3K/OyUa92UY8xNzVa7FO1JnwtpMb+k/pTfSN2cmX+Y4vpGc8wCS00qf1B/k/3vurKsh3kl9aG1Y9PORXenXOumHLNLa2tradaIiN2Z+GZu6RTXm19SH1rnUL2z6D5jiusbzzELLDWp3F5zqq96O8rMe6kOhJxqXkn92ZWJa9v7zs8OY9Zl4if90rgr3/hBtXX7ZMwvqQcRMS8i3g5cxdrLcL8FfG6KoY3nWPn8izSTtija3ZZUTOZO1j4boZxXUn865dAg+dlpjPmpWSMi9gTeVXT/NDOnWo5kfkktEXEEsFFb1zxgM6rVSAex9nmnAJcCR3RZ8dR4jllgqUkbFe2HBpzn4S7zSupPpxwaJD/L3JxsbmnstJYKfZGJG1y8v8tQ80t60ul03sCidB3wkcy8tIdrG88xlwiqSQuL9qoB57HAkupV5iYMlp++AdRsdgYTz9k5LzO/02Wc+SX152bgIuAHPV7feI5ZYKlJGxTtRwec55Eu80rqT6ccGiQ/y9ycbG5prETEu5n47NXtwN/2MNz8kvrzTOBjwB0RcUJEdKtfGs8xCyw1qfw0YdDdj+Z3mVdSfzrl0CD5WebmZHNLYyMiDgU+UXTfD7wmM5f3MIX5JbVk5naZGWv+UH37tD1wMHAmVW6tsZBqSeGFETHVY06N55gFlpq0smgvGHCe8tOEcl5J/emUQ4PkZ6dP+sxPja2IWARcwNo7lD0IHJyZN/U4jfklTSIzH8rM32fmNzLzeKpd+75cXPY64NQppmk8xyyw1KTyh3TDAecpE2DFgPNIqnT6BTJIfnb65WR+aixFxHOptl9vf6P2CHBYZi7pYyrzS+pR6zieNzKxyHpfROw2ybDGc8wCS026r2hvPeA85bhyXkn96ZRDg+TnNj3OLY20iPhz4Epg47bux4E3ZOZVfU5nfkl9yMzVVMchtBc36wDHTjKk8RyzwFKTbi3aO/Y7QURswcQdXcp5JfXnl8Dqoq/v/AR2KNpPAL8aKCKpIRGxO9VBp+3n36wGjszMywaY0vyS+pSZ9wBXFN0vmeTyxnPMAktNuqVo7xQR0eccO/cwr6Q+ZOYq4NdF99MHmKrMz2WZ2WlXJmkkRcSOVMXVU9u6EzgmMy8YZE7zSxrYfxftTu8BRyLHLLDUpPKB4E2oTu7ux/4d+iywpOkr8/N5A8xR5qe5qbEREdsC36Pa0azdezLznGlOb35J/XugaE/1XFWjOWaBpSb9gOrr1nYv7HOO8vobM9M16NL0fa9o7x0RnQ5v7Kh17XOK7u9OOyppBkTEn1F9c1V+6v0PmXlmDbcwv6T+bVm0757i2kZzzAJLjWntDPOjovsNvY6PiI2AQ4ruS6cblyQAvlq05wGH9zH+da0xU80pjZyI2Az4NhNXVHw4M0+v6Tbml9S/5xbt/5vi2kZzzAJLTbukaC+OiD16HPtmJm5wUc4naQCZ+Vsmrnd/Rx9TlLs7Lc3MO6YXlTRcrQ/uvsnET64/mpkn1nUf80vqT0RsA/xF0T3pDp5N55gFlpr2eeCutnYAZ3UbFBFbAh8uui/v46BHSd2Vn9a/ICKO6DYoIo4Ent9lLmmkRMQCqnOuyucuzs7M9w3hluaX5oyImD+NsesCn2biM1df6TK0sRyzwFKjMnMlcErRfVBETPqDHBEbAxcBW7VPBXyg/gilOe0rwH8VfZ+MiPIXz59ExAHAJ4vupZl5cd3BSXWJiHlUB5kuLl46BzhuSLc1vzSXnBoRn5/icOCOImJrqtVJry5euiAzb+wyvLEci8zs53qpdhGxHrAE2Ld46UrgNOC6zHw0Ijaheubqg0xcG/+JzHzP0IOVGtRaPlu+AVxjSyZ+WPFx4LZJrl+amUt7uOd+wDXAem3dDwMfAz4LrFkysSNwDPBeYEHbtY8BL8zMG7rdS2pKRJxNdZBpuzupcqo8T6dX/5GZK6a6wPzSXBERZwDHt5pLqZ5nWgr8DPhjthUkraJqb+Aw4E2sfcA3VHmxf2be2cN9G8kxCyyNhIjYnioBykPdoPrltgLYdJLhVwKvzszHhhSeNBIi4iiqZbV1ODkzT+rxvm9r3bfTOXWrWn8v6PBaAkdl5hcHilCaIRFxNXBgzdPunJm/6eHe5pdmvaLAKiWwnGpn6U2YuLlEu98CizOzPOdqqnvPeI65RFAjITN/BxwAXN/h5XWYvLg6BzjM4koansw8D3g91S/A0gI6/2JaDrzON3/S1MwviaB6n7cFkxdXq4FPAc/up7iCZnLMAksjIzN/T1VkHQP8YqpLqc43eGlmHu3J9dLwZeZXgGcAZ1N9ozyZ5VQb1eyRme7qKfXA/NIccCLwWqplebdQvZfrxe1US3V3z8x3ZWZ52HBPZjrHXCKokRURuwP7ANtRfbqwElgGXJuZd001VtLwRMT6wH7AXjx58OM9wM3A9X6jLA3O/NJc0Hqu/hnATsBTgYXAulQFznKqZ6NuzMz7h3DvoeeYBZYkSZIk1cQlgpIkSZJUEwssSZIkSaqJBZYkSZIk1cQCS5IkSZJqYoElSZIkSTWxwJIkSZKkmlhgSZIkSVJNLLAkSZIkqSYWWJIkSZJUEwssSZIkSaqJBZYkSZIk1cQCS5IkSZJqYoElSdIIioijIiLb/ixqOiZJUncWWJIkSZJUEwssSZIkSaqJBZYkSZIk1cQCS5IkSZJqYoElSZIkSTWxwJIkSZKkmlhgSZIkSVJNLLAkSZIkqSYWWJIkzRIRsW5EnF0cUHx7ROzadGySNFfMazoASZI0fRGxAfCfwKFt3TcAh2Tm3c1EJUlzjwWWJEljLiKeAnwd2L+t+wrg9Zn5UDNRSdLc5BJBSZLGWEQ8HVjC2sXV54BDLa4kaeZZYEmSNKYiYl/gWmC3tu6TMvOvM/OJhsKSpDnNJYKSJI2hiHglcBGwsNX1OHBsZp7bXFSSJL/BkiRpzETE0cBlPFlcPUi1JNDiSpIaZoElSdIYiYiTgc/y5CqUu4BFmfmN5qKSJK3hEkFJksbDvIg4F3h7W9/twF9m5rKGYpIkFSywJEkaD58Bdm5rX091xtU9DcUjSerAJYKSJI2HnYv2hyyuJGn0WGBJkjQevlm0L42IRU0EIkmanAWWJEnj4SPACW3thcA3IuLlDcUjSerAAkuSpDGRmR8Bjm/r2gD4WkQc0lBIkqSCBZYkSWMkM88EjgWy1TUfuCQiXttcVJKkNSywJEkaM5n5Gart2le3utYDLoiINzUXlSQJLLAkSRpLmXke8Bbg8VbXPOBLEXFUY0FJkiywJEkaV5l5AfB64NFW1zrAuRHxjuaikqS5zQJLkqQxlpmXAocDq1pdAXwqIo6ffJQkaVgssCRJGnOZeQXwKuChtu4zIuKESYZIkobEAkuSpFkgM68CXgGsbOs+PSJObCgkSZqTLLAkSZolMvOHwMuAB9q6T46If2ooJEmacyIzu18lSZIkSerKb7AkSZIkqSYWWJIkSZJUEwssSZIkSaqJBZYkSZIk1cQCS5IkSZJqYoElSZIkSTWxwJIkSZKkmlhgSZIkSVJNLLAkSZIkqSYWWJIkSZJUEwssSZIkSaqJBZYkSZIk1cQCS5IkSZJqYoElSZIkSTX5f5/5HURri1A4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot the principal components\n",
    "colors = 200*['green', 'purple', 'orange', 'teal', 'black', 'red']\n",
    "\n",
    "#plt.plot(incpca.explained_variance_ratio_, marker='o', linewidth=1.0, markersize=12, color = 'red', markerfacecolor='None') \n",
    "plt.plot(incpca.explained_variance_ratio_, marker='o', markersize=10) \n",
    "\n",
    "\n",
    "#backround grid details\n",
    "axes = plt.gca()\n",
    "#axes.grid(b = True, which = 'both', axis = 'both', color = 'gray', linestyle = '-', alpha = 0.5, linewidth = 0.5) \n",
    "#axes.set_axis_bgcolor('white')  \n",
    "\n",
    "#font scpecifications\n",
    "title_font = {'family' : 'arial', 'color'  : 'black', 'weight' : 'heavy','size': 20}\n",
    "axis_label_font = {'family' : 'arial', 'color'  : 'black', 'weight' : 'normal','size': 20}                                                   \n",
    "\n",
    "#figure size and tick style\n",
    "plt.rcParams[\"figure.figsize\"] = [12,8]\n",
    "plt.rc('axes',edgecolor='black',linewidth=1)\n",
    "plt.tick_params(which='both', axis='both', color='black', length=4, width=0.5)\n",
    "plt.rcParams['xtick.direction'] = 'in'\n",
    "plt.rcParams['ytick.direction'] = 'in'\n",
    "plt.minorticks_on()\n",
    "\n",
    "#plt.yscale('log')\n",
    "#!!!wp: Why eta? \n",
    "plt.xlabel(r'k', y=3, fontsize=40, fontdict = axis_label_font)\n",
    "plt.ylabel(r'$\\lambda_k$', fontsize=40, fontdict = axis_label_font)\n",
    "plt.xlim([-0.5,30.5]); plt.ylim([1e-3/2,1]);\n",
    "plt.yscale('log')\n",
    "\n",
    "#title and axis labels\n",
    "plt.tick_params(axis='both', labelsize=45)\n",
    "#plt.title('Explained variance for principal components', y=1.05, fontdict = title_font)\n",
    "\n",
    "#saves\n",
    "plt.tight_layout()\n",
    "plt.savefig('figures/OP_variance_4_3N.svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Principal component weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import dot\n",
    "components = dot(incpca.components_, incpca_ig.components_)\n",
    "components_ig = incpca_ig.components_\n",
    "#components = incpca.components_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI4CAYAAAB3HEhGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xuc3HV97/H3Z2Z2k82GiyRETCIhSA1gRBPkInZtAmrx1gasVmywpQZtae0J5/RijyjpA62eh4/T4g1tQWmPUVELRI6XcCpEDJWAkrQil3gJBkiQQECTbJK9zHzOH7OTnf3tzO8y+5uZ38y+no/HPmB+O/P7fXazm/ze8/1+P19zdwEAAAAApi7X7gIAAAAAoFsQsAAAAAAgJQQsAAAAAEgJAQsAAAAAUkLAAgAAAICUELAAAAAAICWFdhfQyczsFElnSlooqU/Sfkk7JG1x96fbWVtQu2o1s5ykl0s6XdLzx659SNIeST+T9IC7H0hwPvYVAAAAQCa4uwWPEbASGgsMl0m6UtJL6jytZGZ3SPqIu29qWXEB7azVzE6VtFbS2yQ9L+SpJTN7QNLtkj7t7o9FnZu92wAAANBuZpOyVfk4N6vxmdkCSV+T9MoEL/tnSe919+HmVFVbu2o1s5mSPirpzyXlE778Ene/KeL8zs8sAAAA2s3MGMGairHAcrekk2p82lWecnd0jc+9W9ICM7vI3UeaV+G4dtVqZvMkfVvS8jpPKUn6laR+STOSnh8AAADIOppcxGBmPZJu1uTAslHSCkl97n6MpGMkXSrpkcDz3qjyqE7TtatWMztW0p2aHK62SnqvpJMlzXD3Oe4+U9IcSa+T9DGV14IBAAAAHY8pgjGY2Z9J+lTg8Mfc/a/rPP8oSV+XtLLqcEnSGe7+YHOqPHLtttRqZl+X9DtVh0Yk/XdJ17l7KeK1JukNkp5y9x9GPJcpggAAAGi7elMECVgRzGyWpJ9LOqHq8Hfc/bURr5sr6UFJ86oOf93dV6Vf5ZFrtqVWM3ubpK9UHSpKeru7/1uswhMgYAEAACAL6gUspghGu0wTA4ur3MAhlLs/I+kDgcO/a2anp1hbUMtrHRsB+0Tg8D82I1wBAAAAWUfAinZx4PEmd98e87VfkhTc4yl4vjS1o9Y/UHl/q4qnJV0T85oAAABAVyFghTCz50l6deDwV2o9t5axDXS/ETjclCmCbaz18sDjL7n7vrjXBQAAALoJASvcb2lyK/v/SHiO4PPPHOu4l7aW12pmL9XkroGxQx0AAADQbQhY4V4SeLxP0kMJz7GlxrHTGisnVDtqHQg8PizpvoTXBAAAALoGGw2HC4aLXzTQwq7WHk+nSbqnsZLqaketZwUeP+LuxeoDYx0KT5P0QpUD/dOSHnX3nySsDQAAAMg8Ala4JYHHjyU9gbs/a2aDkvpDzpuGdtQaDFhHGmqY2fmSrpR0oWr8nJnZEypviPy/3P3JpLUCAAAAWcQUwXDHBR43GgSCrwueNw3tqPWFgce/MrOZZvZpSd+R9CbVD/ELJf03ST83s/c2VCkAAACQMQSscLMDjw82eJ5DEedNQ0trNbOcpKMChwclfUHSFZImbbpWR5+kT5hZcC8tAAAAoOMwRTBcf+Dx4QbP04qA1epaj9bkEPVOSXOrHn9X0nUqdyd8RtIxKk8r/GNJbwm89r1m9qC7/1MDNQMAAKBFdu4d1PWbd2jDtt0aHBpV/4yCVi2br8sHTtaiOcFb0umHgBWuL/B4uMHzDEWcNw2trrVW8KqEK5f0l+7+D4HPPy3pW5K+ZWZvlbReUm/V5681s9ui1mStW7fuyP+vWLFCK1asCHs6AABA00yHsFH9NR4YGpVUfpe90k3twNCobrrvcd18/y5dt3q5Vi6Z17Zas4CAFe6wpFlVj3vrPTHCjBrnTVuraw37Gv6xRriawN2/NtZh8LqqwzMl/Q9Jfxn22uqABQAA4pkOQaDVNm3foyvWb9VIsaTRUjludFvYqPU1SuPhqmK05BotFXXF+q3auHZgWv9MWfJO3tOHmT0lqfq34uPuvraB8/xI0kurDn3F3d8+1foC12hprWY2U5OnE0rSc5IWunusNWA1rveUu58Q8vwGus8DADC91btJLuRMPflcVwSBVqkE1Vu27tLB4WLoc/t68pkMG3HD9s69g7rw2s06NBL+dQb15E2/f9YLuz68m5ncfVLfAUawwh3QxNAyq94TIwSn2e1v8DxhWlqrux82s1FN/hn6ctxwNeZzkq6tevx8MzvV3R9JcA4AABJp12hOO667c++grli/teZNMqMOydQLqvWMFEu6YfOjumbV0hZUF0+SUbfrN+/QSLGU+BojRe+qUbyk6CIY7rnA47ojKxGCrwueNw3tqLXW5+5OeL3NNY69POE5AACIbdP2Pbrw2s266b7HdWBoVK7xG8wLr92sTdv3dMV1d+4d1FUbHtAF//uuyBGIQyNFveYf7tJVGx7Qzr2DqdbRLaqDapxwJZUD7K3bdjW5svjCvobRkuvQSDlsV34GNmzbHftrDaqc77Ibf6DTP7hxWv1sEbDCbQ88XpT0BGZ2nCY3hAieNw3tqLXWKNOjCS/7ixrH5iQ8BwAAsSS9wezU61aHubg3yJVRh2aGzE7W6GjO4PBozeOVALz06tu1+H3f1NKrb68bQpI8d6pfQ2XUTZIGh2rXntTB4eK0+tliimC4hwOPT7Lki4AWxzhvGtpR60OSBgLHkk5/3Ffj2DEJzwEADWHRf33d+r1JcoOZ5rSutK8b9ucjqe6UwCiNTBlM62cl6Xla/TPa6GhOf+/k2+0k0/TqPXf9lse0fstjWnnq8Vr35pfE+prjfA2jJdcXtuyUy9XXm49cZxZXO3+2Wo0mFyHM7CJJtwQOv8TdH0pwjj+T9KnA4ePcPdVpgu2o1cz+XNInA4fPdfd7E1zzWE2eavhedw/WUXk+TS4ApIJF//Vl+Xsz1RuupVfffqTNdJjZMwr68d/9dholJ7puxazevC5evqDm1xW1Dmj+sTO1Z99hjSYfbJkgTqOCJD8rYX92O54ZTPQzF+e6J8/tr9mMIvi9jfsztfh935zUOS9KIWe65OwTJ4Tm7//sGV36+ftUDAk6leYYkmI3mYgTtBr5Gqrbsaeh1vekliz/PVRRr8kFASvE2JS5PZLyVYff4+7/nOAcX5ZU3YVvm7svT6nE6uu0vFYzO1WTR7je6u7/luCaZ0j6r8Dht7v7V+o8n4AFYMridMZK0v2r0Zv+LL47G+d7kzfpC+86R+edMrfuc5rxtU3lhqtSz/otj8W6lpn06Efe2FCdtTRyY5s3qbeQnxRSGunq1qiw723cWmb15nXO4uN0z8/3jo1ijH8n4t68V/8+xgkovfmcciYNF0uq9bTK9/ZPVpysz353R92wWh1akobkYN1S+Wd4zb/8QMWIL7oSQlyeaJqnVA7Yvzo4okPDxUkh9rIbf5Co/maJegMj7b+jm4WA1SAz2yRpRdWhO939gpivnS3pSU1c1/RBd78mvQonXK/ltZrZQ5JOqzp0rbtfGa9iycz+RNJnAodPdfeaa78IWEBny0rntnzOVCx56I2dSfrdl8/XtW9fFnruRm/6s/ru7FUbHoh1Q5fPmW74w1ck/tryOdN5p8zRDx59LnEYbfSGK2nnN6n9I1jVqgNt3D+ftFV/byu/T1++7zE1sCSpIYWc6U1nvEC/PjSiTdufbs1Fq6w89XjlzXTHw3sSBeXqgJY0HBdylsqfcyFnyll5il6Lf2xCrT73RL1h6Qv0rR8/OenfhP2HR/XNHz0Z+vXHHQlrJgJWg8zsvZI+UXXIJZ1WLwAEXvtuSf8UOLzU3R9MscTq67W8VjP7O0kfrDr0lKQXuvtIzJq3SDqn6tAT7v7CkOcTsIAO1YpAUSvAnXXS82q+cx7XjZedVXeKU19vXkOjpVhTfZLuLRPn3dmowJpkr5vK85KGgFrTrZKOsESNkly/eUesUFHrhquRekxSIW8aLXqs71mcoDjVYFQJtO/90raGg9pU9eRNv3nKXN3z8711R4ZQ2/xjZ+qpXx+OHLnKujSnC+YklVR+A6H6+5IkWKb9RkhSBKwGjY3s/FwT95j6jru/NuJ1c1RuAlH9um+4+5tjXHOdpKsDhxe7+y8yWOvcsWseXXU41iidmb1D0hcDhz/k7h8IeQ0BC+hArZju0cgoRRyVuuqtEYlSfdM/1bBQLSqw1pv6FFyjsu62B6c8IlB9zjsefmpKQaKvJ6e3nLkwdF1OmOANVxojPsEA2MibBXGmtUUJ3ogCrdaTN60+d5G+eO9jGp7qIr8UpD2VN/n1CVgNqzEyJEn/y93fV+f5R0naIOn8qsMu6eXu/qMY11unBgJWO2odO8f7JX2o6tCopD909y+FvGaFpNskHVV1eJ/KX+ezIa8jYAEdKM5NbligiDNS06x1KZWpSbc/+FTD5589o6BPvmNZ4rBQGR0Kft1vWPoCvetffzilr7c3n5PkGk7xjr2vJy8zTbnrWN6kQj4nlxq6iVt97olHvmdp/ouRN6mnkNPhkfCaaq27uWL9Vg2NFhn1Qcf7wJtO11d/8Lh+umd/Jn6eV597YtvWrRKwpsDMeiR9X9IrAp/aKOkjkra4+7CZHS3pTZKu0sR1SZL0cXdfG/N669R4wGpprWPX7JO0RdIZgU/dJOnTku5x96KZmaSXSVoj6T2avE3A29z9axHXImABHShp57ZaU9aCU1OqHxdirKmainxOKpWmNjWmsu4rqeB0mULOVHKXe7qdvdKQ1pqRTq+j+s2CNEaumuGkOX167NlDylk2/szQOZYuOFq9edOPntiXqZ+dJO3q00LAmiIze6GkuyWdWOPTJZX3f6q3f9NGSb+TYF3SOjUYsFpda9U1T5J0V8g1f6XyNMJ6e6+9393/PsZ1CFhAB0rSQW3lkuPbsogdSFNl1DJOx7h2yOekYkk6/9R5uvOR7t/4FekySYvn9uuJ5w6p5I2tb22GmYWcPnPpmS1rEETASoGZLZT0b5rYlCHK5yT9mbsPJbjOOk0hYI2doyW1Bq65QNJXJL0qwcv2S1rj7l+NeQ0CFjIjiy22s2oqHdTQedLeN6dTNTpq2Uq9+ZyGW9UKEF3FJF2z6iXa/ssD+sKWne0u54hWtm+vF7ByTb9yF3H3J1QOD++W9EjYUyXdKekCd1/TaGCZinbU6u67JA1IeqekbRFP3y3p71UOjrHCFZAlm7bv0YXXbtZN9z2uA2PrPA4Mjeqm+x7Xhddu1qbt2X1HeOfeQV214QEtvfp2LX7fN7X06tt11YYHtHPvYNOuuWrZfBVyk/4NQpdqdaTI6k9W1sOVJMIVGuaS1t32kFyeqd/BkWJJN2x+tK01MII1BWb2YklnSlooaaakA5J2qLzmKFN3V+2o1cxOlHS2ytMGZ6ncxOJpSf/l7g81eE5GsNB2Wd0AMc6IWjNapVeue8vWXROaG+Rz5RvfDDSaQpdjxKz78GfaOdq95rGWVrVvZ4ogugIBC1kw1Y54cSWZghgnOJ08tz/1YJhGd7Qs/uPcCG4IgXTkTZp/bJ8ef+5Qu0vJtL6evCTXoYiulmF686bf/I3jdd+jz2pweFTdcovVqvbtBCx0BQIWsiBpR7xGxA1MtUaO6onTCS9Jq/Q4m+xOJ52w3gbIut58TjlToo2M8yaVUuiqWfn77/VLT9ANmx/VnTGmeudN+sK7ztF5p8zVzr2D+rvbHor1uqnWWfl34I9v/EHDX/f5p87T1W8+fcIbamnsG5cF7R7BqtfRDQC6TlpNKQZjNmsYHG6sqcPOveXNVWuNNI2WXKOlot514w+U9D3LOEstRkuuL2zZqfVbdqqvN68Fx/Zp168OHQlw1aM0U93rqFuYpJPmzNIv9h5sdylNxyhda+SsfONWnmLbXd/x85fM0/d//szY32UT3zzK50wu6XCCecWFnMlMKk6xVWP16P3OvYO6Z8feyNdUOtadd8pcSdKiOf36/GVnxQ5a/b15nfeiObrzkadVTPDm8SVnn6g1A4u1aE6/+mcUEjUQ6u/N6+LlC4+8PujygZN18/27NFpq/O/3WltLVP5sW7E5cSFnumjZgqZfJ7SGtl4dAFqk1ohQpSnFzffvmjAiFBXA4v6D1t/b2F+x12/eoZGINNTsf6Jc5QD10z0HJh3HRDN78npq/9C0+N7kc6aZPXk6QjbZjEJen/vDV+jbP/6lbt76RNe8mTF7RuFIALlh86O6ddsuDQ6Pqr+3oIuWLdC+wyP65o+eTHTOnnxOhxNsuF3r5r8yGlT5ez7O38Em6cKlJ9Rcs1odtGp9ncFws2n7nljt/GvNMFi1bH6qU9YXzenXdauX60++cL+GGgxDl5x9Yip/to3qyee0ZmBxS65VD1ME0VGYIohGxGlKUZmWUutd1WDzh2avwaKlOZrluFk9eu7gSMNh0Ez67l+uiPx96jStHpnLjU1py5sm3FQH/76pvDF0eKTY8QE+zt+Jjfzdd+NlZ+m9X9oW63WzevN6y/KFkYGnFdPAg+JsSF1rjWyzmi41Ot2x3vekFf+uTaVZU6OYIghg2orzbmS9VsWVKXlXrN965B+oNyx9gb583+Oh56t+B61el71ZvXldvHzBpBGyuFMQka6cpELeNJzFXWFT8uzBRHvIT9LfW4j1+9SImT05jRS95evYak2Z2rR9jy678QdNu+YfnLNIr196gr7941/WvdkPmyocVNlzp9afSmV61qtOmdvWDYXjjCo08nffyiXzYo/ivGX5Ql2zamnkG1/NngZey3mnzNUNf/iKyLW3wZBUGXFK+roowVG4L9/3WOR01bCpeUn+bM9fMq9msKu8EWKS+nrzWvi8Pu167pAOjhTrhuV2IWAB6DhJ11Jt2LZ7yusYKvtqnH/aPF2xfquiRlLPfdFxksK77B0cLurL9z52ZIpi5R23pHPqkQ4bW/fysoVH67+e+HW7y0ms3ihsZbRkqio3T7du25X6uqBCzvTWM1+o80+bF2uqVBrXC3une+WSeZrVm2/K1Lx8znT+afN03ilzdd4pc+ve7CcJsiWVR8NqDXO9+sXHH2lksPh932zaSFicWQBRN75J/+6bPaN8Gxtn3VCSaWPNngZez8ol87Rx7UCsaYVpvC6ORXP6dc2qpVozsDhypCzsexz3exo1jTQrASoKUwTRUZgiiEb2cUrrpmJmT04mi/WOctIpR9VTOLqli1OnSqP1cSsFO0tW35T09aTX6bHyM7riY99N/PsUFfKqf/7jTJVKylQenRwteewbtWb+HsaZspXWlKrqazVjmlbYz1/Sm+Ik3/PglMM09/hr1VYcnWYq3+Nu/Z7Sph1dgYDV/cJGpyQ1NNc8zZuKZq3VMEm/+/L5uvbty5pyg4n4CjnTyXP7teOZwUyF3OqpXpU9a6JuYBsJCWFNAFYumZfo9+nScxdpzcBi7XhmMNGNWb0buUZ//xpdg9Lo2s0ocW4k03pjqPpacX4e6nUwbOTnL6k43/OKeuuR0hj1yOpm8lnQ6Pe4W7+nBCx0BQJW99q5d1DrbntQm7Y/PelzlZuwc190nDb/5JnId8DedMYLNHtm4UhIq+xP1Ak/OWcsPFo/emJfu8voCM0MQrN683JXZho5RLVWrqeRNxcuPXdR6M1To+9EJ70xC3t+vcBWq46pLHqPux9ddZ1x/4mKapCQ5htDlWvFvcmtdDBsx/SsqM3L8yb1FvJNb2SQ5ogYyrrxe0rAQlcgYGVfI3tNbdq+Z0otYTE9VW4E3/WvP0w9CJlJn/+js2LdxDeqEkIq6/qacdORdBQkTle0rLwTHQxgfT3NWfSeNBjG/Z6bSY9+5I11P5/mFMXqa3XCTW7lex5sUd/oGw1TraNT1wFlUbd9TwlY6AoErGxr5B/uJFNCAKl+K+s0g1D1O/7VNwN5S280NLixaTNuOpKMgiRZ/9AJN+ntklaL7zT/bgxeq9tucoF2IWChKxCwsivuzUCwNTkNHaan2TMKOvuk43T3z56ObIuet/JmvmGjErVuGM8+6Th9/+fPpLo+JslNb5K91Zolye9X0lEnbtJrS3Mxfxqj+53YOADoFAQsdAUCVvZUpgQmCUmt3tQT2RG82QvbzDKNIFJrGllUV72ooJFkjchUu6pNVdxAOLOQ02cuPXPajjqlKe0plGG/I8GNiqd6LQDJELDQFQhY2dKMqVlIJiepp5BTscYoSSN/JmY6spZlx9MHlPayuHo3e60cDUljeltW1ojEEfV7ev6p847sk4R0NGMKZb3fkaULjta62x5iuibQBgQsdAUCVnawdiobKl0Tj5rZM+nGa/2WnYlGCoOL7hv5Mw5r8Zylm73pNr1tun29WdDK7zl/vkB7ELDQFQhY2dGpa6dykkqK3vi0k9RbLJ+0zXOt8yQdpWx3i2cAAFqlXsAqtKMYAJ1vw7bdHReuJCmfNy1beIz+8/FfK/aGNTH0FnIablOb+cHh2iFq1bL5sUNwIWe6aNmCScdXLpmnjWsHjrw7XglswXV01SNU550yV+edMpdF9QCAaYkRLHQURrDap9LM4patuyasOekkMws5vf6lL9BPn9qvR365P9WAeONl8fZManRtVJh6I1hJpvglXXTPdCQAwHTHFEF0BQJWe0R1TesElYX8n71rh27Z+kSqmxrX2zOp3sanwcYIYaI6Lka1YE7S8S4La6MAAOgUBCx0BQJW63V6M4u3n/VC3bL1Cb31FS/U1/9zd6I1SXE0ssdMnPVrJuk1p83T3T/bO+V2z53U8Q4AgE5BwEJXIGC1Xqc2s5DKmxqfs/g4bdr+dEOvn9mT00jRp7RnUi1J9snZ8cxg6u2eAQDA1BGw0BUIWK1TWXO1fstj7S6lIbmxv+4azYWVkanzT5vXlICTZJ8c1jwBAJA9BCx0BQLW1FWC04ZtuzU4NKr+GQWtWjZflw+cfORmnQ2EJ45MNSvgEJwAAOhcBCx0BQLW1NQLTtWNFGb15jU0WgqdFpcVvYWcSiWfNAJU+RkpNvAlMPUOAADEQcBCVyBgJdcN7dXrueuvVtQcAUrSpa/CTIwgAQCA2AhY6AoErGS6eapfvb2fJGnx+74Z2to8ybkAAABqqRewcu0oBkDz7dxb7j53aKTYdeGqkDNdtGxB3c/3zyikdi4AAIAkCFhAl7p+8w6NFNPbTLcRhZypryevK1/7G+rryauQm/QmT0N68jmtGVhc9/Orls2Pfa2ocwEAACRBwAK61IZtu9s+cnXJ2Sdq49oB/bcLXqyNawd0ydknavaMgszK0/IuPXeRvrTmHF167qIJx89fMk8zC7lJIakS2K5bvTx0jdTlAyerJx/919vMQi7yXAAAAEmwBgsdhTVY8SVdh5S2qa5rmmoL86j1Z+efOk9Xv/l0whUAAGgITS7QFQhY8S29+nYdGBpty7Urm/Res2ppW65fwT5TAACgWQhY6AoErHH12q/P6s3r4uULtP/wqL75oyfbMk2wepNeAACAbkTAQlcgYJVVpr8NjRZVKz/lTSrkc3JJw6PpNrro783rnMVz9P2fP6PRGpv8skkvAACYDuoFrPi9jAFkQnX79XqKLhVHS6l2sZlZyOkzl555JDgx/Q4AAGAyRrDQURjBkq7a8IBuuu/xlk79oyEEAADARIxgAV2iVe3X+3vzunj5QkakAAAAEiBgAR1msImdAVlDBQAAMDUELCDjKt0CN2zbrcGh0VT3tjJJPvbf/hmsoQIAAJgqAhaQYVGb5U7V6nMXEagAAABSRMACMipOt8CpmD2j0PaNgAEAALpNml2cAaTo+s07NFJMdw+rikLOdNGyBU05NwAAwHRGwAIyqpndAnvyOa0ZWNyUcwMAAExnTBEEMqbS1OJAE7sFXrd6OeuuAAAAmoCABWRIdVOLZqIFOwAAQHMQsICMaHZTi4rZM/i1BwAAaBbWYAEZ8Q///hMdbnK4orkFAABAc/FWNpABm7bv0df/c3fTr0NzCwAAgOZiBAtos8rUwDQVcjbpcV9PnuYWAAAATUbAAtos7f2uZvXmdcnZJ2r2jILMymuuLjn7RG1cO0BzCwAAgCYz9+bsswM0g5l5t/3Mnv7BjTo4nM7aq0LOdMnZJ+qaVUtTOR8AAABqMzO5uwWPM4IFtNGm7XtSC1cSa6wAAADajSYXQItUNhDesG23BodG1deb19BoOlMDCzlTTz7HGisAAIA2Y4ogOkqnThGs3kB4tDT1+k1SIW8aLbn6ewu6aNkCrRlYTLgCAABokXpTBBnBApqsGRsIz+zJa+PaAQIVAABAxhCwgCZLs0sgUwEBAACyjYAFNNmGbbunPC3QTEwFBAAA6AAELKDJBodGG34tbdcBAAA6C23agSbrn9H4+xi0XQcAAOgsBCygyVYtm69CblKDmVhYawUAANBZCFhAk10+cLJ68sl+1UzSqpfP18ol85pTFAAAAJqCgAU02aI5/bpu9XLNKMT/dZvZk9eVr31xE6sCAABAM9DkAmiynXsHdcfDTykfY5pg3qTeQp6pgQAAAB2KgAU00abte3TF+q0aKZYmtGo3ScHG7f29eV28fCFt2AEAADoYAQtokp17B3XF+q06NFKc9LlKuOrryWvj2gECFQAAQJdgDRbQJNdv3qGRYin0OSPFkm7Y/GiLKgIAAECzEbCAJtmwbfeEaYG1jJZct27b1aKKAAAA0GwELKBJBodG4z1vON7zAAAAkH0ELKBJ+mfEW+LY38tSSAAAgG5BwAKa5ILT5imqMXshZ7po2YKW1AMAAIDmI2ABTbBp+x5t/PEvJ7ViD+rJ57RmYHFLagIAAEDzEbCAlFXasw+NhncQnFnIsaEwAABAlyFgASmL057dJF249AStXDKvNUUBAACgJVhdD6Rk595BXb95h9ZveSzyuS7pOw/vaX5RAAAAaCkCFpCCTdv36Ir1WyNHrqrRnh0AAKD7ELCAKaqsuTo0Ukz0OtqzAwAAdB/WYAFTFGfNVRDt2QEAALoTAQuYog3bdmu0FNWQfSLaswMAAHQnAhYwRYNDyddS0Z4dAACgOxGwgCnqn5F8LRXt2QEAALoTAQuYolXL5quQs9jPn91AIAMAAEBnIGABU3T5wMmIOXW9AAAgAElEQVTqycf7VaK5BQAAQHcjYAFTtGhOv65bvVwzCtG/TjS3AAAA6G4ELCAFK5fM0/+78tU6v87aqkLO1NeTp7kFAABAlzP3ZO2lgXYyM8/6z+x5H71DZy06Tnc8skeDw6Pq7y3oomULtGZgMeEKAACgS5iZ3H3SQnxW2wMpcnc9NziiD1/8Un2cZhYAAADTDlMEgRQdGBqVGZ0CAQAApisCFpCip/YN6flHz2x3GQAAAGgTAhaQoj37D+v4o2a0uwwAAAC0CQELSNHT+xnBAgAAmM4IWECKntp3WPMYwQIAAJi2CFhAivbsGyJgAQAATGMELCBFe5giCAAAMK2x0TA6ShY3Gt65d1DXb96hDdt268DQqPp6cnrLmQt1+cDJbCwMAADQpeptNEzAQkfJWsDatH2Prli/VSPFkkZL43UVcqaefE7XrV6ulUvmtbFCAAAANAMBC10hSwHr+z97Rpd+/j4VS/Xr6evJa+PaAUayAAAAuky9gMUaLKABm7bv0aWfuzc0XEnSSLGkGzY/2qKqAAAA0G4ELCChnXsHdcX6rSrGGEgbLblu3bar+UUBAAAgEwhYQELXb96hkWIp9vMHh0ebWA0AAACyhIAFJLRh2+4JDS2i9PcWmlgNAAAAsoSABSQ0OBR/RKqQM120bEETqwEAAECWELCAhPpnxB+R6snntGZgcROrAQAAQJYQsICEVi2br0JuUkfOSfI503Wrl9OiHQAAYBohYAEJXT5wsnry4b86eZO+8Mdns8kwAADANEPAAhJaNKdf161err6e/KSRrELO1NeT1w1/dJbOO2VumyoEAABAu5h7/G5oQLuZmWflZ3bn3kHdsPlRffHenXIvr826aNkCrRlYzLRAAACALmdmcvdJ60YIWOgoWQpYFef8/Xe04c9epRcc09fuUgAAANAi9QIWUwSBKSiWXM8ODmvu7BntLgUAAAAZQMACpuCZA0M6pq83sukFAAAApgfuCoEpeGrfYZ1wDKNXAAAAKCNgAVPw1L4hPf+ome0uAwAAABlBwAKm4Kl9hzXvaAIWAAAAyghYwBTs2XdYzz+aKYIAAAAoI2ABU/DUviE9nxEsAAAAjGEfLHSULOyDtXPvoK7fvEMbtu3WgaFRzezJ6ffOXKjLB05mg2EAAIBpgo2G0RXaHbA2bd+jK9Zv1UixpNHSeB2FnKknn9N1q5dr5ZJ5basPAAAArUHAQldoZ8DauXdQF167WYdGinWf09eT18a1A4xkAQAAdLl6AYs1WEBM12/eoZFiKfQ5I8WSbtj8aIsqAgAAQNYQsICYNmzbPWFaYC2jJdet23a1qCIAAABkDQELiGlwaDTe84bjPQ8AAADdp9DuAjqZmZ0i6UxJCyX1SdovaYekLe7+dDtrC+qkWrOqf0ZBB2KErP5efq0AAACmK+4EEzKznKTLJF0p6SV1nlYyszskfcTdN7WsuICs1Wpmt0t6XY1PLXb3XzTz2mlYtWy+brrv8dBpgoWc6aJlC1pYFQAAALKEKYIJmNkCSXdLukH1A4tU/r6+VtKdZvZPZtbbivqqZa1WM7tUtcNVx7h84GT15MN/ZXryOa0ZWNyiigAAAJA1BKyYqgLLK2t82iXtq/PSd0u6xcx6mlVbUNZqNbO5kv4hzXO2w6I5/bpu9XL19eRVyE3syFnImfp68rpu9XJatAMAAExjBKwYxgLHzZJOCnxqo6QVkvrc/RhJx0i6VNIjgee9UdJHm1tlWUZrvVbS3JTP2RYrl8zTxrUDuuTsE5XPmUzS7BkFXXL2idq4doBNhgEAAKY5NhqOwcz+TNKnAoc/5u5/Xef5R0n6uqSVVYdLks5w9webU+WRa2eqVjO7UNK3qw7dIuniwNNir8Fq50bDQW/+5N360KqletkLj213KQAAAGgxNhpukJnNknRV4PB36gUWSXL3/ZLeJmlP1eGcpA+nX+G4rNVqZv2SPlN16IeSPj3V82bFnv2HNe/oGe0uAwAAABlCwIp2maQTqh67pD+PepG7PyPpA4HDv2tmp6dYW1DWav2QxqcqFiW9R+XRsY5XKrmeHRzWnH4CFgAAAMYRsKIFp7NtcvftMV/7JUkHIs6XpszUamZnSfqLqkOfdPetjZ4va549OKyjZvaot8CvEAAAAMZxdxjCzJ4n6dWBw1+J+3p3PyDpG4HDq6ZaVy1ZqtXMCpKu1/jP1xOSPtjIubJqz74hHT+b0SsAAABMRMAK91uavBnzfyQ8R/D5Z5pZM7oiZKnWv5L0sqrHfzG21qtrPH1giPVXAAAAmISAFS64Qe8+SQ8lPMeWGsdOa6ycUJmo1cxO0cTRqv/r7rcmrCPz9uw7rOOPImABAABgIgJWuGC4+EUDPcJ3xDhvGrJS6z9Lmjn2/4OK0WSjEz19YIiABQAAgEkIWOGWBB4/lvQE7v6sykEj7LxpaHutZvYuTdxP62p3T1xHJ9izb0jzjpoZ/UQAAABMKwSscMcFHj/Z4HmCrwueNw1trdXMni/pY1WH/kvSxxusIfOe3s8IFgAAACYLNkXARLMDjw82eJ5DEedNQ7tr/aSk5439f0nSe9x9tMEaMmnn3kFdv3mHNmzbrQNDo7rzkT2679G9unzgZC2a09/u8gAAAJABBKxwwbvmww2epxUBq221mtmbJb216tA/ufu9DV4/kzZt36Mr1m/VSLGk0VJ5aduhkaJuuu9x3Xz/Ll23erlWLpnX5ioBAADQbkwRDNcXeDzc4HmGIs6bhrbUamZHSbqu6tAvJf1tg9fOpJ17B3XF+q06NFI8Eq4qRkuuQyNFXbF+q3buDS5fAwAAwHTDCFa4w5JmVT3ubfA8wcU6jY4uhWlXrR+RtLDq8Vp3/3WD145l3bp1R/5/xYoVWrFiRTMvp+s379BIsRT6nJFiSTdsflTXrFra1FoAAACQbQSscAc0MbQ02jYuOAp0oMHzhGl5rWb2Skl/WnXodnf/SoPXja06YLXChm27J41cBY2WXLdu20XAAgAAmOaYIhguGC5m1XxWtGBo2d/gecK0tFYz65V0g8Z/hg5JuqLBa2ba4FC8Xh2Dw13V0wMAAAANIGCFey7w+IQGzxN8XfC8aWh1rX8r6fSqxx9y91obFXe8/hnxBnr7exkQBgAAmO4IWOG2Bx4vSnoCMztOkzvxBc+bhpbVamazJP3PqkMPaeIeWF1l1bL5KuQs9DmFnOmiZQtaVBEAAACyioAV7uHA45PMLPxOe7LFMc6bhlbW2quJTTROlzRsZh71IWlTjfM9GnjeLxLW3VSXD5ysnnz4r0pPPqc1A7W+fQAAAJhOCFjhHgw8PlrSaQnPcW6NY80IWJ1Ua0dZNKdf161err6e/KSRrELO1NeT13Wrl7PZMAAAAAhYEe6SVAwc+82E5wg+f5u7N2MNVifV2nFWLpmnjWsHdMnZJ2pGofxrM3tGQZecfaI2rh1gk2EAAABIok17KHd/1sw2S1pRdfj3Jf1znNeb2WxJbwocvjWd6iZqca0HNbE9exIvlnRl4NhVkvZWPW5Gl8UpWzSnX9esWqoFz+vTc4PD+ts3JB0gBAAAQLcjYEW7RRNDy0ozW+LucRpVvEOTm0bcklZhNbSkVncflvTZRgo0sxWaHLC+6O6/aOR87XBwaFR9vfl2lwEAAIAMYopgtBsl7al6bJI+FfUiM5sj6ZrA4W+4e3CtVK3XrqvRIOKkLNY6HR0cLtKSHQAAADURsCK4+wFJHwocfo2ZfbTea8zsKElflVS9MMclvT/9Cqsu0EG1drLB4SIjWAAAAKiJgBXPZyX9MHDsb8zs22b2ajPrlSQzO9rM3iHpXknnB57/CXf/EbV2vkPDo+qfQcACAADAZMxzisHdR8zsYkl3Szqx6lMXjn2UzGy/pGPqnGKjpL9qbpVlnVRrpxocLmoWUwQBAABQAyNYMbn745JepfKIT1BO9QPL5yStcveRZtUW1Em1dqJDw0XNYoogAAAAaiBgJeDuT6gcXN4t6ZGwp0q6U9IF7r7G3YdaUd+EAjqo1k4zODzKCBYAAABqMndvdw0dy8xeLOlMSQslzZR0QNIOSfe4+56w17ZaJ9Uaxsy83T+zv/2P39O1b3+5TnvB0W2tAwAAAO1jZnJ3Cx7nbfgpcPefSPpJu+uIo5NqzbqDI6O0aQcAAEBNqdwlmtkJkuar3Or7WEm/kvS0pN3u/mQa1wCy4uAQbdoBAABQW0MBy8yOlrRK5fbev6WJ3eqCz31c0l0qr/O51d33NXJNICsODhdp0w4AAICaEjW5MLOXm9m/SnpS0o2SLpW0SJKFfJwoabWkz0v6pZn9HzNbntpXALRQseQaGi1qZoGABQAAgMliBSwzO8PMbpN0v8phqU/jAcpjfFSeO1PSH0j6gZndZmZnpPrVAE12aKSomT155XKT1jMCAAAA0QHLzD6rcrB6o2qHqrDRq7Dnv1HS/Wb2mVS/IqCJDtKiHQAAACHi3Cm+W+VgJE0MSJL0uMqb2W5TuanFXkn7VN7I9jiVm14sk3SOyu3BVXUeScqPnf9PG/4KgBY6OMQmwwAAAKgvyVvxlVB1j6QvS7o5SYdAM3uBpLdIukTSKxNcF8iMg8MELAAAANQXt8nFIUkfl7TY3V/l7p9K2n7d3Z8ce92rJC2W9AlJB5OVC7TXweFR9c9giiAAAABqixOwPiRpkbtf6e4707iou+9097WSTpL04TTOCbTCICNYAAAACBH5Vry7f7BZF3f3ZyQ17fxA2g4NjxKwAAAAUFeifbCA6W5wqEgXQQAAANSV2p2imR0tabakA+6+L63zAllycIQpggAAAKivoYBlZjMlvVXSmyQtV7lphVV93iU9qvL+Wd9QuePgoSlXC7TZwSGmCAIAAKC+xFMEzewKSbsl/Yuk35P0orHzVG8snBs7/lZJ/yppl5mx1xU6XrnJBVMEAQAAUFvsgGVmeTP7mqRPSjpW42HKQz4qzzlW0qfM7Gtmxtv/6FiHhkfVP4MfYQAAANSWZATrEypvFBwMVWGCYetilffTAjrS4HBRfYxgAQAAoI5YAcvMXi7pPZoYqizBhzQest5jZi9LqX6gpQ4NF9XPGiwAAADUEfet+HeqHMYqIckl3T728bCkpyQdlDQiqUdSv6TnSzpN0uvGPipykv5I0pVTrh5osUGaXAAAACBE3ID1Go2HqxFJb3L3f4/xuo2S/tHMLpR0m6T82DkuaKBWoO0O0uQCAAAAIeKuwVow9l+X9M2Y4eoId98o6Zsany64IOTpQGYdpMkFAAAAQsQNWLOq/n9Pg9d6qs75gI5xcLiovh5GsAAAAFBb3IBVCVUmaaWZ9SS5iJn1Sjpf4w0ynk7yeiArDg4XGcECAABAXXED1oMab25xiqRvmdkr4rzQzM6S9K2x12nsHD9OWCfQNjv3DuqqDQ9o6dW367FnD+oNH9+sqzY8oJ17B9tdGgAAADLG3KO2spLM7D2SPqOJXQQl6dca7yI4KGlU410ET5C0RNIx1acae+2fuvs/p/MlYDoxM4/zM5uWTdv36Ir1WzVSLGm0NH7dQs7Uk8/putXLtXLJvJbVAwAAgGwwM7m7TToeM2DNUDlILaocqvp02AlqPe8Xkk5396HICwMBrQxYO/cO6sJrN+vQSLHuc/p68tq4dkCL5vS3pCYAAABkQ72AFWuK4FgY+kNJhyqHqj7CNhgOPu+gpHe2MlyZ2evM7IOVj1ZdF53v+s07NFIshT5npFjSDZsfbVFFAAAAyLq4a7Dk7pslvUHSLk0emar3UWGSnpD0Bnf/jynWnNTrJV1d9QHEsmHb7gnTAmsZLblu3barRRUBAAAg62IHLEly9+9JOk3S30l6TOGjV5WPnSoHm9PHQlo7VGoBYhscGo33vOF4zwMAAED3S7yhj7sPqhyw/s7MXiLpTJU7BB6j8v5WB1VufvEzSfe7+4PplQu0Tv+Mgg7ECFn9veyLBQAAgLIp3RmOhScCFLrSqmXzddN9j4dOEyzkTBctW9DCqgAAAJBlmX7r3cx2pHCa4zTeZAOI7fKBk3Xz/bs0WqrfRbAnn9OagcUtrAoAAABZFqtNe7uYWUnphKMj+3e5e37KhaFt2AcLAAAAWdDwPlhmttLdNzWxsBXu/t06n6sErFQuJQJWx2t1wJLK+2HdsPlR3bz1CR0cLmr2jIIuWrZAawYWs/8VAADANDWVgFWSdI+kj7j7N1Is6A2S3i/p3Hqhh4CFoHYErIqHdu/Tf//qf2rj2le35foAAADIjnoBK+4arHMlfX1sTdSXJd3k7g81UMSLJV0y9vEbGt+MOPKlSa8VkN15kOgYJXfljKV8AAAAqC9JkwuT9CKVR53eb2Y/VXlk615J2yQ9LWmvpP2Sjla5ucTxkl4u6SyVQ9ppVeeKY+/YeSTpIUlvTFBvxQclXdbA64AJ3KVcop3jAAAAMN3ECVgPqxyMKqNAlXD0YpVHod4Z81rVoar6XA+HvOZeSW8Y+//fkPRLdx+Keb3yBcz2JXk+UE+RESwAAABEiPN+/BmS1qo8QlWZ0lf5sAQfwdftHTvvy0KufW/V/xdU3tQYaAumCAIAACBKZMBy96K7f0LSyZL+StJjqh2aoj4qr9kl6W8kLXb3T7h7/U2GxgNWZcTr3CRfHJAmd1eOfAUAAIAQsddguftBSf/bzP5R0uskvV3S+ZIWxnj5k5K+K+kmSd+KCFXV7pO0s+rxorj1VvmZpLsaeB0wQbEkRrAAAAAQKkmTC0mSu5ckbRz7kJktlvRKSfMlzZN0rKRfSXpG0m5JW9z9Z40U5+6/krS4kddWnePTkj49lXMA0tgUQYawAAAAECJxwApy90clPZpCLUCmlZgiCAAAgAg0nQZicmeKIAAAAMIRsICYiiW6CAIAACBcSwOWmf25mZ0W/Uwge1iDBQAAgChTXoOV0CckuZn9UtKdku6QdIe7P97iOoDEylME210FAAAAsqzVAUsq74X1AknvGPuQmf1cY2FL0iZ339uGuoBQTBEEAABAlHYErMqmwdV3qqdIepGkd6s8wvWAxgPX99x9sLUlApOVuwgSsAAAAFBfOwJWhQceW9V/XybpDElXSho1s/s0HrjucffRllUJjCkxRRAAAAAR2hmwqm9VXRMDl1V9vkfSeWMfH5B00Mzu1vj6rW0tqBVgBAsAAACR2rUGKxioKscV8TlJ6pf0urEPmdmz7n58E+oEJih3EWx3FQAAAMiyVt8urpZ0vaSfanyUKk6w8sBH9WuPa27JQFmJjYYBAAAQoaUjWO7+JUlfkiQze76k36r6OD34dE2eNhg8zt0uWsaZIggAAIAIbVuD5e5PSfrq2IfMbK6kV2s8cL1UtQNUddACWqbcpr3dVQAAACDL2tnkYgJ3f0bSLWMfMrPnSRpQOWydr3Jnwcr0QKDlSi7lSFgAAAAIkckl+2aWk/QbkpZIerGkxe2tCKCLIAAAAKJlYgTLzHokna3x6YGvVLlb4KSniqmBaJMSUwQBAAAQoS0By8xmqhyiKmuuzpE0s/opgZfUa3ghSb9uRo1AUMmlPAkLAAAAIVoasMzswyqHqrNU3kD4yKcCTw3bB2uvpM2S7pL0PUn/mXKZQE0ldxlTBAEAABCi1SNYf6vajSrCAtWTKgep70n6nrs/2LzygPrKbdrbXQUAAACyrF1rsGqto6rcuu7UeKC6y91/1rKqgBDlNu0kLAAAANTXziYXlYYVlf/eJulD7v7DNtYE1FVyEbAAAAAQqt1t2qvvVt8s6V4ze8jMPmtm7zCzhe0qDAiiTTsAAACitHMEq966q1NV3v/qckkys8qUwbtUXoP185ZVCFQpsQYLAAAAEVodsD6uchfBl2ni6FlYG/aTJC2SdKkkmdmTKncRpOkFWoo27QAAAIjS0oDl7ldKkpkdJek3Nb4P1pma2LY9rKvgfElvG/uQme1193nNqhmooE07AAAAorRliqC775f07bEPmVmfxjcefrWkczVx4+GwwDWneZUC40olpggCAAAgXDvXYB3h7ock3Tn2ITPrkXS2xgPXqyTN1njQqvyX2120DF0EAQAAECUTASvI3UfMbIukoqSSykHqtVVPqbR2B1qm5K4cQ1gAAAAIkZmAZWYzNHmaYF/waa2uC6goj2C1uwoAAABkWdsClpnN1niji1dLeoUmNrrgVhaZUiq5CoV2bx0HAACALGtpwDKzVRoPVMFW7cFAVWsKIFMD0TYld9q0AwAAIFSrR7BuUTkg1bpLjROoajW32JdOaUC4kkv0uAAAAECYdq7Bqheoqj9XK1ANSbpH0h1jHz9oSnVAQMmdLoIAAAAI1a6AFQxOwUBV/bmSpPs1HqjudvfDTa8QCCiVXHkCFgAAAEK0u4tgrUAlSdslfUflQPVdd/9VS6sCamCKIAAAAKK0M2BV36ru0vgI1R3uvrs9JQH1MUUQAAAAUdoRsEzSc5I2aTxQ/aQNdQCJlANWu6sAAABAlrU6YP2NyqFqm7vTbh0dhTbtAAAAiNLSgOXuH2vl9YA0lddgEbAAAABQXy76KQCkchdB1mABAAAgDAELiKk8RbDdVQAAACDLuF0EYmKKIAAAAKIQsICYaNMOAACAKAQsIKbyGqx2VwEAAIAsI2ABMZVctGkHAABAKAIWEFPJnTVYAAAACEXAAmIqr8FqdxUAAADIMgIWEFOpJOUZwQIAAEAIAhYQE1MEAQAAEIWABcTEFEEAAABEIWABMZVc7IMFAACAUAQsIKaSO23aAQAAEIqABcRUcokBLAAAAIQhYAExlUrOFEEAAACEImABMTFFEAAAAFEIWEBM5Tbt7a4CAAAAWUbAAmIqlegiCAAAgHAELCAmpggCAAAgCgELiImNhgEAABCFgAXEVG7TTsICAABAfQQsIKbyCBYBCwAAAPURsICYSu7KE7AAAAAQgoAFxFQqiTbtAAAACEXAAmJiiiAAAACiELCAmGjTDgAAgCgELCCmkos27QAAAAhFwAJiKpWcNu0AAAAIRcACYmKjYQAAAEQhYAExlVyswQIAAEAoAhYQE10EAQAAEIWABcRUXoPV7ioAAACQZQQsICamCAIAACAKAQuIiSmCAAAAiELAAmIqOVMEAQAAEI6ABcRUcilPwgIAAEAIAhYQE1MEAQAAEIWABcRUKhGwAAAAEI6ABcRUcsn4jQEAAEAIbheBmErurMECAABAqEK7C+hkZnaKpDMlLZTUJ2m/pB2Strj70+2sLahVtZrZbEmnS1oiaY6k/rFrPSvpIUk/cvfRtK7XSqzBAgAAQBQCVkJmlpN0maQrJb2kztNKZnaHpI+4+6aWFRfQilrNrCDpfEkXjv33DElhKeSgmX1D0ifd/e6k12unUkm0aQcAAEAopggmYGYLJN0t6QbVDyxS+fv6Wkl3mtk/mVlvK+qr1opazWydpF9Kul3lEPcyhYcrSZol6W2SNpvZV81sTtzrtVvJXfkcCQsAAAD1EbBiqgosr6zxaZe0r85L3y3pFjPraVZtQS2s9UKVpwHWU1J5auBInc+/VdL9ZrYw5vXaiimCAAAAiELAimEscNws6aTApzZKWiGpz92PkXSMpEslPRJ43hslfbS5VZa1udYnJX1S0u9Jmi+px93nuHuvyuuyrpK0N/CaRZK+Y2azGrxmy5RcYgALAAAAYQhY8bxb0jmBYx9z99e7+13uPiRJ7r7P3ddLOltScD3TWjMLm6rXybXeKWmVpBPd/S/c/WZ3f9LdS5UnuPvD7v5hSS+VdG/g9Usk/c8E12u5UsklScYIFgAAAEIQsCKMjaxcFTj8HXf/63qvcff9Kq8z2lN1OCfpw+lXOK4Ntd4j6Xx3v8Ddvx6nO6C7P6nyKNnjgU+tNbOZMa7ZFqy/AgAAQBwErGiXSTqh6rFL+vOoF7n7M5I+EDj8u2Z2eoq1BbW0Vne/spHOg+6+V9K6wOF+SRckPVerMD0QAAAAcRCwol0ceLzJ3bfHfO2XJB2IOF+aOqnWW1UOgNWWNfF6U1JyZ3ogAAAAIhGwQpjZ8yS9OnD4K3Ff7+4HJH0jcHjVVOuqpZNqHbvec5rc8OKEWs/NgpK78gQsAAAARCBghfstTd6M+T8SniP4/DPN7NjGS6qrk2qtCLaDD45oZQZTBAEAABAHAStcsJPePkkPJTzHlhrHTmusnFCdVKvMbK7KreKrPdmMa6WhWGIPLAAAAEQjYIULhotfuHvSUZYdMc6bhk6qVaq9vuuHTbrWlLm7yFcAAACIQsAKtyTw+LGkJ3D3ZyUNRpw3DR1Tq5W7RbwncHi/pO+mfa20lFy0aQcAAEAkAla44wKPG53CFnxd8Lxp6KRaL5O0PHDs8+4+3IRrpaLkTBEEAABANAJWuNmBxwcbPM+hiPOmoSNqNbPFkv4hcPhXavImzFNVKtGmHQAAANEIWOH6A48PN3ieVgSszNdqZn2SvqbJzS2udPen07pOM5SnCLa7CgAAAGRdsK03JuoLPG50CttQxHnTkOlax9Zd/YukMwOf+oq7/0uSc61bt+7I/69YsUIrVqyYWnExMEUQAAAAcRCwwh2WNKvqcW+D55lR47xpy3qtH5P0tsCxrZLWJD1RdcBqFdq0AwAAIA4mPYU7EHg8s8HzBEeBgudNQ2ZrNbP3SfofgcPbJb3e3ZvxvUidu5TjtwUAAAARuGUMF7z5n1XzWdGCoWV/g+cJk8lazexPJX0kcHinpNe6+56pnLuVmCIIAACAOAhY4Z4LPD6hwfMEXxc8bxoyV6uZ/YGkTwUOPynpAnd/vNHztgMBCwAAAHEQsMJtDzxelPQEZnacJnfiC543DZmq1cwulvSvmvgz9oyk17j7zxs5ZzuV3EW+AgAAQBQCVriHA49PsuSbIS2Ocd40ZKZWM3u9pC9Lylcd/rWk33b3h5KeLwtKLuVJWAAAAIhAwAr3YODx0ZJOS3iOc2sca0bAykStZrZC0s2a2KetY0sAABp/SURBVMXwgMoNLbYmrCczmCIIAACAOAhY4e6SVAwc+82E5wg+f5u7N2MNVttrNbNzJf1fTWyUcUjS77j7PQlryZRiiSmCAAAAiEbACuHuz0raHDj8+3Ffb2azJb0pcPjWqdZVS7trNbNlkr6tiWu4hiVd7O6b4p4nq9ylfI6EBQAAgHAErGi3BB6vNLMlMV/7Dk1uGhE8X5raUquZnSbpdknHVh0elfT77r4x5vUzjSmCAAAAiIOAFe1GSdX7NZkmtx6fxMzmSLomcPgb7h5cK1XrtevMzAMfJ2W01sWSviPp+KrDRUmXuvuGyIo7RLHkYgALAAAAUQhYEdz9gKQPBQ6/xsw+Wu81ZnaUpK9Kmld9KknvT7/Cqgu0uFYzmy/pDknzA699l7vfFLfuTlByKXlTRgAAAEw3hXYX0CE+K+mdkl5RdexvzOxlkj4iaYu7D5vZ0SqvY7pKkzv4fcLdf9Rltf67Jrd2v1tSn5n9SQO173f3LzbwuqZzd9ZgAQAAIBIBKwZ3HxnbOPduSSdWferCsY+Sme2XdEydU2yU9FfNrbKsxbWeXuPYwNhHI3ZKymTAKrmYIggAAIBITBGMyd0fl/QqSffW+HRO9QPL5yStcveRZtUW1Em1dopym3YSFgAAAMIRsBJw9ydUDi7vlvRI2FMl3SnpAndf4+5DrahvQgEdVGsncHflCVgAAACIYO7e7ho6lpm9WNKZkhZKminpgKQdku5x9z1hr221Tqo1jJl5O35m7/7pM/rMXT/TF9ec2/JrAwAAIHvMTO4+6R141mBNgbv/RNJP2l1HHJ1UaxYV2QcLAAAAMTBFEIiBjYYBAAAQBwELiMGdjYYBAAAQjYAFxFAqiREsAAAARCJgATEUnTbtAAAAiEbAAmJwd+X5bQEAAEAEbhmBGErOFEEAAABEI2ABMRRLdBEEAABANAIWEEPJXTnaCAIAACACAQuIwV20aQcAAEAkAhYQAxsNAwAAIA4CFhBDseQiXwEAACAKAQuIwV3Kk7AAAAAQgYAFxMAUQQAAAMRBwAJiKLorx28LAAAAInDLCMTARsMAAACIg4AFxOBMEQQAAEAMBCwghmLJ2QcLAAAAkQhYQAwll3IkLAAAAEQgYAExMEUQAAAAcRCwgBjKbdrbXQUAAACyjoAFxFAs0UUQAAAA0QhYQAwld9ZgAQAAIBIBC4jBmSIIAACAGAhYQAxMEQQAAEAcBCwghhJdBAEAABADAQuIgTbtAAAAiIOABcRQ8v/f3r1H63LW9QH/PueSCwkhEAoh3BIuAQTBEyBcU4i6Vm2BGpBVbRsUKaAFVOiyulosqyyLoItaoIBWQCurWrxAQLxgVcAVjQFTEsESUzSQAA33QELI9exf/5g38mb2e5m9z7vPvszns9asnJk9zzyzc+Z33vc788xM3IMFAMBSAhYMcNhTBAEAGEDAggHcgwUAwBACFgxQhggCADCAgAUDHF5zBQsAgOUELBhgzT1YAAAMIGDBAIYIAgAwhIAFAxgiCADAEAIWDNA9RXC79wIAgJ1OwIIB1iruwQIAYCkBCwYo78ECAGAAAQsG6O7B2u69AABgpxOwYIC1iitYAAAsJWDBAIYIAgAwhIAFAxyuyj7VAgDAEr4ywgCGCAIAMISABQOsGSIIAMAAAhYM4B4sAACGELBgAI9pBwBgCAELBlirZJ+EBQDAEgIWDGCIIAAAQwhYMIAhggAADCFgwQCGCAIAMISABQN4TDsAAEMIWDBAVQwRBABgKQELBujuwZKwAABYTMCCAQwRBABgCAELBjBEEACAIQQsGOBwlacIAgCwlIAFAxgiCADAEAIWDLBmiCAAAAMIWDDAmqcIAgAwgIAFA6xVZb9LWAAALCFgwQBrlbiABQDAMgIWDFAecgEAwAACFgxw2D1YAAAMIGDBAN09WNu9FwAA7HS+MsIAVUlzBQsAgCUELBjgsHuwAAAYQMCCAdaqsl/AAgBgCQELBlhb85h2AACWE7BggKrKPi8aBgBgCQELBujuwdruvQAAYKcTsGCAtYp7sAAAWErAggGqymPaAQBYSsCCAQ6vGSIIAMByAhYMsFbJfgkLAIAlBCwYYM0QQQAABhCwYIA1QwQBABhAwIIB1irZ5woWAABLCFgwwFqVe7AAAFhKwIIBqhIXsAAAWEbAggEOVxkiCADAUgIWDLBWlf0CFgAASwhYsERVGSIIAMAgAhYssTYJV96DBQDAMgIWLGF4IAAAQwlYsMSaB1wAADCQgAVLuP8KAIChBCxY4vCaK1gAAAwjYMESa1XZv0/AAgBguVZV270PMFhrrY7WMXvVl2/IWy68Mhd85LO54ZbDOfHYAznv0Gl5wTkPyP1POeGo7AMAADtTay1Vte4svIDFrnK0AtYHrvhCXvQ/PpJbD6/ltrVv9ndgX8vB/fvy5vPPyrkPuceW7wcAADuTgMWecDQC1lVfviHf9boLc+Oth+euc/zB/XnfS89xJQsAYKTmBSz3YEHPWy68MrceXlu4zq2H1/LWCz95lPYIAIDdQsCCnndf+v/uMCxwltvWKhdc+tmjtEcAAOwWAhb03HDzbcPWu2XYegAAjIeABT0nHHtg2HrHDFsPAIDxELCg57xDp+XAkvdeHdjX8sxD9z5KewQAwG4hYEHPC855QA7uX1waB/fvy/PPOeMo7REAALuFgAU99z/lhLz5/LNy/MH9665kHdjXcvzB/Xnz+Wd5RDsAAOt4Dxa7ytF60XDSvQ/rrRd+Mhdc+tnccMttOeGYA3nmoXvn+eecIVwBAIycFw2zJxzNgAUAAPN40TAAAMAWE7AAAABWRMACAABYEQELAABgRQQsAACAFRGwAAAAVkTAAgAAWBEBCwAAYEUELAAAgBURsAAAAFZEwIIBPvjBD273LsCOpDZgPXUB642pLgQsGGBM/yjARqgNWE9dwHpjqgsBCwAAYEUELAAAgBVpVbXd+wCDtdYcsAAA7AhV1frLBCwAAIAVMUQQAABgRQ5s9w7ATtVae1CSRye5T5Ljk1yf5MokF1fVF7dz32C7tNYOJjk7ybckuftk8ZeSfDzJh6vq1i3o86QkT0zyoCR3SXJLkmuSXFZVf73q/tidJsfmmemOzXsmOSnJTUm+kuSqdMfnDVvUr5pgR2qt3SVdXdwvyb2SnJDu+/91Sa5Nd5x+bNXH6djrwhBBmNJa25fkB5O8LMnD56y2luRPkry6qj5wtPYNbtdau0+SxyR57OS/j0lyt95q51bVB1fY56lJ/n2SH0j3xXWWryV5e5JXVdXnV9DnWUl+KsnTkxycs9qnkrwxyRur6uYj7ZPdpbV2KMkzknx7kscnOXbB6oeTXJTkTUneWVW3HWHfaoIdp7X2sHQ18aQkh5Lcd0Czm5P8bpK3VNUfHmH/6iICFvy91tq9k/xWkidsoNkvJfmRqrpla/YKOq21V+SboerUAU1WFrBaa+cl+e/pzggO8dUkz62q92yyv5bklUlenuFD2S9P8j1Vdflm+mR3aa09K8nPpjtTvRl/leT7q+qjm+xfTbAjtdbemOTFR7CJP0jyQ1X16U30rS4m3IMF+ftw9WeZHa4q3aX0WV6Y5F2TS+GwlV6Z7qzkkHC1Mq2185O8K7M/MG9Kd+az7+R0dfEvN9ntf0vyHzL7M+r6dFci+h6W5MLW2iM22Se7y9lZHq6+muTGOT97VJIPt9a+c6Mdqwl2sdvSDQv8WrrROLP84yQXtdYevJENq4s7ErAYvUk4emeS03s/el+SpyY5vqruku4fjeck+Zveek9L8pqt3Us4+lprj0nytiTTj6C9McmrkzwgyZ3S3Z/4wHQ1cNPUevuS/PJkGxvp80eTvKC3+NNJXpTk7lV10qTPs5O8o7feKUneM7nngHG5IcmvJ3l+uvtNjq2qu1bVnZKcluR56c5cTzs2ybtba48c2omaYBf5SrrA87IkT0lyr6o6WFV3q6qTkxyX5NvSBZRrem3vk+T3WmvHDelIXcxQVSbTqKd0l9KrN/3cgvXvnOT9vfUPJ3n4dv8upr079Y63m5P873Rn756fbtx5/xh+6gr6/HBvm9cmeeKC9Z+U7qrBdJuLN9DfPZJ8vdf+o+m+GMxr85IZv/vPbvffl2lrp3Rf0irJFUl+OMlJA9ock+QtM46XCzfQr5ow7egpyT9Jcm6S/Rtoc1K6e7D6x83LB7ZXF/3+tvtAMJm2c0p3VuWaXsH90YB2d0/y+V67d2/372Pau1OSX053Zu6xSY7p/ez0GR8cTz3C/p41Y5vnD2j3/TPanTewzzf02n0jyYMH/r/ptzttu//OTFs3TWrhX23kS+SkXUvy3hnH6JMHtFUTpj07pbui+4necXPFgHbqYsbkIReMWmvtxemeKnO7SvKwqrpiQNsXpruCMO3hVfXxFe4iLNVaOz3JJ3uLj+ghF621v0z3UI3bXVRVTxrY9qLc8X7GD1fV45a0uWuSz6W7ynC7n6mqlw/o765JPpPuhMntfq6qfnLI/jIurbUHpvsiOT2c6bVV9W+XtFMT7GmttZck+a+9xadVVX8I4XQbdTGDe7AYu2f15j8wJFxN/Hq6S9SLtge7TmvtfrnjB2aS/OIGNtE/8XD25NHyizw9d/zAXEv3lM6lquraJL/ZW6wWmamq/i7dUwSnHVrURk0wEn8xY9lp81ZWF/MJWIzW5EzGP+wt/o2h7avq6+nGLE8770j3C3aA7+7N35buZumh3jlpM21ZbfR//qGqumoDffZvYn6Qp6exwN/25pc9nVNNMAbXz1i2f8H66mIOAYsxe0q6t5lP+/MNbqO//qNbaydvfpdgR/iO3vxlVXXD0MaTkw/9KwTLHofd73OjtXhx1j92eMOP4GY0junNL7tfQk0wBveesezqBeurizkELMbs4b3565Js9P6pi2cse9jmdgd2jH5tzDrOl+m3mVsXk/fQ9R+Xu6E+q+prWf8YbrXIPA/szc+9x2RCTTAG/StSH6+qzy1YX13MIWAxZv2C+lRt/KkvVw7YLuwarbVjk5zRWzzrOF+m3+aBrbX+VYPbzaqZVfSpFlmntXZm1n8xvGTB+mqCPW/yHqp/3Vv8ugXrq4sFBCzG7CG9+UWXwWeqqq+ke8Hlou3CbvKgrB9zv+HamNFmf9ZfNbjdrJpZRZ9qkVn6XyKT5HcWrK8m2LNaawdaaz+Y5I9zx6Gzf5jkrQuaqosF+vefwJjcrTe/bIjIPNek+4dm3nZhN5l1/G6mNma1mVcb/eW3VtWXV9DnXTexDfaw1trD0r1Da9plVbVomJGaYFdrrZ2f5MSpRQeSnJzuys13pntx77QL0r3LatGoHnWxgIDFmJ3Ym//GJrdz45Ltwm4y6/jdTG3062LetmctX1UtHmytHVtVN29ye+whk2FHb8/6B1z8xJKmaoLd7jWZ/QCLvovTvRfqggHrqosFDBFkzE7ozd+0ye0IWOwl/bpINlcbG/nQ3KpaXNQn4/O6rH9nz69W1R8taacmGIOPp3tH1J8OXF9dLCBgMWbH9+Zv2eR2+mc8+tuF3WTW8buZ2ph1JnBebWxVLS7qkxFprb046++9+kSSHx3QXE0wBt+S5OeTXN1a+8nW2rKMoC4WELAYs/6Zj3lPrVnm2CXbhd1k1vG7mdro18W8bc9avqpaXNQnI9Fa++4kr+8t/mqSZ1bVdQM2oSbY1arqPlXVbp/SXQm6b5KnJXlDunq43QnphhT+Rmtt0a1E6mIBAYsx+3pv/rhNbqd/1qO/XdhNZh2/m6mNWWcD59XGVtXioj4ZgdbaU5O8I3d82tkNSZ5WVf9n4GbUBHtKVX2jqj5TVb9fVT+W7ql9v9Vb7dlJXrVgM+piAQGLMesX0502uZ1+oV6/ye3ATjDrQ2YztTHrA2xebfT73OxQjX6726rK2fqRaq09Nt3j16e/hN2c5LyqumgDm1IT7GmTV858X9aHrB9vrT14TjN1sYCAxZhd25s/dZPb6bfrbxd2k1nH72Zq414Dtz1r+TGttc287qDfp1ocqdbaI5K8L8mdpxbfluR7q+qPN7g5NcGeV1Vr6V5hMB1u9iX5oTlN1MUCAhZjdkVv/v4b3cCksPtPnulvF3aTv02y1lu24dpIcr/e/OEkfzdn3Vk1s4o+1eIItdbOTPfS1OkvXmtJnlNV79nEJtUEo1BVX0rye73F3z5ndXWxgIDFmF3emz+9tdY2uI0zBmwXdo3JMIlP9hY/YBOb6tfGlQveMTKrZlbRp1ocmdba/dOFq3tOLa4kL6yqd2xmm2qCkflIb37W9xx1sYSAxZj1b3A+Kd1bzTfi8TOW+QBjt+vXxuM2sY1+bcyti6r6TJL+09w21Gdr7aR0jxke1Cd7T2vttCTvT/d0tGkvraq3HeHm1QRj8bXe/KL7qtTFHAIWY/an6S5FT3vyBrfRX//SqjLGnd3u/b35Q621WS+VnGmy7rf1Fv/JBvvcaC0+Ies/05b1yR7RWvsH6a5c9c9m/7uqesMKulATjMUpvfkvLlhXXcwhYDFak6fmXNhb/L1D27fWTkzy9N7iC450v2AHeHdv/kCSZ22g/bMnbRZtc1mfj5sM9xrq+3rzV1bVRzfQnl2qtXZykv+V9SMQfrqqXrOibtQEY/HY3vxnF6yrLuYQsBi7d/Xmz22tPWRg23+R9Q+46G8Pdp2quirrx+H/8AY20X/q1CVVdfWSNu9NcuvU/L4kLxzSWWvtrkn+WW+xWhyByYmuP8j6s+D/uapesap+1ARj0Fq7V5J/1Fs896mb6mI+AYux+5UkX5iab0neuKxRa+2UJD/dW/y7G3hxJex0/TP/T2ytnb+sUWvtOemGYCza1jqTK8pv6S1+2YJ3sEx7be54n8BNSV4/oB27WGvtuHTvuerfw/GmqvrxLehSTbCjtdaOPYK2+5P8Ytbfc/XbS5qqi1mqymQa9ZTkR9I9ZWp6es2C9e+cbrzu9PprSR653b+LaZxTktNnHMNPPcJttiSX9LZ5bZInLGjzpCRf7bX5yw30eWqSG3rt/yrJqQvavGTG7/7a7f47MW3tlG5Y0Xtn/N2/NUnboj7VhGlHT+kCxK8kefAG252a5D0zjpv/OaCtupgxtUmnMFqttYNJLkrymN6P3pfk1UkurqpbJk+eeXqSn8r6sf6vr6qXbvnOMlqToavnzvnxKUn+U2/Zf0nyf+esf0lVXTKgz7OT/FmSg1OLb0zy8+nOIN4+lOP+6YZovCzJcVPr3prkyVX14WV9TfX50sm+T/t0klcl+c2qura1diDdkLB/k+Sf99b9VJJDVfXVoX2y+7TW3pTupajTrklXB/138wz1a1V1/aIV1AQ7WWvtdUl+bDJ7Sbr7lS5JFz4+X1Nf+ltrpyY5lOS8dMfMne+4tVyd5PFVdc2AftVFf/8ELEhaa/dN949D/+VzSfdhfX2Su8xp/r4k/7Sqbp3zczhirbXnpjszuQqvrKr/OLDfH5j0O+sdcTdN/nvcjJ9VkudW1ds3unOttbcled6cH1+fbojH/hk/uzbdlTs38u9xrbUPJnnKijd7RlV9akDfaoIdqRew+irdI84Pp3stTf/hEtOuSnJuVfXfc7Wob3UxxT1YkKSqPp3ukvWHZvx4X+aHq7clOU+4Yq+qql9Nd1Nw/90jSfdhOesD87okz97MB+bEC5L8TGZfibhzZn9gXpHkHF8k2Wpqgl2qpfsuc7fMD1drSX4hyaM2Eq4SddEnYMFEdS+we1K6y9d/s2jVdO9h+I6qen7Nf+M47AlV9dtJHprkTenOCs5zXbqHxDykqjb9ZKaqWquqlyc5O919AbctWP3qJD+R7guBh8xwVKgJdqhXJPmedMPyLk/3fWWIT6QbXntmVb2oqvovGx5EXXyTIYIwR2vtzCSPTnKfdGdevp7kyiR/UVVfWNQW9qrW2jHpPswenm++kPJLST6e5ENbcTV38p6jJyR5cLqhLbemu9/mMmfn2W5qgp1qcu/4Q9M9COmeSU5Id1Xnusl0dZJLt+I+pLHXhYAFAACwIoYIAgAArIiABQAAsCICFgAAwIoIWAAAACsiYAEAAKyIgAUAALAiAhYAAMCKCFgAAAArImABAACsiIAFAACwIgIWAADAighYAAAAKyJgAQAArIiABQAAsCICFgCMWGvtl1prNZlubq0d2O59AtjNBCwAGLdvnfrz5VV127btCcAeIGABwEi11lqSR0wt+th27QvAXiFgAcB4nZ7kxKl5AQvgCAlYADBe39qb/+tt2QuAPaRV1XbvAwAAwJ7gChYAAMCKCFgAAAArImABAACsiIAFACPUWnvG1AuGq7X2tO3eJ4C9QMACgHF6dG/+I9uyFwB7jIAFAON01tSfP1dV12zbngDsIQIWAIzTdMBy9QpgRQQsABiZ1to9ktx7apGABbAiAhYAjM9ZvXkBC2BFBCwAGB8BC2CLCFgAMD7TAevLVXXVtu0JwB4jYAHA+EwHrEu3bS8A9iABCwBGpLV2cpIzphYZHgiwQgIWAIxL//4rV7AAVkjAAoBx8YALgC0kYAHAuEwHrOuTfGK7dgRgLxKwAGBcpgPWZVVV27YnAHuQgAUAI9FaOzHJmVOLDA8EWDEBCwDG41CSNjUvYAGsmIAFAOPhARcAW0zAAoDxmA5YNya5fLt2BGCvErAAYDymA9ZHq+rwtu0JwB4lYAHACLTWjkvy0KlFhgcCbAEBCwDG4VFJDkzNX7pdOwKwlwlYADAOHnABcBQIWAAwDmdP/fmWJB/brh0B2MsELAAYh3Om/vyRqrpl2/YEYA8TsABgj2utPTLJA6cWfWC79gVgrxOwAGDve3lv/ne2ZS8ARqBV1XbvAwCwQq21+yY5Pt1Vq+clefbUjz9aVY/alh0DGIEDy1cBAHaL1tqJSa5K0mb8+MYkLzq6ewQwLoYIAsDe8qjMDldXJXlGVf35Ud4fgFFxBQsA9pZDk/9en+SLSS5L8vtJfq2qbtq2vQIYCfdgAQAArIghggAAACsiYAEAAKyIgAUAALAiAhYAAMCKCFgAAAArImABAACsiIAFAACwIgIWAADAivx/8sEFEj3nFqQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot the principal components\n",
    "#colors = 200*['green', 'purple', 'orange', 'teal', 'black', 'red']\n",
    "for i in range(0,1):\n",
    "    p_data = -components[i]\n",
    "    #plt.plot(range(len(p_data)), p_data, marker='o', linewidth = 1.0, markersize=5, color = colors[i])\n",
    "    plt.plot(range(len(p_data)), p_data, marker='o', linewidth = 1.0, markersize=10)\n",
    "    #p_data = components_ig[i]\n",
    "    #plt.plot(range(len(p_data)), p_data, marker='^', linewidth = 1.0, markersize=5, color = colors[i])  \n",
    "  \n",
    "\n",
    "#backround grid details\n",
    "axes = plt.gca()\n",
    "#axes.grid(b = True, which = 'both', axis = 'both', color = 'gray', linestyle = '-', alpha = 0.5, linewidth = 0.5) \n",
    "#axes.set_axis_bgcolor('white')  \n",
    "\n",
    "#font scpecifications\n",
    "title_font = {'family' : 'arial', 'color'  : 'black', 'weight' : 'heavy','size': 20}\n",
    "axis_label_font = {'family' : 'arial', 'color'  : 'black', 'weight' : 'normal','size': 20}                                                   \n",
    "\n",
    "#figure size and tick style\n",
    "plt.rcParams[\"figure.figsize\"] = [12,8]\n",
    "plt.rc('axes',edgecolor='black',linewidth=1)\n",
    "plt.tick_params(which='both', axis='both', color='black', length=4, width=0.5)\n",
    "plt.rcParams['xtick.direction'] = 'in'\n",
    "plt.rcParams['ytick.direction'] = 'in'\n",
    "#plt.xlim([0,40])\n",
    "\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r'$j$', y=3, fontsize=40, fontdict = axis_label_font)\n",
    "plt.ylabel(r'$\\mathbf{w_1^{(0)}}$', fontsize=40, fontdict = axis_label_font)\n",
    "\n",
    "#title and axis labels\n",
    "plt.tick_params(axis='both', labelsize=45)\n",
    "#plt.title('PCA feature weights', y=1.05, fontdict = title_font)\n",
    "\n",
    "#save figure\n",
    "fig=plt.figure(1)\n",
    "plt.tight_layout()\n",
    "plt.savefig('figures/OP_orientation_kappa_4_weights.svg',dpi=fig.dpi)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Order parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xl8VPXV+PHPSUiAhH1HwEaBgIKSsLsrtkA07qjYqkXx59oqte3TVlsXtOvjggu2PtYNRbQudUEjWkUUhGg1uLGDEQhBdkjCkmTm/P64MzAZEjKZTHLvzJy3r3lN5t4733uSwZx8d1FVjDHGGK9JcTsAY4wxpjaWoIwxxniSJShjjDGeZAnKGGOMJ1mCMsYY40mWoIwxxniSJShjPEpEikVEA4+nQo5nhRxXEZkUcm5S2Lms5o88OiLyVEjcxQ14X9x+z+bQLEEZz6nlF3Doo1JESkTkDRG5sJ5y+ovI30RkkYhsFpEqESkTka9E5EEROUlEpI73nlfLve9tmu/YGFObFm4HYEwDpQGHBR75IvIKcLGqVgcvEJF04M/AL4DwBNQGGBx4/Bw4Aiiu5T6Tazl2qYj8VlWrGvtNNKFPgV+HvN7mViBReB74OvD1TjcDMd5gCcrEg/8CL+AkmyzgMqBt4Nz5wLXAwwAikgLMBCaEvH8f8BrOL78UoB8wDuha281E5DBgfC2nugFnAa805ptpSqr6DfCN23FEQ1XfBt52Ow7jHdbEZ+LBN6p6j6r+r6reQM3kAxDa1Dcp7PwqYLCqXqyqd6nqnap6GdAbuAbYXcv9JgGpga+rgTUh566M9psQkbEi8qKIrBORfSKyS0Q+EZFfiUhGtOWG3aPO/pjwPh4RaSci94rId4F4VovIHSLSspZyJ4rIOyLyfUhTabGIvC0id4lIj1rec7yIzBSRb0Vkr4hUiMiSQPPqkbVcf8g+KBEZEbjfrsD93xORU+r5ebQSkd+ISKGI7BCRahHZJiLLA5/F/0T4ozVuUFV72MNTD5xakoY8ngo7nxl2fkXIuW/Czg1v4L0FJ6kF318A/C7kdTVwWBRl/l9YXOGPL4FuYe8rru1nUMvPZ1LIuUlh57JCzj0VcnwT8EUdsbwDpIa875Z6Ylfg1LDYpwL+Q1xfDpwT9p7Q+IrDzv0IpyYcXo4PmH2I7/md+mJ3+9+7Pep+WBOfiUcnhL0uBRCRnsDRIce/UNX/NrDsU4C+Ia9nAvOBP+IkmlTgpzh9XJH6JfD/Ql4XAAtxmhgvB9oDxwDPAmMbGG+0ugIdgMeALcBFHPi+fwTcADwYeH1jyPv+i5MQAPoAxwLDQgsWkYuBP4QcKsZpos0ArsDpB8wEnheRwaq6+lCBikgrYAaQHjikOP1Vq4B84Mw63jcw8L0E/TsQf1ucGvRx1Pysjde4nSHtYY/wBwfXED4FfoXT+f8wsCvs/M8C7xsRdvz5KO79TMj7K4A2geMLQo6vbEB5KTi1leB7Hwk7nxcWc07IueKQ408d4uczKeTcpLBzWSHnngo7d1nIuU7A9pBzS0LO7Qg5PrqW77EL0D7k9Wch1+8AuoScOy0shgfqiK845PjFYe+5M+RcS2BJbd8zkBNybCeQXkvs/dz+926Puh9WgzLxYHjgUZvXgH/E4iYi0h64IOTQ66paHvj6OeD4wNf9ROQUVZ0XQbEDqDkY4zoRue4Q158ILI405kaowvmeAFDVbSLyBs4AFICjRCRTVSuAecDZgePvikghTu1lBbAIWKSqfoBAX1puyH3eUNUtIfeZG+hfygocCq8N12ZE2OunQ8rbJyLPA3fW8r6lwGacn387oFhEPgvEvhT4UFWXRXB/4xJLUCbeVAFbgc9xajsvaOBPYWB92LUDG1j2JUDrkNfPhXz9L2AaB/6fuRLnF3d9OjUwhlpHFjaBrarqCzv2fdjrDji1yGtxmiFPwWmeOz3wCFolIvmquhzoSM2h/RtrufdGDiSoSH4+HeqJM/w1sD95XYBTMzsS6InTJLifiLyL0xe2J4I4TDOzBGXiwdOqOqm+i1S1VESWcKAfaoiIDFPVzyK8T/jcp9frmMcLMEFEfq6qu+opM3we0ktA4SGuX1hPebHSWURSw5JU97BrdoDzcwVOFZHDgZFAf5ya4Tk4yaMf8HdgDAeaCYM/uING94Udi2Se1o6w192pObIyPO79VPUjEemH01c2BKfP6Vic6QKpOH1Uv8YZ1GE8xhKUSTT3Af8Mef2ciIxX1W9DLxKRNJwO+1dVdZOIHEvdzYi1yQAm4ozOO5TlOIMQugRedwKmacjE4kA8rYGLVHVBA2JojDTgxzi1UESkE84v7aClgeY9RGQI8LWqrgXWhsR8MxBcXWMEgKruFpHFHGjmO0tEugSb+UTkNA7UnsDp26vPp2GvfwrcHigvHedzOEjgXLaqfo0zYvGLkHOvh3y/4U2IxiMsQZlE8yTOqK7zAq+zgSUiEjpRtz8HJuq+E7guvPb0BrXPkcrHGYEWfM8hE5Sq+kXkf4G/Bg6NAb4Skdk4TZWdcP6iPxmnefHpWgtqGo+LyIk4CfRiajalhX5fM4HuIvI+UMKBfp3LQq4JrQn9LweaR9sDnwb6iTKoOY9sH/BQBHG+jtOMF6wp/SFQK1qN83kcVcf72uH8rFcCH+OM9tyF8/mfUUfsxkvcHqVhD3uEP6hnHlQE72+JM0T6UPNw9o/4Cly/NeTYskOU/VjY+wdFEI/U8r565+TQtKP4NgKf1BHHe0CLkPd9HUHs14bF/qd6rt8NnB/2ntD4isPOjaf2eVB+YG5t3zNOrbW+uCuAoW7/m7dH7Q9bScIkHFXdp6o34vxlfQ/OL+KtOJNsK3B+4T6MU2v5DjiXmp31jx2i+PBzta3ZFx6Pqur/A36IM3+nGOeXbRXOX/XvA3fg1KSay16cId9/w/kZVAXiugs4U2s2Qd4KTMf5OZbgxF4JrMNZ9mmcqtYYSamqt+D8fGfhNAtWAntwmjynA8eqasRLRqmzDNLJODXecpzP8UOcYfp11Tp3AtfjNGN+hTPcvxonOS7HqSUOU9XPI43DNC8J/KVhjElw4mzZ8dPAy+9UNcu9aIypn9WgjDHGeJIlKGOMMZ5kCcoYY4wnWR+UMcYYT7IalDHGGE+yBGWMMcaTLEEZY4zxJEtQxhhjPMkSlDHGGE+yBGWMMcaTLEEZY4zxJEtQxhhjPMkSlDHGGE+yBGWMMcaTLEEZY4zxJEtQxhhjPMkSlDHGGE+yBGWMMcaTLEEZY4zxJEtQxhhjPMkSlDHGGE9q4XYAXiUittWwMcY0A1WV2o5bgjoE1Zo5qri4mKysLHeCMQ1mn1f8sc8s/jT2MxOpNTcB1sRnjDHGoyxBGWOM8SRLUMYYYzzJEpQxxhhPsgRljDHGkyxBGWOM8SRLUMYYYzzJkwlKRCaIyEMi8pGI7BIRFZFnoyyrt4g8ISIbRGSfiBSLyDQR6RjruI0xxsSOVyfq/h4YApQD64GB0RQiIn2Bj4FuwGvAMmAkcBMwXkROUNWtMYnYGGNMTHmyBgX8AsgG2gHXNaKcR3CS042qeq6q/lZVxwD3AwOAPzY6UmOMMU3CkwlKVeeq6koNX2uoAQK1p7FAMTA97PTtQAVwmYhkRh2ocZ3P72P2itncNe8uZq+Yjc/vczskY0yMeLWJLxZOCzy/o6r+0BOqWiYiC3AS2GjgveYOzjSez+9j3LPjKCwppKKygsz0TEb1GsWcS+eQmpLqdnjGmEbyZA0qRgYEnlfUcX5l4Dm7GWIxTaBgVQGFJYWUV5ajKOWV5RSWFFKwqsDt0IwxMZDINaj2geeddZwPHu9QVwFTpkzZ//Xo0aNtlWWPeX/p+1RUVtQ4VlFZwdylcxmcPpiNGze6FJmJln1m8acpP7NETlCNNm3atBqvbSsAbxlTOYaHv3yYKn/V/mMiQk5WDl9Xfs3737/PmE5jyOuXZ01+ccT+H4s/TfWZJXKCCtaQ2tdxPnh8RzPEYppA/079UZTMtEx2V+0mMy2Tzq07c83saxAR9lTt4bElj1m/lDFxKpET1PLAc119TP0Dz3X1URmPu2PeHfzh5D8wtOdQFm9cTE6PHHx+Hxe/dDH7qvcB1OiXys/OdzliY0xDJHKCmht4HisiKaEj+USkLXACsBtY5EZwpnE+L/2cecXz+OdZ/yQzPXN/8rlr3l1U+iprXFtRWcHijYstQRkTZ+J+FJ+IpInIwMC8p/1UdTXwDpAF3BD2tjuBTOAZVa3AxJ3fvfc7fn/y78lMrzmNLbdn7kHHWqe1JqdHTnOGZ4yJAU/WoETkXODcwMsegefjROSpwNdbVPVXga97AUuB73CSUajrcZY6elBETg9cNwpnjtQK4NamiN80rffWvMea7Wv4f0P/30Hn8vrlMarXqP1zo9JT0/H5fQzsHNVqWcYYF3kyQQE5wE/Djh0ZeICTjH5FPVR1tYgMB6YC44EzgFLgAeBOVd0es4hNk/P5fby18i2uffNafnLMT0iRgxsAUlNSmXPpHApWFTB36VxOO+o01mxfw/iZ4/ngpx+w+PvFFJUWkdsz10b3GeNxnkxQqnoHcEeE1xYDcojz64ArYhGXcU9w1YgFaxew17eXv3/6dz4v/bzW0XmpKankZ+czOH3w/uGvO/fuZOB0pxa1u2q3rTphTByI+z4okxyCq0bs9e0FoLyqYatG5PTIodJXSUVVha06YUycsARl4kJRaVGtq0Ys3rg4ovcv3riYan911O83xjQ/S1AmLuT2zKVli5Y1jmWmZ0Y8Oq+20X0Neb8xpvlZgjJxIa9fHt0zu5Oeko4gtElvw6heo8jrlxfx+0f1GkWb9DZIoMsyp0dOxO83xjQ/Tw6SMCZcakoqfTv2ZVLOJFqktNifXCId4BA6um/xxsV8+f2X7KnaU+tIQGOMN1iCMnHBr34+K/2MFy58gS4ZXaIqIzi6Lz87n0pfJSMfG8mTi5/kytwrYxytMSYWLEGZuLBi6wo6te4UdXIKl56azjPnPcNpT5+GqrKhbIPNjTLGYyxBmbjwScknjOw1MqZlHt31aDq06sDVs69GVW1ulDEeYw3wJi40RYIqWFXA9+Xf41e/zY0yxoMsQZm40BQJqqi0iIqq6OdWGWOaliUo43n7qvfx9aavye2RG9NybW6UMd5mCcp43hfff0F25+yDkkljhc+NSpGUBs2tMsY0LRskYTyvKZr3oObcqM82fMb9i+7nobyHbICEMR5hCcp43icln3DyD05ukrJD50aVV5bz1OKn+OuP/tok9zLGNIw18RnPa6oaVLjJQyfz9BdPU+WravJ7GWPqZwnKeNqOvTtYv2s9R3c9usnvNbDLQPp16sebK99s8nsZY+pnCcp42n83/JehPYfSIqV5WqMn507m8aLHm+VexphDswRlPK25mveCLhx0IfPXzqdkV0mz3dMYUztLUMbTmjtBtUlvw4VHX8jTXzzdbPc0xtTOEpTxLFWlsKSwWRMUwFVDr+KJoifwq79Z72uMqcmGmRvPKikrwef38YP2P2jW+444bAStWrTiL/P/gs/vs1XOjXGJJSjjWcHmPRFp1vv61U+lr5LbP7gdn99nq5wb4xJr4jOe1dz9T0EFqwrYULaBan+1rXJujIssQRnP+nTDp64kqKLSInZX7a5xzFY5N6b5WYIynuPz+3h9+essWLuAbXu24fP7mvX+tsq5Md5gCcp4is/vY9yz47jk5UvY59vHNbOvYdyz45o1SYWvci4Ifdr1sVXOjWlmlqCMpxSsKqCwpHB/E5sb/T/BVc5nXTCLqadNZfoZ09myewsfrf2o2WIwxtgoPuMxRaVFVFTWvsttfnZ+s8URuso5QHbnbC5+6WI+nPQhA7oMaLY4jElmVoMynuLV/p/TjzydP435E/mz8tmye4ursRiTLCxBGU8J9v+kSiqC0Ca9jWd2uZ08dDIXHHUB571wHvuq97kdjjEJzxKU8ZRg/09WhyyuHnY1sy6Y5akJsn86/U90z+zOVW9chaq6HY4xCc0SlPGc1JRUyirLuP2U28nPzvdMcgJIkRRmnDeD5VuWc9eHd7kdjjEJzRKU8ZxKXyXb92ynW2Y3t0OpVUZaBq9f8jpPFD3Bc18953Y4xiQsS1DGczaUbaB7m+6eqjmF69GmB29c8gZT3p7CgrUL3A7HmIRkCcp4TsmuEnq36+12GPU6pvsxzDhvBhNenMDqbavdDseYhGMJynhOSVkJvdr2cjuMiIzvN57bTr6NM587k+17trsdjjEJxRKU8Zz1u9bHTYICuG7EdeT1y+OCf11Apa/S7XCMSRiWoIznxEsTX6h7xt5D25ZtuW72dTb83JgYsQRlPKekrIRe7eKnBgXO0PiZ58+kaGMRf1vwN7fDMSYh2Fp8xnPW71ofdzUogDbpbXjjkjcY/fhosjpmkZmWSVFpkW0Zb0yULEEZz4mnQRLherXrxb8v+jfHP3E8LVJasLd6r20Zb0yUrInPeIpf/Wwo28BhbQ9zO5SobazYSGpKKnuq99iW8cY0giUo4ylbdm+hbXpbWqe1djuUqBWVFh20mKxtGW9Mw1mCMp5Ssiv+BkiE8+qWIcbEG0tQxlPidYBEqNq2jO/fqb8ntgwxJp54dpCEiPQGpgLjgc5AKfAqcKeqRjxlX0ROBH4NDAF6AJuAr4EHVfXtWMdtGieeB0gEBbcMKVhVwOKNixGE6Z9OZ9uebXTN7Op2eMbEDU8mKBHpC3wMdANeA5YBI4GbgPEicoKqbo2gnOuAR4AK4N/AeqA3cD6QJyK/V9U/Ns13YaJRsiv+ExQcvGV8WWUZV75+Ja9PfB0RcTk6Y+KDV5v4HsFJTjeq6rmq+ltVHQPcDwwA6k0qIpIG/BnYCwxT1ctU9XeqehkwHNgH3CoiLZvsuzANtr4s/pv4ajP1tKlsLN/II58+4nYoxsQNzyWoQO1pLFAMTA87fTtObegyEcnk0DoB7YEVqro89ISqLgVWAK2BNjEI28RIIgySqE16ajrPnf8cd8y7g6++/8rtcIyJC55LUMBpged3VNUfekJVy4AFQAYwup5yNgGbgWwR6R96QkSygf7A4kiaCk3zibeFYhuif+f+3POje7jk5UvYU7XH7XCM8Twv9kENCDyvqOP8SpwaVjbwXl2FqKqKyA3As8BnIvJvYAPQCzgP+AaYeKhApkyZsv/r0aNHk5WVFdl3YKK2fud6/Dv8FO8pbnRZGzdubHxAMXZy+5N5OfNlrnn5GqaOnup2OJ7jxc/MHFpTfmZeTFDtA8876zgfPN6hvoJU9UUR2QDMAi4POfU98CSw5lDvnzZtWo3XxcXFlqSaUNm+Mqq1mmOzj43ZQAIvfl4zes4g99Fcvtz3JWcPONvtcDzHi5+ZObSm+sy82MQXMyJyKfAf4CPgKJymwaNwal4PA8+7F50JV1LmbLOR6KPcOrTqwMzzZ3L1G1ezoWyD2+EY41leTFDBGlL7Os4Hj+84VCGBfqYncJryLlPVZaq6R1WXAZcBnwEXisipjQ/ZxEKiDpCozfF9juf6Eddz+b8vx1+zq9UYE+DFBBUccZddx/nggIe6+qiCxgJpwLxaBlv4gQ8DL4dFE6SJvUQeIFGbW066hX2+ffzvgv9l9orZ3DXvLmavmI3P73M7NGM8wYt9UHMDz2NFJCU0uYhIW+AEYDewqJ5ygvOb6pq6Hzxue3R7RLCJL1m0SGnB0+c8zcDpA21rDmNq4bkalKquBt4BsoAbwk7fCWQCz6hqRfCgiAwUkYFh134UeJ4gIseGnhCRHGACoMD7sYveNEairCLREEu2LCFVbGsOY2rjuQQVcD3OPKYHReRVEfmziLwP/AKnae/WsOuXBh77qeonOCP1WgOfisjzIvJXEXkBKARaAQ+o6jdN/L2YCK0vW580fVBBRaVF7PPZ1hzG1MaLTXyo6moRGc6BxWLPwFks9gEatljsZJy+pknAOKAtsAuYDzymqjaKz0NKdiVXEx8c2JqjvLJ8/7HUlFR+0P4HLkZljDd4MkEBqOo64IoIr611XLKqKvBU4GE8LhFWMm+o4NYchSWFVFRWkJGWQeeMzkx5ewqbKjbx81E/Jz013e0wjXGFZxOUSS5Vviq27t5KjzY93A6lWYVvzZHTI4e8fnms2raKKXOm8M+if/LA+AcY23es26Ea0+wsQRlPKC0vpVtmt6QcuRa+NQfAgC4DeOvHbzF7xWyue/M6ju1+LPeOvZcftP8BBasKKCotIrdnLnn98pLyZ2aSgyUo4wnrdyXfAIn6iAhnDTiLH/X9EfctvI8R/zeCdq3asbliM7urdtuQdJPwvDqKzySZZBwgEalWLVpxy0m3cM/Ye1i/cz0VVRU2JN0kBUtQxhOScYBEQ63ftR6f1lxlwoakm0RmCcp4QrItcxSN4JD0UJnpmeT0yHEpImOaliUo4wnJtsxRNIJD0tukO5tAt0xtyaheo8jrl+dyZMY0DRskYTzBBknUL3RI+stLXmbed/NsgIRJaFaDMp5ggyQiExyS/vg5j+NXP19+/6XbIRnTZCxBGdepKhvKNlgfVAOkSApX5FzB40WPux2KMU3GEpRx3ZbdW8hIy6B1Wmu3Q4krV+RewayvZ7Gnao/boRjTJCxBGdfZAInoHN7+cIYfNpxXl73qdijGNAlLUMZ1ybTVe6xdmXOlNfOZhGUJyrjO5kBF79yB5/LF91/w7fZv3Q7FmJizBGVc5fP7mFs8l9XbVzN7xWx8fl/9bzL7tWzRkh8P/jFPLX7K7VCMiTlLUMY1Pr+Pcc+O4+WlL/NB8Qdc8vIljHt2nCWpBroy90qeXPyk/dxMwrEEZVxTsKqAwpJCqv3VALb4aZSG9BhCt8xu/GfNf9wOxZiYsgRlXFNUWkRFZUWNY7b4aXQm507micVPuB2GMTEVVYISkVQROV1EHhSRYhH5UkSmisjwWAdoElduz1wy0jJqHLPFT6NzyTGXMGfVHLbu3up2KMbETMQJSkTaishFIjIT2Ay8A1wPrA9c8nugUETWi8h0ERkrImmxD9kkirx+eRzZ8UhSJRVBaJPexhY/jVKHVh3Iz87n2S+fdTsUY2Km3sViReQ64BzgVCAd2AO8C7wGzFbVzYHrjgDODVx7NXAtUC4iBcAsVX2tKb4BE79SU1IZc8QYcnvk0r9zf3J65NgW5o0wOXcyN719EzeOuhERcTscYxotktXMp+PUmGbiJKV3VfWgtVVU9VvgfuB+EekE5OMkqzOB7MB7jalhwboF3Df2Pk76wUluhxL3Tsk6hfLKcj4r/Yzhh1lru4l/kSSok4CPVVUjLVRVtwEzgBki0hIYEmV8JoGVV5azZPMSRvQa4XYoCSFFUrgy90oe//xxS1AmIdTbB6WqCxqSnGp5/z5V/STa95vEtWj9InJ75NKqRSu3Q0kYk3Im8cI3L7C7arfboRjTaDbM3Ljmo+8+4qTDrWkvlnq3682o3qN4ZekrbodiTKNZgjKumb9uvvU9NYHJuZNtAVmTECIZxXd4DO+3Q1V3xbA8E6eqfFV8UvIJx/c53u1QEs7ZA87m+jevZ/W21fTt1NftcIyJWiSDJIqBqPugwtwJTI1RWSaOfV76OUd2PJIOrTq4HUrCSU9N5yfH/IQnFz/J3WPudjscY6IWSYKaSuwS1LwYlWPi3Py18zmxz4luh5GwJg+dzPhnx3PnqXfavDITt+pNUKp6RzPEYZLMR2s/YuLgiW6HkbAGdxtMzzY9uevDu0iVVHJ75tokaBN3IqlBGRNTqsr8tfN5+IyH3Q4lYfn8Psoqy7j7w7vxq5/M9ExG9RrFnEvnWJIycaNRo/hE5GURyQh8nR2bkEyiW7ZlGW1btqV3u95uh5KwClYVsH7XenzqQ1HbysTEpcYOMy8FggvCLhORMhH5WEQeEZGrRWSkiNgsTFPDR2tt/lNTKyotOmiyrm1lYuJNo5r4VPVnIS/74CxplBN4/iXQF1ARWamqRzfmXiZxzF87nxMPtwESTSm3Zy6Z6ZmUV5bvP2ZbmZh4E7M+KFUtAUqAt4LHAs1/Q4BjY3UfE/8+WvsRvzvxd26HkdDy+uUxqtcoCksKKa8sJ0VSGNlrpG1lYuJKkw6SUNXdwMLAwxjW71pP2b4yBnYZ6HYoCS01JZU5l86hYFUBRaVFPPPlM1x6zKU2QMLElYZsWHiHiJSKSKWIrBSR39uGhKahgs17tl9R00tNSSU/O58/nPIHZp4/k1vev4Ude3e4HZYxEYsoQYnIlcBtQHecWldfnFUhXmq60EwisgVi3TGi1wjOzj6b2+be5nYoxkQs0hrUtUAl8BOgN/BD4HMgX0QubKLYTAKav84GSLjlT6f/iRe+ecFG8pm4EWmC6gu8pKqzVHWDqr4P/AjYDlzeZNGZhLJj7w7WbF/D0J5D3Q4lKXXO6Mzdp93Nz976GX71ux2OMfWKNEF1BFaFHlDVHcCbgP22MRH5eN3HjDhsBGmp1nXplslDJ1Ppq+SZL55xOxRj6tWQibq1/cm1Fugco1hMgrP+J/elSArTz5jOb9/7rQ2YMJ7X2JUkqjmwkoQxh2QbFHrDiF4jOGfAOTZgwnheQxLUH0TkKxF5PLCMUQ622KyJ0N7qvRSVFjG692i3QzHAH8f80QZMGM+LNEH9B9gJDAKuAP4OfAbcAiAi94jIpSIyWERiso28iPQWkSdEZIOI7BORYhGZJiIdoyhrqIg8JyLrA2V9LyLzRMQGeDSTT0s+5aiuR9EmvY3boRgODJi44a0bbMCE8ayIkomqjlXVLjij+S4G7gE+AHYBAtwMPA18AZSLyCci8mi0QYlIX5wEeAXwCXA/sAa4CVgoIhH3e4nIz4BPgbHAe8C9wL+BVOCMaGM0DWMbFHrP5KGTqfZX24AJ41kNaqJT1W+Bb4EXg8dEpD8wPOSRG3geBlwTZVyPAN2AG1X1oZB73Qf8AvgjztysQxKRscCDwLvABFUtCztv/WfN5KO1H3HV0KvcDsOECA6YOGvWWZwz8Bw6tOrgdkjG1CCqsdrNPVCgs4bNQGC4qjb4T7NA7WkVUAz0VT3Q/iAibXG2+BCgm6pW1FPWF0A/4HBV3drAODT8Z1NcXExWVlZDijE4m+fn+E/RAAAgAElEQVR1/ltnVvx8Bd0yuzXbfe3zisy1s68lVVLJ659HUWmRq7vv2mcWfxr7mYkIqlrr2mcxH+QQ+K2+NPCIxmmB53dCk1Og7DIRWYDTXDcap8muViIyGGcV9VeBbSJyGk6tToHFwNzw8k3T+HrT13Rv071Zk5OJ3NRTp9L7/t48ufhJ9lbvtd13jWfUm6BE5PAY3m+Hqu6q55oBgecVdZxfiZOgsjlEggJGBJ434fSXnRx2/isROV9VV2GalG1Q6G2fbPiEFElhT/UegBq77+Zn57scnUlmkdSginFqHbFwJzC1nmvaB5531nE+eLy+BvPgn+uTcfapOhOYj7Pg7W3ApcCbInKMqlbWVsCUKVP2fz169GhreojSnKVzOKXXKRQXFzfrfTdu3Nis94tX7y99n0pfzf8FKiormLt0LoPTBzdrLPaZxZ+m/MwiSVBTiV2CmhejciIRHKGYCkxU1eCeVLsCw8sH4gzmuACYVVsB06ZNq/Ha2scbTlX5fOvnTDtrGlmdspr9/vZ51W9M5RgeW/LYQbvvnnbUaa78/Owziz9N9ZnVm6BU9Y4muXPdgjWk9nWcDx6vb52W4PmNIckJcPrJROQ1nAQ1kjoSlGm84h3F+NXPkR2PdDsUU4fQ3XcrKp1xRz3b9LTdd43rYjKpNsaWB56z6zjfP/BcVx9VeDl1JbLtgefWEcZlohDsf7INCr0ruPvurAtmMfW0qTx+9uPsq97HU4ufcjs0k+S8uFTR3MDzWBFJqWWY+QnAbmBRPeUsAiqALBHJrGVIerBx/dsYxGzqYAvExofg7rvBQREnHH4Cpz51Km1btuWiQRe5HJ1JVvXWoESkXWNv0pAyVHU18A6QBdwQdvpOIBN4JjThiMhAERkYVs5u4HGgFXC3hPwJLyLHAJNwFru1XYGbkG1QGJ+yO2fz9qVv8/OCn/PWyrfcDsckqUhqUFtEZC7wGvCGqq6LpGARGQScE3ik4vT3ROp64GPgQRE5HWdO1SicOVIrgFvDrg/OuQpvR/oDzvDyKcBxgTlU3YHzcRLXlEBCNE1gc8VmNpRt4Njux7odionCsd2P5bWJr3H2rLN58cIXOSXrFLdDMkkmkj6oGcAQ4GGgWEQ+E5E/iMiQ0IvEcbKI3Csiq4AvgbuBToQsjRSJQNIYDjyFk5h+ibMO4APA6EhXhQjMuToJ+FMgjp8B+TjDzcep6gMNics0zIJ1Cziu93E22TOOje49mucnPM+FL17Ifzf81+1wTJKJZBTfVYHmseOAc4GzcZra7hCRtcAbQDuceUadAm/7DHgSeFVVv4kmsEBN7YoIr62zB15Vy3FqXOG1LtPErP8pMYw5Ygz/PPuf5D+Xz38u/w+DuzXv3CiTvCIaJBFYvujjwON/Av095+Ekqxtw+nI+wFlW6HVVLWmSaE1c+WjtR9wz9h63wzAxcPaAsymvLGf8s+N5//L3WbFthevr9pnEF9UoPlVdBvwZ+LOIdAf2RLCEkUkSPr+Pfy/9N4s3LmZzxWZ8fp/9AksAPz7mx+zcu5Mh/xhCSkoKe6r22Lp9pkk1eh6Uqn5vyckE+fw+xj07jstfvZwqfxWTXpvEuGfH4fP73A7NxECf9n3w42d31W4UrbFunzGx5sWJuiaOFawqoLCksNaFR038KyotospXVeNYRWWFbR1vmoQlKBNTRaVF+5fLCbJfYIkjt2cumemZNY5lpmeS0yPHpYhMIrMEZWIqt2cumWn2CyxRBdfty0jLACAzzemDsnX7TFOwBGViKq9fHoO6DUIC/7VJb2O/wBJIcN2+Fya8wBEdjuDGUTfaAAnTZLy4Fp+JY6kpqUzKmURGWgZjjhhDTo8cG4acYILr9pWWlfLumnftszVNxhKUibnCkkIuPPpCrhtxnduhmCZ0/lHn8+t3f015ZTlt0tu4HY5JQNbEZ2Ju4bqFHNfnOLfDME2sc0Znju9zPG8sf8PtUEyCsgRlYmrbnm2UlJXYcjhJYuLgibzwzQtuh2ESVL1NfCKyJsqyVVX7RvleE6cK1xcy4rARtEix1uNkcM6Ac/h5wc/ZsXcHHVp1cDsck2AiqUGl4GxjEfpoibNfUxbQB2dX2j4hx1pGWLZJMAvXL2R079Fuh2GaSftW7RlzxBheXfaq26GYBFRvElHVLFU9IvjA2XqjBGfH2tOAVqraE2d/pTFAIbAesE2AktDC9Qs5rrf1PyWTiYOsmc80jWhqOX8EOgCnquo8VfUBqKpPVT/ASVqdAteZJOLz+/ik5BOrQSWZ/Ox8Pl73MVt2b3E7FJNgoklQ5wGvqWplbSdVdS/O7rvnNyYwE3+WbllKt8xudM3s6nYophllpmeS1y+PV5a+4nYoJsFEk6A6A2n1XJMWuM4kkYXrrHkvWV086GKe//p5t8MwCSaaBLUamCAi7Ws7KSIdgQlAtKP/TJyy/qfkldc/j6KNRZSWlbodikkg0SSofwCHAZ+IyOUikiUirQPPP8UZJNEDmB7LQI332Qi+5NWqRSvOHnA2Ly15ye1QTAJpcIJS1YeBh4D+wJM4NarywPMTQD/gYVV9JIZxGo/btmcb63et55jux7gdinHJxYMu5vlvrJnPxE5Uc5VU9SbgBJyEVITTnFcEPA6cGDhvkkjh+kKGHzbcJugmsR8e+UOWb1nO2p1r3Q7FJIiof5uo6kJgYQxjMXFs0fpF1v+U5NJT0zlv4Hn865t/8avjf+V2OCYB2GoPJiZsgIQBW5vPxJYlKNNofvVTWFLIqN6j3A7FuOzUrFNZt3Mdq7atcjsUkwCiSlAicoqIzBaRTSJSJSK+Wh7VsQ7WeNOSzUvomtGVbpnd3A7FuCw1JZUJR0/gha+tFmUar8F9UCJyJvAqkAqsBZYDloySmO3/ZEJNHDyR69+8nltPvtXtUEyci2aQxB1AFXCmqr4T23BMPLIBEibU8X2OZ/ve7SzZvISjux7tdjgmjkXTxDcYeMGSkwmyARImVIqkcNHRF1kzn2m0aBJUObAt1oGY+LR9z3bW7VpnE3RNDRMHT+T5b55HVd0OxcSxaBLUe4D9uWwAKCwpZFjPYTZB19Qw/LDhVPurWbxxsduhmDgWTYL6DdBXRH4vIhLrgEx8sRXMTW1EhIsHXWxzokyjRJOgbge+Ae4EVovIKyLyRC2Px2MbqvGiRSWLbASfqVUwQVkzn4lWNO0yk0K+zgo8aqPA5CjKN3HCr34K1xfyzHnPuB2K8aBBXQdR7a/m2tnXctaAs8jrl0dqSqrbYZk4Ek2COiLmUZi4tHTzUrpkdLEJuuYgPr+P8TPHs6l8E//3+f/x3NfPMarXKOZcOseSlIlYgxOUqn7XFIGY+GP7P5m6FKwqoLCkkEp/JQDlleUUlhRSsKqA/Ox8l6Mz8cLW4jNRswESpi5FpUVUVFbUOFZRWWGj+kyDRD02WER6AqcDvYCWtVyiqnpXtOUb71u4fiE3jLzB7TCMB+X2zCUzPZPyyvL9x1qntSanR46LUZl4E1WCEpE7gd+GvV9wBkaEfm0JKkHt2LuDdbvWcWz3Y90OxXhQXr88RvUaRWFJIRWVFaSlptGqRSvG9x3vdmgmjjS4iU9EfgL8AfgImICTjJ4Gfgw8BviB54ExsQvTeE3hepuga+qWmpLKnEvnMOuCWUw9bSovTHiBbhndmL1yttuhmTgSzW+X64D1wHhVrQ7M1S1W1eeB50Xk38CbwKzYhWm8xgZImPqkpqSSn52/f1BEu5btmPz6ZMb2HUtGWobL0Zl4EM0giWOAt1Q1dIuN/eNGVXUOMAf4dSNjMx5mC8SahhpzxBhG9hrJX+b/xe1QTJyIJkGlAVtDXu8B2odd8zUwJNqgjLcFJ+jaChKmoe4dey+PfPqI7bhrIhJNgioFeoa8XguE95Qfhm1imLCWbVlG54zONkHXNFjvdr35nxP+h5vevsmWQDL1iiZBFeHsCRX0PnCSiFwmIpmBHXcnBK4zCcjmP5nGmDJ6Cmu2r+GNFW+4HYrxuGgS1GxgsIgElzz6C7ATeArYBbyOM7Lv97EI0HiP9T+ZxkhPTeehvIeY8vYU9lTtcTsc42ENTlCq+pSqZqjqt4HX64ARwN+Bd4D/A0ao6qLGBCYivQOrom8QkX0iUiwi00SkYyPKPFlEfCKiInJ3Y+JLZjaCzzTWD4/8IcMOG8ZfF/zV7VCMh8VkEksgWf0sFmUBiEhf4GOgG/AasAwYCdwEjBeRE1R16yGKqK3MtjjztXYDbWIVa7LZsXcH3+34zibomka7b+x95D6ay+VDLufIjke6HY7xIK+uxfcITnK6UVXPVdXfquoY4H5gAPDHKMp8AGe04Z9jF2byKVxfyLDDhpGWmuZ2KCbO9Wnfh18e90umvD3F7VCMR3kuQQVqT2OBYmB62OnbgQrgMhHJbECZ5wBXADcCG2ITaXJatH6R9T+ZmLn5uJtZtmUZb6540+1QjAd5LkEBpwWe31FVf+gJVS0DFgAZQESdICLSDWcJpldV9dlYBpqMbICEiaWWLVryUN5D3PT2Teyt3ut2OMZjvJigBgSeV9RxfmXgOTvC8h7D+T6vbUxQyc7n9/H68teZ9908du3bhc/vczskkyDG9RvHMd2O4arXruLBLx5k9orZ9u/LADEaJBFjwVUpdtZxPni8Q30FiciVwNnAxar6fUMDmTLlQNv46NGjycrKamgRCcHn93H5u5fz+ebP2evby7Wzr+Ufi/7BjB/N8PTuqBs3bnQ7BBMBn9/Hum3r+GzzZwjCo18/Sk6XHM//+zKOpvz/zIsJKiZEJAuYBryoqv+Kpoxp06bVeF1cXJyUSWr2itl8ue1L9vqcJpjd1bv5ctuXLK1e6vndUZPx84o3s1fMZvnO5QAoGlf/voyjqf4/82ITX7CGFL6+H2HHd9RTzhM46wReH4ugkpntjmqakv37MnXxYoJaHniuq4+pf+C5rj6qoKE4Q9U3Bybmqogo8GTg/K2BY682LtzEl9szl9ZprWscy0zPtN1RTUwEd98NlZGWYf++TNQ76qYCpwLn4PTx7AJeBV5X1f82Mqa5geexIpISOpIvMNn2BJzJtvWtVDEDZ7RfuP7AycBi4DNszcB65fXLo116O6p91VT5q8hMz2RUr1Hk9ctzOzSTAMJ3302VVLpkdLF/XybyBBVIDnk4SSkPp6lNcRJFO5y1924VkVKc1R9eA+aqalVDAlLV1SLyDs5cqBuAh0JO3wlkAo+q6v42AREZGHjvspBybqzj+5iEk6DeVFVbLzACa7avodJXyYzzZrBy20pyeuSQ1y/POrBNTAR33y1YVcDcpXMZkjWEm+c486MGdRvkdnjGRfUmKBG5DicpnQqk4/TrvIuTgGar6ubAdUcA5wauvRpnWHe5iBQAs1T1tQbEdT3OUkcPisjpwFJgFM4cqRXArWHXLw2G24B7mAj9af6fuHHUjVw8+GK3QzEJKrj77uD0wWRlZbFr3y6uf+t6PvjpBwR27TZJKJI+qOlALjATJwF1CSw/9GQwOYGzHp+q3q+qpwLdcVZu+A9wJs4KEBFT1dXAcJwV0kcBvwT64ixXNLqh6/CZ6K3ZvobXl7/OjaNqrZAa0ySuG34dFZUVzPhihtuhGBdF0sR3EvCxNmB3MVXdhtMHNENEWhLF7rqBVdKviPDaiP/EUtWncBKficCfP/oz1w+/no6to15E3pgGS01J5R/5/yD/uXzOGnAWnVp3cjsk44J6a1CquqAhyamW9+9T1U+ifb9xz3c7vuOVZa8wZbQt5mma3/DDhjPh6An89j+/dTsU4xIvDjM3HvGX+X/h6qFX0zmjs9uhmCR195i7mb1iNgvXLXQ7FOOCRiUoEXlZRDICX0e6Np6JA+t2ruOFb17g5uNudjsUk8Q6tOrAvWPv5do3r6XaX+12OKaZNbYGVQoENwZaJiJlIvKxiDwiIleLyEgRadXIexgX/G3B37hq6FV0zezqdigmyU0cPJGuGV15sPBBt0MxzaxRa/Gpauguun1wBkPkBJ6DI+9URFaq6tGNuZdpPhvKNjDzq5ksvWFp/Rcb08REhEfOfITjHz+eiwZdRO92vd0OyTSTmC0Wq6olQAnwVvBYoPlvCGD7g8eRvy34G5NyJtG9TXe3QzEGgOzO2dww4gamvD2Fly56ye1wTDOJuIlPRO4QkVIRqRSRlSLyexE55L7fqrpbVReq6qOND9U0h43lG5nxxQx+ffyv3Q7FmBp+d9LvWLxxMW+tfKv+i01CiChBBfZVug1nAm4LnKa7OwH7UybB3PPxPVx67KX0bNvT7VCMqaFVi1ZMP2M6P3vrZ+yu2u12OKYZRNrEdy1QiTNxdh4wEPgrkC8iF6rqi00Un2lGmyo28UTRE3x53Zduh2JMrcb1G8eww4ZxxatXMLjbYHJ75tq6kAks0gTVF3hJVWcFXm8QkR8Bq4DLAUtQCeC+hfcxcfBE64Q2nuXz+9hYtpH56+bz4pIX96+sP+fSOZakElCkfVAdcZLRfqq6A3gTZ98lE+e27N7CY58/xm9PtFn7xrsKVhWw+HtnI0NFKa8sp7CkkIJVBS5HZppCQ+ZB+Ws5thawZQYSwLRF05hw1AQOb3+426EYUyfbfTe5NHaYeTUHJuqaOLV9z3b+/t+/89nVn7kdijGHFNx9t7yyfP8x2303cTWkBvUHEflKRB4PrBKRQwznUZnm5/P7mL1iNuf/63yG9hhKn3Z93A7JmEMK7r7bJr0NgpCWkka7lu1s990EFWmC+Q9OX9OgwGNS6EkRuQdnC/XFwJLQbdqNN/n8PsY9O45F6xdRUVVBRloG454dZ53NxtNCd99dvHEx2Z2zuXnOzcxfO59Tsk5xOzwTYxElKFUdC/t3zR0e8hiKs/X7zTjbvwPsE5GvgSJVvSbmEZuYKFhVQGFJIRVVTnv+7qrd+zub87PzXY7OmLoFd98N/jtNS0nj2jevZfE1i2nZoqXL0ZlYatBisYFdc19U1d+o6umq2hEYAPwEmAbMx+mXGg5cFfNoTcwUlRbVaMcH62w28encgeeS3Tmbvy34m9uhmBhrdB+Sqq4EVgKzAEREcCbyDm9s2abp+NR30LHM9EzrbDZxR0R4KO8hhj46lImDJ9K/c3+3QzIxUm8NSkTaNaRAdSxV1WeiLcM0rcUbF/PwJw8zrOew/Z3NbdLbMKrXKOtsNnHp8PaHc8tJt3Ddm9fRiA3AjcdEUoPaIiJzgdeAN1R1XSQFi8gg4JzAIxWrUXlCya4Szp51No+c+QgXHHXB/s7mnB45tmSMiWs3jrqRZ758hplfzeTSYy91OxwTA5EkqBlAPvAj4CERWQy8Cryuql8ELwo07Z3EgaR0BCDAauCfMY7bRKG8spz8WflcP+J6Lhp0EUCNzmZj4lmLlBY8mv8o5zx/Dmf0P4NOrTu5HZJppHqb+FT1KqAncCJwL5CJs5L55yLyrYg8KCJPAZuAucAUYBvO6ufHqGp/Vf1rE8VvIuTz+5j40kSG9RzGb074jdvhGNMkRvYayYSjJvCbd+3feCKIdJi5Ah8HHv8jIgOB84CzgRtwRu59wIGaVUmTRGui9os5v2Cfbx9/P/PvOJVdYxLT3WPuZtAjg5i/dj4nHn6i2+GYRmjQMPMgVV2mqn9W1eOAw4CuqjpOVf9uycl7Hlj0AO99+x4vXfgSaam2MpVJbO1btWfa+GlcM/saKn2VbodjGiGqBBVKVb9X1V2xCMbE3mvLXuOvC/7Kmz9+k/at2rsdjjHN4oKjLuCIDkdwz8f3uB2KaQRbSy/B+Pw+ClYVUFRaRNv0ttz90d289ZO3yOqQ5XZoxjQbEeHhMx5m2KPD6JrRlY3lG21zwzhkCSqBBNfXKywp3L8lweBugxnWc5jLkRnT/Pq060On1p247s3r8KvfNjeMQ41u4jPeEVxfr7yyHA389+2Ob20zN5OUClYVsLF8Iz712eaGccoSVAKxzdyMOaCotGj/YshB9v9DfLEElUAGdBlw0BByW1/PJKvg5oah7P+H+GIJKkHs2reLaYum0bNNT9qk2fp6xgQ3N8xMc5JUmzT7/yHe2CCJBLBr3y7yZuYxpPsQ5k2ax5zVc2x9PZP0Qjc3nPjSRB464yEuO/Yy+/8hjliCinO79u1i/LPjGdJ9CNPPnE6KpNj6esYEBDc3HHbYMA5vf7glpzhjTXxxrLbkZIw5WHanbJZvWe52GKaB7DdanAomp5weOZacjKnHgC4DWLF1hdthmAay32pxaNe+XYx7dhw5PXJ4+IyHLTkZU48BnQewfKvVoOKN9UHFgdDli7I7Z3P/ovsZ2nMo08+YbiuTGxOB7M7ZlqDikCUojwtfvkhE6NmmJ/OvmG/JyZgIHdnxSEp2lbCveh8tW7R0OxwTIWsb8rjw5Yv86mfn3p28vfptt0MzJm6kpaZxePvDWb19tduhmAawBOVhO/fu5LHPHqO8srzG8YoqW67FmIaygRLxx5r4XBLarxS6DUC1v5p3V7/L0188TcGqAgZ1HUSrFq3YW713/3ttuRZjGm5A5wE21DzOWIJyQXi/UmZ6JoO6DuKEPicw6+tZ9Gnfh8uPvZzpZ0ynQ6sOB11ry7UY03DZnbNZtH6R22GYBrAE5YLQfiWA8spyPin5hD7t+vDe5e9xVNejalwfXK7Fli8yJnoDOg9gxhcz3A7DNIAlqGa2t3ovz3zxzEH9SgBDegw5KDnBgeVabPkiY6JnQ83jjyWoGKqrX8mvfhasXcCML2bw8tKXObz94bRMbck+377977V+JWOaVo82PdhXvY9te7bRqXUnt8MxEfBsghKR3sBUYDzQGSgFXgXuVNXtEbw/EzgXOBMYCvQB/MByYBbwkKpWxire2vqVjul2DGOOGMPMr2aSkZbB5cdezpfXfUnPNj2tX8mYZiYi+0fyje492u1wTAQ8maBEpC/wMdANeA1YBowEbgLGi8gJqrq1nmJOAp4FtgFzcZJbR+Bs4B7gfBE5XVX31l1E5GrrV1q4fiFdMrrw8kUvk9sjt8bEWutXMqb5ZXd2Fo21BBUfPJmggEdwktONqvpQ8KCI3Af8AvgjcG09ZWwELgVeDK0picivgA+A44EbgHtjEXBt260LwsheIxnac+hB11u/kjHNb0BnmwsVTzw3UTdQexoLFAPTw07fDlQAlwWa8OqkqotVdWZ4M56qlnEgKZ0ai5jBtpc2Jh7YorHxxXMJCjgt8PyOqvpDTwSSywIgA2hMHb0q8FzdiDJqCG4v3Sbdtls3xqtsJF988WIT34DAc1318JU4Naxs4L0o73Fl4PmQC9pNmTJl/9ejR48mKyvrkIU+euKjfFDyAUu2LeHoTkdzaq9TWbd2XZQhmsbauHGj2yGYBmrqzyy9Kp2VW1ey5ts1tk1NjDTlZ+bFBNU+8LyzjvPB4x2iKVxEfoYzMnAx8MShrp02bVqN18XFxfUmqb5H9o0mLNNE6vu8jPc09WfWJaMLKR1TyOrQtPdJJk31mSXVnxAicj4wDWcAxQWqWlXPW4wxCcYWjY0fXkxQwRpS+zrOB4/vaEihInIu8DywCThVVddEF54xJp7ZorHxw4sJKvgvJ7uO8/0DzxH/CSQiFwIvAt8Dp6iq/es0JknZQIn44cUENTfwPFakZi+miLQFTgB2AxEtSywiP8FZOWIDTnJaGcNYjTFxxuZCxQ/PJShVXQ28A2ThTKQNdSeQCTyjqvtnxYrIQBEZGF6WiPwUmAGsBU62Zj1jjNWg4ocXR/EBXI+z1NGDInI6sBQYhTNHagVwa9j1SwPP+9cSEpHTcEbppeDUyq4IXWooYIeqTgs/aIxJXFkdsthUsYndVbvJSMtwOxxzCJ5MUKq6WkSGc2Cx2DNwFot9gAgXiwV+wIEa4pV1XPMdzqg+Y0ySSE1J5ciOR7Jq2yqO7X6s2+GYQ/BkggJQ1XXAFRFee1DVSFWfAp6KbVTGmEQQXDTWEpS3ea4PyhhjmpoNlIgPlqCMMUnHBkrEB0tQxpikY6uaxwdLUMaYpBNc7khV3Q7FHIIlKGNM0uncujOCsHn3ZrdDMYdgCcoYk3RExBaNjQOWoIwxSckWjfU+S1DGmKRkI/m8zxKUMSYp2Vwo77MEZYxJSlaD8j5LUMaYpNSvUz++3f4t1f5qt0MxdbAEZYxJSq3TWtOzbU+KdxS7HYqpgyUoY0zSCi4aa7zJEpQxJmnZQAlvswRljElaNlDC2yxBGWOSli0a622WoIwxScuWO/I2S1DGmKTVu11vtu/ZTtm+MrdDMbWwBGWMSVopkkL/zv1ZuW2l26GYWliCMsYkNVs01rssQRljkpqN5PMuS1DGmKRmc6G8yxKUMSapWQ3KuyxBGWOSWnCouaq6HYoJYwnKGJPUOrTqQEZaBqXlpW6HYsJYgjLGJD1bNNabLEEZY5KeDZTwJktQxpikZwMlvMkSlDEm6dmisd5kCcoYk/Rs0VhvsgRljEl6R3Y8knU711Hpq3Q7FBPCEpQxJumlp6bTp30f1mxf43YoJoQlKGOMwYaae5ElKGOMwQZKeJElKGOMweZCeZElKGOMweZCeZElKGOMwRlqbn1Q3mIJyhhjgG4Z3SirLOOW925h9orZ+Pw+t0NKei3cDsAYY9zm8/sYP3M8lb5K/jL/L2SmZzKq1yjmXDqH1JRUt8NLWlaDMsYkvYJVBRSWFOJXP4pSXllOYUkhBasKar3e5/cxe8Vs7pp3l9W2mpDVoIwxSa+otIiKyooax8ory7n7w7vZvmc7I3uNpH/n/qRICj6/j3HPjqOwpJCKygqrbTUhS1DGmKSX2zOXzPRMyivL9x9r3aI1g7sN5s2Vb3LbB7exfc92RvQaQcdWHVmwbgF7q/cC1Kht5WfnH1S2z++jYFUBRaVF5PbMJa9fniWyCFmCMsYkvbx+eYzqNeqgWtGj+Y/uTyabKjbxacmn3PPxPfuTU1B5ZTnTFk2jokdURR4AAA2VSURBVLKC/p37079Tf9q2bNvg2pYls5osQRljkl5qSipzLp1DwaoCFm9cTE6PnIOSQ7fMbpyZfSaK8t/S/9aobbVMbUnHVh15ccmLrNi6glXbVtG+VXs6te7Eiq0rqPZXA04iW7h+IW+seINzB55bIwZrOjyYZxOUiPQGpgLjgc5AKfAqcKeqbm9AOZ2A24BzgZ7AVuBt4DZVXR/ruI0x8Sk1JZX87Pxam+lC1VXben7C8/sTiV/9bCjbwG1zb2PJ5iU13r+7ajcX/OsCerTpQZ92fejTvg+92/amrLKMBWsXsNdnTYdBnkxQItIX+BjoBrwGLANGAjcB40XkBFXdGkE5nQPlZAPvA88DA4ErgDNF5DhVteWLjTERi6S2lSIp9G7Xm/OPOp8Xl7xYo7bVJr0Nz573LLk9c1m3cx3rdq1j3c51vLL0lf3JKai8spyb59zMWyvf4ogOR5DVIYsjOh5Bn3Z9+PErP+aTkk/qrW3FcyLzZIICHsFJTjeq6kPBgyJyH/AL4I/AtRGU8yec5HSfqv4ypJwbgQcC9xkfw7iNMUmgsbWt/Ox8UlNSObz94fuvParrUVzy8iUHDdSYOGgiXTO7UryjmIXrF1K8o5gVW1dQUXVg1GF5ZTkffvchNxbcyPh+4+nVrheHtT2Mzq07kzczr0n6wILXvr/0fcZUjmmSxCeqGtMCGytQe1oFFAN9VdUfcq4tTlOfAN1UtaLWQpxr2wCbAD/QU1XLQs6lAGuAHwTucVAtSkQ0/GdTXFxMVlZW1N+baV72ecWfRPzMgr/I66pthV4XaR/UXfPu4vYPbkep+TtqWM9hdG/TnZJdJWwo28C2Pdv2z+0KSktJY9KQSZySdQpdM7vSNaMrXTK60Kl1J855/pyI7h/L/jIRQVWltnNerEGdFnh+JzQ5AahqmYgsAMYCo4H3DlHOaKB1oJyy0BOq6heROcDVgftZM58xpklEWtuKpOkwqLZh8W3S23DHqXfUuM/tc2/nrg/vqvHeKn8VX236irKqMrbs3sLmis1s3r2ZTRWb9g/mAKdWNq94Hmc+dyY5PXLo2KojHVt3pEOrDizfspyP133Mnuo9+689VH9ZtLyYoAYEnuta934lToLK5tAJKpJyCJRjjDGua2zTYV6/vBrXjeg1otZEduvJtx50j6nzpnLHB3fUqG351Ed6ajrtW7Zn656trNq2iu17t1NUWrQ/OQVVVFaweOPihE9Q7QPPO+s4HzzeoanLmTJlyv9v7+5j5KrKOI5/f1Ao24LQBovIW6VQuqS0hQpUC7GlptGYiBEwxPBSY4pREBU0EcVEjC8xoKlo1BIFsUZUVMIfJYFYWkrQaiwWAm0pAoWCVISWN9uipY9/nDN1Op3pzu7e3bmz9/dJJqc75+7pufvs3Wfuueeeu/vfs2bNGnFDDyPd5s2bO90F6yfHrH2Lz1rMiudWsHbLWk4efzJzjprDpmc27bFN76hepo2fxpoX17B953Z6RvUwbfw0ekf1snHjxj22PWq/o+gZ1cO2ndt2v9czqodzjz6XecfM22PbZZuWceXKK/fa9kgduVe7g1HGBFUaixYt2uPrkTg+PtI5Xt3HMWvfpOMn9bnNyoUr2xo2XHDsAm576ra9zsoWzF6w1/b92XYwypigamc2h7aor73/8jC1Y2bWtYbiGlj9tsvXLWdu79whmcVXxgRVe2JYq2tDJ+ayr2czF9WOmVkltJvM6redeuDUITvrLePjNpbncn6eDr5bnmY+G9gGrOqjnVXAdmB2/r76dvYjTbSo///MzKxESpegIuIJ4B5gInB5Q/V1wFhgSf09UJKmSJrS0M7rwJK8/Vcb2rkit3+3V5IwMyunMg7xAXyKtETRjZLmAeuAM0n3LG0Avtyw/bpcNt7s9SVgDnCVpBnAX4Be4FzSTbyNCdDMzEqidGdQsPss6p3Az0iJ6WpgEml5olntrMOX23kJeBdwI3BCbudM4BZgZv5/zMyshMp6BkVEbCIt6trOtk2Xych1W0iLzH6moK6ZmdkwKOUZlJmZmRNUP6xa1dfEQSsTx6v7OGbdZyhj5gTVDz54uovj1X0cs+7jBGVmZpXjBGVmZqVUugcWloUk/2DMzIZBq5nYTlBmZlZKHuIzM7NScoIyM7NScoIyM7NScoIyM7NSqmyCknS0pJsl/UPSG5I2SlokadwA2jpN0i8lPZvb+qek+yRdMhR9r6rBxkzSHEnRxuuYod6XqijqOJN0lqQ78/fvkPSMpLskvW+o+l5VBcbsPEkrJL0iabukRyVdI+nAttuo4iw+SZNIj/OYANwJrAfOID3O4zFgdrsrpku6grTK+lZgKfAcMB6YCjwbERcWvgMVVETMJE0EFrSoPgX4MPBIRJxSSKcrrqjjTNIngR8C/wbuAJ4FjibFawxwbUR8Yyj2oWoKjNk3gWuA14HfAVuAs0lPqVgGvD8i/ttnhyKici/gbiCATze8/938/o/bbGc+sCu3d0iT+gM6va8j5VVUzPbR/m25nSs7va8j5VVEzIADgJdJT8c+qaGuF9hBesL26E7v70h4FRSz0/K2W4Hj694X8KNcd1U7/ancGVT+hPB3YCMwKSJ21dUdAjxP+kFOiLqn9rZo6yHSc6aOjTbPuKz/ioxZi/YPJ30q3wW8PSJeLqLfVVZUzCQdAWwGHo6I6U3qHyad/R7uY3BwCozZ14CvADdExBca6saRzqaejIhJffWpiteg5ubynvoAAETEa8ADpGGDWftqRNJUYBrp8fRbJM2V9HlJV0uaJ6mKP9uhUkjM9uFSYDRwu5NTYYqK2QvAv4DJkk6sr5A0GTgRWOPkVIiiYva2XD7ZWBERW8lnVpLe0VeHqvhH9KRcbmhR/3guJ/fRzum5fAFYAdwLXA/cAPwBWCPphIF30+oUFbNWFuZy8QC/3/ZWSMwiDfFcTvpbtVrSrZK+JennwGrgUeCCAvprxR1nL+ZyrwQk6TCgNtnipMb6RlVMUIfm8pUW9bX3D+ujnQm5/DgwEfhAbnsy8AvSsMPS/sxYsZaKitleJL2HdKA8EhF/HEDfrLnCYhYRtwPnkK5FXQJ8EbiYNGniFpp8UrcBKSpmS3O5ME9MAkCSgPrJLH3OCqxigipK7We3P3BhRNwVEa9GxOOkg+ivpGR1Xqc6aG25LJc3dbQX1pKki0ijEveTJkaMyeUy4AfArzrXO2sUEQ8APyUlsocl3SLpO8CfSR/o1+dNd7VoYrcqJqjap4BDW9TX3u/rWkStfnNE/Km+Ig9L3Jm/PKPfPbRGRcVsD5LGkz5AbAeWDKxr1kIhMcvXmW4mDeVdHBHrI2J7RKwnnUWtBi6QNGfwXa68Io+zhcAnSFPTP5L//SowB3gib/NCX41UMUE9lstW46i1C7GtxmEb22kVrK257GmzX9ZaUTFrVJsc8RtPjihcUTGbT5pqfl+TC/e7gJX5y5kD6aTtobDjLJKbIuL0iBgbEQdHxHsjYhXp8scu4MG+2qliglqey/mNM+3yVMrZpPsq+nqO8SrSGPhESWOb1E/N5VOD6KslRcWsUW1yhIf3ildUzEbn8q0t6mvv/2cgnbQ9DNVxVt/OHOBYYGlEtLrWtVvlElREPEGaGj6RNDuo3nXAWGBJ/Tx/SVMkTWloZxtpnPUg4Ov5AmBt+1NIKxbsBH5b/F5US1ExqyfpbNJ1DE+OGAIFxuz+XJ4vaVp9haQZwPmkGz/vLa731VTkcSbpLU3eOw74CenDxLXt9KlyN+pC0+U81gFnku4D2AC8u/6+CuWn60bDUx9zEO4DZpAuAD4AHEFagqUH+GxEfG+o96cKiopZXf0S4CLSyhHfH9reV1OBx9nNwMdIf9juAJ4m/RH9EHAgsCgiPjfEu1MJBcbsduA40jDeFtKU8w+Shmsvjohft9WhTiynUYYXcAxpiurzpF/8p4FFwLgm2wZ57kOTuoNJUyc3AG+QrkndA8zv9D6OtFeBMRtHmhixDTis0/s1kl9FxIy0esEC0v2GW0kjE1tIs/gu7PQ+jrRXQTG7lPSB/aXcxibgVqC3P32p5BmUmZmVX+WuQZmZWXdwgjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjIzs1JygjLrYpIelbRN0hhJ10haK2mHpOclXS9p/0730WygvNSRWZeS1AO8Rlor7UXSyvr3kxZQ/Shp9WkvWGxdy2dQZt1rOrA/cDzwe2BGRFwREZfx/0fZn9OpzpkN1qhOd8DMBuzUXN4REd9uqKs9rdQfQq1r+ZfXrHudlstmTwSekMsnh6kvZoVzgjLrXqcCb5Iemtloei7/NnzdMSuWE5RZF5J0ADAVWBcR25tsMjOXq4evV2bFcoIy604nA6NpnYBmkp4avHbYemRWMCcos+5Uu/70YGOFpDFAL/BQRLw5rL0yK5ATlFl3qs3ga3YGVZt+7uE962pOUGbd6VRgF7CmSV3t7MoJyrqa74My60IRcfY+qj1BwkYEn0GZjTwzgR14goR1OScosxFE0kGkGX4PRcTOTvfHbDCcoMxGlumkoXsP71nX82rmZmZWSj6DMjOzUnKCMjOzUnKCMjOzUnKCMjOzUnKCMjOzUnKCMjOzUnKCMjOzUvofBnX8ryaSaqwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot the principal components\n",
    "colors = 200*['green', 'purple', 'orange', 'teal', 'black', 'red']\n",
    "for i in range(0,1): #12,13\n",
    "    #p_data = (OPs[:,i]-OPs[0,i])/max(OPs[:,i]-OPs[0,i])\n",
    "    #p_data = -OPs[:,i]\n",
    "    #plt.plot(etas, (p_data), marker='o', linewidth = 1.0, markersize=10, color = colors[i], markerfacecolor='None')\n",
    "    p_data = OPs_std[:,i]/max(OPs_std[:,i])\n",
    "    #plt.plot(etas, (p_data), marker='^', linewidth = 1.0, markersize=12, color = colors[i])\n",
    "    #p_data = OPs_abs[:,i]\n",
    "    plt.plot(etas, (p_data), marker='.', linewidth = 1.0, markersize=10, color = colors[i])\n",
    "\n",
    "#backround grid details\n",
    "axes = plt.gca()\n",
    "axes.grid(b = True, which = 'both', axis = 'both', color = 'gray', linestyle = '-', alpha = 0.5, linewidth = 0.5) \n",
    "\n",
    "\n",
    "#figure size and tick style\n",
    "plt.rcParams[\"figure.figsize\"] = [6,6]\n",
    "plt.rc('axes',edgecolor='black',linewidth=1)\n",
    "plt.tick_params(which='both', axis='both', color='black', length=4, width=0.5)\n",
    "plt.rcParams['xtick.direction'] = 'in'\n",
    "plt.rcParams['ytick.direction'] = 'in'\n",
    "\n",
    "#axis range and labels (also specify if log or not)\n",
    "#plt.xlim(0.60, 0.84)\n",
    "\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r'$\\eta$', y=3, fontsize=20, fontdict = axis_label_font)\n",
    "plt.ylabel(r'$\\langle P_{i} \\rangle$ and $\\langle |P_{i}| \\rangle$', fontsize=20, fontdict = axis_label_font)\n",
    "\n",
    "#title and axis labels\n",
    "plt.tick_params(axis='both', labelsize=20)\n",
    "plt.title('PCA ellipsoids', y=1.05, fontdict = title_font)\n",
    "\n",
    "#save figure\n",
    "fig=plt.figure(1)\n",
    "plt.tight_layout()\n",
    "#plt.savefig('figures/OP_orientation_batch_kappa_4_ave_std.svg',dpi=fig.dpi)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save OPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Append eta ranges\n",
    "etas2=etas.reshape(etas.shape[0],1)\n",
    "\n",
    "#etas2=etas[...,None] #wp: works just as fine\n",
    "OPsave=np.concatenate((etas2,-OPs),1)\n",
    "\n",
    "#save to txt file\n",
    "np.savetxt('OPs_kappa4_30k_300_3N.dat',OPsave)\n",
    "np.savetxt('OPs_kappa4_30k_300_std_3N.dat',OPs_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}