pmf_estimators_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BAdp1h4FSIgQ"
   },
   "source": [
    "# PMF: Analysis of expectation and variance estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The notebook illustrates bias and variance of the estimators used to match Poisson Matrix Factorization priors in pmf_sgd_optimization.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RBAeIwOMrYk8"
   },
   "source": [
    "## Setup \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import absolute_import, division, print_function, unicode_literals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "attributes": {
     "classes": [
      "py"
     ],
     "id": ""
    },
    "colab": {},
    "colab_type": "code",
    "id": "N-2lSGiHqEfq"
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "#tf.enable_eager_execution()\n",
    "import tensorflow_probability as tfp\n",
    "from tensorflow_probability import distributions as tfd\n",
    "#from tensorflow.contrib.distributions import softplus_inverse\n",
    "#tfe = tf.contrib.eager"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(True, '2.1.0', '0.9.0')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.executing_eagerly(), tf.__version__, tfp.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "attributes": {
     "classes": [
      "py"
     ],
     "id": ""
    },
    "colab": {},
    "colab_type": "code",
    "id": "9gsI54pbqEfj"
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from aux import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from aux_plt import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parsing: <-f>\n"
     ]
    }
   ],
   "source": [
    "args = parse_script_args()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "SEED = args.get(\"SEED\", 123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = args.get(\"K\", 25) # factorization dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DESCRIPTION: K=25 SEED=123\n"
     ]
    }
   ],
   "source": [
    "DESCRIPTION = dict2str(globals()).replace(\"NSAMPLES\", \"S\").replace(\"DESIRED_\", \"\").replace(\"EXPECTATION\", \"E\").replace(\"VARIANCE\", \"V\").replace(\"TRAINABLE\", \"TT\")\n",
    "print(\"DESCRIPTION: %s\" % DESCRIPTION)\n",
    "\n",
    "ID = str(args.get(\"ID\", DESCRIPTION))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expectation & variance evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "empirical_moments := empirical_moments_fully_decoupled\n"
     ]
    }
   ],
   "source": [
    "from pmf_model import create_moments_estimator\n",
    "\n",
    "\n",
    "theoretical_moments = create_moments_estimator(K=K, ESTIMATOR_NO=-1, N=1, M=1)\n",
    "empirical_moments = create_moments_estimator(K=K, ESTIMATOR_NO=0, N=1, M=1)\n",
    "\n",
    "print(\"empirical_moments := %s\" % empirical_moments.__name__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def empirical_moments_mean_std(a, b, c, d, nsamples_latent, nsamples_output, nruns):\n",
    "    \"\"\"Returns means and stds of expectation and variance.\"\"\"\n",
    "    ev = []\n",
    "    a = tf.Variable(a, dtype=tf.float64)\n",
    "    b = tf.Variable(b, dtype=tf.float64)\n",
    "    c = tf.Variable(c, dtype=tf.float64)\n",
    "    d = tf.Variable(d, dtype=tf.float64)\n",
    "    for r in range(nruns):\n",
    "      #print(\"%i/%i\" % (r, nruns))\n",
    "      ev.append( empirical_moments(a, b, c, d, nsamples_latent, nsamples_output) )\n",
    "    ev = np.array(ev)\n",
    "    e_mean, var_mean = ev.mean(0);\n",
    "    e_std, var_std = ev.std(0)\n",
    "    return (e_mean, e_std), (var_mean, var_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of number of samples "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ABCD_CFGS = [ (10, 1, 10, 1), (10,2,10,2), (0.001, 0.01, 0.01, 0.1),(0.1, 1, 0.1, 1), \n",
    "              (0.1, 0.1, 0.1, 0.1), (1.0, 1.0, 0.1, 0.1), (1000.0, 1000.0, 1000.0, 1000.0),]\n",
    "LABELS = [\"A\",\"B\",\"C\",\"D\",\"E\",\"F\",\"G\"]\n",
    "\n",
    "LABEL2CFG = dict(zip(LABELS,ABCD_CFGS))\n",
    "\n",
    "NSAMPLES_CFGS = [(100,100),(1000,10),(10,1000)]\n",
    "\n",
    "nruns = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a=10.000\tb=1.000\tc=10.000\td=1.000\t=> E=2500.00 \tV=55000.00\n",
      "a=10.000\tb=2.000\tc=10.000\td=2.000\t=> E=625.00 \tV=3906.25\n",
      "a=0.001\tb=0.010\tc=0.010\td=0.100\t=> E=0.25 \tV=253.00\n",
      "a=0.100\tb=1.000\tc=0.100\td=1.000\t=> E=0.25 \tV=0.55\n",
      "a=0.100\tb=0.100\tc=0.100\td=0.100\t=> E=25.00 \tV=3025.00\n",
      "a=1.000\tb=1.000\tc=0.100\td=0.100\t=> E=25.00 \tV=550.00\n",
      "a=1000.000\tb=1000.000\tc=1000.000\td=1000.000\t=> E=25.00 \tV=25.05\n"
     ]
    }
   ],
   "source": [
    "for a,b,c,d in ABCD_CFGS:\n",
    "    e, v = theoretical_moments(a,b,c,d)\n",
    "    print(\"a=%.3f\\tb=%.3f\\tc=%.3f\\td=%.3f\\t=> E=%.2f \\tV=%.2f\" % (a,b,c,d,e,v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a=10.00\tb=1.00\tc=10.00\td=1.00\tnl=100 no=100 => E=2499.83 +/- 23.15 \tV=54750.87 +/- 161514.53 (0.01s)\n",
      "a=10.00\tb=1.00\tc=10.00\td=1.00\tnl=1000 no=10 => E=2500.44 +/- 7.29 \tV=54497.97 +/- 51442.51 (0.02s)\n",
      "a=10.00\tb=1.00\tc=10.00\td=1.00\tnl=10 no=1000 => E=2501.72 +/- 70.49 \tV=68000.31 +/- 492757.25 (0.01s)\n",
      "a=10.00\tb=2.00\tc=10.00\td=2.00\tnl=100 no=100 => E=625.01 +/- 5.82 \tV=4651.16 +/- 10143.96 (0.01s)\n",
      "a=10.00\tb=2.00\tc=10.00\td=2.00\tnl=1000 no=10 => E=625.07 +/- 1.80 \tV=3858.57 +/- 3320.98 (0.03s)\n",
      "a=10.00\tb=2.00\tc=10.00\td=2.00\tnl=10 no=1000 => E=624.24 +/- 17.75 \tV=2700.05 +/- 32320.39 (0.01s)\n",
      "a=0.00\tb=0.01\tc=0.01\td=0.10\tnl=100 no=100 => E=0.20 +/- 1.14 \tV=82.64 +/- 835.59 (0.01s)\n",
      "a=0.00\tb=0.01\tc=0.01\td=0.10\tnl=1000 no=10 => E=0.25 +/- 0.46 \tV=160.80 +/- 1088.81 (0.02s)\n",
      "a=0.00\tb=0.01\tc=0.01\td=0.10\tnl=10 no=1000 => E=0.17 +/- 2.27 \tV=132.95 +/- 8500.37 (0.00s)\n",
      "a=0.10\tb=1.00\tc=0.10\td=1.00\tnl=100 no=100 => E=0.25 +/- 0.05 \tV=0.56 +/- 0.39 (0.01s)\n",
      "a=0.10\tb=1.00\tc=0.10\td=1.00\tnl=1000 no=10 => E=0.25 +/- 0.02 \tV=0.56 +/- 0.14 (0.02s)\n",
      "a=0.10\tb=1.00\tc=0.10\td=1.00\tnl=10 no=1000 => E=0.26 +/- 0.20 \tV=0.56 +/- 2.03 (0.00s)\n",
      "a=0.10\tb=0.10\tc=0.10\td=0.10\tnl=100 no=100 => E=24.94 +/- 5.55 \tV=2859.72 +/- 3158.15 (0.01s)\n",
      "a=0.10\tb=0.10\tc=0.10\td=0.10\tnl=1000 no=10 => E=25.10 +/- 1.76 \tV=3012.60 +/- 1298.92 (0.02s)\n",
      "a=0.10\tb=0.10\tc=0.10\td=0.10\tnl=10 no=1000 => E=24.96 +/- 19.28 \tV=2440.08 +/- 6949.70 (0.01s)\n",
      "a=1.00\tb=1.00\tc=0.10\td=0.10\tnl=100 no=100 => E=25.02 +/- 2.18 \tV=555.80 +/- 342.22 (0.01s)\n",
      "a=1.00\tb=1.00\tc=0.10\td=0.10\tnl=1000 no=10 => E=24.99 +/- 0.71 \tV=549.23 +/- 102.89 (0.02s)\n",
      "a=1.00\tb=1.00\tc=0.10\td=0.10\tnl=10 no=1000 => E=24.81 +/- 7.04 \tV=525.09 +/- 1136.68 (0.01s)\n",
      "a=1000.00\tb=1000.00\tc=1000.00\td=1000.00\tnl=100 no=100 => E=25.00 +/- 0.05 \tV=25.14 +/- 3.80 (0.01s)\n",
      "a=1000.00\tb=1000.00\tc=1000.00\td=1000.00\tnl=1000 no=10 => E=25.00 +/- 0.05 \tV=25.03 +/- 3.57 (0.02s)\n",
      "a=1000.00\tb=1000.00\tc=1000.00\td=1000.00\tnl=10 no=1000 => E=25.00 +/- 0.09 \tV=24.84 +/- 6.08 (0.01s)\n"
     ]
    }
   ],
   "source": [
    "tf.random.set_seed(SEED)\n",
    "np.random.seed(SEED)\n",
    "results = []\n",
    "for a,b,c,d in ABCD_CFGS:\n",
    "    for nl, no in NSAMPLES_CFGS:\n",
    "        start = time.time()\n",
    "        (e_mean, e_std), (var_mean, var_std) = empirical_moments_mean_std(a,b,c,d, nl, no, nruns)\n",
    "        r = (a,b,c,d,nl,no, e_mean, e_std, var_mean, var_std, ((time.time()-start)/nruns) )\n",
    "        print(\"a=%.2f\\tb=%.2f\\tc=%.2f\\td=%.2f\\tnl=%i no=%i => E=%.2f +/- %.2f \\tV=%.2f +/- %.2f (%.2fs)\" % r)\n",
    "        results.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "path = RESULTS/nsamples_analysis_25_1000.csv\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>c</th>\n",
       "      <th>d</th>\n",
       "      <th>nl</th>\n",
       "      <th>no</th>\n",
       "      <th>Em</th>\n",
       "      <th>Es</th>\n",
       "      <th>Vm</th>\n",
       "      <th>Vs</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>10.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>1.0</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>2499.827486</td>\n",
       "      <td>23.148860</td>\n",
       "      <td>54750.869583</td>\n",
       "      <td>161514.534783</td>\n",
       "      <td>0.007367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>10.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>2500.435710</td>\n",
       "      <td>7.288014</td>\n",
       "      <td>54497.971288</td>\n",
       "      <td>51442.509271</td>\n",
       "      <td>0.017813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>10.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>1.0</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>2501.717178</td>\n",
       "      <td>70.489039</td>\n",
       "      <td>68000.311104</td>\n",
       "      <td>492757.247793</td>\n",
       "      <td>0.006772</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>10.000</td>\n",
       "      <td>2.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>625.007992</td>\n",
       "      <td>5.822493</td>\n",
       "      <td>4651.163930</td>\n",
       "      <td>10143.963276</td>\n",
       "      <td>0.007581</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>4</td>\n",
       "      <td>10.000</td>\n",
       "      <td>2.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>625.070476</td>\n",
       "      <td>1.800209</td>\n",
       "      <td>3858.566589</td>\n",
       "      <td>3320.979347</td>\n",
       "      <td>0.025922</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>5</td>\n",
       "      <td>10.000</td>\n",
       "      <td>2.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>624.241588</td>\n",
       "      <td>17.752174</td>\n",
       "      <td>2700.045504</td>\n",
       "      <td>32320.389505</td>\n",
       "      <td>0.006303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>6</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.1</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>0.203531</td>\n",
       "      <td>1.135151</td>\n",
       "      <td>82.642368</td>\n",
       "      <td>835.587793</td>\n",
       "      <td>0.007328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>7</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.1</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>0.251504</td>\n",
       "      <td>0.455077</td>\n",
       "      <td>160.799738</td>\n",
       "      <td>1088.810950</td>\n",
       "      <td>0.024090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>8</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.1</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>0.171245</td>\n",
       "      <td>2.271866</td>\n",
       "      <td>132.947937</td>\n",
       "      <td>8500.371151</td>\n",
       "      <td>0.004447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>9</td>\n",
       "      <td>0.100</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>0.248745</td>\n",
       "      <td>0.053279</td>\n",
       "      <td>0.559495</td>\n",
       "      <td>0.390089</td>\n",
       "      <td>0.005496</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>10</td>\n",
       "      <td>0.100</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>0.249884</td>\n",
       "      <td>0.017829</td>\n",
       "      <td>0.555999</td>\n",
       "      <td>0.136605</td>\n",
       "      <td>0.022028</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>11</td>\n",
       "      <td>0.100</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>0.256230</td>\n",
       "      <td>0.196088</td>\n",
       "      <td>0.560474</td>\n",
       "      <td>2.031716</td>\n",
       "      <td>0.004026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>12</td>\n",
       "      <td>0.100</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>24.942095</td>\n",
       "      <td>5.546975</td>\n",
       "      <td>2859.720706</td>\n",
       "      <td>3158.149140</td>\n",
       "      <td>0.007417</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>13</td>\n",
       "      <td>0.100</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>25.096340</td>\n",
       "      <td>1.759973</td>\n",
       "      <td>3012.598062</td>\n",
       "      <td>1298.915026</td>\n",
       "      <td>0.022739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>14</td>\n",
       "      <td>0.100</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>24.961940</td>\n",
       "      <td>19.283487</td>\n",
       "      <td>2440.084548</td>\n",
       "      <td>6949.702551</td>\n",
       "      <td>0.005927</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>15</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>25.022691</td>\n",
       "      <td>2.181339</td>\n",
       "      <td>555.795682</td>\n",
       "      <td>342.219168</td>\n",
       "      <td>0.006716</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>16</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>24.989992</td>\n",
       "      <td>0.709236</td>\n",
       "      <td>549.233654</td>\n",
       "      <td>102.885052</td>\n",
       "      <td>0.018581</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>17</td>\n",
       "      <td>1.000</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.1</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>24.807117</td>\n",
       "      <td>7.039204</td>\n",
       "      <td>525.085432</td>\n",
       "      <td>1136.675756</td>\n",
       "      <td>0.006976</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>18</td>\n",
       "      <td>1000.000</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>100</td>\n",
       "      <td>100</td>\n",
       "      <td>25.000428</td>\n",
       "      <td>0.052688</td>\n",
       "      <td>25.140722</td>\n",
       "      <td>3.795858</td>\n",
       "      <td>0.008254</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>19</td>\n",
       "      <td>1000.000</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>1000</td>\n",
       "      <td>10</td>\n",
       "      <td>24.999495</td>\n",
       "      <td>0.051124</td>\n",
       "      <td>25.026953</td>\n",
       "      <td>3.569548</td>\n",
       "      <td>0.022026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>20</td>\n",
       "      <td>1000.000</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.00</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>10</td>\n",
       "      <td>1000</td>\n",
       "      <td>25.000333</td>\n",
       "      <td>0.087162</td>\n",
       "      <td>24.837115</td>\n",
       "      <td>6.079079</td>\n",
       "      <td>0.006753</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           a        b        c       d    nl    no           Em         Es  \\\n",
       "0     10.000     1.00    10.00     1.0   100   100  2499.827486  23.148860   \n",
       "1     10.000     1.00    10.00     1.0  1000    10  2500.435710   7.288014   \n",
       "2     10.000     1.00    10.00     1.0    10  1000  2501.717178  70.489039   \n",
       "3     10.000     2.00    10.00     2.0   100   100   625.007992   5.822493   \n",
       "4     10.000     2.00    10.00     2.0  1000    10   625.070476   1.800209   \n",
       "5     10.000     2.00    10.00     2.0    10  1000   624.241588  17.752174   \n",
       "6      0.001     0.01     0.01     0.1   100   100     0.203531   1.135151   \n",
       "7      0.001     0.01     0.01     0.1  1000    10     0.251504   0.455077   \n",
       "8      0.001     0.01     0.01     0.1    10  1000     0.171245   2.271866   \n",
       "9      0.100     1.00     0.10     1.0   100   100     0.248745   0.053279   \n",
       "10     0.100     1.00     0.10     1.0  1000    10     0.249884   0.017829   \n",
       "11     0.100     1.00     0.10     1.0    10  1000     0.256230   0.196088   \n",
       "12     0.100     0.10     0.10     0.1   100   100    24.942095   5.546975   \n",
       "13     0.100     0.10     0.10     0.1  1000    10    25.096340   1.759973   \n",
       "14     0.100     0.10     0.10     0.1    10  1000    24.961940  19.283487   \n",
       "15     1.000     1.00     0.10     0.1   100   100    25.022691   2.181339   \n",
       "16     1.000     1.00     0.10     0.1  1000    10    24.989992   0.709236   \n",
       "17     1.000     1.00     0.10     0.1    10  1000    24.807117   7.039204   \n",
       "18  1000.000  1000.00  1000.00  1000.0   100   100    25.000428   0.052688   \n",
       "19  1000.000  1000.00  1000.00  1000.0  1000    10    24.999495   0.051124   \n",
       "20  1000.000  1000.00  1000.00  1000.0    10  1000    25.000333   0.087162   \n",
       "\n",
       "              Vm             Vs      time  \n",
       "0   54750.869583  161514.534783  0.007367  \n",
       "1   54497.971288   51442.509271  0.017813  \n",
       "2   68000.311104  492757.247793  0.006772  \n",
       "3    4651.163930   10143.963276  0.007581  \n",
       "4    3858.566589    3320.979347  0.025922  \n",
       "5    2700.045504   32320.389505  0.006303  \n",
       "6      82.642368     835.587793  0.007328  \n",
       "7     160.799738    1088.810950  0.024090  \n",
       "8     132.947937    8500.371151  0.004447  \n",
       "9       0.559495       0.390089  0.005496  \n",
       "10      0.555999       0.136605  0.022028  \n",
       "11      0.560474       2.031716  0.004026  \n",
       "12   2859.720706    3158.149140  0.007417  \n",
       "13   3012.598062    1298.915026  0.022739  \n",
       "14   2440.084548    6949.702551  0.005927  \n",
       "15    555.795682     342.219168  0.006716  \n",
       "16    549.233654     102.885052  0.018581  \n",
       "17    525.085432    1136.675756  0.006976  \n",
       "18     25.140722       3.795858  0.008254  \n",
       "19     25.026953       3.569548  0.022026  \n",
       "20     24.837115       6.079079  0.006753  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "path = \"RESULTS/nsamples_analysis_%s_%s.csv\" % (K, nruns)\n",
    "print(\"path = %s\" % path)\n",
    "pd.DataFrame(results).to_csv(path, index=False, header=False)\n",
    "df = pd.read_csv(path, header=None).rename(columns=dict(enumerate([\"a\", \"b\", \"c\", \"d\", \"nl\", \"no\", \"Em\", \"Es\", \"Vm\", \"Vs\", \"time\"])))\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.\n",
      "findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAFJCAYAAAA/oeJqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3d32/b9r038Lfak4sUWEQrW/cMAR7E1IqDAQVWU/a5X0Slu10jWb3vRKa7PhHt/QG1qabXtZTmD7CpdrenIdvezxKTAAU2oBbjDcjwrJslOgOWiwLVc+FDVrJ+WKIoiZLfL8BATcrkx3LKj76/Pt9Yu91ug4iIKEJemXcARERE5zE5ERFR5DA5ERFR5DA5ERFR5DA5ERFR5DA5ERFR5PzHvAOIqh//+Me4efPmvMMgIlpax8fH+Oc//9n3HJPTADdv3kStVpt3GERES2t9fX3gOXbrERFR5DA5ERFR5DA5ERFR5DA5ERFR5HBCxBLbffxd3+Nba1dmHAkR0XjYciIioshhciIioshhciIiosjhmBMRRd6LFy/w7bff4rvv+o+jUrRcuXIFr7/+Oq5duxb4GkxORBRpL168wN///nfcuHEDV69eRSwWm3dINES73cbLly/x/PlzAAicoNitR0SR9u233+LGjRt47bXXmJgWQCwWw2uvvYYbN27g22+/DXwdJiciirTvvvsOV69enek9K5UKUqkUkskkkskkYrEYYrEYMpnMTOMYRVRjvXr16kTdsExORBR5s2wxlUolqKoKXdfRaDTQaDRgGAYAQJKkmcUxiijHOunfjMmJiKiDpmmQZRmyLPvHstksZFnGxsbGHCPrtUixjovJiYjofzmOAwBoNps95yRJmntrpNMixRpE4Nl6T548CTMOAIAgCNzgj4jmRhRFAIBt27Btu+sBn8/n/fOTcF0XgiBMfJ1ZxDpPgZOTJEmh9wOnUin88Y9/DPWaRETjyGazqFarSKfTePbsmZ9IwmqJaJqGZDKJYrHY97zrukilUjAM48J7TjvWeZponVM8Hh+6k+E4HMdBu90O5VpEREEZhoFkMgnHcZBKpdBoNEK9frlcRi6XQ6lU6klQXmIql8sjJZhpxzpPgZOTKIo4OjoKMxbcvn071OsREQVhmiZSqRQcx0Emk4FpmqFe3zCMngTVmZg6JzjMO9Z5CTwhYhrNxmQyGfo1iYjGJYoivvjiCwCAZVnQNC30exiGAdM0oWla4MQ0q1jnIVKz9ditR0RRIUkSyuUygLP1RJ0sy/IXuqqqCk3ToGkaVlZWEIvFYFnWSPcwTRO2bWN1dTVQYhol1kUVODk9ePAgzDgAALquh35NIqJReFOzOymK4s966zxvmiZarRZM00S5XIau69jY2IDrutB1feQk47ouHMeBKIqwbXsqsS6qwMkpHo+HGcfUrklEdBHXdf2Wx3neA79zavb29nbXdHDHcZDL5ZDNZgfOwut3T68rr16vwzTNkVo948a6qKbWrffkyRPcv38f9+/fx2effYYXL15M61ZERBNxHAfVarXvuVqt1tMS6kxMrusik8lAkiS/dNBF+o0xmaY5UoIaJ1bbtpHJZLrq7HndkFE3leT04MEDv6xGOp3GyckJisUi3n777aks3iUimkStVoPjOFBVteu49/2glgoApNNpNJtNf1LCRYZNfhglQY0TqzcWZVmW39W3CIkJmFJyqlQq+Pzzz/HWW29hbW0NhUIBe3t7+Pjjj/HBBx/gk08+mcZtiYgC8dYHVSoVxGIxpFIppFIpNJtNPHv2bGA3maqqsG0bX3zxxchVHzRNGzr5wTRNHB4eDpxUMW6soihCkqSu1ta8K5aPItaewhS5n//853j8+DF+9KMf9T2/tbWF27dv49atW2HfOjTr6+uo1WrzDmMiu4/7l6vfWrsy40iIgvvTn/6EX/ziF/MOo0elUoGqqjAMA9lsdt7hDFUqlfxWWbVanVm8F/3thj1nA7WcNjY28PTp04Hn0+k00uk0/vWvf/U9v7u7y5l5RLSwbNuGqqooFos9D/pxZt3NiizLftfeokyWCJScbt26BUmSkM/ncXx83HN+a2sLR0dHWF1dxUcffdT3Gv0q6RIRRZ3jOEin05BluedDdrVajeQ0bq9ogmVZC1N3L1D5Il3XcXp6CkEQkM1msbGxAV3X/b3iV1dXcXBwgNu3b6NYLOKDDz7w9xcRRRH7+/twXTfUX4SIaBYymYxfWdybhNBsNuE4DmzbjmwxAVmWoSjKvMMYWeDaerquQ1EU1Go1VCoVSJKEXC6HnZ0dAGdvRK1WQ6FQwOPHj2EYRteA3LLUfyKiy2URi6s6jrNwQykTLcL1PiEoioKjoyOsrKzgjTfewMOHDwGcNSXr9Trq9Tp0XUehUMC9e/dwdHQU6ckQRETLwHEcuK7bs9/TIphoKvn169e7vi8Wizg8PMTh4SHeeOMN/OEPfwAArK2t4d69e9jb28Pu7i5WV1cnuS0REV3Atm2kUikcHBxEfjZhP4G79Z49ewbTNJHP55FIJJBKpbC5uQlBELC3twfHcXD37l188MEHKJVK+NWvfhVm3ERENIQkSWi1WvMOI7DALaetrS1omob9/X2/NZTNZnH//n0AZwu/Hj16hJ2dHRQKBbz99ttDp58TERF5Aicn27ZRKBQAnI0/pdNpPHr0CPF4HO+//77/OlmWcXR0hDt37uBXv/oV8vk8/vKXv0weORERLa3AyandbvddZFsoFBCPx/HZZ591HVcUBc1mEzdv3sTq6ip+97vfBb01EREtucDJSZZl5HK5vtXGt7e3BxZK1HUdzWYT//znP4PemoiIllzg5KRpGh49eoT19XV/6ninYXXpBEHAwcFBoPu6rotcLoeVlRWsrKwgk8kMLReiaZq/W2Umkxl5h0oiIpqfwMnJqwJxdHQERVFw/fp1vPvuu9je3sb6+jrW19fDjBPA2Zz91dVVuK7rV/S1LAupVKpvgkqlUrAsy9+tslwuI5PJDNwLhYiIoiHwVHIAyGazqNVqyOVyePbsmd8aEgRhKhUgVFXFF1980bWYLJfLoVqtolAooF6v+8dLpRJs2+46JooiisUicrkcWq3WyCXuiYhotibez0mSJDQaDTQaDRiGAcMw/IkPYXIcx99tspO38+T5ltPOzg4EQeh5vbePSaVSCTU+IqLLwHXdmdRGDW2zwdXVVdy5cwd37twJ65JdEokEisVi33OiKHa1grySHf26Fr1jw3a2JCKaN2+83Cs0GwWO4yCVSs1kr7vAyen8VPEwDLvmsC44x3Gwubnpf++1ovr9jHcsimXtiYiAHwq1mqYJURQjM5HL21V3FgInp/39/TDjABCsUrllWRAEoavirpd4EolEaLEREU1q1BaQtyGg9/pp1saLSqvsvNC69cIQZANCTdNgGEZXK+miiQ7e+fN/lEql4s80/Mc//jF2LERE/ViWhWQy2TPW7bouVFVFqVSCpmldLSTXdbGzs4ODg4OpzDAOEtMsBZ6tZ1kWvvrqq9AqjHtbCI+jUqkgn8/708o93qeOQcnOS0rnk5iiKP5mXNOYCk9Ek/v666/nHQIA4M033xzpdYPGv4GzijqZTMZ/7iSTSRiGAUmS/B6hfD4PTdNCbT0FjWmWAienVqvVkxQmNc4vb1kWGo1G3w20vDe9X3PVO7Zoe5sQ0WLyPgSf/zBs2zaq1SoePHjgH5NlGTs7O/4sZO/nwn5eTRrTLAROTmtra4jFYmHG4rd4LmLbNkzTHLizoyAIEEWx74wS71jYiZWIaBy1Wg2CIHQliFQqBU3TUK1WUS6Xoaoqms3m0F1sS6USRFHsaVmpqgpN00Z+rl4UEwB/K3rTNLG+vj7VtaKBk1Pn4tZZsm0b+/v7ff9Ymqb5x70/jOM4XX8cb9KFqqqzCZiIqI96vd4zaSuRSMB1XWSz2ZG78YrFop88vJ9RVRWqqo6VmC6KCThrQMxqm/qJKkTMmm3bSKfTftHZTo7jdLWGisUi9vf3oaqqn5Acx0GpVIKu62P/0SLt+V8HnPjZTMMgmoVRx3qiLsxWh67rfoIyTROqqgbqCoxS1ZyFSU6O4yCdTsN13YEzV873idbrdWiahlwuh0QiAcdxYBjGQm5ZTETLJZlM9kzaajabgT8467qOZDIZODFNI6ZJLExyEkUx0JbDw/pqiYjmRZblnklb9Xo9cGJRVRWGYWB/fx/VajXQh/CwY5pEpNY5EREtq/MPfVEUIctyV0+QZVmBPlB7Y0ySJEHXdRweHo60NmqaMU1qYVpORESLqlKpwHEc7O/vQ5Ikf3zcMAx/4pa3NGbcLrRSqdTTleeNQQ0rNzTNmMIQa7fb7ZnfdQGsr6/PpLhhKAZMiNj9tv+EiK21K9OMhihUf/rTn/CLX/xi3mFQABf97YY9Z9mtR0REkcPkREREkTOz5LS9vY2NjQ387ne/w1/+8pdZ3ZaIiBZQqMlpY2MD169fx9tvv40vv/zSP/7++++jVCqhXq9jb28Poiji6dOnYd6aiIiWSKjJqdFooFKp4PPPP8etW7cAAJ988om/62ylUsH333+Pzz//fOCutkRERKEmJ1mWu7ZpPz09haIoiMViKBaL+O1vf+u/LqytNoiIaPmEmpyuX7/e9b1X60kQBOzs7HSdC1LtgYiILofQu/W8yQ6ffPIJKpUKYrFY194gwFmLal67KxIRUfSFWiGiWCxidXUVKysrflkMRVHwzjvv+K958uQJfvvb30Z233oiIpq/UJOTLMuo1Wool8totVrI5/NdY1B37971K9wu1ZYVREQhU1UVlmVBEAQYhhGJZ+YsY5pL+aLj42PcvHlz1rcdC8sXEUXDZSxfZFmWv9OstzGqN+t5kWJauPJFmqbhD3/4wzxuTUQ0N6MOZ8iy7G/8l8lkkEwmL1VMwJSqkh8fH8NxnL7nXNdFvV7H48eP8Zvf/GYatyciihTLsvxtLTrXeLquC03TkEwmcXJygkwm07WjNwAcHh5ie3v7UsTUKdTkdHp6ClmWYdv20Ne1222srKyEeWsiuiwGdGPP3I3/O9LLXNfF+vp633OFQgGZTAaKogA424nWMAx/mwvLspDP50PfPj2KMZ0Xareepmmo1+tYW1tDOp1GPB5HOp3u+orH48jlcpxKTkSXgiAI/lcn27ZRrVaxubnpH5Nl2V8Tatt2135MYc5wjmJM54XacqrVami1WojH4wCABw8eYGNjA2+99Zb/Gtd1oaoq1tbWwrw1EdFCqdVqPQkilUpB0zRYloVcLodEIgHgLJnU6/W+1ymVShBFsWdbdlVV/Q0HZx1TGEJtOYmi6Ccm4Kx5eH42hyAIaLVaePjwYZi3JiJaKPV63X/QexKJBFzXhSzLaLVaaDQaaDQaQ5NAsVjs2ZbdG0sad6p3WDGFIfQxp/Pi8TiePn2KX/7yl/4xURSxt7eH9957L8zbE9FlMOJYT9SFOWbjbcsOAKZp9mzbPo+YJhVqclpbW0M+n0cikUAul8OtW7ewtbUFURTx4YcfIpVKwbIsv6wREdFllUwm0Ww2u455RQqC0HUdyWQycGKaRkyTCDU5bW9v4+bNm37tvG+++QaCIGB3dxeFQqErIZ2fmkhEdJnIstwzoaBerwdOLKqqwjAM7O/vo1qt9oxBzSOmSYQ65hSPx3F8fAxd17G3t+cfVxQFu7u7uHbtGtrtNiRJgmEYYd6aiCjSzj/0RVGELMtdY0WWZUHX9bGv7Y0xSZIEXdd7xqDmEdOkQl+EG4/Hce/evZ7jxWKRGwwS0aVUqVTgOA729/chSZLfc2QYBjRNg+M4aDQa0HV97C60UqnU05XnjUF1TvueZUxhmEttPSD69fVYW48oGi5jbb1lsXC19YD5FzEkIqLoCtStd3x8DMuysLm5iWvXrvnHR6215DgOqtVqz+64REREQMDkJMsynj17hkaj0ZVgTNPE48ePL/z5drvNqeRERDRQoOSUTqdhGAYymUzXcUVRcPfuXciyPHQArdFo4MsvvwxyayIiugQCJadyudx3zCifz8OyLBwcHFx4jbt37wa5NRERXQKhr3Maddwpl8uFeWsiIloioc/WG7XaeDqdDvvWRES0JOY2lZyIiGiQqSanFy9eIJ/P4/r167h+/Tp+/etf46uvvprmLYmIlobrulPd0C/KJkpOL168wNtvv41XX30Vr776Kv7rv/4Lx8fHAIBnz55hZWUF1WrVf4NN04Qsy/j9738fRuxEREvLcRykUqlIVarRNA2ZTAaZTGbqSXOi5CRJEizLQrvdRrvdRq1Ww+3btwGcTXjwjp//0nUdT548CeUXICJaRsPq4s2D4zjQdR2maUIURViWNdX7BU5On376KRzHQTweh6IoKBaLWFtbQ6PRwP3792HbNlKpFOr1Or7//nt8//33qNfrKBQKaLfbqFQqYf4eRESRF8UuulFj8taueq8PsiXHOAInp/39fSSTSTx79gx7e3vY3d1FvV7HO++8A03TsLKygi+++KJr9t7a2hrK5TLS6TQODw9D+QWIiKLOsiwkk8meD+Wu60JVVZRKJWiaNvXWyKQxua6LnZ0dHBwcjLQlxyQCb5lh2za2trYQj8e7juu6jk8//RSKonTV3eukaRry+XzQWxPRJbb7+Lt5hwBg9Or+rutifX2977lCoYBMJgNFUQCc7URrGMbUu/OCxiQIAnRdRz6fh6ZpU209BW45tVotbG5u9hz3mn7nSxt12tjYiGTzlogobIIg+F+dbNtGtVrteo7KsjyTgtiTxiQIwtQTaODk5LruwJaRIAhIJBIDf/Z8a4uI6LKp1Wo9CSKVSvndaI7jwLZtmKY59MN8qVTq28Wmqiocxwktpmq1ikwmg2q1OpPdcUPfCddzPiN3Oj09ndZtiYgWQr1e7/kQn0gk/EQkiiIajcaF1ykWi9A0DcAPkxS8bdvH3cF2WEzZbHbqkyA6BU5O7XYbv//97/v+8q7rolwuI5lM9j1nmmbQ2xLRJbcsOzkP+wA/Lm9bduBs66Lz27bPI6ZJTdRyGtSsi8ViKJVKA3+O+zkR0WWXTCbRbDa7jjWbzbFbOx5d15FMJgMnpmnENImJFuEOWmR70RcR0WUny3LPWFK9Xg+cWFRVhWEYODk5CTzNO+yYJjFRy6larY5chbxTo9HgVHIiulTOP/RFUYQsy6hWq/5YjmVZgYY9vDEmSZIgSVLPGNQ8YppU4OSUzWbxzjvvBPrZ1dVV7udERJdGpVKB4zjY39+HJEmQZRkAYBgGNE2D4zhoNBrQdX3sLrRSqdTTleeNQQ0rgTTNmMIQawfsZ/v0009x586dwDf+4osvIr2n0/r6eqQKLg71/K99D+9++7O+x5dlQJkuhz/96U/4xS9+Me8wKICL/nbDnrOBx5wmSUwANxskIqLBuNkgERFFDpMTERFFDpMTEUUel6Asnkn/ZkxORBRpV65cwcuXL+cdBo3p5cuXuHIl+OSrqdXWozODyvtzxhzRaF5//XU8f/4cN27cwNWrV1ldJuLa7TZevnyJ58+f46c//Wng6zA5EVGkebsf/O1vf8N330VjLyca7sqVK/jpT386cOeKUSx9ctI0DbZtQxRFOI4DTdP8xWZEtBiuXbs20YOOFs9SJ6dUKgXgrDYUcLY/irer4yxLvxMR0XiWNjmVSiXYtu0nJuCsblSxWEQul0Or1YpUeXgiIvrB0s7W29nZ6buVsLd9fKVSmUdYREQ0gqVsOTmOA9d1+44tra+vAwDK5TKKxeKsQ6MZ4SxJosW2lMnJtm0A/Xd19I45jhPeDQcUXj3Tv/jquL7++uuB595cmf5A8aD7v/nmm1O/NxFdPkuZnLzEk0gk5hzJJTBmRXSArZdpY6uRlsFSJqeLJjoIggDXdeG6btdrK5WKPxb15z//2e8CnIZg+1Qu/r3nff95/+7zdJl/d4qm4+PjgeeWMjl5G2M1m82+573dH88nMUVRoCjKdIMjIqILLeVsPa/Fc34L4s5jg3aHJCKi+VvK5CQIAkRR7LvDoneMVSKIiKJrKZMTAKiqCtd1e2blmabpnyciomiKtZd4o5RUKoVEIuEnJK98ka7rXONERBRhS52cgLPCr47jIJFIwHEcqKrKunpERBG39MmJiIgWz9KOORER0eJiciIioshhciIioshhciIioshhciIioshhciIioshhciIioshhciIioshhciIioshhciIiosgZutngkydPQr+hIAi4efNm6NclIqLlMbS23iuvvIJYLBbqDVOpFP74xz+Gek0iIlouF27THo/H/Z1lJ+U4DlhnloiILjI0OYmiiKOjo1BvePv27VCvR0REy2fohAhJkkK/YTKZDP2aRES0XC7s1gvbonTr/fjHP174iRv/79/93+v/81q444jj3HtW9yei6Ds+PsY///nPvueGTog4PT1FPB4PNZhpXHMa1tfXUavV5h3GRHYff9f3+Nbalbnde1b3J6LoG/acHdqtN40ksgiJiYiI5muibr0nT57AsiwAZ5MnZFnGtWvXQgmMiIgur8AVIh48eABN0yDLMtLpNE5OTlAsFvH2229PZfEuERFdHoFbTpVKBYeHh/73a2trKBQKcBwHW1tbuH37Nn7729+GEiQREV0ugVtOrVYL//rXv3qOi6KIg4MDHB0d4csvv5woOCIiupwGJqeNjQ08ffp04A+m02mk0+m+CQoAdnd3oev65BESEdGlMzA53bp1C5IkIZ/P4/j4uOf81tYWjo6OsLq6io8++qjvNZrNZmiBEhHR5TFwzEnXdZyenkIQBGSzWWxsbEDXdX823urqKg4ODnD79m0Ui0V88MEHkGUZGxsbEEUR+/v7cF13Zr8IEREtj6ETInRdh6IoqNVqqFQqkCQJuVwOOzs7AABZllGr1VAoFPD48WMYhoFqter/vGma042eiC6FFy9e4Ntvv8V33w1e3E3RceXKFbz++usTLS0ampzi8bhfbkhRFCiKglKphDfeeANbW1t47733IEkS6vU6Hj9+DMuy0Gg0IAgCVFXF6upq4MCIiICzxPT3v/8dN27cwNWrV0PfxofC1W638fLlSzx//hwAAieoC6eSX79+vev7YrEIRVGwtbWF3d1dlEol/OY3v8Ha2hrW1tYCBUETev7XASd+NtMwiKbh22+/xY0bN/Daa6/NOxQaQSwWw2uvvYYbN27gb3/7W+DkNHQq+bNnz2CaJvL5PN5//3188sknePHiBQRBwN7eHj7//HN8/PHH2NjYwFdffRUoACKiYb777jtcvXp1pvesVCpIpVJIJpNIJpOIxWKIxWLIZDIzjWMUUY316tWrE3XDDk1OW1tb0DQN+/v72N3dxerqKrLZLO7fvw/gbE3To0ePsLOzg0KhgLfffnvo9HMioiBm2ZVXKpWgqip0XUej0UCj0YBhGACms43QJKIc66R/s6HJybZtFAoFAGfjT+l0Go8ePUI8Hsf777/vv06WZRwdHeHOnTv41a9+hXw+j7/85S8TBUZENA9eWTZZlv1j2WzWn40cJYsU67iGJqd2u913kW2hUEA8Hsdnn33WdVxRFDSbTdy8eROrq6v43e9+F260RERT5DgOgP5rNCVJmntrpNMixRrE0OQkyzJyuRxevHjRc257exvlcrnvz+m6jmazOXATKSKiKBJFEcBZr5Ft213n8vm8f34SYa3/nEWs8zQ0OWmahkePHmF9fR0PHz7sOT9sMz5BEHBwcDB5hEREM5TNZgGclWjrTCRhtUQ0TUOpVBp43nVdJJPJnoTTz7RjnaehycmrAnF0dARFUXD9+nW8++672N7exvr6OtbX12cVJxHRTBiGAVEU4bouUqlU6Ncvl8s4PDzsm6C8e5bL5ZESzLRjnacLq5Jns1nUajXcvHkTrVYLBwcH0HUdJycnA7v1iIgWmWmaEAQBjuNMZUq2YRg9CaozMXVOcJh3rPMy0pYZkiR1TVM0DMOf+EBEtGxEUcQXX3wBALAsC5qmhX4PwzBgmiY0TQucmGYV6zyMtZ/T6uoq7ty5gzt37kwrHiKiSJAkye8dOt8FZ1mWv9BVVVVomgZN07CysoJYLAbLska6h2masG0bq6urgRLTKLEuqqHJ6fxU8TBM45pERJPypmZ3UhTFn/XWed40TbRaLZimiXK5DF3XsbGxAdd1oev6yEnGdV04jgNRFEeaABEk1kU1NDnt7++HfkNWKieiqHFdd+AYuvfA75yavb29DUEQ/O8dx0Eul0M2m0WxWBz5nl5XXr1eh2maI7V6xo11UQXepj0obkBIRFHjOE7Xdj+darVaT0uoMzG5rotMJgNJkvzSQRfpN8ZkmuZICWqcWG3bRiaT6Zoo4XVDRt3QquSWZeGrr74KbesLy7KWorlJRMulVqvBcRyoqtrVKlFVFQCGzkxOp9NoNpuo1+sj3WvY5AfTNP1EMqgFNk6s3lhUMpn0uw81TVuImdZDk1Or1Qo8QDfIMiwOI6Ll0mg0AJxV+PY2VgXOuseePXvW1VLqpKoqbNtGvV4f+JrzvOQw6NlqmiZyuRwsy+r7mnFjFUURkiShWq36CW8RppwPTU5ra2uhVwNehr5Qir7dx/1L9W+tXZlxJJfLor7vuq5D1/WxfsZLDoZhjPWhe5RWy7DuwSCx5vN5mKaJYrEI27b9yhJRNjQ5jdpMJSK6TGzbhqqqKBaLPQ9627Yj10MkyzI0TfO79hbBzCdEzJLrusjlclhZWUEymUSlUpl3SES04BzHQTqdhizLPS2YarUayXF1L1lalhW5xDnIhdu0L7JcLgdN0/wK6qqqIpFILESTloiiKZPJwHVdCILgT0JoNptwHAe2baPdbs85wv5kWYaiKPMOY2SRTE6lUgmSJE00GaNarULXdf9TQrlchmVZME2TyYmIAvMmJCwSx3HGHqeat0h161UqFaysrPi1pobRNM0vHZLJZHrKhWSz2Z7mqyAIS1e5l4hoEMdx4LpuJMfBLhKZlpNt29jc3ES9Xr9wbMhLMN6EDcdxkEwmYRjG0FaR4zjY3NwML2giooiybRvpdBq6ri9Ud54nMsnJy+rJZHLo60qlkr+uwCOKIorFInK5HFqtVt/1BpqmwTCMkdciEM1La78AAB84SURBVBEtMkmS0Gq15h1GYJHq1hvFzs4OBEHoaaJ6i8r6tbosy0IymQx9QTEBeP7X/l9ERBNYqOTk9Z/224HXO3Z+gZtt23BddyGbtUREl9XEyen4+Nj/72fPnk11SwyvpHy/rjnvWOcaA9u2/Xn9juPAcRxUKpWxStMTEdHsBUpOx8fHyOfzePXVV7taMaurq2i329jY2MC//vWv0IL0eIknkUhc+FrbtpFKpaBpGpLJpP/VOb2ciIiiaewJEaenp0ilUnBdF+12u2fB2Z07d3BycgJJkvDNN9+EFijQv8V0/rzrunBdF5Ikjb0YzquVBQD/+Mc/AsdJRESTGbvlVCgU0Gq1/CmK/VoxiqKg0Wjgo48+CiVIj1cTatCeUN7aqKAz8hRFQa1WQ61Ww09+8pNgQRIR0cTGTk62bcMwDDx69Aj37t0b+DpRFLG3tzdRcOd5XYj9Fuh6x9hlR0S0+Mbu1ms2m7hz585Irzs9PQ0U1CCCIEAURdRqtZ5z3rFlnS7+9ddfDzz35sq1GUZCRJfZpD1Uoxq75bS+vo4XL14Mfc2DBw/guu5USrOrqgrXdXsq/5qm6Z8nIlp0Xok2r9BsFDiOg1Qq1beBELaxk5OmaSgUCgPPf/LJJ1BVFbFYLFAr5uTkBMDgcaVisQhJkrqSkOM4KJVK0HV9YfYqISIaxCvUapomRFHsqR06L96uurMwdrdeOp3GwcEBrl+/DlmW0Ww2sb29DcdxYFmWn+FFUcTHH3888nUdx0G1WvVny5XLZbiui2w225Nw6vU6NE1DLpdDIpGA4zgX1tUjIpo3b6uNi3jPPO95Os1n26gxzVqg2nrlchmpVApbW1twXbenFLuiKGOXZ/fq43l73F9k0cq/E9HlZVkWVFX1d8/1uK7rr8U8OTlBJpPxe5xc18XOzg4ODg6QyWRCT1BBYpqlwIVfFUWBoih49uyZX1ZIFEWIooh4PB5mjEREC2tQyTXgbGlOJpPxy6t5uytIkgRBEKDrOvL5PDRNCzU5BY1pliauSr66uorV1VUAZwt0W60WkxMRTc2wmauz9Oabb470Oq/L7HzXmW3bqFarePDggX9MlmXs7OzAMIyunw87MUwa0ywEKl909+5dvPrqq3j48GHX8Xg8jt3dXWxvb4cSHBHRsqrVahAEoStBpFIpWJaFarWKTCaDarUKy7KGDmOUSiVUq9We46qq9sxqniQmAP5W9KZpTn0G4dgtpwcPHviTFo6OjnrO7+3t4e7du/j1r3+N//mf/5k8QiIay+7j7/oe31q7MuNIaJh6vd5TYSeRSPgTwUbtxisWi9A0DcAPEye8saRxZy8Piwk4mxswq23qx05O5XIZsixD0zSk0+m+r9nb28Mrr7yChw8f4r333ps4SCKiZRPmDDld1/0EZZomVFUN1BUYpVl7Yycn13VHWoAlCALK5TKTExGFatSxnqhLJpM96zmbzWbgtZq6riOZTAZOTNOIaRJjjzmNEuTp6Slc1+3aSp2IiH4gy3LPuE29Xg+cWFRVhWEYODk56TsGNY+YJjF2chIEAU+fPh36Gq95ySKsRERnzj/0RVGELMtdieSiyQ+DeGNMkiRB13UcHh6OlKCmGdOkxk5Ou7u7yGaz+Oqrr3rOvXjxAu+//z4qlUrg8kVERMumUqnAcRzs7+93lSIyDAOmaaJUKkFV1UAl2Lyf7WwMeAlq2K7f04wpDGOPOYmiiHv37iGdTiOZTPpviDfF0CNJEnZ2dsKLlIhoQXlFC87zxuYnMaiqzkWtnWnGFIZAi3AVRYEoin4f53mapjExERFRYIErRMiyjEaj0VO+aG1tLcz4iIjoEgq1fNF5L168wLVr3AiPiIjGE6h80aiG7ftEREQ0SKCW0/HxMcrlst+d14/jOGPXdSIiIgICJKdPP/0Um5ubAIB2uz30tbFYLFhURER0qY2dnDRNQ7vdhizLyGQyA2sxHR0d4f79+xMHSEREl8/YyanZbCKTyeDzzz+/8LVsORERURBjJydZlnH9+vWRXru7uzt2QES0AJ7/dcCJn800DFpegcoXHR4ejvTaTz75ZOyAiIjorF5eMplEKpWKzOSyWcY0dnL68Y9/jO3tbeTzeRwfHw/8+vLLL+dSLJAumed/7f9FtMC8YquNRgPr6+uReJbOOqaxu/Vu3ryJ09NTAAhclp2I6DJyXXekDf06i2ZnMpmptlKiGBMQoOVUKBTQbrdH+iIiorNWRzKZRKVS6Truui5UVUWpVIKmaV3VwT2Hh4d9C7QuY0ydxm45vfvuu3j8+DEMw0A8Hh/4OsdxcPfu3YmCIyJadK7rYn19ve+5QqGATCbjP+iTySQMw/B3e7AsC/l8PvTt06MY03ljJ6e1tTVomjY0MQFnW2tEoZ903nYff9f3+NbalRlHQrQkojKmeOP/jvQy7yF+/mFu2zaq1SoePHjgH5NlGTs7OzAMA7ZtQxRFfy+lUbvfFjWm8wLV1kun0xe+5vj4mBXKiYgGqNVqEASh6+GeSqVgWRYsy0I6nUYmk/Fnxw1KAqVSqe/4v6qqY48LhRVTGCauSj6Ipml499138Zvf/GZatyAiWlj1eh2JRKLrWCKRgOu6kGUZrVZrpOsUi0VomgYAyGazAH7Ytn3cHWzDiikMgZPT8fHxwKzsui7q9ToeP37M5LQkvv76677H31zhlihEQYTZ6tB13U9Qpmn2bNs+j5gmNXZyOj09hSzLQ/emB86Kwq6srAQOjIiorxHHeqIumUyi2Wx2HWs2m2O3djy6riOZTAZOTNOIaRKBCr/W63VIkoREIoFardYz66NWqyGTyWBrayu0QImIloksyz1bDnnP1iBUVYVhGNjf30e1WvW7+OYZ0yTGTk61Wg2tVsufrffgwQNsbGzgrbfe8l/jzZPnhAgiojPnH/qiKEKW5a5EYlkWTNMc+9reGJMkSZAkqWcMah4xTWrs2XqiKHZNIy8UCiiXy12vEQQBrVYLDx8+nDxCIqIFV6lU4DgO9vf3uxa1GoYB0zRRKpWgqip0XR+7C8372c7Wja7rODw8HDr8Ms2YwhBozOm8eDyOp0+f4pe//KV/TBRF7O3t4b333pssQiKiBacoSt+KCoIg9Hy4H1exWOx7/KJ1ptOMKQyBFuHm83kkEgnkcjncunULW1tbEEURH374oT8nvlKpcD8nIiIKZOzktL297Rd/tSwL33zzDQRBwO7uLgqFQldC6iwUSERENKqxk1M8Hsfx8TEqlUpXH6eiKHBdFzs7Ozg9PYUkSTAMI9RgiYjocgi0CDcej+PevXs9x4vF4sD+TyIiolEFqq03quPj42lenoiIltRUk1MUZnwQEdHiGditd3x8DMuysLm5iWvXfqiftr29PdKFHcdBtVrFzs7O5FESEdGlMjA5ybKMZ8+eodFodCUY0zTx+PHjCy/cbrc5lZyIiAIZmJzS6TQMw0Amk+k6rigK7t69C1mWh64abjQa+PLLL8OLlIh+MHTDvZ/NLAyaLq+8UJSqhc/KwORULpf7jhnl83lYloWDg4MLL85t2omIgnEcB5lMBuVyOTJrRjVN80siGYYx1aQ59oSIVquF//zP/8SLFy8ufG0ulwsUFBHRZSeK4lyqgQ/iOA50XYdpmhBFsase3zSMvc7JG4uKx+P47//+76GvHWU7dyKiy8J13ch10Y0akzeM43U1BtmSYxxjt5xOTk4AoG/BwPNGaV0RES07y7KQTCZRqVS6jnvbC5VKJWiaNvXWyKQxeVWADg4OUK1Wpxrf2C0nwzBw9+7drunlg1QqlQtbV0RE49h9/N28QwAAbK1dGel1ruv2bMjqKRQKyGQy/of9ZDIJwzCm3p0XNCZBEKDrOvL5PDRNm2rraeyWkyzLKJfLyOfzePr06cDXPXv2jItwiejSEwTB/+pk2zaq1So2Nzf9Y7Isz2Rt6KQxCYIw9QQ6dnL6+c9/jtu3b6NarUKSJLz66qt9v37+85/DcZxpxExEtPBqtVpPgvC2HALOJiDYtg3TNHt2rO1UKpX6drGpqjr2M3hYTNVqFZlMBtVqFZZlXbhf1KTG7tZTFAVbW1sjvZaLcImI+qvX60gkEl3HEomEn4hEUUSj0bjwOsVisWdbdm/b9nF3sB0WUzabnfokiE5jJydVVbG7u4svv/wSq6urA193cnIysE+TiCioUcd6oi7MWXu6rvsJyjTNnm3b5xHTpALt56TrOt56662hr2u1WiPX4SMiumySySSazWbXsWazOXZrx6PrOpLJZODENI2YJhGoKnmhULjwNZqmIZlMBrk8EdHSk2W5ZyypXq8HTiyqqsIwDJycnASe5h12TJMItNkgcFa1fNBgm+u6qNfrePz4Md55553AwRERLYvzD31RFCHLMqrVqj+WY1kWTNMc+9reGJMkSZAkqWcMah4xTWrs5HR6egpZlv36SoO0222srKwEDoyIaFlUKhU4joP9/X1IkuTXyjMMA5qmwXEcNBoN6Lo+dhdaqVTq6crzxqCGlUCaZkxhGDs5aZrmN/MSiQRqtVrPxIdarYZMJjPyrD4iomWmKErfqjqCIEy8HrRYLPY9ftFU72nGFIaxk1OtVkOr1UI8HgcAPHjwABsbG10TJLzyF2tra+FFSkREl8bYyUkURT8xAWeTI95//318/PHH/jFBENBqtfDw4UO899574URKRDP19ddfDzz35srF5cuIJhFozOm8eDyOp0+f4pe//KV/TBRF7O3tMTnRxPiQJLp8xk5Oa2tryOfzSCQSyOVyuHXrFra2tiCKIj788EO/1EWlUmGFCCIiCmRocvrss896poJvb2/j5s2bOD09hWVZ+OabbyAIAnZ3d1EoFLoSUlR2bySixdZut/lhd8G02+2Jfn5octI0DZlMBj/60Y/8Y/F4HMfHx6hUKl1TFBVF8ff6OD09hSRJMAxjouDCoKpq12LgQTNbiCiarly5gpcvX+K1116bdyg0hpcvX+LKleClpoZWiGg0GhAEAR999FHXxoHxeBz37t3r2em2WCyi1Wrh+++/R61WG2nPp2lSVRWCIKBYLKJYLOLw8NBfnEZEi+H111/H8+fP8e9//3viT+M0fe12G//+97/x/PlzvP7664Gvc+GY09raGu7du4disYhsNovf//73XRMfpqFUKnUtCgvCcRxUKpWuqr6qqiKTyWB7eztSBQ6JaDDvQ+7f/vY3fPddNDYapOGuXLmCn/70pxM1UIYmJ1mW8ejRIwBnm1Dt7OxgbW0NyWQSd+/eRaFQCLV1VKlUoGkaXNe9sEtQ0zTYtg1RFOE4DjRN60pm3p4onSubvcXCtVqN42FEC+TatWtz74mh2Rraraeqqv/f3hhSq9VCoVDAxx9/jJWVFeTzeXz11VcTB2LbNjY3N7t2YBzEmxFomibK5TLK5bK/CZbH65Ls5H3PTRCJiKJtaHK6c+dOz7F4PI5isYijoyMcHh7i+++/RzqdxhtvvIH79+93jU2Nw9uf/qJK5qVSCbZt48GDB/4xURRRLBaRy+WG7hhJRESLIdCWGR6vNdVsNqEoCvb29rCysoJ33303lNZUPzs7O333r89kMgDOugaBs31Jzicq73tugkhEFG0TJSeP67poNBpwHAftdhuGYUCWZWxsbIRxeZ/jOHBdt29y8Y55BQu9MaXOLrxarQYA4VfYff7XwV9ERDS2iZLTkydPkM/nkUwmu7rZ1tbWcHBwgMPDw4kD7ORt09Fvpt358SRRFKEoStc4lGEYKBaLnKlHRBRxQ2frvXjxou8Mmc8++ww7Ozt+svDWHsiyDE3TetY/hcVLPIlEYqTXl8tlqKqKUqmEk5MTCIIwtIx8pVLxuwX//Oc/T9T9949//AM/+clPBp4Ptk9l9O897/tH/Xef572X+X2nxXR8fDz4ZHuIra0t/79PT0/bH374YTuRSLRfeeWV9iuvvNKOxWLtWCzWVlW17TjOsEuNTNf1NoC2YRg958rlchtAW1GUvj8rCEIbQLvVaoUSyyRSqdSlvPe878/f/fLdm5bT0JZTtVrFxsYG9vf3/e6x9v+2kgRBgKIo2N7e7tpCY5q8saJms9n3vDfhgd12RESLbWhyajQayOVyAH5ISqIoQtM0FAqF6Ud3jtfN1m+6uHds0JbERES0OEbaMqPdbk99PGkUgiBAFEV/1l0n71hUKj/02/74Mtx73vfn73757k1LalifXywWayeTybZlWTPpY2y3h485dZ5vNBpdx4vFYt/jRES0eIZOJRdFEUdHRzNtLZ2cnAAYPK5ULBYhSVJXaSXHcVAqlaDrevhrmIiIaOaGdusNm3YdNsdxUK1W/anc5XIZrusim832JJx6vQ5N05DL5ZBIJOA4DgzDQDabnVm8RPSDVCqFer0+7zBomcy76baoTNNsZ7PZNgD/SxTFtiAIbUEQ2pIktYvFYrter88knkaj0VYUpS1Jkj+l/vxXmAzDaMuy7F/b+50lSWqLotgWRdE/N2jq/yLf31Ov19uKovj36/feh/lvwDCMtqIoXdfv/J0lSWrruj7T5RT1er0NoF0ul2d2T1p+TE4T8h5Kuq53HTcMoy1JUhtAW5blqT4svIeVLMtdD8JGo9GVQKfBu7Zpmn3PF4vFqSaHed7fe98lSep631utlj82Oiy2SXh/1/O/m7cWUBCEmY2/eu+DKIozuR9dDkxOE/IS0KBPjd5EDUEQppKgvPsXi8WBr/EelNO4v9dSGPYAHhbbot7fSw7ZbHbga0zTHDq5ZxJeQuiXeL1/c8NiC0ur1fJbbNNKxHQ5hVL4lQbTdR2yLMN1XX/NWFi87UNEURw6PlgsFiGK4sBJJtM2y7HLWdy/Wq36i9I7a0qeJ8vyXKZYewWXvfJi03RwcABVVbG9vQ3gh8LLRJNicpoB7+FoWVaoGx3u7OwAONsV+CKjvCZslmX5E1zmoVKpTGVjSe+9HKWIsKZpM59B6v3Os1iQXi6XoSgKstksBEFAtVrlnmoUCianGfA2UgTC+zRr27b/EBjl4acoyswfkvP8FO267lTu77qu//AfZUsYURRnWrWkUqn4CXHaLVav1e792/Zaid6HJqJJMDnNiJcYwvok33mdKGyemMvlkEwmkUwmsbKyglgs1rVdyazvv7KyMpVurc7qJIOSveM4UFUVmUzG/+pclxemSqXiv9+xWAyqqiKbzaJer0/9w8jOzo7fnQfA/x3n2Vqm5cHkNGPeIuNlYxgGGo0GGo0GWq0W2u32TNednb9/q9WaeotlUJeeKIool8sQBAGWZaHZbE6tFakoiv9+NxoN6LqOarWKlZWVqX448FqQne+xKIr++OosP5jQcmJymhGvCy6ZTIZyvc6HwjTGVcLQ+al61gRBmEprpbOVelHLzOv2m1XLVhRFFItFGIYB4Kw1Oa1/G954ntda9b68liW79mhSTE4z4DiO/5DY3NwM5ZqdXTaWZYVyzbBJkjT3YqRhd20JguC3mMLe6TksnS3WabVgyuUyWq2W31rtbLUKggDbtmcyW5CWF5PTDHgD096MprAUi0UAizF9d5nGIbyEW6lUIjszzft3No1uZMuyhnaZeh/AFuHfJUUXk9OUefUCRVEcuiYmCF3XIQgCHMe5cKq4pmlz+yQ7j2ns06TrOiRJguu6c9nX7CKdMzkzmUzo19c0Dfl8fuD5zokRUU3eFH1MThPy/uc7/z+hbdtQVRW5XA6yLKNer09lh15vVlapVEIqlYJlWX4sXjHdVCqFjY2NqUwQGPT7e8c0TUOpVJraPlvD7j9N9Xod2WwW1WoVyWSy630Hzv7+XsshlUqFfn9vQfX5hdW2bfu7CGSz2dDfd8uyYNv20H9Lnec49kSBzbtExaIaVPjVK/4qSVJbUZSZFX71CqF65Xy8GKZVBPR84VXvnp3vQefxad9fEIR2sVicacHTdru38KtXyiebzU6lEKppmj2FXzvfc1EU27IsT6WMkGEY/t910D2896MzvmnWVqTlFWu3/3f/dSIioohgtx4REUUOkxMREUUOkxMREUUOkxMREUUOkxMREUUOkxMREUUOkxMREUUOkxMREUUOkxMRRZK3gWIsFkOpVJr4epVKpeuaFG2sEEG0YFRVhWVZME1z6rvdzlMsFsO0Hk/TvDaFgy0nigzLsqCqate247FYDJlMZukqm0/i4OAAjuNEbh+vzpbJysqKvwFh5/edf9dYLIZkMomVlZV5h05RNM/CfkT9mKbpFw3NZrPzDidyTNOcS5HbUSmK0jYMo91ut9uNRqPv39A0zbYkSe1isTiwMHC/x5NhGG1JknoK/jYaja7XdBZllmV5pGtTtPzH/NIiUX+d25ovc7dVULIsT20LkjCkUil/52fHcfr+Db1tZHK53FjbnXjbgHitrc3NTX8zz87XiKKIarUKRVG46eGCYrceRdr169fnHQJNoF/icRzHT15BNuAUBMHfin5Q16amaZBlmYlpgTE5EdFMVatVVKtVAGeJxjTNsa/h7cTrOE5PAiyVSnAcB4ZhTB4szQ2TEy0U27ahaRpWVlbgui4cx0Eul/MH3L2HXj/eazOZDJLJJFKpFCqVStdrXNdFqVRCMpkEcPbJ3Lu292nfO57JZPxzg7qnvN2AM5kMcrkcUqkUVFXt+1rvmt4OyisrKz276Lqui0qlglQq1TNJZJL3xnXdrvfGmxE4Dfv7+13fB+mi9FpOQHfrybIs7OzswDTNqew8TTM070EvovNarZY/mK3run/cNM2u3W+9wXFFUbqOdw6Oe8rlclsUxa6dib2f8XZqLZfLXTv4mqbp73ALwN/ZVlGUdjab9SckGIbhD8533rvVarUFQeiaEOD9bpIkdcVXr9d7Yvd2nu18TedAf7FYDOW9aTQabUEQunbu7fy9vV2VR1Uul/2/m2EYfpze5AhRFLv+roNc9Hjyfi/v/fV+j1F2AeajL/r4F6LIGZSczp87PwvMe1h1PrTb7R9m/51/aHXOCvR4SaLz+q1Wy9/u3ktg551/ULbb7bau6323Kfdmm3UmCkVRBl530L3O/55B3pvOc/3eG0EQ2q1Wa6yZgeeTkxeTd//O88NclEC899eLURTFrgQ7ybVp/titRwuls6vm/CytTCYDAF3db8DZolVRFHu6j2RZRrFY7Bo075xZtr297d+zWCxCEAR/oP08796dXWeSJEEUReRyua7XJhKJnjibzabfZdep3/quQTMYg7w3wFlX2Plrer+j1/04SReZoihoNBool8uoVquhTVLwuvZc10UqlYIsy1AUJZRr0/xxKjktDe8B2jme480M6xyj6HT+Id5JkqSu773Bd8uy/DGpfvf3pk/LsoxGo9H1mkqlglqt1vOz+Xwe1WoVqqqiXq9D13UIghDalPF+7w1wNk4FnCXH8yRJgm3bcByn570Y996iKEJRFGxubmJ1dTXwtTqJoghBEOC6LgRB4My8JcOWEy21fi2FSa/ltQTOf7VaLbRarZ5WiDcpQlVVrK+vd63j8mSzWf9Tf6VSwcrKSk8raho6k9b5xOW18MJcayYIQmitm86YVVUN5ZoUHUxOtNS8h1cYSSrItSqVClZXV5HP51Eul4e2QMrlMgzD8BOGqqpTf+iKougnn4ODg65zzWYTkiSFPuttY2Oj55rjLMT1dMYb5UXJFAyTEy0178Fr2/bQB+AoD0fvWhdNsfaSV6lUgqqq0HV95G6xbDaLVquFYrEI4Cy5TbuuoLceSNM0P3bLsqa2VqizlehJpVJjJygvts4ES8uDyYmWWmdSGPSQV1W175jLed4DsN/EBU+1WvXHcXZ2dgCgbzdevxg66bruj6FMu8CrJEkwTROJRMJfY2WaJp49ezaTh361WoXjOGO30Lz3ha2m5cTkRJF2cnIy8TU6WyHnF8CqquoP2F+kc7xEVdWeBFWpVLCzs9M1iwzoTi6u6/oTIjrjaDabflLzePcaJblNwptB12g0YJomDMPwJ2RMypuwMEilUkGhUBh7wkXne+rNRKQlM++57ETnDatK3mg0uhbJdioWi20AbVEUe67ZWcka/7s2Bn0qVndev9+CVW89Ted1vO9FUexaD9R5T28xbDab9RfSiqLYzmaz/uJUSZK6ft57H87H0W9N1STvTef7cv5LkiS/wviozlclP/8et1otf5FwNpttG4bRsxbMi2uQzoXFQaqz89EXffwLUWSYpukvRu18QMqy3C4Wi+16vd71UOpcdOk9fDuT2vmHVrFY9BOGJEk9CzZ1Xe+6t1fJoN/DT1GUrqTUb1Fpq9Xy4+28X71e95Oa9xD3EpZXUcJLZp0VLbz3p/P39LaLmOS98SphyLLcliSp5/0H0BXHIOVyuStpe1/nPxCc/5BwvjqFp18C0XW96/f0/n3out73w8QgTE7Rx51wiS65XC6H7e3tnq41rz5fuVyGIAhD14RNA3fCvdw45kR0iXnVK/qN+QiCAEmSoOt6oKneRJNgciK6pBzHQbVa9RfbDuJVSyeaJSYnokvu4OBg4MJiy7JgmubA8k9E08IxJ6JLLJfL+cVqZVn2u/ccx4Ft23PdTZZjTpcbW05El5hhGH5ZpVqthlKp5Hf1maY592KqsVgMsVgMpVJp4mtVKhX/ehR9bDkREVHksOVERESRw+RERESRw+RERESRw+RERESRw+RERESRw+RERESR8/8B4Gb70CxNJ3wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 460.8x345.6 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = \"DCEFGBA\"\n",
    "MEANCOL, STDCOL = \"Em\", \"Es\"\n",
    "\n",
    "start_plotting(cmbright=False, font_size=17*1.5)\n",
    "\n",
    "ax1 = plt.subplot(211)  \n",
    "for i, LABEL in enumerate(labels):\n",
    "  a,b,c,d = LABEL2CFG[LABEL]    \n",
    "  shift = i*4    \n",
    "    \n",
    "  d11 = df[(df.a==a)&(df.b==b)&(df.c==c)&(df.d==d)]\n",
    "  exp, var = theoretical_moments(a,b,c,d)\n",
    "  Z, LABEL = 1.0/exp * 100.0, r\"Bias [\\%]\"\n",
    "  #Z, LABEL = 1.0, \"Bias\"    \n",
    "             \n",
    "  d1 = d11[(df.nl==10)&(df.no==1000)]\n",
    "  xs = list(np.arange(len(d1))+shift+1.6)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-exp) * Z, color=\"lightskyblue\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==100)&(df.no==100)]\n",
    "  xs = list(np.arange(len(d1))+shift+0.8)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-exp) * Z, color=\"mistyrose\", width=0.8, align=\"center\")        \n",
    "        \n",
    "  d1 = d11[(df.nl==1000)&(df.no==10)]\n",
    "  xs = list(np.arange(len(d1))+shift)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-exp) * Z, color=\"lightgrey\", width=0.8, align=\"center\")\n",
    "\n",
    "plt.xticks([],[])\n",
    "plt.ylabel(LABEL)\n",
    "plt.axhline(0.0, color=\"k\", zorder=-1, lw=1, ls=\"--\")\n",
    "#plt.yscale(\"log\")\n",
    "#plt.ylim(0.1, 100000)\n",
    "#plt.ylim(0.1, 10)\n",
    "plt.xlim(-1, 43)\n",
    "#plt.yscale(\"log\")\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_lines = [Line2D([0], [0], color=\"lightgrey\", lw=4),\n",
    "                Line2D([0], [0], color=\"mistyrose\", lw=4),\n",
    "                Line2D([0], [0], color=\"lightskyblue\", lw=4),]\n",
    "plt.gca().legend(custom_lines, [r\"$10^3 \\times 10^1$\", r\"$10^2 \\times 10^2$\", r\"$10^1 \\times 10^3$\"], \n",
    "                 title=r\"$S_{Z} \\times S_y$\", fontsize=15, loc=4)\n",
    "\n",
    "ax1 = plt.subplot(212)  \n",
    "for i, LABEL in enumerate(labels):\n",
    "  a,b,c,d = LABEL2CFG[LABEL]    \n",
    "  shift = i*4    \n",
    "    \n",
    "  d11 = df[(df.a==a)&(df.b==b)&(df.c==c)&(df.d==d)]\n",
    "  exp, var = theoretical_moments(a,b,c,d)\n",
    "  Z, LABEL = 1.0/exp * 100.0, r\"Variance [\\%]\"\n",
    "  #Z, LABEL = 1.0, \"Variance\"\n",
    "    \n",
    "     \n",
    "  d1 = d11[(df.nl==10)&(df.no==1000)]\n",
    "  xs = list(np.arange(len(d1))+shift+1.6)\n",
    "  plt.bar(xs, d1[STDCOL]**2 * Z, color=\"lightskyblue\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==100)&(df.no==100)]\n",
    "  xs = list(np.arange(len(d1))+shift+0.8)\n",
    "  plt.bar(xs, d1[STDCOL]**2 * Z, color=\"mistyrose\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==1000)&(df.no==10)]\n",
    "  xs = list(np.arange(len(d1))+shift)\n",
    "  plt.bar(xs, d1[STDCOL]**2 * Z, color=\"lightgrey\", width=0.8, align=\"center\")\n",
    "\n",
    "\n",
    "plt.xticks(np.arange(0, 4*len(labels), 4)+1, labels)  \n",
    "plt.ylabel(LABEL)\n",
    "plt.xlabel(r\"Increasing $\\mathbb{E}[Y]$\")\n",
    "plt.axhline(0.0, color=\"k\", zorder=-1, lw=1, ls=\"--\")\n",
    "plt.yscale(\"log\")\n",
    "#plt.ylim(0.001, 2000.0)\n",
    "plt.ylim(1, 1000.0)\n",
    "plt.xlim(-1, 43)\n",
    "\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_lines = [Line2D([0], [0], color=\"lightgrey\", lw=4),\n",
    "                Line2D([0], [0], color=\"mistyrose\", lw=4),\n",
    "                Line2D([0], [0], color=\"lightskyblue\", lw=4),]\n",
    "plt.gca().legend(custom_lines, [r\"$10^3 \\times 10^1$\", r\"$10^2 \\times 10^2$\", r\"$10^1 \\times 10^3$\"], \n",
    "                 title=r\"$S_{Z} \\times S_y$\", fontsize=15, loc=4)\n",
    "\n",
    "plt.subplots_adjust(left=0.2, right=0.95, top=0.95, bottom=0.2)    \n",
    "plt.savefig(\"RESULTS/nsamples_e.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAFJCAYAAAAR5MP4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3d32/b9r0//qfak4sUOBHtbN13CHAQUysOBhRoTdnnfjGV7HaNFO2+k5judiekvT9gNtXsupac/AE2te72NGSa+1liEqA4G9CI9oZPhs+yWaJdYAE+AcrvhUdWv60fFCXZzwcgoKEk8mXZ5UvvX693zPM8D0RERHPsrWkHQERENC4mMyIimntMZkRENPeYzIiIaO4xmRER0dxjMiMiorn3b9MO4Dz53ve+h+vXr087DCKic+nw8BD/+Mc/uj7HZBai69evo1KpTDsMIqJzaWVlpedz7GYkIqK5x2RGRERzj8mMiIjmHpMZERHNPU4AoWi9/Evv5679R3RxENG5wpYZERHNPSYzIiKae0xmREQ09zhmRkTnzsnJCV69eoU3b95MOxQa0KVLl/Duu+/iypUrI72fyYyIzpWTkxP87W9/w7Vr13D58mXEYrFph0Rn8DwPr1+/xsuXLwFgpITGbkYiOldevXqFa9eu4Z133mEimxOxWAzvvPMOrl27hlevXo10DiYzIjpX3rx5g8uXL0d+3VKphGQyiUQigUQigVgshlgshlQqFXksZ5nVWC9fvjxy1zCTGRGdO1G3yAqFAhRFga7rqNVqqNVqMAwDACBJUqSxnGWWYx3n98ZkRkQ0Jk3TIMsyZFkOjqXTaciyjNXV1SlG1mmeYh0GkxkR0RgcxwEA1Ov1juckSZp6a6fZPMU6rEhnMz579iz0cwqCwA0xiWhqRFEEANi2Ddu2WxJCNpsNnh+H67oQBGHs80QR67REmswkSQq9LzuZTOIPf/hDqOckIhpGOp1GuVzG2toaDg4OgsQTVktH0zQkEgmoqtr1edd1kUwmYRjGmdecdKzTEvk6s3g83ne30GE4jgPP80I5FxHRqAzDQCKRgOM4SCaTqNVqoZ6/WCwik8mgUCh0JDQ/kRWLxYES0qRjnZZIk5koinjx4kWo57x582ao5yMiGoVpmkgmk3AcB6lUCqZphnp+wzA6ElpzImue0DHtWKch0gkgk2jGJhKJ0M9JRDQsURTx+PFjAIBlWdA0LfRrGIYB0zShadrIiSyqWKM297MZ2c1IRLNCkiQUi0UAp+u5mlmWFSxMVhQFmqZB0zQsLCwgFovBsqyBrmGaJmzbxtLS0kiJbJBY51GkyWxnZyf0c+q6Hvo5iYgG4U91b5bP54NZgc3Pm6aJRqMB0zRRLBah6zpWV1fhui50XR84KbmuC8dxIIoibNueSKzzKNJkFo/H5+KcRERncV03aNm08xNE81T3jY2Nlun1juMgk8kgnU73nKXY7Zp+12K1WoVpmgO1qoaNdR7NVDfjs2fPcP/+fdy/fx+ff/45Tk5Oph0SEVFXjuOgXC53fa5SqXS0tJoTmeu6SKVSkCQpKCV1lm5jZKZpDpTQhonVtm2kUqmWOo1+t+gsm5lktrOzE5RZWVtbw9HREVRVxa1btyay2JqIaByVSgWO40BRlJbj/r97tYQAYG1tDfV6PZiEcZZ+kz0GSWjDxOqPpVmWFXQ9znoiA2YomZVKJXzxxRf48MMPsby8jFwuh+3tbXz22Wf4zW9+gwcPHoR6Pdd14bpuqOckoovDX59VKpUQi8WQTCaRTCZRr9dxcHDQs9tOURTYto3Hjx8PXNVD07S+kz1M08T+/n7PSSTDxiqKIiRJamnNTbui/pm8GZFIJLyTk5Oez2ua5j1+/Hikc6uq6gFoeQiC0PO1six7+Xzek2XZM01z4Oskk8mR4rtQ/s+fez+IQvC///u/0w6hp2Kx6AHwDMOYdihn0nXdk2XZ8zwv0nj7/f763WMjWzS9urqKBw8e4IMPPuj6/NraGtbW1vD48WP8+7//e8fzW1tbuHXrFm7cuDH0tUulUsc3mvbmNnBaGgsAqtUqgNN+5kQiAcMwkE6nh74uEZHPtm0oigJVVTvuJ+11EmeBLMvQNC2YOTnrIutmvHHjBiRJQjabxeHhYcfz6+vrePHiBZaWlvDb3/626zm6VXo+S6FQwMbGRtCv7D/a/5gKhQJs225ZPiCKIlRVRSaTYZckEY3McRysra1BluWO5UTlcnkmp8X7ydWyrJlLtN3EPC+6Vcd3796FIAiwLAurq6vQdR1XrlwJnrcsCzdv3kQsFoMgCMH+OqIoYnd3F7Zt4+uvvx7qmgsLC9jY2EA6ne777WJhYQEA0Gg0Wo5bloVUKgVd18+cPruysoJKpTJUfBfOy7/0fu7af0QXB51bf/zjH/HjH/942mG08GshptNpLC4uAjj9cu44DmzbntniD9ModdXv99f3Hht6h2cfrut6d+7c8TzvtO84kUh46+vrLa+pVqueJEleLBbzYrGY99ZbbwWPYcfM/P5pNI2TFYvFjtfVajUPQNA/3KzRaHgAPFEUz7wex8wGwDEzmrBZHjObJ7VazatWq5Ffd9Qxs8gXTXv/+gaSz+fx4sULLCws4L333sPDhw8BnDZtq9UqqtUqdF1HLpfDvXv38OLFi6HHy2RZhmEYUFUVoijCdV0oitIxK8dfRd9tZpF/bBa7AYiIwuY4DlzXnclxvH4in5p/9erVln+rqor9/X3s7+/jvffew+9//3sAwPLyMu7du4ft7W1sbW1haWlp6GuJooh0Og1d11Gr1WCaZtDN2bxuwk9UfvOfiOgism0byWQSe3t7czfpLdJkdnBwANM0kc1m8cknn+DBgwc4OTmBIAjY3t7GF198gc8++wyrq6t48uRJ6NeXZTmYqdi8fuKstR7+890mgZRKJaysrGBlZQV///vfQ4yWiChakiSh0Wggn89PO5ShRZrM1tfXoWkadnd3g9ZWOp3G/fv3AZy2pB49eoTNzU3kcjncunULz58/DzUGURSRz+dbug39iSG9Zkv6Saxb0svn86hUKqhUKvj+978faqxERDSYSJOZbdvI5XIATsfP1tbW8OjRI8TjcXzyySfB62RZxosXL3D79m385Cc/QTabxZ///OfQ4mjfA83f+bpby8s/Nk99x0REF02kyczzPHzzzTcdx3O5HOLxOD7//POW4/l8HvV6HdevX8fS0hJ++ctfhhLH0dFRS3ISBAGiKHad8ukfG3XPICIimrxIk5ksy8hkMl2r4W9sbPQszKnrOur1Ov7xj3+EEke5XO5YuKgoSrBPUDN/jUW3iiEUrq2nb7o+iIjOEmky0zQNjx49wsrKSjAVv1m/BceCIGBvb2/ga9m2jVgshkwm01J8U1EUKIrS0dJSVRWSJLUkLcdxUCgUoOv6XJRzISK6qCKrzQgAS0tL2Nvbw507d5DP56GqKlKpFJaWllAul4OxqzBIkoR0Oh3szSPLMkRRhKZpPRNTtVqFpmnIZDJYXFyE4zisy0hENAciTWYAkE6nUalUkMlkcHBwELS2BEEIvWzKoJveNWvvfiQiotkXeTIDTltNtVoNBwcHQfWN27dvTyMUIiKaoH5Lm8I01c05l5aWcPv2bSYyIqIxaZqGVCqFVCo1M7t8OI6DZDIZSQH2SJNZ+9T7WT0nEdE8cRwHuq7DNE2Iothzx+mo+TtWRyHSZLa7uxv6OaPenoCIKCqDtrD8SW3+6yc5aW1WWn3tptrNGIZRNuwkIppllmUhkUigVCq1HPd3/igUCtA0raUF5rouNjc3sbe311J7dpoxRSnSCSCWZeHJkycjVcDvdT5uzUJE/Xz11VfTDgEA8P777w/0Otd1ey5TyuVySKVSQSHgRCIBwzAgSRIEQYCu68hms9A0LdTW2agxRSnSZNZoNEIvC8WaiUR0nviz/tpn/9m2jXK5jJ2dneCYLMvY3NxsWYYkCELo98VxY4pCpMlseXkZsVgs1HOyMgcRXQSVSgWCILQklGQyCU3TUC6XUSwWoSgK6vV63/WyhUIh2OuxmaIofYtKDBsTcDoxxbZtmKaJlZWViU7PjzSZ+XuJERHRcKrVascGwouLi3BdF+l0euBuRVVVg2Tjv8cv8zds46BfTMBpY6NWqw11zlFNZdE0EVFUBh2rmnVhtmp0XQ8SmmmaUBRlpK7JSS+EHgaTGRHRHEgkEh2zt+v1+shDLbquI5FIjJzIJhHTOOZ+aj4R0UUgy3LHGq9qtTpyIlIUBYZh4OjoaOSp/GHHNA4mMyKiGdSeJERRhCzLLYnHsqyRiqP7Y2SSJEHXdezv7w+U0CYZ07jYzUhENGNKpRIcx8Hu7i4kSQqWNBmGAU3T4DgOarXaSHstFgqFjq5FfwytX/mpScYUhpjneV7kVz2nVlZWIimoOQ96LVR9f+FKz/dsvfph1+Pry5dCiYkuhj/+8Y/48Y9/PO0waET9fn/97rHsZiQiornHZEZERHNvppPZxsYGVldX8ctf/hJ//vOfpx0OERHNqKkns9XVVVy9ehW3bt3Cl19+GRz/5JNPUCgUUK1Wsb29DVEU8fz58ylGSkREs2rqyaxWq6FUKuGLL77AjRs3AAAPHjxAsVgEcDqD5ttvv8UXX3wBVVWnGSoREc2oqSczWZZx+/bt4N/Hx8fI5/OIxWJQVRW/+MUvgteFtXUMERGdL1NPZlevXm35t18vTBAEbG5utjzXaDQii4uIiObH1JNZrVYLJnc8ePAApVIJsVisZX8c4LTFNq0dTImIaLZNvQKIqqpYWlrCwsJCUColn8/jo48+Cl7z7Nkz/OIXv+gopUJE49t6+qbrcS5Wp3ky9WQmyzIqlQqKxSIajQay2WzLGNrdu3eDKszciJOIqDtFUWBZFgRBgGEYM3G/jDKmqSczAJAkKZi92G57exsAcHh4iOvXr0cYFRHRfPCL+wqCAEVRoOt6z3vqeY1p6mNmg9I0Db///e+nHQYRUWQGHVqRZTnYKDOVSiGRSFyomIAZaZkBpy0vx3G6Pue6LqrVKp4+fYqf/exnEUdGRBQty7KCbVqa19e6rgtN05BIJHB0dIRUKhVUr/ft7+9jY2PjQsTUbOrJ7Pj4GLIsw7btvq/zPA8LCwsRRUVE58bLv0w7glPX/mOgl7mui5WVla7P5XI5pFIp5PN5AKc7PRuGEWzbYlkWstls0CIKyyzG1G7q3YyapqFarWJ5eRlra2uIx+NYW1trecTjcWQyGU7NJ6JzTxCE4NHMtm2Uy2XcuXMnOCbLcrAe17btlv3Iwpz9PYsxtZt6y6xSqaDRaCAejwMAdnZ2sLq6ig8//DB4jeu6UBQFy8vL0wqTiGiqKpVKR0JJJpPQNA2WZSGTyWBxcRHAafKpVqtdz1MoFCCKItLpdMtxRVGCDTqjjikMU2+ZiaIYJDLgtMnaPuNFEAQ0Gg08fPgw6vCIiGZCtVoNEoNvcXERrutClmU0Gg3UajXUarW+SUNVVezv76NcLgfH/LGwYafOhxVTGKbeMjs+Pu44Fo/H8fz5c3zwwQfBMVEUsb29jY8//jjK8Iho3g04VjXrwhxz0nU9KB1omiYURQm6AqcV07imnsyWl5eRzWaxuLiITCaDGzduYH19HaIo4tNPP0UymYRlWUGZKyKiiyiRSKBer7cc8wtKjELXdSQSiZET2SRiGsfUk9nGxgauX78e1F78+uuvIQgCtra2kMvlWhJY+3RPIqKLQpbljgkU1Wp15ESkKAoMw8Du7i7K5XLHGNo0YhrH1JNZPB7H4eEhSqVSyweQz+fhui42NzdxfHwMSZJgGMYUIyWisPWqCwmwNmR7khBFEbIstyQey7JgmubQ5/bHyCRJgiRJQZfjWQltkjGNa+rJDDhNaPfu3es4rqoqN+QkogunVCrBcRzs7u5CkqSgV8owDGiaBsdxUKvVoOv60F16hUKho2vRH0NrnkYfZUxhiHme50V+1RHNen3GlZUVVCqVaYcxE7766quux99fuNLzPVuvftj1+EX/hj5p06yaP4mW2R//+Ef8+Mc/HjUkmrJ+v79+99ipT80fxrQLZxIR0WyKrJvx8PAQlmXhzp07uHLlu2/ng9brchwH5XK5Y/dpIiKiyJKZLMs4ODhArVZrSUimaeLp06dnvt/zPE7NJyKiriJLZmtrazAMA6lUquV4Pp/H3bt3Icty30HDWq2GL7/8ctJhEhHRHIosmRWLxa5jXtlsFpZlYW9v78xz3L17dxKhERHRnJv61Px4PD7wuFkmk5lwNETTEcmMwp5boXSfRUo0T2ZiNuOg1fDX1tYmHAkREc2jmUhmRERE45i5ZHZycoJsNourV6/i6tWr+OlPf4onT55MOywiopnnuu5EN8CcZZEns5OTE9y6dQtvv/023n77bfzXf/0XDg8PAQAHBwdYWFhAuVwOfimmaUKWZfz617+OOlQiornhOA6SyeRMVSHSNA2pVAqpVGriSTbyZCZJEizLgud58DwPlUoFN2/eBHA6wcM/3v7QdR3Pnj2LOlwiornQr67iNDiOA13XYZomRFGEZVkTvV6kyex3v/sdHMdBPB5HPp+HqqpYXl5GrVbD/fv3Yds2kskkqtUqvv32W3z77beoVqvI5XLwPA+lUinKcImIpmoWuwwHjclfN+y/fpQtZoYRaTLb3d1FIpHAwcEBtre3sbW1hWq1io8++giapmFhYQGPHz9umd24vLyMYrGItbU17O/vRxkuEdFUWJaFRCLR8QXedV0oioJCoQBN0ybe2hk3Jn8br729PZTL5YnGF+k6M9u2sb6+jng83nJc13X87ne/Qz6fb6nb2EzTNGSz2SjCJKJzpF9l/igNumbQdV2srKx0fS6XyyGVSiGfzwM43enZMIyJdy+OGpMgCNB1HdlsFpqmTbR1FmnLrNFo4M6dOx3H/eZoe6mrZqurqzPZ5CYiCpMgCMGjmW3bKJfLLfdQWZYjKb4+bkyCIEw84UaazFzX7dnyEgQBi4uLPd/b3pojIrpIKpVKR0JJJpNBt57jOLBtG6Zp9v3iXygUunb5KYoCx3FCi6lcLiOVSqFcLsOyLOi6PtS5hzX1clbN2rN+s+Pj4wgjISKaLdVqteML/+LiYpC4RFFErVY78zyqqkLTNADfTcpQFAWKogy9Q3S/mNLp9MQnfTSLNJl5nodf//rXXT8w13VRLBaRSCS6PmeaZhQhEtE5c152Ku/3ZX9Yuq4HCc00TSiKMlI3YJgxjSvyllmvpmYsFkOhUOj5Pu5nRkQXWSKRQL1ebzlWr9eHbk35dF1HIpEYOZFNIqZxRL5outei6LMeREQXmSzLHWNh1Wp15ESkKAoMw8DR0dHI0+bDjmkckbfMyuXywFXym9VqNU7NJ6ILoz1JiKIIWZZRLpeDsSjLskYagvHHyCRJgiRJHWNo04hpXJEms3Q6jY8++mik9y4tLXE/MyK6EEqlEhzHwe7uLiRJgizLAADDMKBpGhzHQa1Wg67rQ3fpFQqFjq5FfwytX0msScYUhkiT2bgtKyYzIroI8vl8sAi5mSAIKBaLY51bVdWux8+aOj/JmMIQ6ZjZ7du3x3o/N+ckIqJuZm4/MyIiomExmRER0dxjMiOic4fLeebTOL83JjMiOlcuXbqE169fTzsMGsHr169x6dJoFVtmqjYj0bn28i99nvxhZGGcd++++y5evnyJa9eu4fLly6wcNAc8z8Pr16/x8uVL/OAHPxjpHExmROi959V5qet3kfg7c/z1r3/FmzezsZcZne3SpUv4wQ9+0HNnlbMwmbXRNA22bUMURTiOA03TgsWBRDQfrly5MvJNkeYTk1mTZDIJ4LS2GHC6P5C/a2qUWxkQEdFwmMz+pVAowLbtIJEBp3XHVFVFJpNBo9GYqe0OiIjoO5zN+C+bm5tdt/ZOpVIATuuSERHRbGLLDKfdia7rdh0bW1lZAQAUi8WeNc1G0mNm29ar7rPaOBGBiKg3JjMAtm0D6L5rqn/McZyhz/vVV1/1fO79hdkbnO41ow84R8m05/T48zE1ftp/c72u//7770/82nSxMZnhu0S1uLg45Ugics5v6NPW84Y+g19gIsO/OZowJjN0b5G1P++6LlzX7XhtqVQKxtP+9Kc/Bd2SYRttH9jzcf2Leu1pX/8i/+w0mw4PD3s+x2QGBBvJ1ev1rs/7u6t2S3q99vghIqLocDYjvpvk0b4lePOxXruvEhHR9DGZ4bTFJYoiKpVKx3P+MVYBISKaXUxm/6IoClzX7Zi1aJpm8DwREc2mmMeNfwLJZBKLi4tBAvPLWem6Hu4aMyIiChWTWRtN0+A4DhYXF+E4DhRFYV1GIqIZx2RGRERzj2NmREQ095jMiIho7jGZERHR3GMyIyKiucdkRkREc4/JjIiI5h6TGRERzT0mMyIimntMZkRENPeYzIiIaO713Zzz2bNnoV9QEARcv3499PMSEdHF1bc241tvvYVYLBbqBZPJJP7whz+Eek4iIrrY+rbMACAejwc7MY/LcRywrjEREYWtbzITRREvXrwI9YI3b94M9XxERER9J4BIkhT6BROJROjnJCKii+3Mbsawneduxu9973uc3DKG//vP7n8b/9874Y7bzqJIfvY3/6/7td9cCu3ar1+/7nr88r+9PdS1R70+nW+Hh4f4xz/+0fW5vhNAjo+PEY/HQw1mEuecFSsrK6hUKtMOY25tPX3T9fj6cu8b3nkRyc/+8i/dr/3qh6Fd+6uvvup6/P2FK0Nde9Tr0/nW7x7bt5txEknnvCYyIiKanrG6GZ89ewbLsgCcThaRZRlXrnT/BkZERDQpI1cA2dnZgaZpkGUZa2trODo6gqqquHXr1kQWWxMREfUycsusVCphf38/+Pfy8jJyuRwcx8H6+jpu3ryJX/ziF6EESfOn59jJ++9HHAkRXQQjt8wajQa++eabjuOiKGJvbw8vXrzAl19+OVZwREREg+iZzFZXV/H8+fOeb1xbW8Pa2lrXhAYAW1tb0HV9/AiJzouXf+n9IKKx9OxmvHHjBiRJQjqdhq7rHeun1tfXkUwmsbS0hI2NDfzqV7/qOEe9Xg89YJpzfW/cvadp02B6de8CvafHE50HPZOZrus4Pj6GIAhIp9NYXV2FruvBbMWlpSXs7e3h5s2bUFUVv/nNbyDLMlZXVyGKInZ3d+G6bmQ/CNE4LvIaN6LzoO8EEF3Xkc/nUalUUCqVIEkSMpkMNjc3AQCyLKNSqSCXy+Hp06cwDAPlcjl4v2mak42eiKiLk5MTvHr1Cm/edP+SQrPn0qVLePfdd0de3tU3mcXj8aD8VD6fRz6fR6FQwHvvvYf19XV8/PHHkCQJ1WoVT58+hWVZqNVqEAQBiqJgaWlppKCIiEZ1cnKCv/3tb7h27RouX74c+jZWFD7P8/D69Wu8fPkSAEZKaGfOZrx69WrLv1VVxf7+Pvb39/Hee+/h97//PYDTqfn37t3D9vY2tra2mMiIaCpevXqFa9eu4Z133mEimxOxWAzvvPMOrl27hlevXo10jr7J7ODgAKZpIpvN4pNPPsGDBw9wcnICQRCwvb2NL774Ap999hlWV1fx5MmTkQIgIgrTmzdvcPny5civWyqVkEwmkUgkkEgkEIvFEIvFkEqlIo/lLLMa6+XLl0fuGu6bzNbX16FpGnZ3d4PWVjqdxv379wGcril79OgRNjc3kcvlcOvWrb7T+YmIohB1i6xQKEBRFOi6jlqthlqtBsMwAExmK61xzHKs4/ze+iYz27aRy+UAnI6fra2t4dGjR4jH4/jkk0+C18myjBcvXuD27dv4yU9+gmw2iz//+c8jB0VENE/80n6yLAfH0ul0MMN7lsxTrMPom8w8z+u6KDqXyyEej+Pzzz9vOZ7P51Gv13H9+nUsLS3hl7/8ZbjREhHNGMdxAHRfVytJ0tRbO83mKdZh9U1msiwjk8ng5OSk47mNjQ0Ui8Wu79N1HfV6vecmakTn2VdffdX1QeeTKIoATnuybNtueS6bzQbPjyOsNbtRxDotfZOZpml49OgRVlZW8PDhw47n+21EKQgC9vb2xo+QiGjGpdNpAKdl/poTT1gtHU3TUCgUej7vui4SiURHgupm0rFOS99k5lf5ePHiBfL5PK5evYqf//zn2NjYwMrKClZWVqKKk4hoZhmGAVEU4boukslk6OcvFovY39/vmtD8axaLxYES0qRjnZYz15ml02lUKhVcv34djUYDe3t70HUdR0dHPbsZiYguGtM0IQgCHMeZyBR3wzA6ElpzImue0DHtWKdhoC1gJElqmcJpGEYw0YOIiE7Hox4/fgwAsCwLmqaFfg3DMGCaJjRNGzmRRRVr1Ibaz2xpaQm3b9/G7du3JxUP0WRxCxaaIEmSgh6r9i5By7KChcmKokDTNGiahoWFBcRiMViWNdA1TNOEbdtYWloaKZENEus86pvM2qfeh2ES5yQimgZ/qnuzfD4fzApsft40TTQaDZimiWKxCF3Xsbq6Ctd1oev6wEnJdV04jgNRFAea8DFKrPOobzLb3d0N/YKspE9E54Hruj3nDfgJonmq+8bGBgRBCP7tOA4ymQzS6TRUVR34mn7XYrVahWmaA7Wqho11Hg3VzRgGbthJROeB4zgtW141q1QqHS2t5kTmui5SqRQkSQpKSZ2l2xiZaZoDJbRhYrVtG6lUqmViiN8tOsv6bgFjWRaePHkSWgV8y7LmvilLRAScJgHHcaAoSkurR1EUAOg723ttbQ31eh3VanWga/Wb7GGaZpB4erXwhonVH0tLJBJBd6amaTM/e71vMms0GiMPLvYy7wvziIgAoFarATitQO9vXgycdtcdHBy0tMSaKYoC27ZRrVZ7vqadn0x63Y9N00Qmk4FlWV1fM2ysoihCkiSUy+UgQc76FP6+yWx5eTn06tPz3i9LRASclu3TdX2o9/jJxDCMob7YD9Iq6tddOUqs2WwWpmlCVVXYth1UDplVfZPZoE1gIiLqz7ZtKIoCVVU7EoNt2zPXayXLMjRNC7oaZ13kE0CIiC4ax3GwtrYGWZY7Wkjlcnkm5xL4ydWyrJlLtN30bZkREdH4UqkUXNeFIAjBpIt6vQ7HcWDbNmJpVLYAACAASURBVDzPm3KE3cmyjHw+P+0wBsJkRkQ0Yf4EjHniOM7Q42zTxG5GIiIKOI4D13VnchyvHyYzIiICcDoRJZlMYm9vb+ZnL7ZjNyMREQE4nfTRaDSmHcZImMzoXPrqq6+6Hn9/4UrEkRBRFNjNSEREc2/sZHZ4eBj898HBAbd4ISKiyI2UzA4PD5HNZvH2229jZWUlOL60tATP87C6uopvvvkmtCCJiIj6GXrM7Pj4GMlkEq7rwvO8jsV+t2/fxtHRESRJwtdffx1aoERERL0MncxyuVxQTT+VSqFUKnW8Jp/P4+7du/jtb3+LX/3qV6EEelFsPX3T9fj68qWIIyEimh9DdzPatg3DMPDo0SPcu3ev5+tEUcT29vZYwREREQ1i6GRWr9dx+/btgV43i8Uz+3FdF5lMBgsLC0gkEl1bnURENDjXdeG67sSvM3QyW1lZwcnJSd/X7OzswHXdudg2oFkmk4GiKHj8+DFkWYaiKD23GicimiWapiGVSgVFjWeB4zhIJpOoVCoTv9bQyUzTNORyuZ7PP3jwAIqiIBaLhb5LdS+FQgGWZY11jnK5DF3XIctysG24KIowTTOkKImIJsMvCmyaJkRRHPt+GBZ/x+ooDJ3M1tbWIAgCrl69imw2i3q9jo2NDWSzWVy9ejXY3kAURXz22WehB9ysVCphYWEBmqad+U3E/9aiKApSqVTHLzudTnd86IIgIJlMhh43EdEgBm1h+b1g/usnWVdxVlp97UZaZ1YsFrG5uQnTNNFoNKDrOgzDQKPRgOd5yOfzE29W2raNO3fu4M6dO2e+NplMwrIsmKaJYrGIYrGIVCp1Zhei4zgDnZ+IKEyWZXUdt3ddF4qioFAoQNO0li/lrutic3MTe3t7ExkeGSWmKI1cmzGfzyOfz+Pg4CDYMkAURYiiiHg8HmaMXfmtqEQi0fd1hUIBtm2jWq0Gx0RRhKqqyGQyaDQaEASh432apsEwjK7PERFNiuu6LcUomuVyOaRSqWDDzEQiAcMwIEkSBEGAruvIZrPQNC3U1tmoMUVp7ELDS0tLWFpaAnC6oLrRaESSzAa1ubkJQRA6PthUKoVCoYBSqQRVVVue87+BRDXmR0ST06vodNTef//9gV7nf4Fu/yJt2zbK5TJ2dnaCY7IsY3NzE4ZhtLw/7EQybkxRGKmb8e7du3j77bfx8OHDluPxeBxbW1vY2NgIJbhx+S3Gbt8o/GPFYrHluG3bcF13brYKJ6KLoVKpQBCEloTiD6GUy+Vg6MSyrL47RBcKha7dkIqiDL2cql9MwOk92LZtmKY58bG2oVtmOzs7QZ/pixcvOp7f3t7G3bt38dOf/hT/8z//M36EY7BtG0Dnt4nmY82/PNu2YVkW0ul0cNyyLKysrMzVjquj6lV9BGAFEqJpq1arWFxcbDm2uLgI13WRTqcH7lZUVRWapgH4bqKIoihQFGXo5VT9YgJOh3RqtdpQ5xzV0MmsWCxClmVomoa1tbWur9ne3sZbb72Fhw8f4uOPPx47yFH5Can9w+7G32EVQPCLBqL9ZRAR9RLm+L2u68F9zjRNKIoy0hf2WZpTMHQyc113oJmKgiCgWCxONZmd9UELghCsTpckqaNo8iBKpVLQUv373/8+UpxENDmDjlXNukQigXq93nKsXq+PXJxC13UkEomRE9kkYhrH0GNmgwR5fHwM13VbZhBOgx9r+4ft85vC43y78JchVCoVfP/73x/5PERE/ciy3DHuVK1WR05EiqLAMAwcHR2NPJU/7JjGMXQyEwQBz58/7/sav/k67XEmf5JHt4FH/9i0YyQi6qb9viWKImRZbkk8Z0326MUfI5MkCbquY39/f6CENsmYxjV0Mtva2kI6ncaTJ086njs5OcEnn3yCUqkUaTmrXgRBgCiKXbtF/WPTjnEqXv6l+4OIZkKpVILjONjd3W1ZhGwYBkzTRKFQgKIo0HV96C49/73NX+T9hOZPmos6pjAMPWYmiiLu3buHtbU1JBKJ4APxp2D6JEnC5uZmeJGOSFEUaJoGx3FaPmC/5qJffus86rW+5v2FKxFHQkTD8ItStPPnIoyjfV2t76zW1CRjCsNI68zy+TwePXqEb7/9FoZhwDAMVKvVYOdpTdOwv78fdqxdHR0dAeg9LqaqKiRJaklajuOgUChM7RsEERGFa+QKILIso1ardZSzWl5eDjO+nhzHQblcDmYSFovFYL1Fe4KqVqvQNA2ZTAaLi4twHAeGYUy0GCcREUUn1HJW7U5OTnDlymS6tPz6ir2azO2mMSDZV88xqh9GGgYR0XkwUjfjoPrte0ZERBSWkVpmh4eHKBaLQfdiN47jDF3ni4iIaBRDJ7Pf/e53wR5fZ1XMiMVio0VFREQ0hKGTmaZp8DwPsiwjlUr1rJ7x4sUL3L9/f+wAiYiIzjJ0MqvX60ilUvjiiy/OfC1bZkREFIWhk5ksy7h69epAr93a2ho6ICIiomGNVM5q0AXRDx48GDogIiIanqIoSCQSSCaTMzP5LsqYhk5m3/ve97CxsYFsNovDw8Oejy+//HL21nYREZ1DfnHfWq2GlZWVmbj3Rh3T0N2M169fx/HxMQCMvG0AERGdzXXdgbaoai6YnkqlJtoKmsWYgBFaZrlcLqjBeNaDiIiGZ1kWEolEUK7P57ouFEVBoVCApmkt1et9+/v7XQsCn8eYmg3dMvv5z3+Op0+fwjAMxOPxnq9zHAd3794dKzgioovGdd1gL8Z2uVwOqVQqSAyJRAKGYQS7l1iWhWw2O9aGw/MSU7uhk9ny8jI0TeubyIDT2omz0G87Tb22YAG4DQtRZGZlr75r/zHQy/ybfvvN37ZtlMtl7OzsBMdkWcbm5iYMw4Bt2xBFMSi0Pmh34LzG1G6k2oxra2tnvubw8DCyCvpEROddpVKBIAgtySCZTMKyLFiWhbW1NaRSqWD2YK+kUSgUus53UBRl6HGtsGIKw9hV83vRNA0///nP8bOf/WxSlyAiujCq1SoWFxdbji0uLsJ1XciyjEajMdB5VFWFpmkAEGyDpSgKFEUZen/HsGIKw8jJ7PDwsGcWd10X1WoVT58+ZTIjIgpBmK0aXdeDhGaaJhRFCca4phXTuIZOZsfHx5BlGbZt932d53lYWFgYOTAiolAMOFY16xKJBOr1esuxer0+dGvKp+s6EonEyIlsEjGNY6RCw9VqFZIkYXFxEZVKpWOWS6VSQSqVwvr6emiBEhFdZLIsd2y55d+LR6EoCgzDwO7uLsrlctDlOM2YxjF0MqtUKmg0GsFsxp2dHayuruLDDz8MXuOvO+AEECKi0bQnCVEUIctyS+KxLAumaQ59bn+MTJIkSJLUMYY2jZjGNfRsRlEUW6bl53I5FIvFltcIgoBGo4GHDx+OHyER0QVTKpXgOA52d3dbFiEbhgHTNFEoFKAoCnRdH7pLz39vc+tJ13Xs7+/3HT6aZExhGGnMrF08Hsfz58/xwQcfBMdEUcT29jY+/vjj8SIkIrpg8vl814oZgiB0NB6Gpapq1+NnrQueZExhGGnRdDabxeLiIjKZDG7cuIH19XWIoohPP/00WGNQKpW4nxkREUVi6GS2sbERFBu2LAtff/01BEHA1tYWcrlcSwJrLjRJREQ0KUMns3g8jsPDQ5RKpZY+13w+D9d1sbm5iePjY0iSBMMwQg2WiIiom5EWTcfjcdy7d6/juKqqPftjiYiIJmWk2oyDOjw8nOTpiYiIAEw4mc3CDBciIjr/enYzHh4ewrIs3LlzB1eufLddycbGxkAndhwH5XIZm5ub40dJRETUR89kJssyDg4OUKvVWhKSaZp4+vTpmSf2PI9T84mIKBI9k9na2hoMw0AqlWo5ns/ncffuXciy3HeVd61Ww5dffhlepERE1JdfbmqWqtlHpWcyKxaLXce8stksLMvC3t7emSe/e/fueNEREdFAHMdBKpVCsVicmTW+mqYFJbIMw5hokh16Akij0cB//ud/4uTk5MzXZjKZkYIiIqLhiKI4lWr1vTiOA13XYZomRFFsqec4CUOvM/PH0uLxOP77v/+772vX1tZGDoyI6KJzXXfmugwHjckfhvK7PkfZYmYYQ7fMjo6OAKBrwcl2g7TeiIiolWVZSCQSKJVKLcf97bUKhQI0TZt4a2fcmPyqUHt7eyiXyxONb+iWmWEYuHv3bst0/V5KpdKZrTcioknaevpm2iEAANaXLw30Otd1OzY89uVyOaRSqaAxkUgkYBjGxLsXR41JEATouo5sNgtN0ybaOhu6ZSbLMorFIrLZLJ4/f97zdQcHB1w0TUQ0JEEQgkcz27ZRLpdx586d4Jgsy5Gs5R03JkEQJp5wh05mP/rRj3Dz5k2Uy2VIkoS333676+NHP/oRHMeZRMxERBdOpVLpSCj+llvA6YQL27ZhmmbHjtDNCoVC1y4/RVGGvmf3i6lcLiOVSqFcLsOyrDP3SxvX0N2M+Xwe6+vrA72Wi6aJiMJRrVaxuLjYcmxxcTFIXKIoolarnXkeVVWhaRqA7yZlKIoCRVGG3iG6X0zpdHrikz6aDZ3MFEXB1tYWvvzySywtLfV83dHRUc8+ViKiqAw6VjXrwpzVqOt6kNBM04SiKCN1A87STMuR9jPTdR0ffvhh39c1Go2B6zgSEVF/iUQC9Xq95Vi9Xh+6NeXTdR2JRGLkRDaJmMYxUtX8XC535ms0TUMikRjl9ERE1EaW5Y6xsGq1OnIiUhQFhmHg6Oho5GnzYcc0jpE25wROq+r3Gix0XRfVahVPnz7FRx99NHJwREQXVXuSEEURsiyjXC4HY1GWZcE0zaHP7Y+RSZIESZI6xtCmEdO4hk5mx8fHkGU5qLfVi+d5WFhYGDkwIqKLqlQqwXEc7O7uQpKkoNaiYRjQNA2O46BWq0HX9aG79AqFQkfXoj+G1q8k1iRjCsPQyUzTtKAZubi4iEql0jHRo1KpIJVKDTzrkYiIvpPP57tWWRIEYez1u6qqdj1+1tT5ScYUhqGTWaVSQaPRQDweBwDs7OxgdXW1ZUKIX95keXk5vEiJiIh6GHoCiCiKQSIDTieDtGdlQRDQaDTw8OHD8SMkIiI6w9DJ7Pj4uONYPB7vKG0liiK2t7dHj4yIiGhAQ3czLi8vI5vNYnFxEZlMBjdu3MD6+jpEUcSnn34alDIplUqsAEJERJHom8w+//zzjqn1GxsbuH79Oo6Pj2FZFr7++msIgoCtrS3kcrmWBDYru50S0cXieR6/TM8hz/NGfm/fbkZN0/DNN9+0HIvH4zg8PISu6y3diPl8HltbW7hy5Qo8z4MkSTAMY+TApsF1XWQyGSwsLGBhYaFj3x4imn2XLl3C69evpx0GjeD169e4dGm08mN9k1mtVoMgCPjtb3/bstFmPB7HvXv3OnaSVlUVjUYD3377LSqVykB7ns0STdOg63pQimuUKtJENF3vvvsuXr58iX/+859jfdOn6Hieh3/+8594+fIl3n333ZHOceaY2fLyMu7duwdVVZFOp/HrX/8aH3zwwUgXG1ShUGhZlBcF13Wh63pQONOvLN1eEZqIZpv/Jfqvf/0r3ryZjY056WyXLl3CD37wg5EbQX2TmSzLePToEYDTTdg2NzexvLyMRCKBu3fvIpfLhdr6KpVK0DQNruue2UWpaRps24YoinAcB5qmjZX82qs/l8tlFIvFmaoKTUSDuXLlytz1DNF4+nYzKooS/Lc/BtZoNJDL5fDZZ59hYWEB2WwWT548GTsQ27Zx586dlh1Le/FnTJqmiWKxiGKxGGwCF0YcqVQKmUwG1Wp17PMREdHk9U1mt2/f7jgWj8ehqipevHiB/f19fPvtt1hbW8N7772H+/fvt4ytDUOSJAiCcGal/UKhANu2sbOzExwTRRGqqiKTyfTdYXXQOAzDgK7rQUuRiIhm28hV84Hvbvyu62JnZwfb29vQNA2ZTAaKouAnP/lJWHEGNjc3IQhCRzHMVCqFQqGAUqkU1B4rl8tnVm9OJBIdtcoEQYCqqqjVasGW5ERENLvGSmY+13VRq9WCmX+GYcAwDEiShP39/TAuAQBwHAeu63YdG/OLHReLxSA5jbttdzKZHPm9REQUnZE25/Q9e/YM2WwWiUSipdtveXkZe3t7oSYyAMG2M90mZfjHxplK394KM02T3YxERHOgb8vs5OSk64ygzz//HJubm0Fy8ddyyLIMTdM61p+FxU9Uk5gu71f6r9fr2NjYgCAI2NnZOXM2Y6lUChZX/+lPf+rYDmdQf//73/H973+/5/PjT22Z3evP8s/Oz/38Xp/mz+HhYe8nvT7W19eD/z4+PvY+/fRTb3Fx0Xvrrbe8t956y4vFYl4sFvMURfEcx+l3qoHpuu4B8AzD6HiuWCx6ALx8Pt/1vYIgeAC8RqMRSixRSiaTF/b6/Nkv3rVn4fp0vvRtmZXLZayurmJ3dzeY9u79qxUmCALy+Tw2NjZatoSZJH/30nq93vV5fyYj14YREV0sfZNZrVZDJpMB8F0SE0URmqYhl8tNPro2fhdet+n3/rFeW34TEdH5NdBsRs/zJj4eNghBECCKIiqVSsdz/rF5rdTfbTvyi3J9/uwX79qzcH06Z/r1QcZiMS+RSHiWZUXS5+l5/cfMmp+v1Wotx1VV7XqciIjOv75T80VRxIsXLyJtjR0dHQHoPS6mqiokSWopteU4DgqFAnRdD8bViIjo4ujbzajrelRxwHEclMvlYJp7sViE67pIp9MdCaparQaVRhYXF+E4DgzDGGuBNBFFJ5lMsvYphWvaTcOLxDRNL51OewCChyiKniAIniAIniRJnqqqXrVanXgs1WrVy+fzniiKHoBgWUPzY5Jx1Go1L5/Pe5Ikdb32JP40DcPwZFkOzu9/5pIkeaIoBp8F+iz/mKfr+tfO5/Mtn2vzdSVJ8nRdj3Q5S7Va9QB4xWIxsmvS+cdkNgX+zUvX9ZbjhmF4kiR5ADxZlid2g/FvbpIktSSsRqMRjEkC8EzTnOj1ZVluuX6tVmtJ9pNy1s+nqmroSWWa1/U8L/hc28/vr90UBCGy8Wb/9y+KYiTXo4uByWwK/ITV65upP5lFEITQE5p/U0un0z1fY5pm30k44/B/dlVVe77GT6iTSuZ+S7Bfsu4X37xd1/O+SyDdkqX/99bvbyIsjUYjaBFO8gsTXTxj1WakydB1HbIsw3XdYJ1fGMrlcrD4vbmWZjtZlicybdrfvkcUxb7jsaqqQhTFnpOAohDlePG0r7u6ugrgu9qnk7S3twdFUbCxsQHgdGycKAxMZjPKv6lZljVW8eRmftFkVVXPrJKiaVroM0M3Nzdb4jjr+tNQKpVC+7yHYVlWMPkpav7PG0XBgWKxiHw+j3Q6DUEQUC6Xx96DkAhgMptZ/malQDjfmF3XDW5a/jfxfkRRDPXmZtt2cNMaJEnm8/nIl1m4rju1lsK0rutvQHtWazkMfqvc/7v2W//+lxyicTCZzTD/Zh5GS6G5akqvJOE4DhRFQSqVCh7N6/nG0fwzjLqzQNgymQwSiQQSiQQWFhawsLAQSVdb+3VjsVjQ/TtppVIpuGYsFoOiKEin06hWqxP/8rC5uRl0LwII/ram1SKl84XJbA74C8nD0quLURRFFItFCIIAy7JQr9fP9ZiGYRio1Wqo1WpoNBpoNBqRdLW1X9fzvMjWSObz+eCatVoNuq6jXC5jYWFhognV7xlo/nxFUQzGhqNK5nR+MZnNML9bLpFIjH2u5tbQWa0PvxsyzBZU801sGmNSgxAEIbSW6LCaWyxREUURqqrCMAwApy3GSf1u/LFIv0XqP/weA3Y10riYzGaU4zjBjeXOnTtjn08QhKBFFvYO4INo7sJq39F7lkxjrA44TfbTKrzb3CqcVAupWCyi0WgELdLmlqkgCLBtO5IuXjq/mMxmlD8Y78/6CoN/syyVSlOZQaaqKgBOxz7LNMaQ/L+xsLu0gdMvL/26b/0va/y7oHEwmc0gv0alKIp914MNS9d1SJIE13Wnsh+drusQBAGO45w59V7TtAv5TX0aSxKaZ5qmUqnQz69pGrLZbM/nmyeCcJo+jYrJbAr8/2Hb/8e1bRuKoiCTyUCWZVSr1dB3za5Wq0in0yiXy0gkErAsqyUO27aDb8jJZDLUa/vXF0URhUIByWSy5fp+selkMonV1dWJTcbo9flPWr/ruq4LTdNQKBQmsiefvwC9fSG6bdvBrhjpdDr0a1uWBdu2+/4um5/j2BmNbNolSC6SXoWG/WLDkiR5+Xx+KoWG/RJD6XQ6kgKwfvFdv8ST//NPsuhtt4K/qqpOvMhu+3X9azf/7puPh8k0zY5Cw83XFUXRk2V5ImWlDMMIfrZe1/D/Dpvjm1R9SjrfYp7neRHnTyIiolCxm5GIiOYekxkREc09JjMiIpp7TGZERDT3mMyIiGjuMZkREdHcYzIjIqK5x2RGRERzj8mMiM4Nf9PRWCyGQqEQyjlLpVLLeWk2sQII0TmnKAosy4JpmlPZ3iZKsVgMk7ylTfr8NDq2zGhuWZYFRVGwsLDQ8s05lUpNpfr8rNrb24PjODO3j1wqlUIsFsPCwkKwWaf/u1xYWOjYKLRcLgfPNb+OlfYJYMuMzgHLsoKtS9LpdLBzMp3yW2UbGxuh78IwLk3T4LpusFOD/yXE38+vXblcxu7ubtDCbH9dt5ZTuVzG5uZmsKWQIAjI5/NQFCU4j39ef3NSWZZhmmbH9dkym2HTqnBMFJZGoxFUXFdVddrh0JBEUfRqtZrneaeV9vtVza/VakHV/2563dKa/0Z6nb9arZ5ZtZ+3zNnFbkY6V65evTrtEGhImqYFLSxZlvt2hzqOA8dxgg09ByUIAtLpNAD0PL+maZBlmTtezykmMyKaqnw+H2zSKggCRFHsucu4fzyfzw99HX+3a8dxOsbZCoUCHMdhF/UcYzKjc822bWiaFkwUcBwHmUwmmHTgj5F04782lUohkUggmUyiVCq1vMZ1XRQKBSQSCQCn3/r9czdPYPDH9fznMplM3x2nU6kUMpkMkskkFEXp+lr/nP7u5AsLCx27g7uui1KphGQy2TEpZpzPxnXdls/GnzE5KkVRgl2mNU3rueP07u4u8vn8SGN/fssMaG2dWZaFzc1NmKY5c2OKNIRp93MSjat5PETX9eC4aZotOzwbhhHs5t183B+vaVYsFj1RFFt2/fbf44+pFIvFll2iTdMMdu4GEOzYnc/nvXQ6HexobRhGsKt087UbjYYnCIKXTqc7fjZJklri88d3mt/v7+zc/Jrmnc2bxxPH+WxqtZonCELLjuTNP7e/a/iwBEEIPiNRFDt2APc/i24x+c66pfk/m/8Z+z/LoDtt85Y5u/ibobnXK5m1P9ecJDzvuxtb+6QR0zSD5NTtePMNzU8qzedvNBqeruteo9EIEl679puq53merutdJyBIktRxE8/n8z3P2+ta7T/nKJ9N83PdPhs/IbUnokHouh78/pr/u/n5bj9fs7OSjf8Z+3GKotiSlM/CZDa7+JuhudcvmXme17OV4d/Y2m/k/WbLqaracvNrvnZzK87X3tJqv3bzzdFv2bUnUT95NB/3W1ztN+JuLYx8Pt8zMQ372fjv6fb5+OcaJZF53nctU1/7Nbp9Nt1i6KdWqwVxiqLYd+biKOen6eGYGV1Y/vhI83iUP1tOkqSu79F1vefkg/b3+BMNLMsKFgX7j2KxCEEQIAhCMLYmyzJqtRpkWQ7OUSqVUKlUOq7lT2ZQFKVlTK35vePo9tkA303AqNfrHe/xf/72xc7DXDOfzwdlqNLpdPDf/hjXuD+fKIrBzyYIAmcuniNMZkRNRr0R9ztXPp9HrVbreDQaDTQajY4SU/4kEEVRsLKygpWVlY5zp9PpIKmWSiUsLCx0TE6ZhOYk157oFhcXAWCsklkbGxtBgmn+b13XQ6nq0hz3sNP7abYxmRE18W90YSS1Uc5VKpWwtLSEbDaLYrHYs4UIAMViEYZhBAnGb6VNkiiKQbLa29trea5er0OSpLFmBAqCAFmWUSqVgv/WNA2VSmWk6fjtmmMOqxVLs4HJjKiJf6O2bbtvzb9B6gH65zpryrqf7AqFAhRFga7rfZNYs3Q6jUajAVVVAZwmw0nXpfTXYmmaFsRuWVZo67SaF1FrmoZCoRBKIgO+i705KdP5wGRG1KQ5ifRKCoqidB0zauffLP21Xt2Uy+VgHMpfW9WtW7FbDM10XQ+65CZdUFiSJJimicXFxWCNm2maODg4CCVBiKIISZJQLpeDpBNWizOssTeaPf827QCIwnR0dDT2OVRVRaFQCBKQrustXXl+lYqz+BMaSqVScDNubmGUSiUUi0VUq1UA37X2LMsKkqrrusEEkObWYL1eh23bLcnXL547SDIch1+Ut1arTewauq4jk8lAEARIkhRKkmxO8n5hajpHpj2dkmhczeu/2qeSN0/Fbp/Wrapqz2nm/tou/+Evjm5f59R8/m6Lef21TM3n8f/dvjC4+Zr+4uV0Oh1MwxdF0Uun016tVvPS6bQnSVLL+/3PoT2Obmvaxvlsmj+X9ockSZ5hGB3vGYUsyx0L18/S75bWvBh81OUDvGXOLv5maG6ZphksHm6+ocqy7Kmq6lWr1ZYbWPMCWf9m3ZwE229wqqoGCUaSpI41Xbqut1xbFMVgsXS7fD7fksS6rYdrNBpBvM3Xq1arQRL0E4Wf4Px1bH7ya77x+59P88+pqqpXq9XG+mz89XCyLHuSJHV8/uix5m5YfozD6JZs/MXW7X8juq73rSYy6PlpNnA/MyIaSiaTwcbGRsckFb++o7+GrteeZJPEnaYvLk4AIaKBZTIZAJ0LxAEE41u6rnP3Z4ockxkRDcRxHJTL5WBxdC/NO38TRYXJjIiGsre313MhuGVZME2zZbsVoihwzIyIBpbJZIJ9zmRZbqnHaNv21Hdq5pjZxcWWGRENnrxr+wAAAFZJREFUzDCMoMxWpVJBoVAIuh5N05yJwr2xWAyxWCwoUjyuUqkUnJNmF1tmREQ099gyIyKiucdkRkREc4/JjIiI5h6TGRERzT0mMyIimntMZkRENPf+f+5qvd7sZIONAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 460.8x345.6 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = \"DGCFEBA\"\n",
    "MEANCOL, STDCOL = \"Vm\", \"Vs\"\n",
    "\n",
    "\n",
    "start_plotting(cmbright=False, font_size=17*1.5)\n",
    "\n",
    "ax1 = plt.subplot(211)  \n",
    "for i, LABEL in enumerate(labels):\n",
    "  a,b,c,d = LABEL2CFG[LABEL]    \n",
    "  shift = i*4    \n",
    "      \n",
    "\n",
    "  d11 = df[(df.a==a)&(df.b==b)&(df.c==c)&(df.d==d)]\n",
    "  exp, var = theoretical_moments(a,b,c,d)\n",
    "  Z, LABEL = 1.0/var * 100.0, r\"Bias [\\%]\"\n",
    "  #Z, LABEL = 1.0, \"Bias\"\n",
    "             \n",
    "  d1 = d11[(df.nl==10)&(df.no==1000)]\n",
    "  xs = list(np.arange(len(d1))+shift+1.6)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-var)*Z, color=\"lightskyblue\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==100)&(df.no==100)]\n",
    "  xs = list(np.arange(len(d1))+shift+0.8)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-var)*Z, color=\"mistyrose\", width=0.8, align=\"center\")        \n",
    "        \n",
    "  d1 = d11[(df.nl==1000)&(df.no==10)]\n",
    "  xs = list(np.arange(len(d1))+shift)\n",
    "  plt.bar(xs, abs(d1[MEANCOL]-var)*Z, color=\"lightgrey\", width=0.8, align=\"center\")\n",
    "\n",
    "plt.xticks([],[])\n",
    "plt.ylabel(LABEL)\n",
    "plt.axhline(0.0, color=\"k\", zorder=-1, lw=1, ls=\"--\")\n",
    "#plt.yscale(\"log\")\n",
    "#plt.ylim(0.1, 100000)\n",
    "plt.xlim(-1, 43)\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_lines = [Line2D([0], [0], color=\"lightgrey\", lw=4),\n",
    "                Line2D([0], [0], color=\"mistyrose\", lw=4),\n",
    "                Line2D([0], [0], color=\"lightskyblue\", lw=4),]\n",
    "plt.gca().legend(custom_lines, [r\"$10^3 \\times 10^1$\", r\"$10^2 \\times 10^2$\", r\"$10^1 \\times 10^3$\"], \n",
    "                 title=r\"$S_{Z} \\times S_y$\", fontsize=15, loc=4)\n",
    "\n",
    "ax1 = plt.subplot(212)  \n",
    "for i, LABEL in enumerate(labels):\n",
    "  a,b,c,d = LABEL2CFG[LABEL]    \n",
    "  shift = i*4    \n",
    "    \n",
    "  d11 = df[(df.a==a)&(df.b==b)&(df.c==c)&(df.d==d)]\n",
    "  exp, var = theoretical_moments(a,b,c,d)\n",
    "  Z, LABEL = 1.0/var * 100.0, r\"Variance [\\%]\"\n",
    "  #Z, LABEL = 1.0, \"Variance\"    \n",
    "     \n",
    "  d1 = d11[(df.nl==10)&(df.no==1000)]\n",
    "  xs = list(np.arange(len(d1))+shift+1.6)\n",
    "  plt.bar(xs, d1[STDCOL]**2*Z, color=\"lightskyblue\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==100)&(df.no==100)]\n",
    "  xs = list(np.arange(len(d1))+shift+0.8)\n",
    "  plt.bar(xs, d1[STDCOL]**2*Z, color=\"mistyrose\", width=0.8, align=\"center\")\n",
    "\n",
    "  d1 = d11[(df.nl==1000)&(df.no==10)]\n",
    "  xs = list(np.arange(len(d1))+shift)\n",
    "  plt.bar(xs, d1[STDCOL]**2*Z, color=\"lightgrey\", width=0.8, align=\"center\")\n",
    "\n",
    "\n",
    "plt.xticks(np.arange(0, 4*len(labels), 4)+1, labels)  \n",
    "plt.ylabel(LABEL)\n",
    "plt.xlabel(r\"Increasing $\\mathbb{V}[Y]$\")\n",
    "plt.axhline(0.0, color=\"k\", zorder=-1, lw=1, ls=\"--\")\n",
    "plt.yscale(\"log\")\n",
    "plt.ylim(0.001, 1e6)\n",
    "plt.xlim(-1, 43)\n",
    "\n",
    "\n",
    "from matplotlib.lines import Line2D\n",
    "custom_lines = [Line2D([0], [0], color=\"lightgrey\", lw=4),\n",
    "                Line2D([0], [0], color=\"mistyrose\", lw=4),\n",
    "                Line2D([0], [0], color=\"lightskyblue\", lw=4),]\n",
    "plt.gca().legend(custom_lines, [r\"$10^3 \\times 10^1$\", r\"$10^2 \\times 10^2$\", r\"$10^1 \\times 10^3$\"], \n",
    "                 title=r\"$S_{Z} \\times S_y$\", fontsize=15, loc=4)\n",
    "\n",
    "plt.subplots_adjust(left=0.2, right=0.95, top=0.95, bottom=0.2)    \n",
    "plt.savefig(\"RESULTS/nsamples_v.pdf\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "HPF Matching priors",
   "private_outputs": true,
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}