{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import sqrt\n",
    "from numpy.linalg import inv\n",
    "from numpy import log\n",
    "from numpy.linalg import norm\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "os.chdir('/Users/sifanliu/Dropbox/Dual Sketch/RidgeRegression')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a uniformly distributed orthogonal matrix\n",
    "def generate_haar_matrix(n, p):\n",
    "    if n <= p:\n",
    "        return np.linalg.qr(np.random.randn(p, n))[0].T\n",
    "    else:\n",
    "        return np.linalg.qr(np.random.randn(n, p))[0]\n",
    "\n",
    "\n",
    "def theta_1(lbd, gamma):\n",
    "    z_2 = -0.5 * (sqrt(gamma) + (1 + lbd) / sqrt(gamma) + sqrt((sqrt(gamma) + (1 + lbd) / sqrt(gamma)) ** 2 - 4))\n",
    "    return -(1 + lbd) / (gamma * lbd) - z_2 / (sqrt(gamma) * lbd)\n",
    "\n",
    "\n",
    "def theta_2(lbd, gamma):\n",
    "    delta = (sqrt(gamma) + (1 + lbd) / sqrt(gamma)) ** 2 - 4\n",
    "    return -1 / (gamma * lbd ** 2) + (gamma + 1) / (2 * gamma * lbd ** 2) - 1 / (2 * sqrt(gamma)) * (\n",
    "            (lbd + 1) / gamma + 1) / (lbd * sqrt(delta)) + 1 / (2 * sqrt(gamma)) * sqrt(delta) / (lbd ** 2)\n",
    "\n",
    "\n",
    "# MSE of ridge regression, no sketching\n",
    "def MSE_original(lbd, gamma, alpha=1, sigma=1, verbose=0):\n",
    "    the_1 = theta_1(lbd, gamma)\n",
    "    the_2 = theta_2(lbd, gamma)\n",
    "    bias = alpha ** 2 * lbd ** 2 * the_2\n",
    "    variance = gamma * sigma ** 2 * (the_1 - lbd * the_2)\n",
    "    if verbose == 0:\n",
    "        return bias + variance\n",
    "    else:\n",
    "        return bias + variance, bias, variance\n",
    "    \n",
    "\n",
    "# MSE of orthogonal dual sketch\n",
    "def MSE_dual(lbd, gamma, zeta, alpha=1.0, sigma=1.0, verbose=0):\n",
    "    the_1 = theta_1(lbd, zeta) * zeta + (1 - zeta) / lbd\n",
    "    the_2 = theta_2(lbd, zeta) * zeta + (1 - zeta) / lbd ** 2\n",
    "    bias_1 = alpha ** 2 / gamma * (\n",
    "            1 + (lbd + zeta - gamma) ** 2 * the_2 + 2 * (gamma - zeta - lbd) * the_1 + (\n",
    "            gamma - zeta) * the_1 ** 2)\n",
    "    bias_2 = -2 * alpha ** 2 / gamma * (1 - (lbd + zeta - gamma) * the_1)\n",
    "    bias = bias_1 + bias_2 + alpha ** 2\n",
    "    variance = sigma ** 2 * (the_1 - (lbd + zeta - gamma) * the_2)\n",
    "    if verbose == 0:\n",
    "        return bias + variance\n",
    "    else:\n",
    "        return bias + variance, bias, variance\n",
    "\n",
    "\n",
    "# compute the Stieltjes transform of MP\n",
    "def MP_Stieltjes(z, gamma):\n",
    "    return (1 - gamma - z - sqrt((1 + gamma - z) ** 2 - 4 * gamma)) / (2 * gamma * z)\n",
    "\n",
    "\n",
    "# compute the Stieltjes transform of inv-MP\n",
    "def inv_MP_Stieltjes(z, gamma):\n",
    "    return -1 / z - 1 / z ** 2 * MP_Stieltjes(1 / z, gamma)\n",
    "\n",
    "\n",
    "# compute the R-transform of MP\n",
    "def MP_R(z, gamma):\n",
    "    return 1 / (1 - gamma * z)\n",
    "\n",
    "\n",
    "# compute the R-transform of inv-MP\n",
    "def inv_MP_R(z, gamma):\n",
    "    a = (z + 1) * gamma\n",
    "    b = -z * gamma * (1 + gamma)\n",
    "    c = (z * gamma) ** 2\n",
    "    w = (-b + sqrt(b ** 2 - 4 * a * c)) / (2 * a)\n",
    "    w2 = (-b - sqrt(b ** 2 - 4 * a * c)) / (2 * a)\n",
    "    # print(w2)\n",
    "    # w2 = (z * (gamma + 1) - sqrt(z ** 2 * (gamma + 1) ** 2 - 4 * gamma * (z + 1) * z ** 2)) / (2 * (z + 1) * gamma)\n",
    "    num = 8 * (z + 1) ** 2 * gamma ** 2\n",
    "    den = z * (gamma + 1 + sqrt((gamma + 1) ** 2 - 4 * gamma * (z + 1)))\n",
    "    # return num / den - 1 / z\n",
    "    return 1 / w2 - 1 / z\n",
    "\n",
    "\n",
    "# compute the inverse of the Stieltjes transform of inv-MP\n",
    "def inv_MP_Stieltjes_inv(z, gamma, xi, lbd):\n",
    "    return 1 / (1 + z / xi) - (gamma - 1 - sqrt((gamma - 1) ** 2 + 4 * lbd * z)) / (2 * z) - 1 / z\n",
    "\n",
    "\n",
    "# compute the Stieltjes transform of inv-MP evaluated at zero\n",
    "def inv_MP_Stieltjes_0(gamma, xi, lbd):\n",
    "    maxiter = 100\n",
    "    tol = 1e-8\n",
    "    low = 0.01\n",
    "    high = 2\n",
    "    if inv_MP_Stieltjes_inv(low, gamma, xi, lbd) * inv_MP_Stieltjes_inv(high, gamma, xi, lbd) > 0:\n",
    "        print(\"no root\")\n",
    "        return\n",
    "    err = 1\n",
    "    old = 0\n",
    "    for i in range(maxiter):\n",
    "        mid = (low + high) / 2\n",
    "        if inv_MP_Stieltjes_inv(low, gamma, xi, lbd) * inv_MP_Stieltjes_inv(mid, gamma, xi, lbd) < 0:\n",
    "            high = mid\n",
    "        else:\n",
    "            low = mid\n",
    "        if abs(old - inv_MP_Stieltjes_inv(mid, gamma, xi, lbd)) < tol:\n",
    "            print(\"converged\")\n",
    "            break\n",
    "        old = inv_MP_Stieltjes_inv(mid, gamma, xi, lbd)\n",
    "    return mid\n",
    "\n",
    "\n",
    "# compute the first order derivative of the inverse function of the Stieltjes transform of inv-MP\n",
    "def diff_MP(z, gamma, xi, lbd):\n",
    "    epsilon = 1e-5\n",
    "    tol = 1e-8\n",
    "    maxiter = 100\n",
    "    ans = (inv_MP_Stieltjes_inv(z + epsilon, gamma, xi, lbd) - inv_MP_Stieltjes_inv(z - epsilon, gamma, xi, lbd)) / (2 * epsilon)\n",
    "    old = ans\n",
    "    for i in range(maxiter):\n",
    "        ans = (inv_MP_Stieltjes_inv(z + epsilon, gamma, xi, lbd) \n",
    "               - inv_MP_Stieltjes_inv(z - epsilon, gamma, xi, lbd)) / (2 * epsilon) \n",
    "        if abs(old - ans) < tol:\n",
    "            print('converged')\n",
    "            return ans\n",
    "        old = ans\n",
    "        epsilon = epsilon / 2\n",
    "    return ans\n",
    "        \n",
    "\n",
    "# compute the first order derivative evaluated at 0 of the Stieltjes transform of inv-MP\n",
    "def inv_MP_Stieltjes_derivative_0(gamma, xi, lbd):\n",
    "    m_0 = inv_MP_Stieltjes_0(gamma, xi, lbd)\n",
    "    print(m_0)\n",
    "    return 1 / diff_MP(m_0, gamma, xi, lbd)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 500\n",
    "gamma = 0.9\n",
    "p = int(n * gamma)\n",
    "alpha = 1\n",
    "sigma = 1e-7\n",
    "num_steps = 20\n",
    "zeta_seq = np.linspace(1, 2, num_steps)\n",
    "dual_simu = np.zeros((num_steps, 3))\n",
    "dual_theo = np.zeros((100, 3))\n",
    "rep = 10\n",
    "lbd_seq = [0.1, 0.3, 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# orthogonal dual projection\n",
    "for j in range(1):\n",
    "    lbd = gamma * sigma ** 2 / alpha ** 2\n",
    "    for i in range(num_steps):\n",
    "        zeta = zeta_seq[i]\n",
    "        d = int(n * zeta)\n",
    "        for k in range(rep):\n",
    "            X = np.random.randn(n, p)\n",
    "            beta = np.random.randn(p, 1) * alpha / sqrt(p)\n",
    "            epsilon = np.random.randn(n, 1) * sigma\n",
    "            Y = X @ beta + epsilon\n",
    "            R = generate_haar_matrix(p, d)\n",
    "            beta_dual = X.T / n @ inv(X @ R @ R.T @ X.T / n + lbd * np.identity(n)) @ Y\n",
    "            dual_simu[i, j] = dual_simu[i, j] + np.linalg.norm(beta_dual - beta) ** 2 / rep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "zeta_seq_2 = np.linspace(0.1, 0.5, 100) \n",
    "for j in range(1):\n",
    "    lbd = gamma * sigma ** 2 / alpha ** 2\n",
    "    for i in range(100):\n",
    "        zeta = zeta_seq_2[i]\n",
    "        dual_theo[i, j] = MSE_dual(lbd, gamma, zeta, alpha, sigma, verbose=1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 500\n",
    "p = 700\n",
    "d = 600\n",
    "gamma = p / n\n",
    "xi = d / n\n",
    "alpha = 1\n",
    "sigma = 1e-8\n",
    "num_steps = 10\n",
    "zeta_seq = np.linspace(0.1, 2, num_steps)\n",
    "rep = 50\n",
    "dual_simu_gaus = np.zeros((rep, num_steps))\n",
    "lbd = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gaussian dual projection\n",
    "for i in range(num_steps):\n",
    "    zeta = zeta_seq[i]\n",
    "    d = int(n * zeta)\n",
    "    for k in range(rep):\n",
    "        X = np.random.randn(n, p)\n",
    "        beta = np.random.randn(p, 1) * alpha / sqrt(p)\n",
    "        epsilon = np.random.randn(n, 1) * sigma\n",
    "        Y = X @ beta \n",
    "        R = np.random.randn(p, d) / sqrt(d)\n",
    "        # R = generate_haar_matrix(p, d)\n",
    "        beta_dual = X.T / n @ inv(X @ R @ R.T @ X.T / n + lbd * np.identity(n)) @ Y\n",
    "        dual_simu_gaus[k, i] = norm(beta_dual - beta) ** 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged\nconverged\n1.1871523558348418\nconverged\nconverged\nconverged\n1.152420945242047\nconverged\nconverged\nconverged\n1.1195826909691098\nconverged\nconverged\nconverged\n1.0885865003615618\nconverged\nconverged\nconverged\n1.0593732743710285\nconverged\nconverged\nconverged\n1.0318772419542075\nconverged\nconverged\nconverged\n1.0060270275920624\nconverged\nconverged\nconverged\n0.9817471339553592\nconverged\nconverged\nconverged\n0.9589591873437164\nconverged\nconverged\nconverged\n0.9375833610445257\nconverged\nconverged\nconverged\n0.9175394428521395\nconverged\nconverged\nconverged\n0.8987479322403668\nconverged\nconverged\nconverged\n0.8811308335885404\nconverged\nconverged\nconverged\n0.864612575434148\nconverged\nconverged\nconverged\n0.8491203514859078\nconverged\nconverged\nconverged\n0.8345846988633274\nconverged\nconverged\nconverged\n0.8209396686032413\nconverged\nconverged\nconverged\n0.8081230480596424\nconverged\nconverged\nconverged\n0.7960763757303355\nconverged\nconverged\nconverged\n0.7847450153902171\nconverged\nconverged\nconverged\n0.7740779633447528\nconverged\nconverged\nconverged\n0.7640278336033224\nconverged\nconverged\nconverged\n0.7545506651327012\nconverged\nconverged\nconverged\n0.7456057142838834\nconverged\nconverged\nconverged\n0.7371553510054945\nconverged\nconverged\nconverged\n0.7291647771373391\nconverged\nconverged\nconverged\n0.7216019077971578\nconverged\nconverged\nconverged\n0.7144371193274854\nconverged\nconverged\nconverged\n0.707643101029098\nconverged\nconverged\nconverged\n0.7011946920678018\nconverged\nconverged\nconverged\n0.6950686590746045\nconverged\nconverged\nconverged\n0.6892435923591256\nconverged\nconverged\nconverged\n0.6836997576430441\nconverged\nconverged\nconverged\n0.6784189033135772\nconverged\nconverged\nconverged\n0.673384215943515\nconverged\nconverged\nconverged\n0.6685801571980117\nconverged\nconverged\nconverged\n0.6639923303946853\nconverged\nconverged\nconverged\n0.6596074545569717\nconverged\nconverged\nconverged\n0.6554132161475716\nconverged\nconverged\nconverged\n0.6513982245884835\nconverged\nconverged\nconverged\n0.6475519010610877\nconverged\nconverged\nconverged\n0.6438644377328455\nconverged\nconverged\nconverged\n0.6403267310373484\nconverged\nconverged\nconverged\n0.6369303075410426\nconverged\nconverged\nconverged\n0.6336673017032446\nconverged\nconverged\nconverged\n0.630530374329537\nconverged\nconverged\nconverged\n0.6275127051584422\nconverged\nconverged\nconverged\n0.6246079187281428\nconverged\nconverged\nconverged\n0.6218100695498288\nconverged\nconverged\nconverged\n0.6191136050410568\nconverged\nconverged\nconverged\n0.6165133358724415\nconverged\nconverged\nconverged\n0.6140043989010153\nconverged\nconverged\nconverged\n0.6115822497569023\nconverged\nconverged\nconverged\n0.6092426109500226\nconverged\nconverged\nconverged\n0.6069815015234052\nconverged\nconverged\nconverged\n0.6047951555065809\nconverged\nconverged\nconverged\n0.6026800515688955\nconverged\nconverged\nconverged\n0.6006328759528694\nconverged\nconverged\nconverged\n0.5986504928208889\nconverged\nconverged\nconverged\n0.5967299813218416\nconverged\nconverged\nconverged\n0.5948685688711703\nconverged\nconverged\nconverged\n0.5930636385641992\nconverged\nconverged\nconverged\n0.5913127365894615\nconverged\nconverged\nconverged\n0.5896135351620613\nconverged\nconverged\nconverged\n0.5879638251103461\nconverged\nconverged\nconverged\n0.5863615455292166\nconverged\nconverged\nconverged\n0.5848047170601784\nconverged\nconverged\nconverged\n0.5832914789579808\nconverged\nconverged\nconverged\n0.5818200742639601\nconverged\nconverged\nconverged\n0.5803888275660574\nconverged\nconverged\nconverged\n0.5789961524121463\nconverged\nconverged\nconverged\n0.5776405438967049\nconverged\nconverged\nconverged\n0.5763205712474881\nconverged\nconverged\nconverged\n0.5750348704122006\nconverged\nconverged\nconverged\n0.573782158885151\nconverged\nconverged\nconverged\n0.5725612060539425\nconverged\nconverged\nconverged\n0.5713708331994711\nconverged\nconverged\nconverged\n0.5702099357359112\nconverged\nconverged\nconverged\n0.569077446144074\nconverged\nconverged\nconverged\n0.5679723562113941\nconverged\nconverged\nconverged\n0.5668936947919427\nconverged\nconverged\nconverged\n0.5658405426330864\nconverged\nconverged\nconverged\n0.5648120249621569\nconverged\nconverged\nconverged\n0.5638072892464696\nconverged\nconverged\nconverged\n0.562825549673289\nconverged\nconverged\nconverged\n0.5618660278432068\nconverged\nconverged\nconverged\n0.5609279824234543\nconverged\nconverged\nconverged\n0.5600107239745555\nconverged\nconverged\nconverged\n0.5591135704703629\nconverged\nconverged\nconverged\n0.5582358843646942\nconverged\nconverged\nconverged\n0.5573770355246961\nconverged\nconverged\nconverged\n0.5565364457108077\nconverged\nconverged\nconverged\n0.5557135366834698\nconverged\nconverged\nconverged\n0.5549077672697604\nconverged\nconverged\nconverged\n0.5541186111234129\nconverged\nconverged\nconverged\n0.5533455641381442\nconverged\nconverged\nconverged\n0.5525881518609821\nconverged\nconverged\nconverged\n0.5518459072522819\nconverged\nconverged\nconverged\n0.5511183855123818\nconverged\nconverged\nconverged\n0.5504051492549478\nconverged\n"
     ]
    }
   ],
   "source": [
    "# theoretical part\n",
    "zeta_seq_2 = np.linspace(0.1, 2, 100)\n",
    "theo_gaus = np.zeros(100)\n",
    "for j in range(1):\n",
    "    lbd = gamma * sigma ** 2 / alpha ** 2\n",
    "    lbd = 1\n",
    "    for i in range(100):\n",
    "        zeta = zeta_seq_2[i]\n",
    "        bias_1 = alpha ** 2\n",
    "        bias_2 = 2 * alpha ** 2 / gamma * inv_MP_Stieltjes_0(gamma, zeta, lbd)\n",
    "        bias_3 = alpha ** 2 / gamma * inv_MP_Stieltjes_derivative_0(gamma, zeta, lbd)\n",
    "        theo_gaus[i] = bias_1 - bias_2 + bias_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEfCAYAAABf1YHgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXl8FOX9gJ9vDhLCDYGAXEFOQTwA8ayCiOKNB6g9FI+qtWrVqpXWerZKFbWt0nq0ivbXqqhYqaKWQ8SDCiESIyEQAiEkxIRAMIRsjs2+vz9md9lsNskmszv7Bt7n89lPdmbeeed5Zyb7znvM+4pSCoPBYDAYmiMu1gIGg8Fg0BuTURgMBoOhRUxGYTAYDIYWMRmFwWAwGFrEZBQGg8FgaBGTURgMBoOhRUxGYTAYDIYWMRmFw4jITBH5r4jsEZE6ESkWkTdE5NRYu4VCRB4SkfIYOxwtIkpEpkQh7ineuI+2EcccbxxdbbpMFpGH2rGfrTSIyEARqRKRI9uzf0dERApEZH6E47xHRFZEMk5dMBmFg4jIM8A7QDFwA3AWcB/QDfhcRIbHUK85/gacE2uJw4TJwIMxOO79wH+UUtticOxYcQnw5wjH+TwwIRoPNLEmIdYChwsicjFwB3CtUmph0OZ/iMiFgMtxsVZQShUBRbH2MEQHEekOXANc7MCxOiultLjHlVJfRyHO/SLyDnAbsCrS8ccSU6JwjjuAdSEyCQCUUv9RSu3yLYvIL0VknYh8LyKlIvIfERkRuE+o4nNwNYiIJIrIfBEpFJFaEdklIu+KSCfv9p4i8jfv+hpvuJcC4mtU9SQiXUTkORHZLCLVIrJdRBZ4f3ACPZSI/EJEHhOR3SJS5g2X1NqJEpFbRGSniBwQkf8AA4K2p3vjvyBo/UIRyQhYHuOt1tvpdd0oIneISJvu+9bOYTP73OM9nxcFrDtNRD71uuwRkZdEpJt32xzgWe935f2sCtj3GO89sM9bTbRWRKYHHTZVRN7ybt8mIreEkbzZWA8oKwOOtTDAQTXj01tEXvDemzUi8qWInBh0DpSI3CUifxSR3UB2wLZbRSTPez63isidrYmKyCoReVtEfujdp1JEPhSRQUHhUkXkVe85rvbuNykoTKP/HREZJyIfiche7323SUR+HrTPxSKS4U3vdyLyhIgkBmm+A1wgIr1bS09HwpQoHEBEEoCTgbbUiQ4CngN2AN2Bm4EvRGSUUur7NsQzF/gRVhXXdqA/cB4Q793+NHAKcCfwHTAYOL2F+FK8+/4G2O0N/xvgLZpWUf0S6wfox8AxwOPe9DzRXORilbwWYBXj/w2cAbwcTkJDMBDYDPwT2A8cBzwMdPa6hEtr57ARIvJb7z4XK6U+9q47FViBlabLgT7APKCXd/kD4Cmsc3ayN6pK775jgC+8abkZ2ANMwjr3gbwEvAq8CFwFLBCRDKXU2hbSNg1Yq5RqCFj3KNb599EfeB3Y4vVJApYDPYF7gDLgZ8ByERmplPouYN97gNXAT/A+mIrIT7EyxaeBj4GpwFMikqSUmteCK8CJwBFY56kz8Cdves8LCPNvYARwN1DudfhERI5XSm1tJt4lQC7WvVoLjMb6v8PrPNt7Dl4Afg0Mx7qH4rzH8fElkAj8AHivlbR0HJRS5hPlD5AGKOCmoPWClVn7PtLM/vFY/xT7gasD1hcA84PCzvEeq6t3+X3gqRbcvgVua2H7Q0B5C9sTgFO9xxwSsF4Bq4PC/hv4Xyvnai3wYdC6l7zxTfEup3uXLwgKtxDIaCZe37n+NbAtYP0Ub1xHt+DU2jn0n3PgMe91mhIU5jPgk6B1ZwYeG7jV+pdsEv/rWNV/nZs5vi8NjwSsS8TKyOe1cr63AE+2sD0R+Nx7n3TxrrseqANGBt0H+YFxeZ2+DoovDquN7pWg9X8BvgeSW3BZ5Q3TK2DdHd7jdPYuz/AunxEQpov3XLwQ6n8HSPXuM76Fe2dHCOfrsEpjfYLWFwC/b+m8d7SPqXpyBvH+DR6q95dAfcDHX9QVkZNEZJmI7AHcQDXWD9GoNh57AzBHRO71Vl9IiO33eKt7wopbRH4iIl+LSJXX+3PvpuD9/xu0nINVUmou3njgeJo+iS0OxytEfMki8rCIbMV6SqwHfg8M85bywqW1c+jjaeAW4Byl1KoAjxSsUsIiEUnwfbDOWz0wsZXjnwm8qVqv3/efb6VUPZBHC+fbS3+sp+7m+DNwNHCJUuqAd91ZwHpge0BaAD7FKukE8kHQ8iCsEsFbQevfxHqCH9+K7zqlVEXAco7370Dv38nAbqXUp74AXu/3gdOaiXMvsBN4XkSuEJF+QdtHAUNoev1WAslY5yeQcqzzeshgMgpnKMf6oQr+p/0HcIL340dEhmD90wtwE9YT+wlYRfzkNh77d1hVObcAWcBOEflFwPZbsZ70HwA2e+uNr2wuMhG5BHgNWAPMAk7C6kFCCLd9Qct1rfj3xXoyLQtaH7wcLn/AqhbwVU2cgHU+aMUjmNbOoY/LsH5Ag6t6emGVCv9C4weDWqwn9uAqpGD6ACVheLb1fOPdXhtqg4hcj3X//UQplRewKRXrutcHfa6laVpKg5YHNLPet9xa3X6oNMLBdA4IEbcv/pBxK6U8wNlYVa8vA9+JyGcicrw3SKr371Iap3e7d31wmmtp+/+p1pg2CgdQSrlFZA3WzfhAwPpSvDd10EPqDKy2gIt9T3HeJ5jgG70GCG5QbRRGKVXjPeYDIjISq477jyKyWSn1kVJqH3A7cLuIHAPcC/xTRL5RSuXQlFnAV0opf0OpiJwRxmkIh91YpafgJ7rg5Rrv3xbTjuX6rFLK3yYiIue3Vaq1cxgQ9AKsJ9fXROTH3h8gsH7cFFY13tIQh9gVYl0gewhq0I8ge7HaGhrhbZheADyqlPpPiH0ysNolggnOdIJL0b4ML/iapgXEbYeSEHH74m82bqVULnCZt3H6B1gPGR94G8p9+90IhOottT1ouWdLx+qImBKFc/wROFFEfhJG2M6AB+tH08dsmmbsRcBRQeuCe8L48T4V3o31zzw2xPZvsBr+4oAxLbgF/xj8qLljtgVlNahuoGlXzUuDlsuwnuj8aRerl9fJQeEauXqrtpotLYXp2NI5zAbOxcowng/Y5wDwP2C0UiojxMeXUdR5PYOfRlcAs0OsjwSbgWGBK0SkP1bvneVYmVswK7AaiwtDpCU7RPhAirAyxllB62djNd63tn9rfAX0ExF/hwxv1d/5HKwibRalVL1SaiVWNeIArB/9zVjtKunNXL89AceKw6qm2mIzHVphShQOoZR6T0T+CCwUkanAf7CqpPpw8Me9yvt3JVZVxSsi8ndgHNaPU3Cx+13gWRH5NbAO6wd1XGAAEXkXqzrka6yGt8uxrvtq7/bPvfF8i/X091PgAE2rT3wsw+pN8xusf8rzsHrORIrHgMUi8lev1xlYJSw/SimPiLwH3CkiO7DOyy9p+h7KMuDn3jaKvVhtQK12zw2mtXMY5LZWrG67H4lIpVLK1yPmXmCFiHiAt7EavIdg/YD9Rim1BavXDcAvRGQlUKmU2ozVU2sdsFpEnsIqYRwP7FFKtbdHmI8vgIuC1r2G9RLoc1gPN771ld5S5mtYpapV3i6m27Du48nAd0qpZ5o7mPfaPQS84G1/W4Z1jX8G/Npbems3SqmPReQL4E0RuQ/rXN2N9dDwZKh9vCXp+VjtJNuwqgp/BWQppfZ6w/wS632n7sCHWJn6kcBM4HKlVLU3utFYbYlf2EmHdsS6Nf1w+2DV5y/D+uGqx3q6egc4Nyjc1Vi9SFxYT6MnEtTLCat++2msutUKrK6CN9K419M9WNUE32P9OH2FVaXli+NJrKe4/Vg/uJ8APwjY/hABvZ6wMrD5WE/1lV73EwnqheRdvjUoTY3iauEc3Yr15FmNVVVzNgG9nrxh0rAavSuxeqTcSFCvJ2+Yd71hSrG65f406PxMofVeT62dwzmBcXrXzcAqdTwYsO5E4COvzwGshtingR7e7eJ13IVVolwVsO8x3nOxP8BhWktpwOol9HYr53qS91iBPdYKvPEFfwJ9emDdbzuxfjSLsDodnNrSPRB0jbd6990G3BnGfdEkPaHSjtXW9RrW/4QLq5H9hKD9CjjY66kfVnvhNqxqze+wepoNCdrnXKzeawe813ADVvtVQkCYO73xhOzB2FE/4k2cwWA4TBGRDcA/lVIhn7gN4eNti/xAKfW7VgN3IExGYTAc5ojILKyS5QillLu18IbQeDsAfAQMU1YnkUMG00ZhMBjexqpvH4hVjWdoH72Baw61TAJMicJgMBgMreBoiUJEXsbqOlimlGoydr53nJ9HOdg19A6lVKtd2lJTU1V6enqEbQ0Gg+HQZv369eVKqb6thXO0ROHt21wFvNZMRtEVOKCUUt4ua4uUUs315/czadIklZGR0VqwmJCfn8/w4TpOM2Ghux/o72j87GH87GHHT0TWK6WCh11pgqMv3CmlVtPy25FV6mDO1YWmb3V2OHr31nu0Yd39QH9H42cP42cPJ/y0ezNbRC4RkVyswcSuayHcjd6x4TNKSkooLy+npKSE4uJiKioqyM/Px+VykZOTg8fjITMzE4D169cDkJmZicfjIScnB5fLRX5+PhUVFRQXF+OLr6CggKqqKnJzc3G73WRlZTWKw/c3Ozub2tpa8vLyqKyspLCwkLKyMsrKyigsLKSyspK8vDxqa2vJzs4OGUdWVhZut5vc3FyqqqooKChwJE2bNm3SPk0HDhzQ+jqVlpZqfe/l5+dree/50vTtt99qe+95PB62bNmi7b1XXl7Orl272n2dwsXxxmwRSQfeD1X1FBTudOABpdRZrcWpc9VTSUkJAwZEa5ge++juB/o7Gj97GD972PHTsuqpLXirqYaLSGqrgTUmMTF4Aiy90N0P9Hc0fvYwfvZwwk+rjEJERvjG+heRCVijg+5peS+9qaqqaj1QDNHdD/R3NH72MH72cMLP6e6xr2ONzZIqIkXAg1jjFaGUeh5rPP+rRaQea4yWK1QHf9EjNVXvApHufqC/o/Gzh/GzhxN+jmYUSqmrWtn+B6xx4A8ZioqKGDOm1R6+MUN3P9DfURe/yspKysrKqK+vb7S+vr5e6+oT42eP5vwSExPp168f3bt3D7FX2zBDeESZESNGxFqhRXT3A/0ddfCrrKyktLSUgQMH0rlz50YTYSmlgifG0grjZ49QfkopXC4XxcXFALYzC63aKJzimWVbSL/vgyafZ5ZFfq6RjRs3RjzOSKK7H+jvqINfWVkZAwcOJCUlpcmPRlu6QcYC42ePUH4iQkpKCgMHDqSsrL0zCQfE18GbAID2d49Nv+/gvO8F89o8Q6bBoA2bNm1izJgxWj/5GpxHKUVubi5HHRU8EaZFh+8ee6jge9FFV3T3A/0ddfFrLpM4cOCAwyZtw/jZoyW/SD04mIwiykycODHWCi2iux/o76i7X5cuXWKt0CLGzx5O+JmMIsro8rTZHLr7gf6Ouvt15CdiHTB+JqOIOro/beruB/o76u7XEZ6IRaTVz6pVq1i4cCEi4uhLcB3h/EUbk1FEGd8AX7qiux/o76i7X3V1dawVWqS6upo1a9b4PytXrgTg/vvvb7R+woQJMfPTGSf8zHsUUWbUqFGxVmgR3f1Af0fd/ZKTk2Ot0CLJycmcdNJJ/mVfaWH48OGN1seK4PPncrno3LlzjGya4sT1NSWKKFNYWBhrhRbR3Q/0d9Tdr66uLtYKLdIev+3btzN9+nS6dOnCmDFjWLx4cZMw7733HpMmTSI5OZn+/ftz7733NnlrfeXKlZx44okkJyeTlpbGLbfc0qhaa9WqVcTHx/Pxxx9z0UUX0bVrV2699VZmzZrF1KlTmxzzwQcfJC0trclxookT19dkFFEmLS0t1gotorsf6O+ou5/Ow09A+/x++MMfctFFF/Huu+8ycuRIrrzySoqKivzbFy1axKWXXsrkyZNZsmQJDz74IC+++CJz5871h8nJyWHGjBmkpqbyzjvv8PDDD/Ovf/2Lyy+/vMnxrr/+eo499liWLFnC9ddfzw033MCnn37K9u3b/WGUUrz22mv8+Mc/dvScO3EsU/UUZfbt2xeRsVaihe5+oL+j1n6fPE78p/PCCzvhGrjoz43XLbkdMl8Nb/8z7oOpc1sPF4Tb7SY+Pr5N+9x5551cd501r9nEiRNJS0vj/fff5+abb0YpxT333MPVV1/NX/7yF/8+SUlJ/PznP2fu3Ln06dOHRx55hKFDh7JkyRL/8Xv37s0VV1zBmjVrOPnkk/37zpo1i0cffdS/7PF4GDx4MAsXLuThhx8G4JNPPqGgoIBrr722zefADu05f23FlCiiTEeoH9Yd3R1199OduLi2/wydffbZ/u99+vShX79+/hLFli1bKCwsZPbs2bjdbv/nzDPPpKamxj+j3tq1a7nkkksa/chedtllJCQk8Pnnnzc63vnnNx65IS4ujjlz5vDaa6/hG91i4cKFTJo0iaOPbnFOtojTnvPX5mNE/QgGg8EQYXr27NlouVOnTtTU1ABQXl4OwHnnnUdiYqL/M2zYMAB27twJWDPDBVcbxsfH06dPH/bu3dtofajqxWuvvZYdO3bwySefsH//ft555x1/KedQw1Q9RRnfzasruvuB/o5a+02dS+0pd5GUlNS+/S/6c9PqqAjj8XgiGl/v3r0BePHFFzn++OObbPdlGAMGDGgyYF5DQwN79uzxx+Ej1FAY6enpnHXWWSxcuJDt27fj8Xi46qoWZ1KICpE+f6EwGUWUCX7y0Q3d/UB/R939EhL0/jePtN/o0aMZOHAgBQUF/PSnP2023Iknnsi7777LY4895q9+Wrx4MW63m9NOOy2sY11//fVcd911bNy4kZkzZ8bkXnDi+jpa9SQiL4tImYh828z2H4nIN97PlyJyrJN+0aC0tDTWCi2iux/o76i7n5NdNdtDpP3i4uJ46qmnePLJJ7nttttYunQpy5cv58UXX+S8887zv6B2//33U1BQwMyZM1m6dCkvvvgiN954I+ecc06jhuyWmDlzJsnJyWRmZjreiO3Dievr9KPGQuA54LVmtm8HzlBKVYjIucCLwIkOuUWFIUOGxFqhRXT3A/0ddffr1KlTrBVaJBp+V1xxBd27d+exxx7j5ZdfJj4+niOPPJILLrjAf7xx48bx4Ycf8utf/5pLL72U7t27c9VVV/HEE0+EfZykpCTOPfdcVq9ezVlnnRXxdISDE9fX6alQV4tIegvbvwxY/B8wKNpO0WbLli2MHz8+1hrNorsf6O+ou19NTQ0pKSmx1miWYL+uXbvS3Dw5c+bMYc6cOU3WFxQUNFl37rnncu6557Z47GnTpvHVV181u33KlCkcOHCg2fPndrv55JNPuO666xzpfRQKJ66vzr2ergc+bG6jiNwoIhkiklFSUkJ5eTklJSUUFxdTUVFBfn4+LpeLnJwcPB4PmZmZwMGRPn3LPvLz86moqKC4uBhffAUFBVRVVZGbm4vb7SYrK6tRHL6/2dnZ1NbWkpeXR2VlJYWFhZSVlVFWVkaPHj2orKwkLy+P2tpa/7hAwXFkZWXhdrvJzc2lqqqKgoKCdqXJ4/GQk5ODy+UKK02+hrC2pKmwsNDRNI0bN65NaWrPdbKTpgEDBkT9OrWWpvr6ehoaGvzf6+vrqa2tpaGhgbi4ODwej7/KxTfaqO9vdXW1f+rMwDjq6uqoq6vD7XZTU1ODx+PB5XKhlGoSx4EDB/xxeDweampqcLvd/jgCfXxx+Hx8mUKgjy+O5tLk83EiTSLSJE1VVVWsWbOGu+66iz179nD11Vc3SlMon2ilKSEhodU0NXfvhY1SytEPkA5820qYqcAmoE84cU6cOFG1h6G/et//iRYZGRlRizsS6O6nlP6OOvjl5OQ0u62qqspBk7bTEf22b9+uANWvXz+1cOHCGFgdpLXz19K9AWSoMH5jtesOISLHAH8DzlVK7Ym1j110H4Jadz/Q31F3PzNMtj1C+aWnpzdbPeY0h90w4yIyBFgM/EQptSXWPpFA90ltdPcD/R119zMT79jD+DncmC0irwNTgFQRKQIeBBIBlFLPAw8AfYC/eF9wcaswJv7WGd2fNnX3A/0ddffriE/sOmH8HC5RKKWuUkoNUEolKqUGKaX+rpR63ptJoJS6QSnVSyl1nPfToTMJwN8IqSu6+4H+jrr7mYl37GH8NKt6OhQZN25crBVaRHc/0N9Rdz+dJtkJhfGzhxN+h21GkVe635HjbN261ZHjtBfd/UB/R939tB6LCuNnFyf8tOv15BR7DzSeFaqmvoHkxMiP6T5okN7vDOruB/o76u7Xljd3n1m2hT+tyGuy/hfTRnLn9OhM+Xo4vjkeSZzwO2xLFJOHNR4d8sNvS6JyHN+Qx7qiux/o76i7n9vtDjvsndNHUTCv8dwLBfPOj1omAY39Fi5cyMSJE+nWrRu9evXi+OOP56677jroUlCAiPD+++9HzSeYqVOnhpz1riW2bNnCQw89xL59+xqtX7hwISLSaLpVu7Tl+raXwzajEBE6cbBUsWjdzqgcp2vXrlGJN1Lo7gf6O+ruF6uhJcLF5/f4449zww03cM4557B48WJee+01Lr74YpYsWeIPO2DAANasWRP26K6RINQQ462xZcsWHn744SYZxfnnn8+aNWsiOuSGE9f3sK16QikWdXqES+sewUMca7btZceeAwztE9muZofbyJ3RQHdH3f10eTGsOXx+zz33HDfddBOPPfaYf9uFF17Igw8+6F9OSkripJNOctwxUvTt25e+fftGNE4nrq/ejxrRRITNniFMidvgX/X2+qIWdmgfTkwqYgfd/UB/R939Ogr79u2jf//+TdYHPtGHqnpKT0/n7rvvZt68eQwYMIAePXrwy1/+EqUUS5cuZdy4cXTr1o2ZM2dSUVHh36+5aiBffM2Rm5vLlVdeyeDBg0lJSWHcuHH88Y9/9N8Hq1at4sILLwSsSZJEhPT09GaPWV5ezjXXXEOfPn1ISUlhypQpZGRkhHR65plnGDRoEL169eLKK69sUmKJFodviQJY7pnA7PhVrPRMAKyM4o6zRhEf1/aiZnPoPGon6O8H+jvq7tdRqp4mTJjAs88+y5AhQ7jgggvo06dP2HG88cYbTJ48mVdeeYX169dz//334/F4WL16NY8++igul4tbb72VuXPn8vzzz9vyLS4uZvTo0fzoRz+iW7dubNiwgQcffBCXy8XcuXOZMGEC8+fP5+6772bx4sUMGDCgxRkGZ86cydatW5k/fz6pqak8+eSTTJ06la+//poRI0b4wy1atIhjjjmGF198kaKiIu666y5+/etf8+c/R3cGQjjMM4rPPUfzTMJfSOV7yulByfc1fJa3mymj+0XsGHv37qVXr14Riy/S6O4H+jvq6pd+3wcxjye4YTwUbrebhIQEFixYwMyZM5kzZw4iwlFHHcVll13G3XffTffu3VuMIzk5mbfeeov4+HhmzJjBe++9x7PPPkteXp5/6tOsrCxeffXVNmcUwVU706ZNY9q0af5tp512GtXV1bz00kvMnTuX7t27M3r0aACOP/54f2kiFB999BFffPEFq1at4owzzgDgzDPPJD09nSeffJIXXnjBHzYxMZF///vf/hntcnJyeOONN3j66aejPsud3o8aUcZFMmvVGC6J/8y/blFGZBu1jzjiiIjGF2l09wP9HXX30x1f985jjjmGTZs2sWTJEm655RaUUjz66KNMmjSp1V5CU6ZM8U9nCjBixAjS09P9mYRv3e7du6mrqwsVRbMEl8hqamp48MEHGTFiBElJSSQmJvKb3/yG7du3t7kH0tq1a+nbt68/kwBrSI4LLriAzz//vFHYqVOnNsoQxo4d22TO72hxWGcUACs8E5gd/6l/eVlOKXuqaiMW//bt2yMWVzTQ3Q/0d9TdT3dqaw/+vyUlJXHhhRfy3HPPkZOTw9/+9jfy8vL4+9//3mIcwXNVd+rUKeQ6pVSbM4rgNqhf/epXzJ8/nxtvvJGlS5eybt067r//fqDtL7+VlJSQlpbWZH1aWhp79+5ttK659OzfH/2Xhw/rqieAFQ3H8/vElzle8vhajaS+QfHvDbu4/rRhre8cBmPGjIlIPNFCdz/Q31FXP1+1j/JOvtMWAqubwqk+skNycnKz266//nruvfdecnNzo3bc4IwjsMEbmpYo3nrrLW677Tbuvfde/7oPPmhf9dyAAQNClgpKS0vp3bt3iD2a0tL5ixSHfYniO/rwrSedK+I/8a9btG5nxLqcbdiwofVAMUR3P9DfUXe/jjKoXagfzN27d/P999+HfOq2i++N+k2bNvnXffXVV1RWVjYK19DQ0GjZ5XI1apxuaGjgjTfeaBTGV53WWgnjxBNPpKysjNWrV/vXVVdX88EHH4T9rogT1/ewL1GAVf10XfyHPOy+GhfJbC7dzzdF33Ps4J6t79wKEyZMiIBh9NDdD/R31N2vowyTPX78eC6++GLOPvts+vXrx44dO5g/fz4pKSlcc801ET/u5MmTGThwILfffjuPPvooe/fu5YknnmjScB7Y9gEwffp0FixYwIgRI+jduzcLFixoVH0G+BuzX3jhBa688kpSUlJCzqt+zjnncOqpp3LFFVcwb948+vTpw/z583G5XNxzzz1hpeOQG2ZcV5Y3TKCbuDg//uAk629GqFFb90ltdPcD/R119+soE+888MADFBQUcPvtt3P22Wfz29/+lnHjxrF27dpGjdKRolOnTrz77rvExcVx+eWX89RTT/HXv/61SQ+24BLFs88+yw9+8AN+/vOfc91113H00Uczd+7cRmGGDh3K/PnzWbx4Maeeeqr/vYpQvPvuu0yfPp077riDWbNmoZRi5cqVjbrGtoQT11d0f2szHCZNmqSCX1AJB189rOBh+0kfsjZlCrNXWrlzt6QE1v7mLDp3ivxAgQZDpNm0aRNHHXVUxOJzso3CEF1aujdEZH048/44WqIQkZdFpExEvm1m+xgRWSMitSLS/KuREUYRBzP/wgnTZ3FkqpVR7K91R2SgwMzMTNtxRBPd/UB/R939OkqPyJZfAAAgAElEQVSJQleMn/NVTwuBGS1s3wvcDsx3xCYIEWHWpMH+5TcjMFDgcccdZzuOaKK7H+jvqLtfW94cf2bZliYv2KXf9wHPLIveFPa6v9lu/BxuzFZKrRaR9Ba2lwFlIhKzsu5lEwYy/7+bafAovtq+l4LyA6Sntr+xKDc3l7Fjx0bQMLLo7gf6O+ruV1NTE/YsaHdOHxXVIcVD0Ra/WGD8OnBjtojcKCIZIpJRUlJCeXk5JSUlFBcXU1FRQX5+Pi6Xi5ycHDwej796wNfwGFxdkJ+fT0VFBfXfl3FK/4Mv2PxzzTZyc3Nxu93+uZF9cfj+ZmdnU1tbS15eHpWVlRQWFlJWVkZZWRmdOnWisrKSvLw8amtryc7ODhlHVlYWbreb3NxcqqqqKCgoaFeaPB4POTk5uFwuf5qKi4vxnaOCggKqqqr8aXK5XG1OU2FhoaNpGjp0aJvS1J7rZCdNPXr0iPp1ai1N9fX1NDQ0+L/X19dTW1vrb4j1eDz+bpS+qgrf3+rqapRSuFyuRnHU1dVRV1eH2+2mpqYGj8eDy+VCKdUkjgMHDvjj8Hg81NTU4Ha7/XEE+vji8Pn4XmgL9PHF0VyafD5OpEkp1eY0hfKJVpri4uJaTVNz9164ON6Y7S1RvK+UOrqFMA8BVUqpsKqg7DZmg7fBruQb+L9L+bhyKDfVW5OlpHVP4sv7prV7oMD8/HyGDx/ern2dQHc/0N9RB79NmzYxZsyYkC/W1dTUOPJSVnsxfvZoyU8pRW5ubsdqzNaePiOgtooz474mle8BKK2sZfWW3e2OMty3K2OF7n6gv6MOfomJic0+IUZ7wDi7GD97tOTncrlITEy0fQyTUQTSKQVGTCNRGrg0QgMFdpS3YnVGd0cd/Pr160dxcbG/eiIQ3efLMH72COXnqwYrLi6mXz/7o2E7mlWKyOvAFCBVRIqAB4FEAKXU8yLSH8gAugMeEbkDGKuUqmwmyshz1IWQ+z6z41fxYsMFACzfZA0U2Kdr82PKN0dHmQtAZ3R31MHP9zbxrl27msy419DQ0OTtYp0wfvZozi8xMZG0tLRWh2gPB6d7PV3VyvbvgEEO6YRm1DkQl8AIdjFBtpCpRlHfoHj362Ju+MGRbY4uEsW+aKK7H+jvqItf9+7dQ/4olJeXk5qaGgOj8DB+9nDCL/aPQrrRuRcMOx2AK+JX+Ve/2c6BAlsbRz/W6O4H+jsaP3sYP3s44WcyilCMnQnA+fH/I0WsIYjzyqrYsLPt89Pq/CQC+vuB/o7Gzx7Gzx5O+JmMIhRjLgCJp6vUcH7cl/7VizKK2hxVUVHb93ES3f1Af0fjZw/jZw8n/ExGEYoufeBIa2rCwOqn/2TtorqubVMdhjsCZKzQ3Q/0dzR+9jB+9nDC77DMKMIaz8Zb/TRRtnBkojUlYVWtmw+zv2vTsTZu3GhPNsro7gf6Oxo/exg/ezjhd1gPM94i1Xvhn5fDURfy/P7TmPepNfvW5GG9WXTTyZE9lsFgMMQA82a2XVJ6w09Xwml3culp4/1DeKzdvpft5eEP66v7pDa6+4H+jsbPHsbPHk74mRJFmNzwagbLN5UCcMuU4dw7Y0xUj2cwGAzRxpQoIszsSQffA3x7fRHuhvBe6zdPI/bR3dH42cP42cOUKMLEiRJF/XebOfmFfMpd1vl6ec4kzhyTFtVjGgwGQzQxJYpIUZYLz/+AxOcnc1niwXcqwp39zjevga7o7gf6Oxo/exg/ezjhZzKK1uh+BJRb3WZn1bzjX71iUxnlVbWt7j5qlLOzhbUV3f1Af0fjZw/jZw8n/ExG0RrJ3WH0uQCMiNvFpB77AXB7FO9mFre6e2FhYVT17KK7H+jvaPzsYfzs4YSfySjCYfws/9fZ6mP/90UZrQ8UmJamdzuG7n6gv6Pxs4fxs4cTfiajCIcRZ0FyDwDOr/2AFO/g7HllVXzdykCB+/a1fSBBJ9HdD/R3NH72MH72cMLPZBThkJAEYy8GoIvUckGfXf5Ni1pp1NZ5rl3Q3w/0dzR+9jB+9nDCz9GMQkReFpEyEfm2me0iIn8Wka0i8o2ITHDSr0XGz/Z/vaL6Tf/39gwUaDAYDB0Jp0sUC4EZLWw/Fxjp/dwI/NUBp/AYegp0OwKACXXrONKqieJAXQMffFPS7G41NTVO2LUb3f1Af0fjZw/jZw8n/BzNKJRSq4G9LQS5GHhNWfwP6CkiA5yxa4W4eBh/OQAicEXXDf5Nb7UwT0XPnj2jrmYH3f1Af0fjZw/jZw8n/HRroxgIBFb6F3nX6cFxP/R/vbTiZeKtcQJZW7CXbbtDT0dYWlrqhFm70d0P9Hc0fvYwfvZwwk+3jEJCrAvZ/1REbhSRDBHJKCkpoby8nJKSEoqLi6moqCA/Px+Xy0VOTg4ej4fMzEzg4LgomZmZeDwecnJycLlc5OfnU1FRQXFxMb74CgoKqKqqIjc3F7fbTVZJHYw8m9Jhl9L3p4uZMCDJ7/PGVzvIy8ujsrKSwsJCysrKKCsrIz4+nsrKSvLy8qitrfW/Renz8P3NysrC7XaTm5tLVVUVBQUFjqRp//79IX2ys7Opra0NmabCwkJH0zRo0KC2XaesLEfT1LVr1+jfezbS5Ha7tbz3fGmqqKjQ9t7zeDy4XC5t773y8nKSk5PbfZ3CxfGxnkQkHXhfKXV0iG0vAKuUUq97lzcDU5RSzTcC4MxYT6FYnlPKDa9Zx+3XLYkv7zuThPjGeW92djbjx4933C1cdPcD/R2Nnz2Mnz3s+HXUsZ6WAFd7ez+dBHzfWiYRS6aM7kvfblapomx/Las2724SRucbDPT3A/0djZ89jJ89nPBzunvs68AaYLSIFInI9SJys4jc7A2yFNgGbAVeAm5x0q+tJMTHcdmEg8OPL8po+k6FGaLYPro7Gj97GD97mGHGwyRWVU8cKCf/f0uYtsx6hT4hTlgzd5q/lGEwGAw601GrnjoOS26Hp0Yz/LM7OcE71Irbo3j368ZdZc3TiH10dzR+9jB+9jAlijCJSYniv/fDl88CsKjvbdy782QAhvftwvK7zkAkVAcug8Fg0AdToog2E67xfz1/99/p0sk6lfm7D5BZeHCQLl93OV3R3Q/0dzR+9jB+9nDCz2QU7SV1JAw9DYAuVHNh2h7/psCBAseNG+e4WlvQ3Q/0dzR+9jB+9nDCz2QUdpg4x/911oHX/d/f/2YXB2qtgQK3bt3qtFWb0N0P9Hc0fvYwfvZwws9kFHY46kLo3AuACQc+Z0RPq13iQF0DH2Rbr38MGjSo2d11QHc/0N/R+NnD+NnDCT+TUdghMRmOvQqwBgqcnZLp3+SrfiovL4+JWrjo7gf6Oxo/exg/ezjhZzIKuwQ0al9S8XcS4qxSRcaOCvJ3V9G1a9dYmYWF7n6gv6Pxs4fxs4cTfiajsEu/MTD4JAD6qr2c2fd7/6ZFGTupr6+PlVlY6O4H+jsaP3sYP3s44WcyikhwwvX+r1ckrPZ/f2d9MXX1DbEwChuPxxNrhVbR3dH42cP42cMJP5NRRIKxF8O4S+GHizjj5j/TzzuER3lVLZnf1cZYrmVSUlJirdAqujsaP3sYP3s44WcyikiQkASzXoFR55CQmMhlEw/2QnhrffOz3+nA3r0tTTioB7o7Gj97GD97OOFnMoooMHvSYP/3rwoPULZf3zl3jzjiiFgrtIrujsbPHsbPHk74mYwiCgxL7cLk9N4ANCjF4sziGBs1z/bt22Ot0Cq6Oxo/exg/ezjh12pGISJbROSYgGURkZdFZEhQuMkiUhcNyQ5FQz18+w6z697xr1qUsRNdB18cM2ZMrBVaRXdH42cP42cPJ/zCKVGMAJKD9rkGSA0KJ0B8hLw6Lv+cBW9fx3l7XqNrgtUbYdvuA6zfURFjsdBs2LAh1gqtoruj8bOH8bOHE37trXoyY2g3x7iZAKRILRcmrPOvDjX7nQ5MmDAh1gqtoruj8bOH8bOHE36Ot1GIyAwR2SwiW0XkvhDbh4rIChH5RkRWiYjeA60Ec8wVkGIVtmZ5PvCvfv+bEqq8AwXqhO6TsoD+jsbPHsbPHk74OT1ndjywADgXGAtcJSJjg4LNB15TSh0DPAI87qSjbRI7w+SfAnC8bGVkojUOS3VdA0u/KYmlWUgmTpwYa4VW0d3R+NnD+NnDCb9wM4rLROQWEbkF+BmggFm+dd71l4URz2Rgq1Jqm1KqDngDuDgozFhghff7JyG2688JN0BCsjVQoPrIv/pNDaufMjMzWw8UY3R3NH72MH72cMIv3IziHuA57+fPWG0UvwpY9xxwdxjxDAQCfy2LvOsCyeJgpnMJ0E1E+gRHJCI3ikiGiGSUlJRQXl5OSUkJxcXFVFRUkJ+fj8vlIicnB4/H4z+ZvmJaZmYmHo+HnJwcXC4X+fn5VFRUUFxcjC++goICqqqqyM3Nxe12+2eS8sXh+5udnU1tbS15eXlUVlZSuKea6tGXWAmI/5wErEbt9TsqWPp5Zsg4srKycLvd5ObmUlVVRUFBgSNpiouLCy9NhYWUlZVRVlZGYWEhlZWV5OXlUVtbS3Z2dlTTdMwxx0TnOkUoTYMHD9bn3guRpq5duzpyndqbJt8QFDreex6Ph6SkJG3vvfLycvr379/u6xQujs6ZLSKzgHOUUjd4l38CTFZK3RYQ5gisjGcYsBor0xinlPo+RJRAjObMbo3yrfDcJEBxc90dfOSZDMBNpx/J3POOiq1bADk5OYwdG1z7pxe6Oxo/exg/e9jx03XO7CJgcMDyIGBXYACl1C6l1KVKqeOB33jXNZtJaEvqCBh9HgCz41f5V7+TWUR9gz6DjA0bNizWCq2iu6Pxs4fxs4cTfu3OKEQkRURuE5EFIvKAiAwNY7d1wEgRGSYinYArgSVB8aaKiM9rLvByex1jzqm/AOD0uG9IwxqPpbyqjk9yy2Jp1Yhdu3a1HijG6O5o/Oxh/OzhhF84b2Y/JSJbgtZ1AzKBPwJXAL8FskRkVEtxKaXcwK3Ax8AmYJFSaqOIPCIiF3mDTQE2e4+ZBvy+bUnSiCEnUj/oZBISO3HZ4Cr/ap3eqejdu3esFVpFd0fjZw/jZw8n/BLCCDMV+L+gdXcDo4AblFIvi0hfYBlWhvGTliJTSi0FlgateyDg+9vA22F4dQj2nnw/aUNHMbumC3+ZvwqATzbvpqyyhn9+VcifVuQ12ecX00Zy5/QW89yIUV1dTa9evRw5VnvR3dH42cP42cMJv3CqntKB4Dc6LgNylFIvAyildgNPAadG1O4QwNN7OHTtR3pqFyYP8w4U6FG8k1nMndNHUTDv/EbhC+ad71gmAfh7PemM7o7Gzx7Gzx5O+IVzhATAP062iPQGjgJWBoUrAPpHzOwQITEx0f/9ioDhx9/SZKDAQD9d0d3R+NnD+NnDCb9wMootWO0GPi7w/v04KFw/QO8ZPmJAVdXBtolzx/enayfrlG8rP0CGBgMFBvrpiu6Oxs8exs8eTviF00bxHPCSiPQASoHbge3Af4PCnQ18G1m9jk9qqneQ3bpqUlY8zIWeBl5nKgBvrtvJCemxbSjz+2mM7o7Gzx7Gzx5O+LVaolBKLQQeAC7F6q66GbhEKVXvC+NtzL4YeC86mh2XoiLvVKgJyVDwBVfEHayx+0CDgQL9fhqju6Pxs4fxs4cTfmG1giilHldKDVJKdVVKna6Uyg7avlsp1V8p9dfoaHZcRowYYX2Ji4Mz7+dYyWeUWN1jXfUNvJ8V2z7afj+N0d3R+NnD+NnDCb9w3qN4oA2f30bduIOxcePGgwujzkEGn9DoTe1Yv1PRyE9TdHc0fvYwfvZwwq/VsZ5ExAO4gAO0PmGRUkr1i5Bb2Gg51lNzbPuUPa/+iJNqF1AfookouLuswWAwRItIjvW0DUjEepfibmC4UqpvMx/HMwndaTKpyJFn0GfYcZwVp8dkKLpPygL6Oxo/exg/ezjhF9bosSIyCWtcptlYc2V/BLwOvK+UCn+s2ijRoUoUADvX8smL93Bt/b0AJODG7S1dmBKFwWBwioiOHquUylBK3a2UGgLMAL7D6jZbJiL/FJHT7ekeuoTM7QdP5vSj0+nPHgBvJhGbl+90f1oC/R2Nnz2Mnz20KVGE3NEa/fX3wJ3AEqXUpZEUawsdrkQBsHcb85/+A8+5rbEQ+1LBbnqZEoXBYHCMqM1HISKnisizwA6saVHfBv7UdsXDA99sVU3ofSSzJh6c3G83PR0yakyzfhqhu6Pxs4fxs4cTfmFlFCIyQUSeEJEdWPNZD8YqSfRTSl2plPo0mpIdmVGjmh/gb+iM2zipn++9xdY6lEWHlvx0QXdH42cP42cPJ/zCeY9iM/A/4BjgQazMYaZS6g2lVHW0BTs6hYWFzW/s3IvZUxqX+j7PK4+yUWNa9NME3R2Nnz2Mnz2c8AunRDEScAMTgSeArSJS1twnqrYdkLS0tBa3X3DMEY2Wf/Z/69lSuj+aSo1ozU8HdHc0fvYwfvZwwi+cQQEfjuQBRWQGVptGPPA3pdS8oO1DgFeBnt4w93knO+qQ7Nu3j+7duze7vVNC47x6f62ba1/6nHd/MZV+3ZKjrdeqnw7o7mj87GH87OGEX6sZhVIqYhmFiMQDC4DpQBGwTkSWKKVyAoLdjzVF6l9FZCzWbHjpkXJwmuTk8H7sB1LGHnpQQxLFVR5uWLiON286hc6d4rXwiyW6Oxo/exg/ezjh5/TUTZOBrUqpbUqpOuANrFFnA1GAL3vsAeg9s3mE8BDPHxMWEIcHgG+KK/nFG1/T4In95EYGg+HwxumMYiAQOApekXddIA8BPxaRIqzSxG2hIhKRG0UkQ0QySkpKKC8vp6SkhOLiYioqKsjPz8flcpGTk4PH4yEzMxM4+HJKZmYmHo+HnJwcXC4X+fn5VFRUUFxcjC++goICqqqqyM3Nxe12k5WV1SgO39/s7Gxqa2vJy8ujsrKSwsJCysrKKCsro7i4mMrKSvLy8qitrfV3ZQuOo4Q+ZKkRPJyw0J/G/+aU8pu3MqKapry8vDanqbCwMKw0ZWVl4Xa7yc3NpaqqioKCgnZdJ9/2aF4np9PkxL3nS9OOHTu0TlNubq7W12nbtm1a33ulpaXtvk7h0u4X7tqDiMwCzlFK3eBd/gkwWSl1W0CYu7xeT4nIycDfgaOVUp7m4tX5hbvKyspW6w/T7/sAgETc5A36Pb8vmchLDRf4tz9y8TiuPjk9Zn6xRndH42cP42cPO35Re+HOJkVY72D4GETTqqXrgUUASqk1QDLW+FIdktLS0rDD1pMAFz/H3IQ3mBG31r/+oSUbWZkbfjxtoS1+sUJ3R+NnD+NnDyf8nM4o1gEjRWSYdwiQK4ElQWEKgWkAInIUVkax21HLCDJkyJC27TB4MnEn3cgziX/hWNkKgEfBrf/6mm+Lv4+9XwzQ3dH42cP42cMJP0czCqWUG7gV+BjYhNW7aaOIPCIiF3mD/RL4qYhkYY1QO0c5WT8WYbZs2dL2nc78LZ179udvneYzyPtqSnVdA9e/uo6S7yM7WG+7/BxGd0fjZw/jZw8n/Bxto4gWOrdRhIOvjQIChhnfuhz+7zLyPAO5tO4h9tMFgDH9u/H2z06ha1I4r8AYDAZD8+jaRnHY0e4hgEecBcf/hJFxxbyQ+AwJNACQ+91+fv7PTNwNzbbtO+PnILo7Gj97GD97aD3MuE4ckiUKgNr98PwPYMRZvJX6M+55N9e/6UcnDuF3M49GJDaDCRoMho6PKVFogq3cPqkb3PwZnD+fWScO57YzR/g3/fOrQl76bFts/RxCd0fjZw/jZw9TogiTQ7ZEEYRSijve3MB7Gw72KP7rjyZw7vgBUfUzGAyHJqZEoQm+NzUjgYjwxGXjOSG13r/ujjc38HVhRbvjjKRftNDd0fjZw/jZwwk/U6KIMm63m4SE0D2Unlm2hT+tyGuy/hfTRnLn9BCTkdRUwn9+QcW3/+VS5rO91nobs0+XTvz756cyuHdKRP10QXdH42cP42cPO36mRKEJW7dubXbbndNHUTDv/CafkJkEwI4vYeNiekkVr6gH6dXJ6gm150Adc15Zy/fV9aH3a6efLujuaPzsYfzs4YSfySiizKBBgyIX2egZMPkmANLjSnlJHsM3Cnn+7gPc/H/rqXO3rdtsRP2ihO6Oxs8exs8eTviZjCLKlJdHeGrTsx+FIyYAMEk2MT/5Ff+mNdv2MHdxNm2pToy4XxTQ3dH42cP42cMJP5NRRJmuXbtGNsKEJLji/6BLPwAualjGPT1X+Te/k1nEsyvDL4pG3C8K6O5o/Oxh/OzhhJ/JKKJMfX3b2w1apcdAuOIfEJcIwC2uF7miz8F3Kp5etoV/f10cO78Io7uj8bOH8bOHE34mo4gyHk9khtpowpCT4Pz5AIjA76oe5LRe+/yb7337G9Zu3xs7vwiiu6Pxs4fxs4cTfiajiDIpKW3vsho2E+f4G7cTpYG/VP+Skd2sp4u6Bg83/iODbburYucXIXR3NH72MH72cMLPZBRRZu/e1p/qbTHjcRhtvc3dvWdfXrlyBKldkwDYV13PtQvXsaeqNnZ+EUB3R+NnD+NnDyf8TEYRZY444ojoHiAuHi77G0y4Bm5YxqDhR/P3ayaRnGhd2h17qrnxH+upqW+IjV8E0N3R+NnD+NnDCT+TUUSZ7du3R/8gnVLgoj9Dt/4AHDu4J3+68nh8A8uu31HB3W9l4fE07TbriJ9NdHc0fvYwfvZwws/xITxEZAbwJyAe+JtSal7Q9meAqd7FFKCfUqpnS3HqPISHx+MhLi42+fHf3lvB79bU+JdvmTKce2eMaRQmln7horuj8bOH8bOHHT8th/AQkXhgAXAuMBa4SkTGBoZRSt2plDpOKXUc8Cyw2EnHSLNhw4bYHDhvOdd/80OuTlrtX/WXVfm8ua6wUbCY+bUB3R2Nnz2Mnz2c8HO0RCEiJwMPKaXO8S7PBVBKPd5M+C+BB5VSy1qKV+cSRUyoOwB/OhYO7Mat4rjRcx8r648GICFOWHjtZE4bmeoP3ubBCQ0GwyGBliUKYCCwM2C5yLuuCSIyFBgGrGxm+40ikiEiGSUlJZSXl1NSUkJxcTEVFRXk5+fjcrnIycnB4/GQmZkJHJzkIzMzE4/HQ05ODi6Xi/z8fCoqKiguLsYXX0FBAVVVVeTm5uJ2u/3D+fri8P3Nzs6mtraWvLw8KisrKSwspKysjLKyMlavXk1lZSV5eXnU1taSnZ0dMo6srCzcbje5ublUVVVRUFDQ/jRt3UHNJQtpSOxKgnh4Nu4pxsVZJQm3R3Hz/2Xw0RrreMuXL+fO6aN4Z1b/Ruf3Pz8awi2nDw2ZpsLCQkfTlJGREfXrZCdNX375pZb3ni9Nq1atcu7ea0eali1b5sh1am+aVq5cqe29V15ezueff97u6xQuTpcoZgHnKKVu8C7/BJislLotRNhfAYNCbQvGlCiaoXg9/OMSqPme71QvZtb9ju9ULwAG9uzMu7ecQr/uyf7g4U6gZDAYDg10LVEUAYMDlgcBu5oJeyXwetSNoowvR48JAyfCNe9D5970lwpeTvwDXbCeIor3ubjhtQy+XKt/BhvTcxgGxs8exs8eTvg5nVGsA0aKyDAR6YSVGSwJDiQio4FewBqH/SLOcccdF1uBAcfAnA+g2wDGxhWyIPFPxGO9U/FN0fe8sknREKLbrE7E/By2gvGzh/GzhxN+jmYUSik3cCvwMbAJWKSU2igij4jIRQFBrwLeUIfA9Hu5ubmxVoC0sXD9f6HPSKbEf8MjCQeHJl+2qYzff7AphnKto8U5bAHjZw/jZw8n/Byf308ptRRYGrTugaDlh5x0iibDhg2LtYJFzyFWZvGv2fyoaCU7VH9ebLgAgJe/2M7QPvqOZ6PNOWwG42cP42cPJ/z0fYvkEGHXruaaYGJASm+4egmMvZj7Tk/l3KMP9nR6+D8bYyjWMlqdwxAYP3sYP3s44afvjOGHCL179461QmM6pcDlC4lD8UwDFO35nOySKqxmCgVIjAWbot05DML42cP42cMJP1OiiDLV1dWxVmhKXBzExZOcGM/j5w1lUK/O3g1CV6qJpyHkuFCxQstzGIDxs4fxs4cTfiajiDI6jxED0KdrEgvnnED3+DoAqkihgXjOeWo5b2XspM4d+0lbdD+Hxs8exs8eTvjpfQYOARITE2Ot0CKJiYmM6J3A8+mrSeHgAIJ5e+q45+1vOP2JT3hxdT77a2I3HWRHOIc6Y/zsYfxMRhF1qqpanmEu1lRVVUFiZ0654WkuiFvDdfFL6crBoux3lTU8tjSXUx5fyeMfbqK0sqaF2KLoqDHGzx7Gzx5O+JmMIsqkpqa2HiiG+P1EWOSZypeeo3k1cR6/SnidvlT4w+2vdfPCp9s47Q8rufftLLaW7XfeUVOMnz2Mnz2c8DMZRZQpKiqKtUKLBPvlqiH8sP5+fnbqID5PuoM/JLzIkXKw+119g2JRRhFnPb2aG17NIKMg+tMwdrRzqBvGzx7GLwYTF0UDnQcFdLvdJCTo2ws50K/JoIA71sB7P8ezZxvLPRN4wX0B69XoJnFMGNKTm84YzvSj0oiLi3z32o50DnXE+NnjUPbTdVDAw46NG/V9kQ1a8Rt6Mtz8OXGn3sbZCRt4Z+znvH3TSUwfm9YoWGbhPm76x3rOeuZT3lhb2Oz83FFx1ADjZw/jZw8n/EyJwuCnxWHGSzdCpy7QKx2ArWVVvLR6G+9m7qQuqAdtatckrj01nR+fOJQeKeH3yDATKBkMzmJKFJrgmzBEV8L2SxvnzyQARvTryh9mjubz1Mf5Wfx7dIur9W8rr6rlyY83c8q8Ffzu/Rx27QtvgpQ7p49qkkEVzDuf03s713DeHg6ZaxwjjJ89nPAzJQqDnzZPXPTFn2CZNZ7jftWZNxqm8ne5lO/qGw8wmBAnXHTcEdx4+pGM6d898h4Gg1cJoR8AABnRSURBVKFdmBKFJhzSTyPjLoXjfgwSRzdx8dOEpayOu4n5iX9lVNLBrrVuj2JxZjEz/vgZc15Zy5r8PbTlAeWQPocOYPzsYfxMicJABNoGdm+GlY/Cpv/4V3mUsMpzLC8k/JCvXIOa7HLsoB7cdMZwzhnXn/ignlKmRGEwOEO4JQqTUUSZ7Oxsxo8fH2uNZomo33ffwuonIec9rJFoLb72DOfFfvfzUXESwbfb0D4p3PCDI5k1cRDJifFA04zisDqHUcD42eNQ9jMZhSbU1taSlJQUa41miYpfWa7VfpH9FnjqoVNXuGsT26vieemzbby9vqjJYIN9unTimlPS+clJQzn+0WX+9QXzzj88z2EEMX72OJT9tG2jEJEZIrJZRLaKyH3NhJktIjkislFE/uW0YyQpLCyMtUKLRMWv3xi45K/wiyw45XY48SZI7s6w1C48dsl4vvjVmdw6IZnucQfHjdpzoI6nl23hlHkrnXGMIMbPHsbPHk74Ofq6oYjEAwuA6UARsE5EliilcgLCjATmAqcqpSpEpJ+TjpEmLS2t9UAxJKp+PQbC2Y82Wd23WxJ3J77FzxIX82bDVP7ecB7FyhqvxhX0st6TH+dy8pCuDPeoqLz1HQkO62scAYyfPZzwc7pEMRnYqpTappSqA94ALg4K81NggVKqAkApVeawY0TZt29frBVaJCZ+7lrY/CFdpJbrEj5iVac7+WPiAsbIjiZBF3ySz49fzWLyYyu49+0sPt74HdV1buedW8BcY3sYP3s44ed0RjEQ2BmwXORdF8goYJSIfCEi/xORGaEiEpEbRSRDRDJKSkooLy+npKSE4uJiKioqyM/Px+VykZOTg8fjITMzEzjYlSwzMxOPx0NOTg4ul4v8/HwqKiooLi7GF19BQQFVVVXk5ubidrvJyspqFIfvb3Z2NrW1teTl5VFZWUlhYSFlZWWUlZVRWVlJZWUleXl51NbWkp2dHTKOrKws3G43ubm5VFVVUVBQ4EiafPPttiVNhYWF9tK0ey+7Zn9I9ZRHqOkzjkRpYGb8F3zYaS6vJs5jatzXdKLx/BflVbUsyijipn+s57iH/8tVf/2MBR9/w8btxRG5ToFpmv9RDun3fdDk88yyLSHTVFdXp+W950vTnj17tLz3fGnauXOnc/deO9JUWlrqyHVqb5qqq6vbfZ3CxdHGbBGZBZyjlLrBu/wTYLJS6raAMO8D9cBsYBDwGXC0UqrZbFPnxuyysjL69dO39kwLv92b4Zs34dt3oKIAgAMqic8841numcSKxClU1DQ/097RA7szbUwa08emMe6I7ohEpooq3G66WpzDFjB+9jiU/cJtzHZ6SMQiYHDA8iBgV4gw/1NK1QPbRWQzMBJY54xiZKmpcX6in7aghV/f0TDtATjzt7Ark5f+Op9z49cyIz6DGcM7s33KLeyV7izfVMqKTaVsKW08Ucu3xZV8W1zJn1bk0b97Mmce1Y+zjurHKcNT/V1uo4kW57AFjJ89jJ/zGcU6YKSIDAOKgSuBHwaF+TdwFbBQRFKxqqK2OWoZQXr27BlrhRbRyk8EBk7k9+4f83v3jxgnBXwwbTJ9evRiWPfuTBzai1/NGEPhksdYsXYDKzwT+J/nKNwBt/F3lTX866tC/vVVIZ0T4zltZCpnHdWPM8ek0bdbdLo4anUOQ2D87GH8HM4olFJuEbkV+BiIB15WSm0UkUeADKXUEu+2s0UkB2gA7lFK7XHSM5KUlpbSvXvr4xvFCn39hI1qGAyeTGleXiPHIaUruDYhg2v5mErVmdWeY1jRMIFPPMexj27+cK76BpbllLIspxSRbI4d1JOzjurHtKPSGNO/W8SqqPQ9hxbGzx7Gz7xwF3UO5Zd1okVw20ATR3ct7PgStq2C7Z/Crg2Awq3iWK9GsaJhAss9E9imjmj2GAN7JHHW2P5MOyqNE4/sTVJC0yqqcNsoonUOIzXsuo7XOBDjZw8nXrjTd9qmQ4QtW7Zo/fq/7n4QwjEhCYZPtT4A1XuhcA0JBV9wYsFnnPjd6/yaf7HN058VJ73K8p2KjB0VNHgOPhQVf1/Lq2t28OqaHXSNb+D0wfFMGzeIqcePoXfXtv3TResc3jl9FHdOH2V77Cvdr3E4frGcq+RQOH92MSUKg3bYHhSwphKK1kHxevjB3RAXx77qOlZt3s3yrG18mlvCfrqE3DUODxM7f8e0flXk7Cxjqeck3CTEdHDCWA6SqNtkUmbAyMhiShSasH79eiZOnBhrjWbR3Q/a4ZjcHUZMsz5eeqZ0YubxA5mZtpv62odYV1LPsvpjWOGZQKE6+GarhzjWuY5g3Q6w+lFYD1JXvLCG/j2S6d+pjrSqjfTvm0pa2hGkDTqSoh3bOOmEVv/XYoadaxypUk1L6H4PGj9TojBohKNPr+46KN+M2vUNW7dtZfmOepbv7UdmwzBUO95DTY2vJi2pnrQuQlr3ZPr37kH/fmmk9etrLXdPpmdKYrsa0HV4itbBQSePWBOp/xVTotAE8zQSPr6n12CiMjFLQifoPx7pP56RE6wXdX6mFHu+28EnX29mRd4+lpckUx/mv0h5Qwrl1bCxGv6/vXMPjqs6D/jv0z4k7coPSX5ibMs0xA4G8wwlhJLwmOHVQjLQJDRNQiYJgZmEPtIUSKZhhqZNZviDNoU2DQwzDSWQliQ0CaTBxSEmgJMYg7CxjW3ZsvySZElrS/vwalf79Y97ZV+tpdVKq717Md9v5s7ee87Z3d89e3Q/nfs4h8NARwbodBeH+nAdC2OwKJplwax6FjU3sWheCwtb57JoTiMLZzWwYHb9jD/7EaTfeDzMb+r40dPzYoGiygStgRUTdD/w0VGE1sVt3LK4jVsY+9/r45+7mO6jx+jp3Eb3wb10D+XpPVZHdz7OYZ1TVi8kmy/QNQhdRKCvAHsGgcGTyjWHcyyMFVjUPJuF88c+cbt+x2Fa4lHmxiK0xKM0RkKT9lKC/hubX/CxQFFl2tvbOffcc2utMSFB94NgOP7RmfOdlYuWjs0YybHlledpnV1Pd/cBevr66Ukk6R7K0TNnDd2RpfQMHqNnMEsyW95ghol8hMQgbB/Mwt59Y/I+/djvxmxHpUBzdITmBmhuCNEcj9Lc1EjzrDhz58yhpamBRM9+Lli90g0wUWY3hGfsGZKZIAi/bymC7ucHFiiqzOrVq2utUJKg+0HAHUMRVn3gGsLhMIsnOZYks3m6f/M4Pfs66D6aoSc1Qk8mRPdIE93aQo82c5i5jFD+qadhraMnW0dPFjgKkAeG3KX7RMFfv3JCmQLNkTxz65WWBmFuQ4jmWITmpgY3yDTRPGe2E2zi9bTEo+XXR5moKqrOrQIrV72P/EgBBTdNj8+EWLxdC7ztL2h3gfmFBYoqs2vXLlatWlVrjQkJuh8E37Fcv6b6MO+56rO8x5uoCtkhGDoEQ4cYOdpDf38v3QODdC+5mp7QIv7umS3Hi18a2cFAvp4j2sQAsxhm6gfxEeroy0XpywHHh80qAGl3mWggBKWRYc7+2jMogrq9EkXcxVnHDQCjac5u6vFAMFPc+r0NtM2Lsbw1zvIW97U1Rrx+Zg9r3t/X72sDQcECRZU5/fTTa61QkqD7QfAdK/ITcW7nbZgN81cSAha4yxq3iDdQ/OCrH4fUYUgdRpN9ZAYPkhg8SuJokkQ6SyKVI3FshERWSKz8OAltom8ow2B2hERqmMSRI6RpmK4sGeqdmFIW1e0GvLq7n1d3nxzU5jVFaWuNs7w1TltrjGWtMdpa47S1xpkTi0z5e4Le/vzAAkWV6evro6mpqdYaExJ0Pwi+o69+s09zFkCAmLsUT+ripbOzk7a2NhjJwbafkk12c2QoycBQmkQqy5F0joHMCEeyMDAc4kguTGKknoFFl3HkmBNgBo9VZ7IocaOOwPF+ibjbhKPOuoAgiI4g+Qxp6kvePNCXHKYvOczGvYmT8ubUZWiLHGV5dJC2hiTL6jO0xY6xPDbMvAZFwlEI1cOic+D8TzqfN/r7HngNet6CUJQ/rttCjpBzV9yuBghFIRRxlroIxOfD7MVjvzyXAS04+aGIs2PvECxQVJkgH+Ag+H4QfMd3jF8oAmffTD2w0F3Kpe2eZ2nkGO+X7fzLzWciw0NINoXkkshwEhl217NJyKUBkM884xzgxQ0Ee9Yjj3/ECQaTHSNbzoC7Xh+btvlp+NHnyGqY/TqfLl1Ipy5kr2fZpwtK3tJ8tNBIe7aR9uwi5zKOhzgZlkkPbdLD8nkHWZ7rYnlrjNkSoVBQ6rb9DH7zIAAPec/4/ec4X3Tx7XD9A2PTfv5X6BtPkqGeFPWkiJOSOKm6JlISIy1xUsRISSOpZVeQmn8e6WyeZHaE9HCe5IFtpDPHSGmUVCFCuhAhTGTM6MnVwgJFlcnlcpMXqiFB94PgO75b/DI0sF7PY85F0zwnv+IyuHs3DKec/65zacil6e/ZT2tTg7udcZboOEOszFrE2pELaSBLveS4YqlArgvyO5yBIvMZ8rkch4Zj7F10FZ0X3Mve/hR7+9POcvgoxwoT90RSNLJN29imbdAL/Hjz8bxo+E2WRVfTNvwVlksPy6WHuBwjrfWkaCCljaSoJ00DSW0kvWUlyX2vkh7Ok86OkMzmSaeuJVW4obwHOt8C6ChKbC6jkquDBYoqUyiUfUK3JgTdD/x3HO/OlrZ7np3wzpag12Fg/EJhiLU4i4dj4QOwpNTJM5e2y/hC7ivHNztvPzlghYGlqiwtjHBZaOzhTYd66T20j85Ehr0DWfYeGabzaIG9gwX2DglDuYm7OcP5ArvyMXZR5jMVA8DAQFHi1K+PBAULFFUmFovVWqEkQfcD/x0nekJ8IoJeh+86PxEnKBUnz1rAwlkLWAj8YVGeqpJI5+jsT9HVn6bzeE8kxZ6+JIn0zF2jaYjU0VQfJhYNEYvU0RQNEYsKTZE6YhGIRyAeixOPx4lFQ8Trw8SjYeLpfcQLKWLhEZpCSiyU54tPbWMTK2fMbSIsUFSZgYEBmptr12WcjKD7QfAdq+U31Z7NRLxb628qiAgt8Sgt8SgXLBvr0tHRwYIly46fwursT/HAL98+nv/Fy88gFg0Tr3cO6rFoyA0EJ9Li7nosGiZUN92L2ItPStn0lD+9Rd8DhYhcC/wzzgx3j6rqt4vybwMewJkqFeAhVX3UV8kZ5LTTJp48JwgE3Q+C71gtv6n2bCbi3Vp/M8Vpp51GY0OEs5fM4ewlcwDGBIp7r39frdR8Y+rDZFaAiISAh4HrgLOAW0XkrHGK/lBVz3OXd2yQANizZ0+tFUoSdD8IvqP5VYb5BR9fAwVwMbBLVXer6jDwFHCTzw6+EuQniiH4fhB8x1PZ78G1O8Y8hQzO6a8H1+6oVOs4p3L9nSr4HSiWAN5RzvYz/rNCN4vImyLytIgsHScfEbldRDaKyMZDhw7R19fHoUOHOHDgAIlEgo6ODjKZDFu3bqVQKLBp0ybgxJDVmzZtolAosHXrVjKZDB0dHSQSCQ4cOMDo53V2dpJMJtm+fTv5fJ729vYxnzH6unnzZrLZLDt37mRwcJCuri56e3vp7e3l5ZdfZnBwkJ07d5LNZtm8efO4n9He3k4+n2f79u0kk0k6Ozt92ad169ZNeZ+6urp83afXX3+96r9TJfu0YcOGQLa90X1av379tH+nj50V5/W7L+WFL6xi231X8tynV7D7H6/jQ63JGdunF154YdJ9+uYzr48bsO594qWq/z29+OKLJ+2Tl1q2PS/TaXvl4uvERSLyp8A1qvp5d/tTwMWq+mVPmVYgqapZEbkD+JiqXlnqc23iIsMw/CQoYz1V6lHuxEV+9yj2A94ewunAQW8BVe1X1ay7+QiUe+NyMKnKpDszSND9IPiO5lcZ5hd8/O5RhIEdwFU4dzX9HvgzVX3LU2axqh5y1z8K3K2ql5T6XOtRGIbhB0EZZtzvqVDdceH9W4DrcYJFB/B1N+1+4EZ3/Vs4D7C3A78CVk32mRdeeKEGlddee63WCiUJup9q8B3NrzLMrzIq8QM2ahnHbV97FNUiyD2KQqFAXZ3fZ/jKJ+h+EHxH86sM86uMSvyCeo3iXcf27dtrrVCSoPtB8B3NrzLMrzL88LNAUWVWrFhRa4WSBN0Pgu9ofpVhfpXhh58Fiipz8ODByQvVkKD7QfAdza8yzK8y/PCzQFFlWlpaJi9UQ4LuB8F3NL/KML/K8MPPAkWVSafTtVYoSdD9IPiO5lcZ5lcZfvhZoKgyQb5bAoLvB8F3NL/KML/K8MMv2DVwChCJBHtWq6D7QfAdza8yzK8y/PA7JZ6jEJHDwN5ae0zAPKCv1hIlCLofBN/R/CrD/CqjEr/lqjp/skKnRKAIMiKysZwHWmpF0P0g+I7mVxnmVxl++NmpJ8MwDKMkFigMwzCMkligqD7fq7XAJATdD4LvaH6VYX6VUXU/u0ZhGIZhlMR6FIZhGEZJLFAYhmEYJbFAUQEicq2IvC0iu0TknnHy/1pEtorImyLygogs9+SNiMgb7vLTGvndJiKHPR6f9+R9RkR2ustnauT3oMdth4gc8eT5UX+PiUiviGyZIF9E5Duu/5sicoEnz4/6m8zvk67XmyLyioic68nrFJHNbv1VZTKXMvw+LCJHPb/jNzx5JduGT35f9bhtcdtci5tX1foTkaUi8isR2SYib4nIX4xTxr/2V87sRraMO1NfCGeWvjOAKM6MfGcVlbkCiLnrdwI/9OQlA+B3G/DQOO9tAXa7r83uerPffkXlvww85lf9ud9xOXABsGWC/OuBXwACXAL81q/6K9Pv0tHvBa4b9XO3O4F5Na6/DwM/r7RtVMuvqOyfAOv8qj9gMXCBuz4LZ1bQ4r9f39qf9Simz8XALlXdrarDwFPATd4CqvorVR0dsWsDcHqQ/EpwDbBWVQdUNQGsBa6tsd+twJMz7FASVV0PDJQochPwfXXYAMwVkcX4U3+T+qnqK+73g//tr5z6m4hK2m7ZTNHP1/anqodUdZO7PgRsA5YUFfOt/VmgmD5LgH2e7f2c/EN6+RxO9B+lQUQ2isgGEflIDf1udrutT4vI0im+1w8/3FN2K4B1nuRq1185TLQPftTfVClufwo8LyKvicjtNXIC+ICItIvIL0RktZsWqPoTkRjOgfZHnmTf6k9E2oDzgd8WZfnW/sKVvPldjoyTNu69xiLy58BFwIc8yctU9aCInAGsE5HNqtrhs9/PgCdVNSsidwD/AVxZ5nv98BvlE8DTqjriSat2/ZXDRPvgR/2VjYhcgRMoLvMkf9CtvwXAWhHZ7v6H7SebcMYaSorI9cAzwJkErP5wTju9rKre3ocv9SciTTgB6i9VdbA4e5y3VKX9WY9i+uwHlnq2TwdOmmpKRK4Gvg7cqKrZ0XRVPei+7gZexPmPwVc/Ve33OD0CXFjue/3w8/AJirr9PtRfOUy0D37UX1mIyBrgUeAmVe0fTffUXy/wE5zTPb6iqoOqmnTXnwMiIjKPANWfS6n2V7X6E5EITpB4QlV/PE4R/9pftS7GnOoLTm9sN84pkdELbquLypyPc1HuzKL0ZqDeXZ8H7GSGL9aV6bfYs/5RYIOeuBi2x/Vsdtdb/PZzy63EuXAoftaf57vamPhi7A2MvZj4O7/qr0y/ZcAu4NKi9Dgwy7P+CnBtDfwWjf6uOAfaLrcuy2ob1fZz8+fgXMeI+1l/bj18H/inEmV8a3926mmaqGpeRL4E/BLnLo3HVPUtEbkf2KiqPwUeAJqA/xYRgC5VvRF4H/DvIlLA6dV9W1W31sDvLhG5Ecjj/DHc5r53QET+Hvi9+3H369hut19+4FxEfErdvwCXqtcfgIg8iXNnzjwR2Q/cB0Rc/+8Cz+HcebILSAOfdfOqXn9l+n0DaAX+1W1/eXVGGV0I/MRNCwM/UNX/rYHfLcCdIpIHMsAn3N953LZRAz9w/oF6XlVTnrf6UX8fBD4FbBaRN9y0r+EEf9/bnw3hYRiGYZTErlEYhmEYJbFAYRiGYZTEAoVhGIZREgsUhmEYRkksUBiGYRglsUBhGIZhlMQChWEYhlESCxSGMcOISIeI3OfZ/piIdIv7hJZhvNOwQGEYM4iIzMYZeqLdk3wD8Jza063GOxQLFIYxs5yHM/ZOO4CI1OEMUf1sLaUMoxIsUBjGNBGR2SLyXREZEJE+EflbnIEgB1V1j1vs/TgDs61131MnIikRuUNEviki+93pQB9xg4phBA4bFNAwpoGIRHEGrZuFM81tCvgW0MjJp51e0hNzCZwBxIC/Af4HZyC3S4D7cQZ5+4kf/oYxFSxQGMb0uAd4L/Beded5EJEjwEuMnUnuBuAJz/Ya9/VhVX3QXV8rInfiTNpjGIHDurqGMUVEJATcBfybeiYDwpl7BOANt9xinFNR3usT5wCDwMOezxNgLtBXRW3DmDYWKAxj6pyDM8/D80Xpbe7r6Kmn64Hdqvp20XtfUtVhT9of4JyymvE5FwxjJrBAYRhTZ7H72lWUfg0wwokD/g2cfLfTGsZewwA4FygAW2bQ0TBmDAsUhjF1Rk83Hb+mICKtwJeAt1U1417svhpPoBCRRpzewxuMZQ1OzyOFYQQQCxSGMXXacSaw/46I3CgiNwPrcOZPHu0tXI7z9/Vrz/tWu2nFPYo1wJtVNTaMCrBAYRhTRFWzOPM9Z4H/Av4BeAinpzEaBG4A/s8tO8o5OHMb7yr6SAsURqCxObMNowqIyA7gAVV9pNYuhlEpFigMwzCMktipJ8MwDKMkFigMwzCMkligMAzDMEpigcIwDMMoiQUKwzAMoyQWKAzDMIySWKAwDMMwSvL/KkUTozyJ2s4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.figure(1, figsize=(12, 8))\n",
    "xx = np.mean(dual_simu_gaus, axis=0)\n",
    "yerr = np.std(dual_simu_gaus, axis=0)\n",
    "plt.errorbar(zeta_seq, xx, yerr, capsize=3, lw=3, label='Simulation')\n",
    "plt.plot(zeta_seq_2, theo_gaus, lw=3, ls='--', label='Theory')\n",
    "# plt.plot(zeta_seq_2, np.ones(100) * MSE_original(lbd, gamma, alpha, sigma), lw=3, ls='-.', label='No sketching')\n",
    "plt.grid(linestyle='dotted')\n",
    "plt.legend(fontsize=15)\n",
    "plt.xlabel(r'$d/n$', fontsize=15)\n",
    "plt.ylabel(r'MSE', fontsize=15)\n",
    "plt.title('Gaussian dual sketch (zero noise)', fontsize=15)\n",
    "plt.savefig('./Plots/dual_gaus_gamma_{}_alpha_{}_lbd_{}.png'.format(gamma, alpha, lbd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged\nconverged\n0.7403707632794974\nconverged\nconverged\nconverged\n0.6342223634943365\nconverged\nconverged\nconverged\n0.5659191832132635\nconverged\nconverged\nconverged\n0.5212173954583703\nconverged\nconverged\nconverged\n0.49076707040891054\nconverged\nconverged\nconverged\n0.4691004266776143\nconverged\nconverged\nconverged\n0.45306392362341275\nconverged\nconverged\nconverged\n0.440790129173547\nconverged\nconverged\nconverged\n0.43112973104231045\nconverged\nconverged\nconverged\n0.4233468637894837\nconverged\nconverged\nconverged\n0.4169527648668737\nconverged\nconverged\nconverged\n0.4116119737830012\nconverged\nconverged\nconverged\n0.4070875198114663\nconverged\nconverged\nconverged\n0.40320769916288557\nconverged\nconverged\nconverged\n0.3998453065473585\nconverged\nconverged\nconverged\n0.39690423558466137\nconverged\nconverged\nconverged\n0.3943105939310044\nconverged\nconverged\nconverged\n0.3920066910702735\nconverged\nconverged\nconverged\n0.3899468609038741\nconverged\nconverged\nconverged\n0.3880944964196533\nconverged\nconverged\nconverged\n0.38641991836018863\nconverged\nconverged\nconverged\n0.3848988258372992\nconverged\nconverged\nconverged\n0.38351113985292606\nconverged\nconverged\nconverged\n0.3822401136998087\nconverged\nconverged\nconverged\n0.3810716917086392\nconverged\nconverged\nconverged\n0.3799939680751413\nconverged\nconverged\nconverged\n0.37899681619368486\nconverged\nconverged\nconverged\n0.3780715476442128\nconverged\nconverged\nconverged\n0.3772106860857456\nconverged\nconverged\nconverged\n0.376407752269879\nconverged\nconverged\nconverged\n0.3756570972409099\nconverged\nconverged\nconverged\n0.37495379484258584\nconverged\nconverged\nconverged\n0.37429350827820596\nconverged\nconverged\nconverged\n0.3736724159773439\nconverged\nconverged\nconverged\n0.3730871337559074\nconverged\nconverged\nconverged\n0.37253465921618045\nconverged\nconverged\nconverged\n0.3720123161468655\nconverged\nconverged\nconverged\n0.37151770633645353\nconverged\nconverged\nconverged\n0.3710486799199135\nconverged\nconverged\nconverged\n0.37060330943204467\nconverged\nconverged\nconverged\n0.37017984532751147\nconverged\nconverged\nconverged\n0.36977671968750647\nconverged\nconverged\nconverged\n0.36939249803312124\nconverged\nconverged\nconverged\n0.36902588673867276\nconverged\nconverged\nconverged\n0.3686756996717303\nconverged\nconverged\nconverged\n0.3683408618997782\nconverged\nconverged\nconverged\n0.36802038374356927\nconverged\nconverged\nconverged\n0.3677133681904523\nconverged\nconverged\nconverged\n0.36741897753439845\nconverged\nconverged\nconverged\n0.36713645190931854\nconverged\nconverged\nconverged\n0.3668650907557457\nconverged\nconverged\nconverged\n0.3666042454075068\nconverged\nconverged\nconverged\n0.36635331538505855\nconverged\nconverged\nconverged\n0.3661117446888238\nconverged\nconverged\nconverged\n0.3658790255058556\nconverged\nconverged\nconverged\n0.36565467226319015\nconverged\nconverged\nconverged\n0.36543825128115714\nconverged\nconverged\nconverged\n0.36522934370674187\nconverged\nconverged\nconverged\n0.3650275640469044\nconverged\nconverged\nconverged\n0.3648325601685791\nconverged\nconverged\nconverged\n0.3646439910586923\nconverged\nconverged\nconverged\n0.3644615416508168\nconverged\nconverged\nconverged\n0.3642849265318363\nconverged\nconverged\nconverged\n0.36411386770196263\nconverged\nconverged\nconverged\n0.36394810198806216\nconverged\nconverged\nconverged\n0.3637873921636492\nconverged\nconverged\nconverged\n0.3636315121222288\nconverged\nconverged\nconverged\n0.3634802431706339\nconverged\nconverged\nconverged\n0.3633333888556808\nconverged\nconverged\nconverged\n0.36319075272418555\nconverged\nconverged\nconverged\n0.3630521605629474\nconverged\nconverged\nconverged\n0.3629174381587654\nconverged\nconverged\nconverged\n0.36278642983175813\nconverged\nconverged\nconverged\n0.36265898360870774\nconverged\nconverged\nconverged\n0.3625349586363882\nconverged\nconverged\nconverged\n0.3624142103549093\nconverged\nconverged\nconverged\n0.3622966201510279\nconverged\nconverged\nconverged\n0.3621820619981735\nconverged\nconverged\nconverged\n0.36207041728310274\nconverged\nconverged\nconverged\n0.36196158221922814\nconverged\nconverged\nconverged\n0.3618554456066341\nconverged\nconverged\nconverged\n0.36175191477872426\nconverged\nconverged\nconverged\n0.3616508896555751\nconverged\nconverged\nconverged\n0.361552284983918\nconverged\nconverged\nconverged\n0.36145600809715683\nconverged\nconverged\nconverged\n0.36136198486201465\nconverged\nconverged\nconverged\n0.36127013743855063\nconverged\nconverged\nconverged\n0.36118038428015997\nconverged\nconverged\nconverged\n0.3610926586668938\nconverged\nconverged\nconverged\n0.36100689017213883\nconverged\nconverged\nconverged\n0.3609230194892733\nconverged\nconverged\nconverged\n0.3608409798983484\nconverged\nconverged\nconverged\n0.3607607120927423\nconverged\nconverged\nconverged\n0.36068216047249724\nconverged\nconverged\nconverged\n0.36060527314431956\nconverged\nconverged\nconverged\n0.36052999080158776\nconverged\nconverged\nconverged\n0.3604562726709991\nconverged\nconverged\nconverged\n0.3603840631525963\nconverged\nconverged\nconverged\n0.36031331776641307\nconverged\nconverged\nconverged\n0.3602439957391469\nconverged\n"
     ]
    }
   ],
   "source": [
    "zeta_seq_2 = np.linspace(0.1, 10, 100)\n",
    "theo_gaus = np.zeros(100)\n",
    "for j in range(1):\n",
    "    lbd = gamma * sigma ** 2 / alpha ** 2\n",
    "    lbd = 1\n",
    "    for i in range(100):\n",
    "        zeta = zeta_seq_2[i]\n",
    "        bias_1 = alpha ** 2\n",
    "        bias_2 = 2 * alpha ** 2 / gamma * inv_MP_Stieltjes_0(gamma, zeta, lbd)\n",
    "        bias_3 = alpha ** 2 / gamma * inv_MP_Stieltjes_derivative_0(gamma, zeta, lbd)\n",
    "        theo_gaus[i] = bias_1 - bias_2 + bias_3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1218e1080>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAGopJREFUeJzt3X9w3PV95/Hna3e1lgSWMbGMjW2wQwyYQA8XlaahzYRQwOn1gJvkcpDpHck0ZTpTkjTp9AZu5pIbMjfD3LQN7dTXxsm5l5kmcTO0k6g59wghJck1gbGc+EhsYzAmYMUYC8vY8k/9et8f+13rq9VKWtmS1ny/r8dko/1+9vPd7/sL6LVffb6f/X4VEZiZWT4Uml2AmZnNH4e+mVmOOPTNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHHPpmZjni0Dczy5FSI50kbQD+HCgCX4qIR2te/zxwa7LYDiyNiEuS10aAnyavvRoRd021rSVLlsTq1asb3gEzM4Pt27e/ERGd0/WbNvQlFYGNwO1AL7BNUndE7Kr2iYhPpfp/HFifeotTEXFjo4WvXr2anp6eRrubmRkg6ZVG+jUyvHMzsDci9kXEILAFuHuK/vcBX2tk42ZmNr8aCf0VwP7Ucm/SNoGkK4E1wHdTza2SeiQ9I+meSdZ7IOnT09fX12DpZmY2U42Evuq0TXZpznuBxyNiJNV2RUR0AR8GHpN01YQ3i9gUEV0R0dXZOe2QlJmZnaNGQr8XWJVaXgkcmKTvvdQM7UTEgeTnPuBpxo/3m5nZPGok9LcBayWtkVSmEuzdtZ0kXQMsBn6UalssaUHyfAlwC7Crdl0zM5sf087eiYhhSQ8CT1CZsrk5InZKegToiYjqB8B9wJYYf1eWdcAXJI1S+YB5ND3rx8zM5pcutDtndXV1hadsmpnNjKTtyfnTKWXmG7nHzwzzZ0++wI79bza7FDOzC1ZmQn9oeJS/eOpFfvLqkWaXYmZ2wcpM6LeViwCcGhqZpqeZWX5lJvQXlApIcGrQoW9mNpnMhL4k2lqKnHTom5lNKjOhD9BeLnp4x8xsCpkK/daWood3zMymkKnQby879M3MppKp0G8rlzjp4R0zs0llK/RbCpz2kb6Z2aQyFfrt5RInh4abXYaZ2QUrU6HvKZtmZlPLVuiXix7eMTObQqZCv71c9IlcM7MpZCr02zxP38xsStkK/XKRM8OjjIxeWPcIMDO7UGQr9Ft8pU0zs6lkKvTbq5dX9hCPmVldmQr9tnLllr8OfTOz+rIV+h7eMTObUqZCvzq8c3LQ38o1M6snU6Hf2uIxfTOzqTQU+pI2SNojaa+kh+q8/nlJO5LHC5LeTL12v6QXk8f9s1l8rXbfJ9fMbEql6TpIKgIbgduBXmCbpO6I2FXtExGfSvX/OLA+eX4p8FmgCwhge7LukVndi0Tb2eEdh76ZWT2NHOnfDOyNiH0RMQhsAe6eov99wNeS53cCT0ZEfxL0TwIbzqfgqfhErpnZ1BoJ/RXA/tRyb9I2gaQrgTXAd2e67mzwPH0zs6k1Evqq0zbZdQ7uBR6PiGrqNrSupAck9Ujq6evra6Ck+to8pm9mNqVGQr8XWJVaXgkcmKTvvYwN7TS8bkRsioiuiOjq7OxsoKT6Wkse0zczm0ojob8NWCtpjaQylWDvru0k6RpgMfCjVPMTwB2SFktaDNyRtM2JQkG0thQ45Xn6ZmZ1TTt7JyKGJT1IJayLwOaI2CnpEaAnIqofAPcBWyIiUuv2S/oclQ8OgEcion92d2G89nLJwztmZpOYNvQBImIrsLWm7TM1y/91knU3A5vPsb4Z8y0Tzcwml6lv5EJyy0Qf6ZuZ1ZW50G8v+0jfzGwymQv9Vg/vmJlNKnOh3+7hHTOzSWUu9H0i18xsctkL/XLRl2EwM5tE5kK/vVz0PH0zs0lkLvQrwzv+Rq6ZWT3ZC/1yidNDo4yOTnZNODOz/Mpe6CfX1D897CEeM7NamQt9X1PfzGxymQt93zLRzGxy2Qv96vCOZ/CYmU2QudBv95G+mdmkMhf61SN9h76Z2UTZC/2yh3fMzCaTudBvL1fuC+MjfTOziTIX+tXhHV+KwcxsouyF/tl5+r4Ug5lZrcyGvod3zMwmyl7oe3jHzGxSmQv9YkGUSwVfhsHMrI7MhT74mvpmZpNpKPQlbZC0R9JeSQ9N0udDknZJ2inpq6n2EUk7kkf3bBU+lXbfMtHMrK7SdB0kFYGNwO1AL7BNUndE7Er1WQs8DNwSEUckLU29xamIuHGW655Sq2+ZaGZWVyNH+jcDeyNiX0QMAluAu2v6/B6wMSKOAETEodktc2Y8vGNmVl8job8C2J9a7k3a0q4Grpb0L5KekbQh9VqrpJ6k/Z56G5D0QNKnp6+vb0Y7UI9vmWhmVt+0wzuA6rTV3ouwBKwF3gusBH4g6fqIeBO4IiIOSHo78F1JP42Il8a9WcQmYBNAV1fXed/nsK1c4uipofN9GzOzzGnkSL8XWJVaXgkcqNPnmxExFBEvA3uofAgQEQeSn/uAp4H151nztNpbiv5GrplZHY2E/jZgraQ1ksrAvUDtLJxvALcCSFpCZbhnn6TFkhak2m8BdjHH2sqevWNmVs+0wzsRMSzpQeAJoAhsjoidkh4BeiKiO3ntDkm7gBHgjyPisKR3A1+QNErlA+bR9KyfudJWLvrSymZmdTQypk9EbAW21rR9JvU8gE8nj3SfHwI3nH+ZM9PmefpmZnVl+hu5lc8iMzOrymTot7YUiYAzw6PNLsXM7IKSydBvP3tNfQ/xmJmlZTr0T/pkrpnZOJkM/dYW3z3LzKyeTIZ+9ebopwY9pm9mlpbJ0K/ePcvX3zEzGy+boV/2LRPNzOrJZOh79o6ZWX2ZDP2x4R2HvplZWiZDv93DO2ZmdWUy9Fs9vGNmVlcmQ7/dwztmZnVlMvRLxQLt5SLHTvvuWWZmaZkMfYCO1haO+ZaJZmbjZDb0F7W1+D65ZmY1Mh36Ht4xMxsvs6Hf0Vbi6ClfhsHMLC3Doe8xfTOzWtkNfZ/INTObILOhv6ithYEzw4yM+j65ZmZVmQ59gAGfzDUzO6uh0Je0QdIeSXslPTRJnw9J2iVpp6Svptrvl/Ri8rh/tgqfTkcS+p62aWY2pjRdB0lFYCNwO9ALbJPUHRG7Un3WAg8Dt0TEEUlLk/ZLgc8CXUAA25N1j8z+roxXPdI/5hk8ZmZnNXKkfzOwNyL2RcQgsAW4u6bP7wEbq2EeEYeS9juBJyOiP3ntSWDD7JQ+tY7WyueZj/TNzMY0EvorgP2p5d6kLe1q4GpJ/yLpGUkbZrDunFjUnhzpe0zfzOysaYd3ANVpq50SUwLWAu8FVgI/kHR9g+si6QHgAYArrriigZKmt8hj+mZmEzRypN8LrEotrwQO1OnzzYgYioiXgT1UPgQaWZeI2BQRXRHR1dnZOZP6J9XR6tA3M6vVSOhvA9ZKWiOpDNwLdNf0+QZwK4CkJVSGe/YBTwB3SFosaTFwR9I259rLRUoF+QtaZmYp0w7vRMSwpAephHUR2BwROyU9AvRERDdj4b4LGAH+OCIOA0j6HJUPDoBHIqJ/LnakliQ6fKVNM7NxGhnTJyK2Altr2j6Teh7Ap5NH7bqbgc3nV+a5qVxp01M2zcyqMvuNXKhM2/SRvpnZmGyHvq+0aWY2TqZDf5FD38xsnEyHvk/kmpmNl+nQr94ysXKe2czMMh36Ha0tDI0Ep4ZGml2KmdkFIdOh7yttmpmNl4vQ97i+mVlFpkO/o63y3TNfadPMrCLToX/2SP+kQ9/MDDIe+r7SppnZeJkO/bMncj28Y2YGZDz0F/qWiWZm42Q69EvFAhcvKHnKpplZItOhD5UhHh/pm5lVZD70F7aWPKZvZpbIfOj7SN/MbEwuQt+XVzYzq8h86PtGKmZmYzIf+h7eMTMbk/nQ72ht4cTgCMMjo80uxcys6TIf+ovOXnTNc/XNzBoKfUkbJO2RtFfSQ3Ve/4ikPkk7ksfHUq+NpNq7Z7P4Rixqr15T30M8Zmal6TpIKgIbgduBXmCbpO6I2FXT9e8i4sE6b3EqIm48/1LPjS+6ZmY2ppEj/ZuBvRGxLyIGgS3A3XNb1uzxjVTMzMY0EvorgP2p5d6krdYHJD0n6XFJq1LtrZJ6JD0j6Z7zKfZcdPhKm2ZmZzUS+qrTFjXL/wisjohfAr4DfDn12hUR0QV8GHhM0lUTNiA9kHww9PT19TVYemN8pG9mNqaR0O8F0kfuK4ED6Q4RcTgiziSLXwRuSr12IPm5D3gaWF+7gYjYFBFdEdHV2dk5ox2Yjm+ObmY2ppHQ3waslbRGUhm4Fxg3C0fS8tTiXcDupH2xpAXJ8yXALUDtCeA5taBUoFws+EjfzIwGZu9ExLCkB4EngCKwOSJ2SnoE6ImIbuATku4ChoF+4CPJ6uuAL0gapfIB82idWT9zShKXXlTm8PEz03c2M8u4aUMfICK2Altr2j6Tev4w8HCd9X4I3HCeNZ63ZYtaOXjsdLPLMDNrusx/Ixdg+aJWXjvq0Dczy0XoL1vUykGHvplZPkJ/+aJWjp8ZZsBz9c0s53IR+pd1tALwusf1zSznchH6yxe1AXhc38xyLyehXznSd+ibWd7lIvSXdiwA8MlcM8u9XIT+glKRJReXfaRvZrmXi9CHyrRNn8g1s7zLT+h3+AtaZmb5Cf1FrRw8eqrZZZiZNVVuQn/5ojaOnBzi9NBIs0sxM2ua3IT+suQLWp7BY2Z5lpvQ91x9M7Mchf6yJPQPHvO4vpnlV/5C/6hvpmJm+ZWb0G8vl+hoLXkGj5nlWm5CHyozeDymb2Z5lqvQ920TzSzvchX6vm2imeVdrkL/so5W3jh+hqGR0WaXYmbWFLkK/eWLWomAQwOewWNm+ZSr0B+btukZPGaWTw2FvqQNkvZI2ivpoTqvf0RSn6QdyeNjqdful/Ri8rh/NoufKd820czyrjRdB0lFYCNwO9ALbJPUHRG7arr+XUQ8WLPupcBngS4ggO3JukdmpfoZGjvSd+ibWT41cqR/M7A3IvZFxCCwBbi7wfe/E3gyIvqToH8S2HBupZ6/jtYS7eWij/TNLLcaCf0VwP7Ucm/SVusDkp6T9LikVTNZV9IDknok9fT19TVY+sxJ4vJL2ni1/+ScbcPM7ELWSOirTlvULP8jsDoifgn4DvDlGaxLRGyKiK6I6Ors7GygpHN3zbKFPH/w2Jxuw8zsQtVI6PcCq1LLK4ED6Q4RcTgiqvMgvwjc1Oi68+265R3s7z/FwOmhZpZhZtYUjYT+NmCtpDWSysC9QHe6g6TlqcW7gN3J8yeAOyQtlrQYuCNpa5p1yxcC8PzBgWaWYWbWFNPO3omIYUkPUgnrIrA5InZKegToiYhu4BOS7gKGgX7gI8m6/ZI+R+WDA+CRiOifg/1o2LrlHQDsfu0Yv7L60maWYmY276YNfYCI2ApsrWn7TOr5w8DDk6y7Gdh8HjXOqmUdrVzS3sLu1zyub2b5k6tv5EJlBs+6ZR3ses3DO2aWP7kLfYBrly/khYMDjIxOmEhkZpZpuQz9dcs7ODU0wiuHTzS7FDOzeZXL0L/u7MlcD/GYWb7kMvTfsfRiigX5ZK6Z5U4uQ7+1pchVnRc59M0sd3IZ+lAZ13fom1ne5Dr0Dxw9zZsnB5tdipnZvMlt6F+7rHI5Bp/MNbM8yW3oX5e6HIOZWV7kNvQ7Fy7gbReVHfpmliu5DX1J3LjqEp59uZ8IfzPXzPIht6EPcOu1S3m1/yR7Dx1vdilmZvMi16F/27qlAHxn96EmV2JmNj9yHfrLF7Xxzss7eGr3680uxcxsXuQ69AFuW3cZP371CP0nPF/fzLIv96H/m+uWMhrwz897iMfMsi/3oX/95YtYunABTz3vIR4zy77ch36hIG5bt5Tvv/AGg8OjzS7HzGxO5T70Ad537WUcPzPMsy8fbnYpZmZzyqEP/Po7lrCgVOApT900s4xz6ANt5SLvubqTbz13gNNDI80ux8xszjQU+pI2SNojaa+kh6bo90FJIakrWV4t6ZSkHcnjr2er8Nn2u7++hjeOD/L1nv3NLsXMbM5MG/qSisBG4P3AdcB9kq6r028h8Ang2ZqXXoqIG5PH789CzXPiV9dcyk1XLuYL39vH0IhP6JpZNjVypH8zsDci9kXEILAFuLtOv88B/x04PYv1zRtJ/MGtV/GLN0/RveNAs8sxM5sTjYT+CiA95tGbtJ0laT2wKiK+VWf9NZJ+Iul7kn7j3Eude7des5Rrly3kfzy9l9FRX3nTzLKnkdBXnbaziSipAHwe+KM6/V4DroiI9cCnga9K6piwAekBST2Sevr6+hqrfA5UjvbfwUt9J/j2roNNq8PMbK40Evq9wKrU8kogPf6xELgeeFrSz4F3Ad2SuiLiTEQcBoiI7cBLwNW1G4iITRHRFRFdnZ2d57Yns+S3bljO6re189h3XvSXtcwscxoJ/W3AWklrJJWBe4Hu6osRcTQilkTE6ohYDTwD3BURPZI6kxPBSHo7sBbYN+t7MYuKBfHwb63j+YMD/Mm39zS7HDOzWTVt6EfEMPAg8ASwG/h6ROyU9Iiku6ZZ/T3Ac5L+H/A48PsR0X++Rc+1O9+5jN951xVs+v4+vvdC84abzMxmmy60WwV2dXVFT09Ps8vg9NAId/3l/6X/xCBbP/kbLF3Y2uySzMwmJWl7RHRN18/fyJ1Ea0uRv/zwLzNwepg/3LLD39Q1s0xw6E/h6ssW8t/+7Q388KXDfPRvtjFweqjZJZmZnReH/jQ+eNNKHvv3N7Lt5/18+IvPcvj4mWaXZGZ2zhz6Dbhn/Qo2/cebeOH1AT7wVz9k+ysX/LloM7O6HPoNet+1l/GVj/0qQyPBB//6R/yXb/zMwz1m9pbj0J+BrtWX8u1PvYePvnsNf/vsK9z2p99j0/df4pjD38zeIjxl8xzt2P8mj/7Tbp7Z18/FC0p8qGsV96y/nBtWLEKqd+UKM7O50+iUTYf+efrZL47yxR/s438/9xrDo8GKS9q4853L+LWr3savrF7MJe3lZpdoZjmQz9D/p4fg4E9nt6AGDY2OcuTEIP0nBjl6aujsFenaWoq0l4u0l0u0l4u0thRZUCpQ8F8DZrMuiJrl2id1F5ksB2vXnyotgxjXYdK+Mflrw0vfycJ7/nSKrUyu0dAvndO72wQthQJLF7aydGEroxEcPzPMwOkhBk4PM3B6mMMnBsf3LxYoFwuUS6KlWKBULNBSEMWCKBUKFJPnhQIUJIoSUuW5BKp78dMLU/UXMWp+caq/aJH837TtycJYvynePxpoTxVTd9updSHq1F+zhzW/zLXbm/b12lpTjfXWG98vSOdWOqzG9anpMNn7nnPf1D5MqCW14oTQnWJb9WqZ2Ccb+o708W/meBvZCv33P9rsCoDK2fGO5FF17PQQL74+wKv9J9nff4r9/Sc5NHCGQwNn6Bs4zZGTQ4zM4Br+5WKBUlGUCqJULCQfFqKgsQ8KUflJ5X91zzVEjA+s6nIkYTSaXGh0NILRqATLaFT6VdpItafaRsee55EERensP//i2Q9rKBR09t9HsSAKAqj8LNR+uKvSr5D8HL9+tV+1vXIwMG4dag4URNI+fjuM6ze2LaXrrqkp3S9dH+PeP11T5cWxbVReq65Lqm1cnzrvqzp9keq2V/dvym0w9s+h2nfstfH/HKod662bbqv+3qXf7+z2J9nm4vaWmf/HNkPZCv0LWEdrCzddeSk3XXlp3dcjgoEzw7x5YoiBM0McPz3M8TPDnBoa4dTgCKeGRjgzNMrgyChnhkcZGhlleGSUoZFgaGSU0QhGRoPh0Up6VwM3qDwff4QU4/9SmPCLPvZLXv1FK6YCpVDQ2V/mYmEsUIqFsRAqiOS1yntW/3Kh+jwVboWz76k6yxPfpxo+xSQ0q6Fab91xz6VxoVosMOF9qdlGdV1V/9oqMKF/uk86wMwuRA79C4QkOlpb6Gid+096M8svz9M3M8sRh76ZWY449M3McsShb2aWIw59M7McceibmeWIQ9/MLEcc+mZmOXLBXXBNUh/wygxXWwK8MQflXMjyuM+Qz/3O4z5DPvf7fPb5yojonK7TBRf650JSTyNXl8uSPO4z5HO/87jPkM/9no999vCOmVmOOPTNzHIkK6G/qdkFNEEe9xnyud953GfI537P+T5nYkzfzMwak5UjfTMza8BbOvQlbZC0R9JeSQ81u575IGmVpH+WtFvSTkmfbHZN80VSUdJPJH2r2bXMF0mXSHpc0vPJv/Nfa3ZNc03Sp5L/tn8m6WuSWptd01yQtFnSIUk/S7VdKulJSS8mPxfP9nbfsqEvqQhsBN4PXAfcJ+m65lY1L4aBP4qIdcC7gD/IyX4DfBLY3ewi5tmfA/8nIq4F/hUZ339JK4BPAF0RcT1QBO5tblVz5n8BG2raHgKeioi1wFPJ8qx6y4Y+cDOwNyL2RcQgsAW4u8k1zbmIeC0ifpw8H6ASAiuaW9Xck7QS+NfAl5pdy3yR1AG8B/ifABExGBFvNreqeVEC2iSVgHbgQJPrmRMR8X2gv6b5buDLyfMvA/fM9nbfyqG/AtifWu4lB+GXJmk1sB54trmVzIvHgP8EjDa7kHn0dqAP+JtkWOtLki5qdlFzKSJ+AfwJ8CrwGnA0Ir7d3Krm1WUR8RpUDvCApbO9gbdy6Ne7+3RupiJJuhj4e+API+JYs+uZS5J+GzgUEdubXcs8KwG/DPxVRKwHTjAHf+5fSJIx7LuBNcDlwEWSfqe5VWXLWzn0e4FVqeWVZPTPwFqSWqgE/lci4h+aXc88uAW4S9LPqQzjvU/S3za3pHnRC/RGRPUvucepfAhk2W8CL0dEX0QMAf8AvLvJNc2n1yUtB0h+HprtDbyVQ38bsFbSGkllKid7uptc05yTJCpjvLsj4s+aXc98iIiHI2JlRKym8u/5uxGR+aO/iDgI7Jd0TdJ0G7CriSXNh1eBd0lqT/5bv42Mn7yu0Q3cnzy/H/jmbG+gNNtvOF8iYljSg8ATVM7wb46InU0uaz7cAvwH4KeSdiRt/zkitjaxJps7Hwe+khzY7AM+2uR65lREPCvpceDHVGaq/YSMfjNX0teA9wJLJPUCnwUeBb4u6XepfAD+u1nfrr+Ra2aWH2/l4R0zM5shh76ZWY449M3McsShb2aWIw59M7McceibmeWIQ9/MLEcc+mZmOfL/AexlGA/hpGlvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(zeta_seq_2, theo_gaus)\n",
    "plt.plot(zeta_seq_2, np.ones(100) * MSE_original(lbd, gamma, alpha, sigma))"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}