2D_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PDE 1 - 2D\n",
    "\n",
    "\n",
    "#### Problem Setup\n",
    "\n",
    "$\\phi u + u_{x} + u_{y,y} = f(x,y)$\n",
    "\n",
    "For the generation of our initial data samples we use:\n",
    "\n",
    "$\\phi = 2$ <br>\n",
    "$u: \\mathbb{R}^2 \\rightarrow \\mathbb{R}, \\; u(x,y) = x^2 + y$ <br>\n",
    "$f: \\mathbb{R}^2 \\rightarrow \\mathbb{R}, \\;f(x,y) = 2(x^2 + x + y)$ <br>\n",
    "$X_i := (x_i, y_i) \\in [0,1] \\times [0,1] \\subseteq \\mathbb{R}^2$ for $i \\in \\{1, \\dotsc, n\\}$ \n",
    "\n",
    "and our known function values will be $\\{u(x_i,y_i), f(x_i,y_i)\\}_{i \\in \\{1, \\dotsc, n\\}}$.\n",
    "\n",
    "We assume that $u$ can be represented as a Gaussian process with SE kernel.\n",
    "\n",
    "$u \\sim \\mathcal{GP}(0, k_{uu}(X_i, X_j; \\theta))$, where $\\theta = \\{\\sigma, l_x, l_y\\}$.\n",
    "\n",
    "Set the linear operator to:\n",
    "\n",
    "$\\mathcal{L}_X^{\\phi} := \\phi + \\partial_x + \\partial_{y,y}$\n",
    "\n",
    "so that\n",
    "\n",
    "$\\mathcal{L}_X^{\\phi} u = f$\n",
    "\n",
    "Problem at hand: Estimate $\\phi$ (we expect $\\phi = 2$).\n",
    "\n",
    "\n",
    "#### Step 1: Simulate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "from scipy.linalg import solve_triangular\n",
    "import scipy.optimize as opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Global variables: x, y, n, y_u, y_f, s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Parameters, that can be modified:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of data samples:\n",
    "n = 5\n",
    "\n",
    "# Noise of our data:\n",
    "s = 1e-7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_data():\n",
    "    x = np.random.rand(n)\n",
    "    y = np.random.rand(n)\n",
    "    y_u = np.multiply(x,x) + y\n",
    "    y_f = 2*(np.multiply(x,x) + x + y)\n",
    "    return (x,y,y_u,y_f)\n",
    "(x,y,y_u,y_f) = simulate_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Evaluate kernels\n",
    "\n",
    "$k_{uu}(X_i, X_j; \\theta) = \\sigma exp(-\\frac{1}{2l_x}(x_i-x_j)^2 - \\frac{1}{2l_y}(y_i-y_j)^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi = sp.symbols('x_i x_j y_i y_j sigma l_x l_y phi')\n",
    "kuu_sym = sigma*sp.exp(-1/(2*l_x)*((x_i - x_j)**2) - 1/(2*l_y)*((y_i - y_j)**2))\n",
    "kuu_fn = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y), kuu_sym, \"numpy\")\n",
    "def kuu(x, y, sigma, l_x, l_y):\n",
    "    k = np.zeros((x.size, x.size))\n",
    "    for i in range(x.size):\n",
    "        for j in range(x.size):\n",
    "            k[i,j] = kuu_fn(x[i], x[j], y[i], y[j], sigma, l_x, l_y)\n",
    "    return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$k_{ff}(X_i,X_j;\\theta,\\phi) \\\\\n",
    "= \\mathcal{L}_{X_i}^{\\phi} \\mathcal{L}_{X_j}^{\\phi} k_{uu}(X_i, X_j; \\theta) \\\\\n",
    "= \\phi^2k_{uu} + \\phi \\frac{\\partial}{\\partial x_i}k_{uu} + \\phi \\frac{\\partial^2}{\\partial y_i^2}k_{uu} + \\phi \\frac{\\partial}{\\partial x_j}k_{uu} + \\frac{\\partial^2}{\\partial x_i, x_j}k_{uu} + \\frac{\\partial^3}{\\partial y_i^2 \\partial x_j}k_{uu} + \\phi \\frac{\\partial^2}{\\partial y_j^2}k_{uu} + \\frac{\\partial^3}{\\partial x_i \\partial y_j^2}k_{uu} + \\frac{\\partial^4}{\\partial y_i^2 \\partial y_j^2}k_{uu}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "kff_sym = phi**2*kuu_sym \\\n",
    "        + phi*sp.diff(kuu_sym, x_i) \\\n",
    "        + phi*sp.diff(kuu_sym, y_i, y_i) \\\n",
    "        + phi*sp.diff(kuu_sym, x_j) \\\n",
    "        + sp.diff(kuu_sym, x_i, x_j) \\\n",
    "        + sp.diff(kuu_sym, y_i, y_i, x_j) \\\n",
    "        + phi*sp.diff(kuu_sym, y_j, y_j) \\\n",
    "        + sp.diff(kuu_sym, x_i, y_j, y_j) \\\n",
    "        + sp.diff(kuu_sym, y_i, y_i, y_j, y_j)\n",
    "kff_fn = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kff_sym, \"numpy\")\n",
    "def kff(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros((x.size, x.size))\n",
    "    for i in range(x.size):\n",
    "        for j in range(x.size):\n",
    "            k[i,j] = kff_fn(x[i], x[j], y[i], y[j], sigma, l_x, l_y, phi)\n",
    "    return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$k_{fu}(X_i,X_j;\\theta,\\phi) \\\\\n",
    "= \\mathcal{L}_{X_i}^{\\phi} k_{uu}(X_i, X_j; \\theta) \\\\\n",
    "= \\phi k_{uu} + \\frac{\\partial}{\\partial x_i}k_{uu} + \\frac{\\partial^2}{\\partial y_i^2}k_{uu}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "kfu_sym = phi*kuu_sym \\\n",
    "        + sp.diff(kuu_sym, x_i) \\\n",
    "        + sp.diff(kuu_sym, y_i, y_i)\n",
    "kfu_fn = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kfu_sym, \"numpy\")\n",
    "def kfu(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros((x.size, x.size))\n",
    "    for i in range(x.size):\n",
    "        for j in range(x.size):\n",
    "            k[i,j] = kfu_fn(x[i], x[j], y[i], y[j], sigma, l_x, l_y, phi)\n",
    "    return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kuf(x, y, sigma, l_x, l_y, phi):\n",
    "    return kfu(x, y, sigma, l_x, l_y, phi).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Computing the negative log-likelihood \n",
    "(with block matrix inversion, Cholesky decomposition, potentially SVD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the block-inversion technique: Let\n",
    "$ K = \\begin{pmatrix} K_{uu} & K_{uf} \\\\ K_{fu} & K_{ff} \\end{pmatrix} = \\begin{pmatrix} A & B \\\\ B^T & C \\end{pmatrix}$. \n",
    "\n",
    "Then $det(K) = det(A) det(C-B^T A^{-1} B)$.\n",
    "\n",
    "$K^{-1} = \\begin{pmatrix} A^{-1} + A^{-1} B(C-B^T A^{-1} B)^{-1}B^T A^{-1} & -A^{-1}B(C-B^T A^{-1} B)^{-1} \\\\\n",
    "            -(C - B^T A^{-1}B)^{-1}B^T A^{-1} & (C-B^T A^{-1} B)^{-1} \\end{pmatrix}$\n",
    "            \n",
    "So it suffices to invert $A$ and $C-B^T A^{-1} B$.\n",
    "\n",
    "A theorem about Schur-complements ensures that $K$ positive-definite implies the positive-definiteness of $K/A = C-B^T A^{-1} B$ as well, so Cholesky should work as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nlml(params):\n",
    "    \n",
    "    sigma_exp = np.exp(params[0])\n",
    "    l_x_exp = np.exp(params[1])\n",
    "    l_y_exp = np.exp(params[2])\n",
    "    # phi = params[3]\n",
    "    \n",
    "    A = kuu(x, y, sigma_exp, l_x_exp, l_y_exp) + s*np.eye(n)\n",
    "    B = kfu(x, y, sigma_exp, l_x_exp, l_y_exp, params[3]).T\n",
    "    C = kff(x, y, sigma_exp, l_x_exp, l_y_exp, params[3]) + s*np.eye(n)\n",
    "    \n",
    "    # Inversion of A\n",
    "    A_inv = np.zeros((n, n))\n",
    "    \n",
    "    try:\n",
    "        L = np.linalg.cholesky(A)\n",
    "        L_inv = solve_triangular(L, np.identity(n), lower=True) # Slight performance boost \n",
    "                                                                # over np.linalg.inv\n",
    "        A_inv = L_inv.T @ L_inv\n",
    "        logdet_A = 2*np.log(np.abs(np.diag(L))).sum()\n",
    "    except np.linalg.LinAlgError:\n",
    "        # Inverse of K via SVD\n",
    "        u, s_mat, vt = np.linalg.svd(A)\n",
    "        A_inv = vt.T @ np.linalg.inv(np.diag(s_mat)) @ u.T\n",
    "        logdet_A = np.log(s_mat).sum()\n",
    "        \n",
    "    # Inversion of $C-B^T A^{-1} B$\n",
    "    KA_inv = np.zeros((n, n))\n",
    "    KA = C - B.T @ A_inv @ B\n",
    "    \n",
    "    try:\n",
    "        L = np.linalg.cholesky(KA)\n",
    "        L_inv = solve_triangular(L, np.identity(n), lower=True) # Slight performance boost \n",
    "                                                                # over np.linalg.inv\n",
    "        KA_inv = L_inv.T @ L_inv\n",
    "        logdet_KA = 2*np.log(np.abs(np.diag(L))).sum()\n",
    "    except np.linalg.LinAlgError:\n",
    "        # Inverse of K via SVD\n",
    "        u, s_mat, vt = np.linalg.svd(KA)\n",
    "        KA_inv = vt.T @ np.linalg.inv(np.diag(s_mat)) @ u.T\n",
    "        logdet_KA = np.log(s_mat).sum()\n",
    "        \n",
    "    # Piecing it together\n",
    "    T = A_inv @ B @ KA_inv\n",
    "    yKy = y_u @ (A_inv + T @ B.T @ A_inv) @ y_u - 2*y_u @ T @ y_f + y_f @ KA_inv @ y_f\n",
    "    \n",
    "    return (yKy + logdet_A + logdet_KA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4: Optimize hyperparameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1. Nelder-Mead**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nfeval = 1\n",
    "def callbackF(Xi):\n",
    "    global Nfeval\n",
    "    print('{0:4d}   {1: 3.6f}   {2: 3.6f}   {3: 3.6f}   {4: 3.6f}'.\n",
    "          format(Nfeval, Xi[0], Xi[1], Xi[2], Xi[3]))\n",
    "    Nfeval += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = time.time()\n",
    "m_n = opt.minimize(nlml, np.random.rand(4), method=\"Nelder-Mead\", callback = callbackF,\n",
    "                                        options={'maxfev':5000, 'fatol':0.001, 'xatol':0.001})\n",
    "t_Nelder = time.time() - t0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. Nonlinear Conjugate Gradient**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The partial derivatives for $L(\\theta, \\phi) := \\log(\\lvert K(\\theta, \\phi) \\rvert) + y^T [K(\\theta, \\phi)]^{-1} y$ are given by: <br>\n",
    "\\begin{equation}\n",
    "\\frac{\\partial L}{\\partial \\theta_i}(\\theta, \\phi) = tr(K^{-1}\\frac{\\partial K}{\\partial \\theta_i}) - y^T K^{-1}\\frac{\\partial K}{\\partial \\theta_i} K^{-1}y\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial L}{\\partial \\phi}(\\theta, \\phi) = tr(K^{-1}\\frac{\\partial K}{\\partial \\phi}) - y^T K^{-1}\\frac{\\partial K}{\\partial \\phi} K^{-1}y\n",
    "\\end{equation}\n",
    "\n",
    "Since $nlml(\\sigma, l_x, l_y, \\phi) = L(e^{\\sigma}, e^{l_x}, e^{l_y}, \\phi)$,  it follows:\n",
    "\n",
    "\\begin{equation}\n",
    "\\nabla nlml(\\sigma, l_x, l_y, \\phi) = (\\nabla_{(e^{\\sigma}, e^{l_x}, e^{l_y}, \\phi)}L)^T\n",
    "\\begin{pmatrix}\n",
    "e^{\\sigma} &0&0&0\\\\\n",
    "0&e^{l_x}&0&0 \\\\\n",
    "0&0&e^{l_y}&0 \\\\\n",
    "0&0&0&1\n",
    "\\end{pmatrix}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Building the gradient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "kuu_sym_sigma, kuu_sym_l_x, kuu_sym_l_y, kuu_sym_phi = [sp.diff(kuu_sym, l) for l in [sigma, l_x, l_y, phi]]\n",
    "kfu_sym_sigma, kfu_sym_l_x, kfu_sym_l_y, kfu_sym_phi = [sp.diff(kfu_sym, l) for l in [sigma, l_x, l_y, phi]]\n",
    "kff_sym_sigma, kff_sym_l_x, kff_sym_l_y, kff_sym_phi = [sp.diff(kff_sym, l) for l in [sigma, l_x, l_y, phi]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "kuu_sigma = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y), kuu_sym_sigma, \"numpy\")\n",
    "kff_sigma = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kff_sym_sigma, \"numpy\")\n",
    "kfu_sigma = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kfu_sym_sigma, \"numpy\")\n",
    "\n",
    "kuu_l_x = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y), kuu_sym_l_x, \"numpy\")\n",
    "kff_l_x = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kff_sym_l_x, \"numpy\")\n",
    "kfu_l_x = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kfu_sym_l_x, \"numpy\")\n",
    "\n",
    "kuu_l_y = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y), kuu_sym_l_y, \"numpy\")\n",
    "kff_l_y = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kff_sym_l_y, \"numpy\")\n",
    "kfu_l_y = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kfu_sym_l_y, \"numpy\")\n",
    "\n",
    "kuu_phi = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y), kuu_sym_phi, \"numpy\")\n",
    "kff_phi = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kff_sym_phi, \"numpy\")\n",
    "kfu_phi = sp.lambdify((x_i, x_j, y_i, y_j, sigma, l_x, l_y, phi), kfu_sym_phi, \"numpy\")\n",
    "\n",
    "# Derivative of K w.r.t. sigma:\n",
    "def K_sigma(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros([2*n, 2*n])\n",
    "    for i in range(2*n):\n",
    "        for j in range(2*n):\n",
    "            if i<n and j<n:\n",
    "                k[i,j] = kuu_sigma(x[i], x[j], y[i], y[j], sigma, l_x, l_y)\n",
    "            if i<n and j>=n:\n",
    "                k[i,j] = kfu_sigma(x[j-n], x[i], y[j-n], y[i], sigma, l_x, l_y, phi) \n",
    "                # Transpose of kfu\n",
    "            if i>=n and j<n:\n",
    "                k[i,j] = kfu_sigma(x[i-n], x[j], y[i-n], y[j], sigma, l_x, l_y, phi)\n",
    "            if i>=n and j>=n:\n",
    "                k[i,j] = kff_sigma(x[i-n], x[j-n], y[i-n], y[j-n], sigma, l_x, l_y, phi)\n",
    "    return k\n",
    "\n",
    "# Derivative of K w.r.t. l_x:\n",
    "def K_l_x(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros([2*n, 2*n])\n",
    "    for i in range(2*n):\n",
    "        for j in range(2*n):\n",
    "            if i<n and j<n:\n",
    "                k[i,j] = kuu_l_x(x[i], x[j], y[i], y[j], sigma, l_x, l_y)\n",
    "            if i<n and j>=n:\n",
    "                k[i,j] = kfu_l_x(x[j-n], x[i], y[j-n], y[i], sigma, l_x, l_y, phi) # Transpose of kfu\n",
    "            if i>=n and j<n:\n",
    "                k[i,j] = kfu_l_x(x[i-n], x[j], y[i-n], y[j], sigma, l_x, l_y, phi)\n",
    "            if i>=n and j>=n:\n",
    "                k[i,j] = kff_l_x(x[i-n], x[j-n], y[i-n], y[j-n], sigma, l_x, l_y, phi)\n",
    "    return k\n",
    "\n",
    "# Derivative of K w.r.t. l_y:\n",
    "def K_l_y(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros([2*n, 2*n])\n",
    "    for i in range(2*n):\n",
    "        for j in range(2*n):\n",
    "            if i<n and j<n:\n",
    "                k[i,j] = kuu_l_y(x[i], x[j], y[i], y[j], sigma, l_x, l_y)\n",
    "            if i<n and j>=n:\n",
    "                k[i,j] = kfu_l_y(x[j-n], x[i], y[j-n], y[i], sigma, l_x, l_y, phi) # Transpose of kfu\n",
    "            if i>=n and j<n:\n",
    "                k[i,j] = kfu_l_y(x[i-n], x[j], y[i-n], y[j], sigma, l_x, l_y, phi)\n",
    "            if i>=n and j>=n:\n",
    "                k[i,j] = kff_l_y(x[i-n], x[j-n], y[i-n], y[j-n], sigma, l_x, l_y, phi)\n",
    "    return k\n",
    "\n",
    "# Derivative of K w.r.t. phi:\n",
    "def K_phi(x, y, sigma, l_x, l_y, phi):\n",
    "    k = np.zeros([2*n, 2*n])\n",
    "    for i in range(2*n):\n",
    "        for j in range(2*n):\n",
    "            if i<n and j<n:\n",
    "                k[i,j] = kuu_phi(x[i], x[j], y[i], y[j], sigma, l_x, l_y)\n",
    "            if i<n and j>=n:\n",
    "                k[i,j] = kfu_phi(x[j-n], x[i], y[j-n], y[i], sigma, l_x, l_y, phi) # Transpose of kfu\n",
    "            if i>=n and j<n:\n",
    "                k[i,j] = kfu_phi(x[i-n], x[j], y[i-n], y[j], sigma, l_x, l_y, phi)\n",
    "            if i>=n and j>=n:\n",
    "                k[i,j] = kff_phi(x[i-n], x[j-n], y[i-n], y[j-n], sigma, l_x, l_y, phi)\n",
    "    return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def K_inverse(sigma, l_x, l_y, phi):\n",
    "    \n",
    "    A = kuu(x, y, sigma, l_x, l_y) + s*np.eye(n)\n",
    "    B = kfu(x, y, sigma, l_x, l_y, phi).T\n",
    "    C = kff(x, y, sigma, l_x, l_y, phi) + s*np.eye(n)\n",
    "    \n",
    "    # Inversion of A\n",
    "    A_inv = np.zeros((n, n))\n",
    "    \n",
    "    try:\n",
    "        L = np.linalg.cholesky(A)\n",
    "        L_inv = solve_triangular(L, np.identity(n), lower=True) # Slight performance boost \n",
    "                                                                # over np.linalg.inv\n",
    "        A_inv = L_inv.T @ L_inv\n",
    "    except np.linalg.LinAlgError:\n",
    "        # Inverse of K via SVD\n",
    "        u, s_mat, vt = np.linalg.svd(A)\n",
    "        A_inv = vt.T @ np.linalg.inv(np.diag(s_mat)) @ u.T\n",
    "        \n",
    "    # Inversion of $C-B^T A^{-1} B$\n",
    "    KA_inv = np.zeros((n, n))\n",
    "    KA = C - B.T @ A_inv @ B\n",
    "    \n",
    "    try:\n",
    "        L = np.linalg.cholesky(KA)\n",
    "        L_inv = solve_triangular(L, np.identity(n), lower=True) # Slight performance boost \n",
    "                                                                # over np.linalg.inv\n",
    "        KA_inv = L_inv.T @ L_inv\n",
    "    except np.linalg.LinAlgError:\n",
    "        # Inverse of K via SVD\n",
    "        u, s_mat, vt = np.linalg.svd(KA)\n",
    "        KA_inv = vt.T @ np.linalg.inv(np.diag(s_mat)) @ u.T\n",
    "        \n",
    "    # Piecing it together\n",
    "    T = A_inv @ B @ KA_inv\n",
    "        \n",
    "    K_inv = np.block([\n",
    "        [A_inv + T @ B.T @ A_inv, -T],\n",
    "        [-T.T, KA_inv]\n",
    "    ])\n",
    "    \n",
    "    return K_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting the gradient:\n",
    "def grad_L(sigma, l_x, l_y, phi, s):\n",
    "    K_inv = K_inverse(sigma, l_x, l_y, phi)\n",
    "#     print(K_inverse(sigma, l_x, l_y, phi))\n",
    "    \n",
    "    y_con = np.concatenate((y_u, y_f))\n",
    "    comp_1 = np.trace(K_inv.dot(K_sigma(x, y, sigma, l_x, l_y, phi))) - \\\n",
    "             (y_con.T).dot(K_inv).dot(K_sigma(x, y, sigma, l_x, l_y, phi)).dot(K_inv).dot(y_con)\n",
    "    comp_2 = np.trace(K_inv.dot(K_l_x(x, y, sigma, l_x, l_y, phi))) - \\\n",
    "             (y_con.T).dot(K_inv).dot(K_l_x(x, y, sigma, l_x, l_y, phi)).dot(K_inv).dot(y_con)\n",
    "    comp_3 = np.trace(K_inv.dot(K_l_y(x, y, sigma, l_x, l_y, phi))) - \\\n",
    "             (y_con.T).dot(K_inv).dot(K_l_y(x, y, sigma, l_x, l_y, phi)).dot(K_inv).dot(y_con)\n",
    "    comp_4 = np.trace(K_inv.dot(K_phi(x, y, sigma, l_x, l_y, phi))) - \\\n",
    "             (y_con.T).dot(K_inv).dot(K_phi(x, y, sigma, l_x, l_y, phi)).dot(K_inv).dot(y_con)\n",
    "    return [comp_1, comp_2, comp_3, comp_4]\n",
    "\n",
    "def grad_nlml(par):\n",
    "    # sigma = par[0]\n",
    "    # l_x = par[1]\n",
    "    # l_y = par[2]\n",
    "    # phi = par[3]\n",
    "    M = np.diag([np.exp(par[0]), np.exp(par[1]), np.exp(par[2]), 1])\n",
    "    grad = np.array([grad_L(np.exp(par[0]), np.exp(par[1]), np.exp(par[2]), par[3], s)]).dot(M)\n",
    "    \n",
    "    # Scipy.minimize wants an array as the gradient, not a matrix. \n",
    "    # We therefore remove any one-dimensional entries of the grad by squeezing.\n",
    "    return np.squeeze(np.asarray(grad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the conjugate gradient method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nfeval = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import conjugate_gradient._minimize as mi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   1    0.132043    0.332153    0.154280    1.628256\n",
      "   2    0.133594    0.329288    0.147850    1.933769\n",
      "   3    0.133708    0.329171    0.147714    1.995426\n",
      "   4    0.133759    0.329136    0.147876    2.000100\n",
      "   5    0.133853    0.329068    0.148194    2.000475\n",
      "   6    0.295579    0.210490    0.698028    1.991637\n",
      "   7    0.620095   -0.027310    1.801335    1.997384\n",
      "   8    0.620123   -0.027223    1.801400    1.999851\n",
      "   9    0.852281    0.686640    2.325283    1.999854\n",
      "  10    0.852281    0.686640    2.325283    1.999854\n"
     ]
    }
   ],
   "source": [
    "t1 = time.time()\n",
    "m = mi.minimize(nlml, np.random.rand(4), jac=grad_nlml, method='CG', callback = callbackF)\n",
    "t_CG = time.time() - t1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plotting nlml close to the optimum in direction of the gradient**\n",
    "\n",
    "*Gradient-based optimization is difficult as we're dealing with numerical instabilities:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'$\\\\mathcal{NLL}(x_{opt} + \\\\alpha \\\\cdot \\\\nabla \\\\mathcal{NLL}(x_{opt}))$')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAERCAYAAABVU/GxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl83FW9//HXJ3vStOmkSdekaYHS0hba0ljWiwWVtnpZVNAWUVZxA/S6XbgoKojgggIqKkpl+V0KinipCJZ9b6Ep0J0udN+XdG+TNMnn98d8U8aSNDOTWbK8n4/HPGbmfLdPvknmM+d7zvccc3dERERilZHuAEREpGNSAhERkbgogYiISFyUQEREJC5KICIiEhclEBERiUunTyBmNtXMtpjZgijWHWhmL5jZ22Y2z8w+nooYRUQ6ok6fQID7gIlRrvs94C/uPgaYDNydrKBERDq6Tp9A3P1loDqyzMyONrN/mdkcM3vFzIY1rQ70CF4XARtSGKqISIeSle4A0uQe4MvuvszMTiJc0zgL+CHwtJldA3QDPpq+EEVE2rcul0DMrBA4FfirmTUV5wbPU4D73P12MzsFeNDMRrp7YxpCFRFp17pcAiF82W6nu49uZtkVBO0l7j7TzPKAEmBLCuMTEekQOn0byOHcfTew0swuBLCwUcHiNcBHgvLjgDxga1oCFRFp56yzj8ZrZtOA8YRrEpuBHwDPA78D+gHZwMPufpOZDQf+CBQSblD/rrs/nY64RUTau06fQEREJDm63CUsERFJjE7diF5SUuKDBg1KdxgiIh3KnDlztrl7aWvrdeoEMmjQIKqqqtIdhohIh2Jmq6NZT5ewREQkLkogIiISFyUQERGJixKIiIjERQlERETiogQiIiJxUQIREZG4dOr7QOL1wpItvL16R/jN+0O+0/Qqwwyz8HszsOB9VoaRmZFBVoaRkWHBeyMvO5PC3Ey65WTRLbfpkUmoIIfsTOVwEemYlECa8eqybUx9bSXJHiYsO9M4urSQY/t0Z2jf7gwNnstC+UTMVSIi0i516sEUKysrPdF3ors77tDojgPu4ITLGhqd+kYPnhvDzw1OzcEG9tU1sK+2nr219eyvq2dvbQMbdh5gyaY9LNm0h/U7Dxw6Ru/uuXx7wlAuOLGMjAwlEhFJLTOb4+6Vra2nGkiMmi5XZZDYD/Y9NQdZunkvSzbt4a9z1vLdR+fx4MzV/PDc4YytKE7osUREEkE1kHbI3Xn8nQ3c+tRiNu+u5dxR/blu0jD698xPd2gi0gVEWwNRC247ZGacP2YAz39rPNecdQz/WriJs25/kbueW0ZDY+dN+CLSsSiBtGPdcrP41tlDee6bH+bMob355TNLuf/1VekOS0QEUALpEMqLC7j7cydyxrGl/PKZpWzeXZPukERElEA6CjPj5vNGUNfQyE1PLEp3OCIiSiAdSUWvblx95jH8c95GXlq6Nd3hiEgXpwTSwXzpw0dxVEk3bnx8ATUHG9Idjoh0YUogHUxuViY/Pn8kq7fv5+4Xlqc7HBHpwpRAOqBTjynh/NH9+f1LK1ixdW+6wxGRLkoJpIO64RPDyc3O4PuPL6Az3wwqIu2XEkgHVdo9l+9OHMZry7czfe6GdIcjIl2QEkgHdtG4gYwqK+LmJxaz68DBdIcjIl1MShOImU01sy1mtqCF5ePNbJeZvRM8boxYNtHMlpjZcjO7LnVRt1+ZGcYtnzyebXtr+WvV2nSHIyJdTKprIPcBE1tZ5xV3Hx08bgIws0zgt8AkYDgwxcyGJzXSDmLkgCIG9SrgzZXV6Q5FRLqYlCYQd38ZiOeTbhyw3N1XuHsd8DBwXkKD68AqBxUzZ/UONaaLSEq1xzaQU8xsrpk9ZWYjgrIBQOQ1mnVB2QeY2VVmVmVmVVu3do27tSsrQmzfV8fKbfvSHYqIdCHtLYG8BVS4+yjg18D/xboDd7/H3SvdvbK0tDThAbZHlYPCE05VrdqR5khEpCtpVwnE3Xe7+97g9ZNAtpmVAOuB8ohVy4IyAY4u7UaoIJvZq9QOIiKp064SiJn1NTMLXo8jHN92YDYwxMwGm1kOMBmYnr5I2xczY2xFuB1ERCRVUjonuplNA8YDJWa2DvgBkA3g7r8HLgC+Ymb1wAFgsodbhuvN7GpgBpAJTHX3hamMvb370KAQzy7ezLa9tZQU5qY7HBHpAlKaQNx9SivLfwP8poVlTwJPJiOuzqByUAiAOat3MGFE3zRHIyJdQbu6hCXxGzmgiJysDKrUDiIiKaIE0knkZmUyqqyI2eqJJSIpogTSiVQOKmbhhl0cqNNEUyKSfEognUhlRYiDDc7cdTvTHYqIdAFKIJ3I2IpwQ7raQUQkFZRAOpGeBTkM6V1Ile4HEZEUUALpZJoGVmxs1MCKIpJcSiCdzIcGhdhTU8/SLXvSHYqIdHJKIJ1MZUV4YEV15xWRZFMC6WTKi/Pp3T1XDekiknRKIJ2MmVE5KKSh3UUk6ZRAOqHKimLW7zzAxl0H0h2KiHRiMScQM+sWzFEu7VTTwIqqhYhIMrWaQMwsw8wuMrN/mtkW4F1go5ktMrOfm9kxyQ9TYjG8Xw8KcjLVDiIiSRVNDeQF4GjgeqCvu5e7e2/gdGAW8FMzuziJMUqMsjIzGDOwp3piiUhSRTMfyEfd/eDhhe5eDfwN+JuZZSc8MmmTsRXF/Ob5ZeypOUj3PP16RCTxWq2BHJ48mmsDaS7BSHp9aFCIRoe312hgRRFJDrWBdFJjBobIMA2sKCLJozaQTqowN4ujSwtZvElDmohIcqgNpBMrLy5g3Q7dCyIiyRFzGwiAmX0r4vXQaNtAzGyqmW0xswUtLP+cmc0zs/lm9rqZjYpYtioof8fMqqI5XldXFspn3Y796Q5DRDqpaGogh5hZT+BXwDAzOwDMA64ALotyF/cBvwEeaGH5SuDD7r7DzCYB9wAnRSw/0923xRJzV1YeKmBPTT27DhykKF+VRBFJrJgSiLvvBC4LPty3AMcDj8Ww/ctmNugIy1+PeDsLKIslPvl3ZaF8ANbt2E9RflGaoxGRzibesbCuDB5nA3sTF86/uQJ4KuK9A0+b2Rwzu6qljczsKjOrMrOqrVu3Jim0jqEsVADA2mq1g4hI4sVUA4nwrrvfAGBmvyXcUythzOxMwgnk9Iji0919vZn1Bp4xs3fd/eXDt3X3ewhf+qKysrJLT8sXWQMREUm0eBPIRDOrBuYCCf10MrMTgD8Bk9x9e1O5u68PnreY2d+BccAHEoi8r2dBNoW5WeqJJSJJEe8lrInAYuA0YICZ3Z+IYMxsIOE2lc+7+9KI8m5m1r3pNeFLZ8325JL3mZl6YolI0sRVA3H3rcCTwSNqZjYNGA+UmNk64AdAdrDP3wM3Ar2Au80MoN7dK4E+wN+DsizgIXf/VzyxdzXhBKIaiIgkXlwJxMxuBoYB+4Db3X1+NNu5+5RWljc1zh9evgIY9cEtpDVloQJmrajG3QkSsIhIQsR7CSvf3S8ErgK+msB4JMHKQvnsrQ3fCyIikkjxJpA8MzvR3esAfa1tx9SVV0SSJd4E8h3gLDObCjyewHgkwcqL1ZVXRJIj6jYQM7sT+IaHHQB+kbywJFGaaiBqSBeRRIulBrIHmB50o8XMJpjZa8kJSxKlKD+b7nlZrFUNREQSLOoaiLt/z8wuAl40szrCQ5hcl7TIJGHKQhrWXUQSL5ZLWB8Bvki4624/4HJ3X5KswCRxykP5rNq+L91hiEgnE8slrBuA77v7eOAC4BEzOyspUUlCNdVA3Lv00GAikmDRzIluAO5+lru/GryeD0wCbo5cR9qnslA+++saqN5Xl+5QRKQTiWpOdDO7JhinKtJ24LZgHKxLEh+aJEp5sXpiiUjiRZNAJgINwDQz22Bmi8xsBbAMOBe4w93vS2KM0kbvD+uuBCIiidNqI7q71wB3Ex7gMBsoAQ4EsxNKBzAgSCDqyisiiRTTnejuftDdN7r7TjMblqygJLF65GVTlJ+tu9FFJKGiSiBm9h0zm2lmx0UUrzezLycpLkmw8mIN6y4iiRVtDeQY4BvAoUnG3X0PcE4ygpLEK+tZwNpq1UBEJHGiTSDPE56f/FA/UDMrITwjoXQATRNL6V4QEUmUqBKIuz8SrPuemc02s1uAUwHdid5BlBcXUFvfyLa9uhdERBIj6kZ0d/85MJDwNLSZwLcJD7AoHcD7XXl1GUtEEiOmKW2DYdwPzYVuZuOTEJMkwaGJpXYcYMzAUJqjEZHOIN4JpTCzYe7+YgJjkSRSDUREEi2l3XjNbKqZbTGzBS0sNzO7y8yWm9k8MzsxYtklZrYseGjolBh1y82iuFuOuvKKSMKkuhvvfYSHRmnJJGBI8LgK+B2AmRUTbns5CRgH/MDMdB0mRmWhfHXlFZGESWk3Xnd/Gag+wirnAQ8E0+bOAnqaWT9gAvCMu1e7+w7gGY6ciKQZZaF81qsGIiIJ0t668Q4A1ka8XxeUtVT+AWZ2lZlVmVnV1q1bm1ulyyoPFbBu5wEaG3UviIi0Xafrxuvu97h7pbtXlpaWpjucdqUslE9dfSPb9tamOxQR6QRa7cZrZl8H3nD3WSnoxrseKI94XxaUrQfGH1b+YoKP3em935V3P7175KU5GhHp6KKpgVwJbGxuQRK68U4HvhD0xjoZ2OXuG4EZwNlmFgoaz88OyiQG5cWaF0REEieaGwk/A0w0s2rg2aAROy5mNo1wTaLEzNYRvhyWDeDuvydcs/k4sBzYD1wWLKs2s5uB2cGubnL3IzXGSzMG9NTMhCKSONFMKLUYWGxmucB4M+sJ1AKz3H1TLAdz9ymtLHfgay0smwpMjeV48u/yczIpKcxRV14RSYiohzJx91qCy0Zm1g34u5lVufv/JCs4SbwBoQLVQEQkIWIaysTMupnZNcBbhJPJDUmJSpKmPJSv4UxEJCGiHcrkeDO7G1hG+K70j7j77a7JJTqcslAB63UviIgkQLSXsBqA+cBmYDcwwszWK4F0PGWhfA42OJv31NCvKD/d4YhIBxZVAnH3RcCipvfBoIrfM7OXguFJpIMoL36/J5YSiIi0RVzDuQc9sx4gPDiidCAa1l1EEiXu+UCA7sAvExWIpMaAnuEEsrZaPbFEpG2ibUS/2cwqIsvcfYG7/yY5YUmy5GVn0rt7rmogItJm0dZAvg48bWYzzOxCM4tpKlxpX8pC+boXRETaLNoEstndhwK3AecCy83s52Y2NHmhSbKUhQpYqxqIiLRRtAnEAdz9BXf/PHAC4bGqFpjZ8GQFJ8lRXpzPhp011Dc0pjsUEenAYr0TPcPMzgUeBD4LfB9YmYzAJHnKQwU0NDobd9WkOxQR6cCibcsoNLPbgMnA68Cd7v588sKSZGq6F2Ttjv2HXouIxCraBFJNeE6QEzWMesdXHkwsta76AByd5mBEpMOK9k70kckORFKnX888Mgw1pItIm7TlRkLpoLIzM+hXlK95QUSkTWJOIGY20MwsGcFI6pQX57NW94KISBvE2gsrH3gD6J2ccCRVykMFqoGISJvEdEe5ux8A+iUpFkmh8uICtuyppeZgA3nZmekOR0Q6ILWBdFHlxU2j8uoylojEJ+oEYmZj23owM5toZkvMbLmZXdfM8l+Z2TvBY6mZ7YxY1hCxbHpbY+nqmrryqieWiMQrlhrIJWb2qJmd3FRgZlEP525mmcBvgUnAcGDK4cOguPt/uftodx8N/Bp4LGLxgaZl7n5uDHFLMw5NLKV2EBGJUywJZAtQCTwW1CJWAn1j2H4csNzdV7h7HfAwcN4R1p8CTIth/xKD0sJccrIy1BNLROIWSwK5GBjq7v2BM4GXgTdj2H4AsDbi/bqg7AOCuUcGA5HDpeSZWZWZzTKz81s6iJldFaxXtXXr1hjC61oyMoyyUD5rtqsGIiLxiSWBrCX8oY67b3D3S4AvJSWq8Jhbj7p7Q0RZhbtXAhcBd5hZs4NwuPs97l7p7pWlpaVJCq9zKNew7iLSBrF04/068Dczewt4m3DtYV8M268HyiPelwVlzZkMfC2ywN3XB88rzOxFYAzwXgzHl8OUF+fz9pod6Q5DRDqoqGsg7r4IOJFw20U+sIkjt2EcbjYwxMwGm1kO4STxgd5UZjYMCAEzI8pCZpYbvC4BTgMWxXBsaUZ5qIDdNfXsOnAw3aGISAcU642EtcA/g0dM3L3ezK4GZgCZwFR3X2hmNwFV7t6UTCYDD7u7R2x+HPAHM2sknPRuCxKatMGhYd2r91M0oCjN0YhIR5PSuc3d/UngycPKbjzs/Q+b2e514PikBtcFHRrWfcd+RiqBiEiMdCd6F9Z0N/raanXlFZHYxTMa78eSEYikXlF+Nt1zs9QTS0TiEk8N5KcJj0LSwswoL9aovCISH13C6uI0L4iIxCuqRnQz+zPggAEDzWxq0zJ3vzxJsUkKlIcKeGnpVtwdzRMmIrGIthfWfRGvTwfuT3wokg7lxQXUHGxk695aenfPS3c4ItKBRJVA3P2lptdmtifyvXRskT2xlEBEJBbxtIHUJTwKSZvIe0FERGIRcwJx95NbX0s6irLQ+3eji4jEQr2wurj8nExKCnN1M6GIxEwJRIKuvKqBiEhs4kogZvatiNdDExeOpIPmBRGReMSUQMysZ3BPyAVm9lUzOx24LjmhSaqUF+ezYWcN9Q2N6Q5FRDqQWIdz3wlcZmaTCM+RfjzwWDICk9QpDxXQ0Ohs3FVzaIh3EZHWxNsGcmXwOBvYm7hwJB0OzQuiy1giEoN4E8i77v4Vd78IuCCRAUnqHboXRD2xRCQG8U4oNdHMqoG5gL62dnD9euaRYaqBiEhs4q2BTAQWA6cCA8xMY2N1YNmZGfQrytfNhCISk3hrINcQnqd8P3Cru89PXEiSDhrWXURiFW8NpMDdLwS+CHw1gfFImpSHNLGUiMQm3gSSZ2Ynunsd4TlComJmE81siZktN7MP3D9iZpea2VYzeyd4XBmx7BIzWxY8LokzbmlBeXEBW/bUUnOwId2hiEgHEW8C+Q5wVjCx1OPRbGBmmcBvgUnAcGCKmQ1vZtVH3H108PhTsG0x8APgJGAc8AMzC8UZuzSjaVh3jcorItGKOoGY2Z0WTFnn7gfc/Rfufrm7PxXlLsYBy919RVBzeRg4L8ptJwDPuHu1u+8AniHckC8JUn5oVF61g4hIdGKpgewBpptZNwAzm2Bmr8Ww/QBgbcT7dUHZ4T5tZvPM7FEzK49xW8zsKjOrMrOqrVu3xhBe16abCUUkVlEnEHf/HjANeDFIHN8k8eNg/QMY5O4nEK5lxNw92N3vcfdKd68sLS1NcHidV2lhLjlZGWpIF5GoxXIJ6yOEe13tA0qAa939lRiOtR4oj3hfFpQd4u7b3b02ePsnYGy020rbZGQYZaF8XcISkajFcgnrBuD77j6e8PAlj5jZWTFsPxsYYmaDzSwHmAxMj1zBzPpFvD2X8M2KADOAs80sFDSenx2USQJpWHcRiUXUNxK6+1kRr+cHI/L+jfDd6NFsX29mVxP+4M8Eprr7QjO7Cahy9+nAtWZ2LlAPVAOXBttWm9nNhJMQwE3uXh1t7BKd8uJ83l6zI91hiEgH0WoCMTNzdz+83N03Bpe1WlynmW2eBJ48rOzGiNfXA9e3sO1UYGprx5D4VRR3Y3dNPdv31tKrMDfd4YhIOxfNJawXzOwaMxsYWRhchjolGAdLN/Z1AqPKewLw9pqdaY5ERDqCaBLIRKABmGZmG8xskZmtAJYBU4A73P2+JMYoKXJCWRHZmcYcXcYSkSi0egnL3WuAu4G7zSybcA+sA8HshNKJ5GVnMqJ/EXNWK4GISOtiGsrE3Q+6+0Z332lmw5IVlKTP2IoQc9fupK5e86OLyJFFlUDM7DtmNtPMjosoXm9mX05SXJImYytC1NY3smjj7nSHErete2rZsFP3s4gkW7Q1kGOAbwCHxgZx9z3AOckIStJnbEV4jMqOehlr064azvn1q5z20+e54r7ZvLBkC42NrXYQFJE4RJtAngdOB+qaCsysBDgtGUFJ+vTpkUdZKJ+3OmAC2V9Xz5UPzGZPzUGuOG0wc9ft4rI/z2b8L17kDy+9x459da3vRESiFlUCcfdHgnXfM7PZZnYL4RsIlyQzOEmPsRUhqlZXE8WtPe1GY6PzX4+8w6INu/n1RWP43n8O5/XrzuLXU8bQtyiPW596l5NufY7rH5vPnpqD6Q5XpFOIZTDFnwMDCc/LkQl8m/AIvdLJjK0IsXl3Les7UDvCz2YsYcbCzdzwieGcNawPADlZGZwzqj9/+dIpzPjGGVw4toy/VK3lP3/9KvPWqROhSFvF2gvrgLs/6e7XufsZwI+TFJekUUdrB/lL1Vp+/9J7fO6kgVx+2qBm1xnatzu3fPJ4Hr7qZA7WN/Lp373Ova+u7FC1LJH2Jt4ZCTGzYe7+YgJjkXZiaJ/udMvJ7BAJZNaK7dzw9/n8x5ASfnjuCII5z1r0oUHFPPn1/2D80N7c/MQirry/Sm0jInFSN175gKzMDEYP7NnuE8jKbfv48v+bw8DiAn5z0YlkZ0b3fahnQQ73fH4sPzxnOK8s28akO1/hzZUam1MkVurGK80aW1HM4o272Vdbn+5QmlXf0MgXH6jCgKmXfoii/OyYtjczLj1tMI999VTysjOY8sdZvLRUM1iKxELdeKVZYytCNDrMXds+G5ufXrSZ5Vv28pNPHk9Fr25x72fkgCL+cc3pDOldyDUPvcXq7fsSGKVI56ZuvNKs0eU9MYOqdnoZ695XVzKwuICzR/Rt876652Xzh8+Pxcz40oNz2F/XPmtdIu1NqwnEzL5uZierG2/XUpSfzbG9u7fLdpB31u5kzuodXHbaIDIzjtxoHq2KXt24c/Jolmzew3cfnafeWSJRiKYGciWwEdSNt6sZOyjEW2t2tLuhQO59dSXdc7O4sLI8ofsdP7Q33z57KE/M28gfX1mR0H2LdEbRJJDPABPN7MJgPvJD1I23cxs7MMSemnqWbdmb7lAO2bDzAE/O38jkceUU5kY9I3PUvjr+aCaN7MttT73Lq8u2JXz/Ip1JqwnE3Re7+x+A6cA4M/usmZ1vZm2/+CztWnu8ofD+matwdy45dVBS9m9m/PzCURxdWsg1095ibfX+pBxHpDOIZSiTWnefETSoPwM8YGY/SV5okm4VvQooKcxpNwlkX209095Yw6SR/SgLFSTtOIW5WdzzhUrqG50vPTiHA3UNSTuWSEcW053oZtbNzK4B3gJmADfEuP1EM1tiZsvN7Lpmln8zmDJ3npk9Z2YVEcsazOyd4DE9luNKfMyMEweG20Hag7+9tY7dNfVcfvrgpB9rcEm4UX3Rxt3c/rQ6G4o0J9o70Y83s7sJz4N+DPARd7/dY+iqYmaZwG+BScBwYIqZDT9stbeBSnc/AXgU+FnEsgPuPjp4nBvtcaVtxlaEWLltH9v21qY1jsZG58+vrWJ0ec9Dl9aS7axhfZgybiB/fn0VCzfsSskxRTqSaGsgDcB84A/AKmCEtTbo0AeNA5a7+wp3rwMeBs6LXMHdX3D3povOs4CyGI8hCdb0YZ3u+UGef3cLK7ft44oU1D4iXTdxGKGCbP7n7wtoaGe90UTSLdobCRe5++/c/Ufu/itgDfA9MzsjhmMNANZGvF8XlLXkCuCpiPd5ZlZlZrPM7PyWNjKzq4L1qrZu1dAUbTVyQBE5mRnMSfNlrHtfXUn/ojwmjUxt342igmy+94nhzF27k4feWJ3SY4u0d3GNxuvui4EHgPsSGk3AzC4GKoGfRxRXuHslcBFwh5kd3UJs97h7pbtXlpaWJiO8LiUvO5ORA3owZ1X6EsjCDbuYuWI7l5w6iKwoB0xMpPNG9+f0Y0r42b+WsGV3TcqPL9JeteW/sTvwyxjWXw9E3vlVFpT9GzP7KOHG+XPd/dCFd3dfHzyvAF4ExsQessRjbEWIeet3UVufnt5IU19dRUFOJpPHDUzL8c2Mm88fSW1DIzc9sSgtMYi0R3EnEHdf4O6/iWGT2cAQMxtsZjnAZML3lhxiZmMIt7Oc6+5bIspDZpYbvG4axFH/ySkytiJEXX1jWoY837DzAP+Yu4ELx5bFPOJuIg0u6cbXxh/DE/M2atRekUDKrge4ez1wNeHuv4uBv7j7QjO7ycyaelX9HCgE/npYd93jgCozmwu8ANzm7kogKXLGsaX0K8rjln8upr6hMWXHdXe+938LyMwwrvyPo1J23JZ8efxRHFXaje//3wJqDureEJGUXlAOxtE61t2PdvdbgrIb3X168Pqj7t7n8O667v66ux/v7qOC53tTGXdXV5CTxQ/OGc67m/Zw3+urUnbc6XM38Py7W/jW2cdSXpy8GwejlZuVyY/PH8ma6v38+vll6Q5HJO1S3yIpHdKEEX0ZP7SUXz2zlE27kt+QXL2vjh/9YxGjynty2Wmp7bp7JKceXcKnThzAPS+vYNlmDUYtXZsSiETFzLjp3JHUNzo3p6Ah+aZ/LGRPzUF+9ukTEjZke6Lc8PHj6JabpQZ16fKUQCRqA3sV8LUzj+Gf85PbkPzCki383zsb+Mr4Yxjat3vSjhOvXoW5XH3mMbyybBuvL9eIvdJ1KYFITL704aMYXNKNHzyenIbkvbX13PDYfIb0LuRrZzZ7q0+7cPHJFfQvyuOn/3pXk09Jl6UEIjHJzcrk5vNGsmr7fn734nsJ3//P/vUuG3fXcNunTyA3KzPh+0+UvOxMvvGxY5m7bhczFm5KdzgiaaEEIjE7fUgJ54zqz+9eeo+V2/YlbL9Vq6p5cNZqLjllUMoGTGyLT40ZwDG9C/nZjCUp7d4s0l4ogUhcvv+J48jJzODGxxck5BLO9r21/Pff5tG/KJ/vTBiagAiTLyszg2+fPZQVW/fxt7fWpTsckZRL/Jyg0iX07pHHt84+lh/9YxHfeXQe548ewElHFZMdxVhVW/fUsmDDLhas2xV+Xr+b9TsPAHD/5ePoloSpapNlwog+jC7vyR3PLuO80QPIy26/l91EEq3j/KdKu/P5kyt4d+Meps/dwKNz1tGzIJuPHdeHScf35bRjSsjJzGDz7lrmr9/FgvW7WBgki00RAxLowAidAAAQIUlEQVQOLunGiRUhvnBKBScf1YtR5T3T+BPFzsz474nDmPLHWTwwcxVXnZG4hv85q6s5rl8PCnL0byrtk3XmHiSVlZVeVVWV7jA6vQN1Dby0dCszFm7i2cWb2VNTT/fcLHKzM9i2tw4AMzi6tJCR/XswckARIwcUMaJ/D7rnpW98q0T6wtQ3mbt2Jy9/98yEjNk1e1U1F/5+JteedQzfPLtjXNKTzsPM5gSjnx+RvtpIm+XnZDJxZF8mjuxLXX0jr723jacXbqa+oTFIFj06/Tfp704Yyn/++lX++PIKvt3GNpzGRufHwU2KTy3YpAQi7Vbn/Y+WtMjJyuDMob05c2jvdIeSUiMHFHHOqP7c++pKvnBKBb175MW9r3/M28DcdbuorAhRtXoHK7bu5ajSwgRGK5IY6oUlkiDf+tixHGxo5I7n4h9oseZgAz996l1GDujBHZNHAzBj4eZEhSiSUEogIgkyqKQbnz+lgofeWMMry+Ib6uXeV1eyYVcNN3x8OGWhAk4oK+JfulFR2iklEJEE+u6EYQzpXcg3/zKX7XtrW98gwtY9tdz9wnI+NrwPpxzdCwiPgjx37U427jqQjHBF2kQJRCSB8nMyuWvKGHYdOMh3Hp0X002Wv3p2KbX1jVw/adihsgkj+gLwtC5jSTukBCKSYMf168H/TBrG8+9u4YGZq6PaZsmmPTz85houPrni3xrMj+ldyNGl3brMeFt19Y386pmlcV8ClNRSAhFJgktOHcSZQ0u55cnFvLtpd6vr/+TJxRTmZvH1jwz5wLKJI/vyxspqduyrS0aoCeHu1NW3bTywvbX1XHH/bO58bhmX3zebp7tI0uzIlEBEksDM+PmFo+iRl821094+4tD3Ly3dyktLt3LtR4YQ6pbzgeUTRvSlodF5dvGRL2PVHGxg5/7UJ5kF63dx/m9f45Rbn+PVZfHNj7JlTw2f/cNMXn9vOz86dwQj+hfx1f99i38tUBJpz5RARJKkpDCXX35mFEs37+WWfy5udp0tu2v4yT8XM7C4gM+fUtHsOscPKKJ/Ud4RL2O5O198oIrKHz/L1x9+m7lrdybkZziSfbX1/PiJRZz7m1dZv7OGngXZfGHqG/z2heU0Nkbf9rNi614+dffrrNy2jz9dUsklpw7igSvGcXxZEVc/9BZPzd+YxJ9C2iLlNxKa2UTgTiAT+JO733bY8lzgAWAssB34rLuvCpZdD1wBNADXuvuMFIYuErMzji3li/8xmD++spIxA3sS6pbD/HW7mLduF/PX72Tz7nBPrd9ffGKL85+YGWeP6MtDb65hX219s4NNPvTmGl5Zto0zh5by3OItPP7OBsZWhLj8tMFMGNGHrCgGuYzFs4s2c+PjC9iwq4aLThrIf08YRlamcf1j8/n5jCW8vWYHt39mdKvDury1ZgdX3DebDDOmffHkQ2Oh9cjL5oHLx3HJ1De5etrb3OXwiRP6JfRnkLZL6VhYZpYJLAU+BqwDZgNT3H1RxDpfBU5w9y+b2WTgk+7+WTMbDkwDxgH9gWeBY929xWsDGgtL2oPa+gY+dffrLNwQbgsxg6NKunFCWU+OH1DEhwYVc3xZ0RH3MWvFdibfM4vfXnTiBz5I11bvZ+IdLzNmYIgHrxjH3tp6Hp2zjvteX8Xq7fvpX5TH5acP5tJTB0WVSBobnX/M20D1vjryszPJz8k89JydmcGfX1vJjIWbObZPIbd+6njGVhQf2tbdeWDmam5+YhH9e+bzu4tPZET/D/5sjcEluWsffps+PfK4/7JxDCrp9oH19tbWc+nUN3l77U7unDya/zyhf6vxS9tFOxZWqhPIKcAP3X1C8P56AHe/NWKdGcE6M80sC9gElALXRa4buV5Lx1MCkfZi/c4DPLNwE8P69YhrEMmGRudDtzzL6ceUcNeUMYfK3Z2L732Dd9bsZMZ/nUFZqODftnn+3S1MfXUlM1dsp7IixF1TxtC/Z36Lx9m+t5Zv/mXuEee8z83K4OsfHcKVpx9FTlbzCWnO6h187X/fYsf+Or519rHkZGawuno/a7bvZ3X1ftZW76e2vpETyoqYeumHKCnMbfF4e2vruezPb/LWmp184ZQK8rMzaWh06hud+oZG6hudpitmGRZO0IYFz13X5HEDOa5fj7i2ba+DKQ4A1ka8Xwec1NI67l5vZruAXkH5rMO2HXD4AczsKuAqgIEDByYscJG2GNAzn0tPGxz39pkZxseO68OT8zdSW99w6HLXtDfX8try7dzyyZH/ljwObTO8Dx8b3ofH31nP/zw2n4/f9Qq3XziKjxzX5wPHmLViO19/+G127D/IzeeP5JwT+rG/roEDBxs4UNfA/roG9tfVc2yf7kdMQgBjK0I8ce3pXDvtbX7y5LsA5GdnUtGrgKNKunHm0FIGlxRy/pj+rQ6yWZibxX2XjeMr//sW972+iqwMIzPDyMrIIDPDyM40zIzwd+FwMnF3HHAPJ5SuaPzQ3nEnkGh1usEU3f0e4B4I10DSHI5IwkwY2YdHqtby+nvbOXNob9bt2M8t/1zEqUf34qJxR/6ydN7oAZxQ1pOv/e9bXHF/FVeePpjvThxGTlYGDY3Or59fxl3PLWNQr278+dJxDO8f/uDpWXDE3R5RSWEuD15xEks27aGkew6lhblYnJ/m3XKzeODycfEHI0mR6gSyHiiPeF8WlDW3zrrgElYR4cb0aLYV6bROPbqEbjmZPL1wE+OPLeX6x+bjwE8/fUJUH8yDS7rx2FdP5dYnF/OnV1cye1U1N54znF/MWMrMFdv51JgB3Hz+yITOCJmZYYeSkXQ+qU4gs4EhZjaY8If/ZOCiw9aZDlwCzAQuAJ53dzez6cBDZvZLwo3oQ4A3Uxa5SJrlZWdy5rDePL1wMyMHFPHKsm3cfP5IyoujrybkZWfyo/NGcvJRvfju3+bx6d/NJD87k19cOIoLxpYlMXrpjFKaQII2jauBGYS78U5194VmdhNQ5e7TgXuBB81sOVBNOMkQrPcXYBFQD3ztSD2wRDqjCSP68sS8jdz4+EJOPboXn2vl0lVLJh3fj5EDirj31ZV87qSBDOnTPcGRSlegKW1FOpC9tfWceNMzZGUaM75xRky1D5FotddeWCLSBoW5WXz/nOH06Z6r5CFppwQi0sF8/uTmhzwRSTWNhSUiInFRAhERkbgogYiISFyUQEREJC5KICIiEhclEBERiYsSiIiIxEUJRERE4tKphzIxs63A6jg3LwG2JTCcRFFcsVFcsVFcsWuvsbUlrgp3L21tpU6dQNrCzKqiGQsm1RRXbBRXbBRX7NprbKmIS5ewREQkLkogIiISFyWQlt2T7gBaoLhio7hio7hi115jS3pcagMREZG4qAYiIiJxUQIREZG4dOkEYmYXmtlCM2s0sxa7u5nZRDNbYmbLzey6iPLBZvZGUP6ImeUkKK5iM3vGzJYFz6Fm1jnTzN6JeNSY2fnBsvvMbGXEstGpiitYryHi2NMjytN5vkab2czg9z3PzD4bsSyh56ulv5eI5bnBz788OB+DIpZdH5QvMbMJbYkjjri+aWaLgvPznJlVRCxr9neaorguNbOtEce/MmLZJcHvfZmZXZLiuH4VEdNSM9sZsSyZ52uqmW0xswUtLDczuyuIe56ZnRixLLHny9277AM4DhgKvAhUtrBOJvAecBSQA8wFhgfL/gJMDl7/HvhKguL6GXBd8Po64KetrF8MVAMFwfv7gAuScL6iigvY20J52s4XcCwwJHjdH9gI9Ez0+TrS30vEOl8Ffh+8ngw8ErweHqyfCwwO9pOZwrjOjPgb+kpTXEf6naYorkuB3zSzbTGwIngOBa9DqYrrsPWvAaYm+3wF+z4DOBFY0MLyjwNPAQacDLyRrPPVpWsg7r7Y3Ze0sto4YLm7r3D3OuBh4DwzM+As4NFgvfuB8xMU2nnB/qLd7wXAU+6+P0HHb0mscR2S7vPl7kvdfVnwegOwBWj1Tts4NPv3coR4HwU+Epyf84CH3b3W3VcCy4P9pSQud38h4m9oFlCWoGO3Ka4jmAA84+7V7r4DeAaYmKa4pgDTEnTsI3L3lwl/YWzJecADHjYL6Glm/UjC+erSCSRKA4C1Ee/XBWW9gJ3uXn9YeSL0cfeNwetNQJ9W1p/MB/94bwmqr78ys9wUx5VnZlVmNqvpshrt6HyZ2TjC3yrfiyhO1Plq6e+l2XWC87GL8PmJZttkxhXpCsLfYps09ztNZVyfDn4/j5pZeYzbJjMugkt9g4HnI4qTdb6i0VLsCT9fWW3ZuCMws2eBvs0susHdH091PE2OFFfkG3d3M2uxr3XwzeJ4YEZE8fWEP0hzCPcF/2/gphTGVeHu683sKOB5M5tP+EMybgk+Xw8Cl7h7Y1Ac9/nqjMzsYqAS+HBE8Qd+p+7+XvN7SLh/ANPcvdbMvkS49nZWio4djcnAo+7eEFGWzvOVMp0+gbj7R9u4i/VAecT7sqBsO+GqYVbwLbKpvM1xmdlmM+vn7huDD7wtR9jVZ4C/u/vBiH03fRuvNbM/A99OZVzuvj54XmFmLwJjgL+R5vNlZj2AfxL+8jArYt9xn69mtPT30tw668wsCygi/PcUzbbJjAsz+yjhpPxhd69tKm/hd5qID8RW43L37RFv/0S4zatp2/GHbftiAmKKKq4Ik4GvRRYk8XxFo6XYE36+dAmrdbOBIRbuQZRD+I9luodbpV4g3P4AcAmQqBrN9GB/0ez3A9degw/RpnaH84Fme2skIy4zCzVdAjKzEuA0YFG6z1fwu/s74WvDjx62LJHnq9m/lyPEewHwfHB+pgOTLdxLazAwBHizDbHEFJeZjQH+AJzr7lsiypv9naYwrn4Rb88FFgevZwBnB/GFgLP595p4UuMKYhtGuEF6ZkRZMs9XNKYDXwh6Y50M7Aq+JCX+fCW6h0BHegCfJHwdsBbYDMwIyvsDT0as93FgKeFvEDdElB9F+B98OfBXIDdBcfUCngOWAc8CxUF5JfCniPUGEf5WkXHY9s8D8wl/EP4/oDBVcQGnBseeGzxf0R7OF3AxcBB4J+IxOhnnq7m/F8KXxM4NXucFP//y4HwcFbHtDcF2S4BJCf57by2uZ4P/g6bzM72132mK4roVWBgc/wVgWMS2lwfncTlwWSrjCt7/ELjtsO2Sfb6mEe5FeJDw59cVwJeBLwfLDfhtEPd8InqYJvp8aSgTERGJiy5hiYhIXJRAREQkLkogIiISFyUQERGJixKIiIjERQlERETiogQiIiJx6fRDmYi0J2Y2ArgTGEh4TK7ehO+On53WwETioBsJRVLEzPKAt4ALCc/F8C4wx90/ldbAROKkGohI6nwUeNvdF8Kh8bluT29IIvFTG4hI6owG3gYws/6EZ617Lb0hicRPCUQkdep4fwKfWwnPPyLSYSmBiKTOQ8AZZraE8EitM83sjjTHJBI3NaKLiEhcVAMREZG4KIGIiEhclEBERCQuSiAiIhIXJRAREYmLEoiIiMRFCUREROLy/wEdi/acqP4L/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lin = np.linspace(-1, 1)\n",
    "\n",
    "def plotting(alpha):\n",
    "    return nlml(np.array(m.x) - alpha*np.array(m.jac))\n",
    "\n",
    "out = [plotting(l) for l in lin]\n",
    "\n",
    "plt.plot(lin, out)\n",
    "plt.xlabel(r'$\\alpha$')\n",
    "plt.ylabel(r'$\\mathcal{NLL}(x_{opt} + \\alpha \\cdot \\nabla \\mathcal{NLL}(x_{opt}))$')\n",
    "#plt.savefig('Long_range.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'$\\\\alpha$')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXmcHHd55/9++pyes+fQMTOSLMmSLMvGGFs+CMHGYLAhCQ5HgklCyOlwZcnxSwLLZsmGlzfXJvwISw4n61yQOMTAoiQGY3OGgI1kI9s6LGssWZqRRnNoru7p6fu7f1RVT01PH9V3T/f3/Xr1Sz3VR33V1V1PPdfnEaUUGo1Go2lfXI1egEaj0WgaizYEGo1G0+ZoQ6DRaDRtjjYEGo1G0+ZoQ6DRaDRtjjYEGo1G0+ZoQ6DRaDRtjjYEGo1G0+ZoQ6DRaDRtjqfRC3DC0NCQ2rlzZ6OXodFoNBuKp556alYptanY8zaEIdi5cydHjhxp9DI0Go1mQyEi55w8T4eGNBqNps3RhkCj0WjaHG0INBqNps3RhkCj0WjaHG0INBqNps3RhkCj0WjaHG0INBqNps3RhkCjyUE6rfjs4XESqXSjl6LR1JyKDIGI/JiIHBeRtIgczHrswyIyJiKnROQu2/a7zW1jIvKhSvav0dSKI+fm+c3PPcu3x2YbvRSNpuZU6hEcA94KfMu+UUQOAPcC1wB3A38mIm4RcQOfAt4IHADeaT5Xo2kqppaiACytJBq8Eo2m9lQkMaGUOgkgItkP3QM8pJSKAWdFZAy42XxsTCl1xnzdQ+ZzT1SyDo2m2syGYwAsx1INXolGU3tqlSMYBcZtf0+Y2/JtX4eI3CciR0TkyMzMTI2WqdHkxjIEkXiywSvRaGpPUY9ARB4HtuZ46CNKqS9Wf0kGSqkHgAcADh48qGq1H40mF7OhOADhmDYEmtanqCFQSt1ZxvteALbb/t5mbqPAdo2maVgNDWlDoGl9ahUaOgTcKyJ+EdkF7AW+BxwG9orILhHxYSSUD9VoDRpN2ViGIKxzBJo2oKJksYi8BfgksAn4dxE5qpS6Syl1XEQ+i5EETgLvV0qlzNd8AHgUcAMPKqWOV/Q/0GhqwGzYCA1pj0DTDlRaNfQF4At5HrsfuD/H9keARyrZr0ZTS5RSzOjQkKaN0J3FGk0WoViSeNLoKNbJYk07oA2BRpPFbCiWub+sy0c1bYA2BBpNFlZ+oL/TqxvKNG2BNgQaTRZWxdAVg106R6BpC7Qh0GiysAzBzsFObQg0bYE2BBpNFrOhGC6Bbf2dLMdTpNO6sV3T2mhDoNFkMROOMdDlozdgVFdHEjpPoGlttCHQaLKYCcUZ6vbT5TcMgQ4PaVodbQg0mixmwzGGuv10m4ZA9xJoWh1tCDSaLAxD4KPTpz0CTXugDYFGY0MpxWw4xqYeP11+N6A9Ak3row2BRmNjOZ4imkivCQ3ppjJNq6MNgUZjw5KXsCeL9ZQyTaujDYFGY8NqJhvq0cliTeOJJlIoVfs+Fm0INBobGUPQ7dPlo5qG8/7PPM2Pfuo/a74fbQg0GhszpuDcpm4/nV4rWaxzBJrGMBOOEez01Xw/2hBoNDZmQzFEYKDLh8sldPnc2iPQNIyZkFHBVmu0IdBobMyGY/R3+vC4jZ9Gl9+jDYGmIdhLmWuNNgQajQ2rmcyiy+/RyWJNQ1hcSZBIKYa6m9wQiMgficjzIvKsiHxBRIK2xz4sImMickpE7rJtv9vcNiYiH6pk/xpNtZkJxdb88Lr8OjSkaQwzZinzRvAIHgOuVUpdB7wAfBhARA4A9wLXAHcDfyYibhFxA58C3ggcAN5pPlejaQpmw/G1hsDn0Q1lmoYwY1awbWp2j0Ap9RWllHW59ASwzbx/D/CQUiqmlDoLjAE3m7cxpdQZpVQceMh8rkbTFGTHZLv9Hj23WNMQVj2CjVU19HPAl8z7o8C47bEJc1u+7esQkftE5IiIHJmZmaniMjWa3ETiSSLxVFZoSCeLNY1hNlPK3FHzfXmKPUFEHge25njoI0qpL5rP+QiQBD5TrYUppR4AHgA4ePCgHhGlqTmzIeOHtz5ZrENDmvozE4rhc7syA5JqSdE9KKXuLPS4iPwM8MPA69RqL/QFYLvtadvMbRTYrtE0lBmbvIRFt04WaxqEUbjgQ0Rqvq9Kq4buBn4TeLNSKmJ76BBwr4j4RWQXsBf4HnAY2Csiu0TEh5FQPlTJGjSaajGbIznX5fewkkiR0nOLNXWmXj0E4MAjKML/BvzAY6bVekIp9R6l1HER+SxwAiNk9H6lVApARD4APAq4gQeVUscrXINGUxVWdYbWJosBluNJeju8DVmXpj2ZCcUY7qt9fgAqNARKqT0FHrsfuD/H9keARyrZr0ZTC6wcwaAtR2CfUqYNgaaezIRjXLetry770p3Fmrpx/OIiT52bq+k+vnj0At87W94+ZsMxgp1evO7Vn4U1pUznCTT1JJVWzC3HN0xoSKNxzH/7v8dIpNL82y+/umb7+MMvn+KqrT3cvGug5NdaQ+vtrM4k0JVDmvoxH4mTStdHXgK0R6CpE9FEimMXFonU+IS6FE1wcWGlrNdm6wwBeiaBpiFkChfq5BG0rSFQSpHWlSB145nxBRIpxUqidoYgnVaEY8kKDEE8r0egDYGmnsyE1hcu1JK2NQRv/4vvcv8jJxu9jLbhyLl5gJoagnA8iVKwFE0SiiZKfn0u7fcuW9WQRlMv6ik4B21qCFJpxbMTC3z+6QmSqXSjl9MWHHnJSOCuxGtnCELR1ZP15GK0pNdGEynCseS6KzArWaxzBJp6okNDdWA6FCWRUsxHEmVXmGick04rnjI9glgyXbOQnN0LuFBieChzBaZDQ5omYCYUo8Prosvnrsv+2tIQjM+tniS+dOxSA1fSHozNhFmKJtm7uRuAaLI2V9d2j6DUPEGmmSxL6THgdeMSbQg09cUKU9ZDXgLa1BBMzBtqGPu39vDo8Us6aVxjjrxkeAOv3rsJqF14yO4RlG4ILMG5tR6BiNDl01PKNPVlNhyvyxwCi7Y0BJZH8HM/uIvpUIzvj883eEWtzZGX5hjq9rF/aw9Qu4Sx5RGIwMWF0nIEueQlLDq18JymzmRPyqs1bWkIJuYjbOn1c/e1W/G6hS89p8NDteTIuXluvKKfgBnvjNbIECyZhmDnYBcX5kv0CMwcwWD3+iEgxkwCnSzW1I96Cs5BmxqC8fkI2/s76e3w8oN7hvjy8UusKmhrqsn0UpTzcxFu2jlAwGsYgpV4bSq1rNDQVVt6Sk4Wz4Zj9HZ48HvWJ+e69QB7TQmcnV3mfZ95quwLnkQqzVykfvIS0KaGYGJ+hW39AQDeeO0wE/MrHL+41OBVtSZW/4DdI6hlaMjjEnZv6uLSUrQk6ejZcHzNHAI7XT4PEd1HoHHIl45N8shzlxibDpf1+rnlOErVr5kM2tAQJFNpJhejbOvvBODOA1twu4QvHZts8MpakyMvzeP3uLhmpI8Ob60NQYKeDg+j/QFSacV0yHmeYCaHzpBFu00pW4jEtYdcASfMi0or71Qq9W4mgzY0BJOLxpXi9gHDIxjo8nHLrgG+rMtIa8JT5+Z4+fYgPo+LDq/xdatd1VCSng4vI0Hj2JZSOTQbjuWt0minKWWXwzFu+Z9f5RGdNyubE5OGIZhbjpf1+pkChQu1ou0MwYSZRLQ8AoA3XruVF2eWOT0VatSyWpJIPMnxi0vctLMfIJMjqFWy2DAEHkZNQ3ChhMqh2dB6wTmLdhpg/+LMMrFkmucv6VBpOUTiSc7OLgNwOVymITA9gs3aI6gd42YPwXabIXjDNVsB3VxWbY6OL5BMKw5eYUhC1z5HYISGrKlOTj2CWDLFUjSZ1xVvp2Tx+Jzx+yi1/FZjcOpSCCuqNrtcXmioUClzrWg7QzAxv4JLYDi4OgJuS28HN17Rr8NDVeYps5Hshh1rPYJah4Z6Orz0dngcG4LLeZrJLLr8HmLJdFvoUp03DcHkYnkKru2OFRbyuV0VeQTdfk/mwqketJ8hmIsw3BdYM4UKjPDQicklzl+ONGhlrceRc/Ps29JNX6cx4rH2yWIjNAQwEgw4NgTFJH87fdaUstZPGI9nDIH2CMrhxMUlejo87N3SzeUKksX1TBRDhYZARD4mIs+KyFER+YqIjJjbRUT+VETGzMdvsL3m3SJy2ry9u9L/QKlMzK8wapaO2rnLDA99+biuHqoGqbTi6fPzHNy5OinM73EhUsuGskRmrvBoMOA4R7CqM5Q/NASGzHU5jM9F+M6Ls2W9tt6cz4SGVnTlUBmcmFzi6uFehrr9ZSeLCxUu1IpKPYI/Ukpdp5S6Hvg34L+b298I7DVv9wF/DiAiA8BHgVuAm4GPikh/hWsoifH5SKaHwM72gU6uHe3VeYIq8cJUiFA0ycErVg+viBDwumtiCKyhNHaP4MK8M+9uNSabP1kM5QvP/cljL/DeTz9d1mvrzfm5CCKGSmy5J7J2JZVWnLoU4sBwL4Pdvox+VanMhGLrxA9rTUWGQCllLy3oAqxLiHuAv1cGTwBBERkG7gIeU0rNKaXmgceAuytZQynEk2kuLUXXJIrtvPHaYb5/fqEt46OPn5ji1z/7TNXez2oksxLFFgGvuyahoWVzKI3dEDgdUJNPcM6iUinq4xcXWVxJ1MwTqhYr8RTToRgHhnsBHR4qlXOXl4nEUxwYMTyCy8uxsryqegvOQRVyBCJyv4iMAz/JqkcwCozbnjZhbsu3Pdf73iciR0TkyMzMTKXLBIwEmFLk9AhgNTz0leNTVdnfRuKxE1N87umJsptgsnnqpTk29/gz/RoWHV53TSQmLMG5His0ZB5jJyezmVCMHr8nk8PIZtUjKP1EHk2keHHGLCds8itsS5X31t2DQOkKru2OlSg+MNzLQJePaCJNpMTCiFgyxeJKovlyBCLyuIgcy3G7B0Ap9RGl1HbgM8AHqrUwpdQDSqmDSqmDmzZtqsp7Wqqj2wdyewR7Nnezd3M3jx5vv/CQ1YVbLamNI+fmObizf52eesDnPDQUiSf5k8deIOZgfsGqITBO2qNmVZgTzaHZcCxvfgDsU8pK9whemAplpC7KTR7WCys/cMsuw4vTHkFpnJxcwuMS9m7pZrDLCO2UGl4r5p3WiqKGQCl1p1Lq2hy3L2Y99TPA28z7F4Dttse2mdvyba8L1hVPPo8A4ODOfl5ow8ayabNy5tiFxYrf69JilIn5FW7MCgtBaaGhJ85c5k+/ejoz3awQVgjI8ghK6S6eDedvJoPKQkMnJ1cNa7W8rVphGYLrtwfxuV1cbMMQaSWcuLjEns3d+D3uzIm81GM+2wB5Cai8amiv7c97gOfN+4eAnzarh24FFpVSk8CjwBtEpN9MEr/B3FYXxucjeFzC1t6OvM8Z7gswG447ugptJaaWjC/giSp4BEfOGeM/7Ylii4DX7biPwArFLESKx/mzPYLNPR24XeLQEMQLXoFVMsDe/nnOhpo7NHR+LkLA62ZTj5+tfR26qaxETkwuZfIrlpx5qb0ExUqZa0WlOYLfN8NEz2Kc1D9obn8EOAOMAX8FvA9AKTUHfAw4bN5+19xWFybmVxgOduBx5/9vW12pl9rILU6m0lw2uyCPXazcIzjy0jwBr5sDI73rHuvwOfcILMVPJ+71kukR9JqGwG0afCcns9kCgnNgKx8twyM4MbnENebnUG6nab0Yn1thx0AnIsJIsINJnSNwzOVwjKmlWOY7P2CGhi6XeMxn6jy03sJTyYuVUm/Ls10B78/z2IPAg5Xst1zG5yJsC+bOD1ishhSiXDHYVY9lNZzLpuztSF8H5y5H1tTjl8P3zs5x/fbguqY9gIDXxfSSU0NgeQTFDUF2shisXoLCJ7NEKs1CpHByzu9xlTW3OJ1WnJwM8dYbRnlpdtmxR6CU4l+fneQNB7bkTWDXgvG5SCZ/NtIX4MmzdbtG2/CcnDTCyVdbHkGX8X0qtUCg0ICkWtJWncUT8yvrqliyyXgES+1zNTRthoXu2L8ZqCw89PVT05yYXOLOA1tyPl5KjsAyBPNlhIYARoIdRUNDxeQlwJxbXMaUsvH5COFYkgPDvQz1+B3Hi09Nhfgv//R9HvzPsyXtrxKUUpyfi7DDNATDwY6SZzq0MycmDU/aMgQBn5sun7v00FA4Rl/Am3NAUi1pG0MQTRg10tvy9BBYDPetegTtwtSS8X99rWkIyk0Yx5NpPvavJ9g91MW7br0i53MCvlJyBMbJfd6RR5DA7ZKMnhEY3t2lxcIns2LNZBblCM9ZBvXq4V4Gu3yOwwST5nfvc09N1K27dzYcZyWRYod5oTTcZ8x0sGLWmsKcuLjEcF9HJiQEMNjtL7lSrN4jKi3axhBYIYJiHkHA5ybY6W2rpjKrYujq4V629PrL9gj+9jtnOTO7zG//yAF8ntxfrY4yPAKnyeKeDs+actWRYIBkkZPZTBF5CYtypKhPTi7hErhqaw9D3X7HoSGrlPfFmWWOji+UtM9ysSqGdgyaoaESym81axPFFoPdvpJDQzOh+stLQBsZAktMq5hHAMbV0GQbeQTWiWeo28+1I31lJYynl6J84vHTvG7/Zu64anPe55UiMVFKstiSoLazOpcgv9REZhpUkR9fl9/DconNQScml7hyUzcdXndJoSGrgsvvcfG5pydK2me5WL8PKzRk5cpqcUG0HEvye4+cbBlpb6tp8OpsQ9DlK6tqqNhFSS1oG0NgDaTJJy9hZ6Svg4ttVDU0HYox2OXD53FxzUgvY9PhkqWi/+DLp0ikFL/9wwcKPi/gdZNIKRIOJJ1LTRb3+NcmuEccDKg5cXEJv8fF5t7CP75yppSduLiUqSIZ6vIxF4k7irlPh6L0d3q5+9qt/Oszk3WRpjifdaFkhUhrcUF0+KU5/vJbZ3jk2dYQeDw9FSaVVuuq5Aa7/CVXDTVCXgLazBD43C5HU3+Ggx1caqfQ0FI0E5e8ZrSPtKKkCVVPn5/nc09P8POv3sXOocKVVpbGupOTW6nJ4myPwApvFEoYf+uFGV555WDR5FyXr7TQ0PxynIuL0Uy4YKjHj1LOvJuppRibezp42w3bWFxJ8NWT0473Wy7n5yJs6fVnqpR6Ozx0+dw1aSqzPIFvnq6OdEyjsRLFOUNDYefznyPxJOFYsu6Cc9BGhmB8PsJofwCXS4o+d7gvwHwkUbMBKs3GdCjGZrPJzqp5dyo1kU4rfufQcbb0+vnAHXuKPr+UmQTWiXcpmih6Jb0UTawpHQXMITX5B9ScvxzhzOwyt+8rLmFSarLY6ijOeAQldJpOL0XZ3OvnVXuG2NrbUZfwkL1iCIxKqeFgbUKk1nH99unZlqhKOnFxiS6fe83nB0ayOJlWLK04+95YOSTtEdSQifmVgtISdqwS0nZJGE8vxTKe0mgwQLDTy3GHeYKHn5rg2YlFPvzGqzMduIXIzC12IDxnGQulYHGlsFcQiiYzzWR2RgsMqPnGC8aV9msK5DQsSk0WWwJkq3XlzjtNp0MxtvQandFvuWGUb74wk8nj1Ap7D4HFcF9HTTwCq9R3cSXBMxP1SYbXkpOTIfYP9667yLSOudNGwkY1k0E7GYK53HMIcpGJj7ZBniCdVsyEY2wxY+QiwjUjvRy7UNwjWFxJ8Adffp4br+jnnutHHO2vlLnFy7EkVhFQsRLSXMliKDyg5hunZrhisJNdRcJZAJ1+d0l9BCcuLrGl15/xBKwEYDGPIG1WOVmG+W03bCOVVnzx+xcd77tUYskUl5ai665oDSNa/d+A5VmJGKG5jUw6rXJWDMFqU5hT4blGyUtAmxiCSDzJ5eW4o4ohcBZbbhUuLxsJzM09q/pL1470cepSqGhC9xOPn2YuEud/vPmadSqj+QiUEBqKxFOZE+J8gR+TUtZQmvXd0PlGVkYTKb7z4iyvcRAWAuj2eYin0sSTziS0rUlVFk5DQ3OROMm0YosZqtuzuZvrtwd5uIY9BRfmDXn27EIKQ3crVnXdrXA0SZfPzcu3BfnmBjcEE/MrRtNgDjmVTHexw2oxyyNwksesNm1hCKyKIacegfUjbAe9ISvkYP/yXTPaRzyV5vRUOO/rxuci/N13X+Lem3Zw7Wif4/11lDDAPhJPZUpACyWMl+Mp0oqcHsFIMMDiSmJdfP97Z+eIJtKOwkJQ2pSyWDLF2HR4zVVib4cHn9tVdGqV1dxnPx5vu3Ebp6ZCVZMIzya7h8Bi2LwgmlqsblNZOJaky+/htn2beGZ8gUUHxQDNin0GQTZWk6LTSWWzoRgirGlKqxdtYgiML3q+OQTZdHjdDHb52qKE1Gom22xTZF1NGOfPEzx0+DxKKT7w2uIJYjulVQ0lGTWvUguFhrIlqO1Y3l22gNo3Ts3g87gyQ1iKUYrw3OmpMMmsckIRMccXFj6p5joeb75uBJ/bxcNP1SZpnN1DYDFiddlXOU8QiiXp7vBw+74h0gq+PbYx5jnn4oStaTCb/hLyQmB4BAOdvoKimLWiLQyBNZDGqUcAxtVQOySLZ5bWu6O7Brvo8rnzXoHGk2n++fAEr92/OXPF7hSnoaF4Mk0ipTLvX6iXIJfOkMVqU9naY/nNF6a5ZddAxjAVw/IInEycyneVONRdvKls2vQIttj6Gvo6vbz+wBYOPXPRcWiqFM7PRfB7XOuqVYZrFCJdjiXp8Xt4+bYgvR0evvlC7ctja8WJi0vsNpsGs/G6XfQFvMw5TBbPhhojLwFtYggm5nN/0QvRLt3FVijC/gV0uYSrh3vzegSPnZhiNhzjJ2/JrSdUiA6v8ZUrFhqyHt/U48frloKhoVWPIHdoCNZqR43PRXhxZtlxWAhKm1J24uISnT73OvVaq668EJYAYPYJ4W03jjK3HOfrp6p/0jxvVgxlV72M1KhoIhw1PAKP28UP7h3iWy/M1k1TqdqczJMothjs9jHrNFncIJ0haBtDYJSOOk1ogtVd3PoewXTIUDvMvqK5drSP4xeXSOeo8/7Mk+cYDQa4zWGi1Y5Tj8AaAtPlcxPs9BVMFi/lkKC22Nzjx+2SNTIT3zATlK+5yvn6S5lSdmJyif1be3BnnVideARTZldxdoPbbXs3MdTt53M1CA+dN+cQZBPwuenv9FbdIwjHkpnP8/Z9m7i0FOX0dP58VLOyGElwYWElZ6LYYqjLufDcTKjwXIxa0haGYHw+4rhiyGI4GCAUTbaMHko+pkPRNWEIiwMjvUTiKc5eXl6z/cWZMN958TI/ccuOdSc6J3Q4zBFYIRjrZFQ4R2Aco1x9BB63a92Amm+emmb7QIDdDspGLZwmi5VSxlVijpPDULe/aKep1VWcjcft4i2vGOFrz09XdfaxUorxrGYyO8N9gap7BKFoMvN5WhcT3zy18aqHCiWKLQYc6g0ppRqmPAptYgiczCHIZnVSWWt7BflOPNeOGJVA2XmCf3ryPB6X8GMHt5W1v4DDqqFIxiPw0N/pcxgayj1MZyTYkckRxJIpvvPiZV6zb3NJHqLTZPHE/AqhaJIDw+srqYa6fcRT6YwHkwujyzv3yeBtN24jmVb86zPV6ymYjxgVVfkKKZzMdCiVsJkjAMPQ7N3czbc2oNxEdtNgLga7fY76CMKxJNFEuqgceq1oeUMQiiZYiCRK9wjaZC6BvXnJzt4t3fjcLo7bZhNEEykefnqCu67ZmtN4OMHrduFxCdEitelW81an301/p89Rsrg3kLuz2d5LcPjsPJF4qqSwEECn6ckU8whOZElL2HHSSzC9FM372e7f2svmHn9Vy0jP56kYshjuy9+ZXQ5KKZbNqiGL2/dt4smzcxtK0kUpwyDvGOgseBU/2O13JDY406Ch9RZVMQQi8usiokRkyPxbRORPRWRMRJ4VkRtsz323iJw2b++uxv4LUYrqqJ12kJlQyuhi3ZTjCtTrdnHV1p41J51HnptkIZLgJ2/ZUdF+jQH2hatfVhLGCbfT56G/y1vUI8geSmPHPqDmmy9M43O7eOWVzspGLVYH2Bc+WZ24aJYTbllfTpgxBHnmI1hdxblCdRaj/YGq5q6KGYKRYIClaLJk5dV8xJJpkmlFt00p9rZ9m4gn0zxx9nJV9lEPvnpymqPjC7zvNVcWfN5Qtw+linfGW70Gm7rLu8CqlIoNgYhsxxhcf962+Y3AXvN2H/Dn5nMHgI8CtwA3Ax8Vkf5K11CI1TkEpYWGtvZ1IFKZR/CdF2ebullmIZIgnkqzJc8V6LWjvRy7uJiJaX/myfPsHuoq+SSajZMB9pZHYE8W54ut5xpKY8c+oOYbp2a4edcAnb7SxnX7PYYnUyw0dGJyiV1DXTnLUi3JgXzDSrK7inMxUuVqNuv3kS90munDqJLxsbw3u0dw864BOryuDSM3kU4r/tdXTrFzsJO33Vg4RDrgsJcgIy/RAOVRqI5H8HHgNwH7r/Qe4O+VwRNAUESGgbuAx5RSc0qpeeAx4O4qrCEvGY/AYTOZhddtlJuW2128uJLgp/76Sd7z6aeaVmFxyuoqznMFemCkjwWzMuL5S0s8dW6en7hlR0mx9Vw4GU5j5Qg6/R76O70k0yrvSTiXBLWdbWYJ6eGX5jg9HS45LAT2ucVFDMHFJQ6M5O60LhYams7R05GNle+oVrnl+FyEoW5/XsNY7RCpdQx7bAKFHV43t+wa3DByE//+3CTPXwrxK3fuw1uk+Wt1iH3hBL/1nWiE8ihUaAhE5B7gglLqmayHRoFx298T5rZ822vG+HyETrPypFQqUV8cmw6TVvDdM5d54FtnynqPWrN64snjEdgkqf/xyfP4PC7edkN5SWI7RmjIWdVQp9fwCCD/yMqllcS6oTR2rF6Cf3zScFrLMQRgJIwLCc9lygnzJA8HunyI5A8NrRrm/B7BcF+AWDLtaEaDEwz56fzecrVDpGHTI8hWqr1t3ybOzCxnVACalWQqzccfe4F9W7r5kZcXF1q0kr9OPAK3S+jvbFKPQEQeF5FjOW73AP8V+O+1WJiI3CciR0TkyMxM+VcKE/MrbO/vLOsqtpLSubHpEAA3XtHPH3/lFM9NlDcQvpZk5AzyXIHu39qLS+Dw2Tk+//QFfvhlw5m2+UpwEhrKGAK/mwGlguHMAAAgAElEQVTzx5EvzlrMI7DCG989c5nRYIArN3WXs2y6ikwpO3kpf6IYwO0SBjrzNxhN59AZyma1Qa46J+bsOQTZWCHSQlPeSiEUMwxYd5YhsGZCfOuF5pab+ML3L3Bmdplfe/1VjsqnB7udCc/NmFMCncxLqQVFDYFS6k6l1LXZN+AMsAt4RkReArYBT4vIVuACsN32NtvMbfm259rvA0qpg0qpg5s2lXcFB6XNIchmONjBZJlu+OmpMB1eF3/10wfZ1OPngw99PxPuyMdyLFlz3Xk700VCQwGfmz2bu/mHJ84RjiX5yVsrSxJn3tfrcmAIknhcgs/tor/LuNrPdxWcayiNHWtADRjeQLmhLWNucf5jeOJi8bpyY4h9kdBQoWRxHsmMckik0lxcWCkYNvWaU/2ytZrKxfKosg33lZu6GA0GmlpuIp5M84mvnuZlo33cdc0WR68JBry4JH9eyKKRPQRQQWhIKfWcUmqzUmqnUmonRpjnBqXUJeAQ8NNm9dCtwKJSahJ4FHiDiPSbSeI3mNtqglKqpDkE2Yz0BViOpwrWfefj9HSY3UPdDHT5+JMfv56zl5f53X89kff5T565zJ1/8k1+5JPfrltOYXopRo/fUzBxes1IH7Fkmv1be7hhR3Xy+k5yBMuxFAGfGxGxhYbyewS5msnsWCfQUmQlsik2pezE5BKbevwFf9BDPfmF56ZCUYI5uortDOcR0SuHiwsrpFXx/Fk1m8rCeTwCEeG2fZv4zthlR/OsG8E/Hz7PxPwKv/6GfY4vJlwuYaDLV1SBdCbcuK5iqF0fwSMYHsMY8FfA+wCUUnPAx4DD5u13zW01YWklSahAs0wxtvaVL0c9Nh1m7xYjBPHKKwd5z+1X8tDhcb58bO3A7lRa8YnHT/POv3qCUDTJ1FLM8dQmpRRfPzVd9g9nOhTNWTpqx1Ii/ckqJIktAj4nOYIkXaaBsuKm+Rpz8g2lsTMSDOB1S0UVT8XmFp+4uFSwuQisgeb5QkOxvBVcq6/34fO4qqKMW6x01GIkWDhX9nffeYk/eewFR/sM56gasrh93xChWJK//+45Dr80x9h0iNlwjGQdDEMylebJM5fz/tZX4ik++bUxbtrZ72i0qZ3BLn9R4bmZBgrOAZRWQ1cA0yuw7ivg/Xme9yDwYLX2WwiXC37nRw5wi0Op4WwyA2oWV3LKzOZjOZbkwsIK79y8GgX71Tv38Z9js3zo889x/fZ+tvZ1cGkxygcf+j5Pnp3jnutH+I27ruK2P/w633h+2tHV93+OXeZn/+Ywf/muG7nrmq0l///sIyrz8UPXDfPiTJi3ViFJbNHhdZYjsJq4+gJeRHKHhgoNpbHzrluv4FV7htZdiZZCV4FkcTyZ5vR0iNuLJKILhYamCnQVW4iIoYNVBY/AqSEY7gvw9ednUEqtuxiIxJP8r0dP0d/l49dev6/oPkOmIc11HH5gzxBdPjcf+7f1nnNvh4f3vmYP7y1St18q5y4v89kj4zz81ARTSzG8buHNLx/lvtt2r/nNf/qJc0yHYnzyna8o+YKomNhgOJbMOSGunlTNEDQjPR1efuZVu8p+fWZkZYmJshdnDAGtPZtXk5I+j4v//x3X80N/+m1+7bNH+blX7eI3Hn6GaCLNH739Ot5+4zZEhBt29PP1UzP82huuKrqfLx83vAsnA9FzMR2K8YodwYLPGe4L8Htvva6s98+Hs/LRFJ2m4qfbJfQFvDlDQ4WG0ti5Y/9m7ih/yYCRLM4XGhqbDpNIqYL5ATBCQ8vxFCvx1Lpeg5mlKHs2DRVdR76pa6Vyfi6Cz+0q2LcARuXQSiLFQiSxrljg0NGLhGJJx0nOcNTI/fg964MRvR1evv1br+XCwgrzkTjzkQTzy3HmI3H+7dlJPvf0RFUMQTSR4svHLvHQ4fM8cWYOlxghw4/80ChPn5vnnw+P87mnJ7h93yZ+6bbdXLc9yJ9/80VevXeorIvKwW4/xy7kLxY5fmERpeBlJQx4qjYtbQgqZXOPH5eUXjpnTfbas3mtF7F7Uze/8+YD/NbnnuM7L17m6uFePvnOV6wxGHfs38wfPXqK6VB+qQEwmlq+cnwKyF9WWQilFFNL0YaMxXNSProcS67JXeTTGyqmM1RNrD6CXFfGhaQl7Ax1rfYS2EOW6bQyh9YXPx7DfQG+82Ll1TXjZv6sWPWLlV+5uLiyxhAopfiHJ84BRsI+nVZFDULYlJfId1Xd3+XLWZnmdbv4o0dPMb8cr6hybXopyg998tvMhGJcMdjJb9x1FW+7YVsmDPzml4/wwdft5dNPnOPvvvsSP/HXTzLU7WduOc6vO7g4y8VgV+GBRM+ZRqKUSX/VpuW1hirBY14tldpMMzYTxusWrhhc7+r9+MHt/OKrd3Hfbbv5wvt+YI0RgNUa92JqjEcnFjLln4V0ePKxFE0SS6bL1gyqhIBZPlqoGmslsRoaAgh25vYICg2lqTbdfg/JtCKWYzjMiYtLBLxudg4WVjS1OkezTwxWV7ETwzwa7GBqKVpxUnV8rnDFkMVwMLdnfHR8geMXl9i9qQulVsM+hbBLUJfCTTsHAHjq3HzJr7XzH6dnmQnF+NRP3MDXf/01vP+OPRkjYNHf5eOXX7eXb//Wa/m9t76MgS4vb7thG9dvL+w952Owy0comsw7VOjYhUWG+zo2ZtVQu7C1r4NLS6V7BLuGunJ2HYoIH/mhA/zXN12dc6rRgWFDWOwbRbosHz1+CY9L6O3wlOURzBQpHa0lHV43aQXxAiey5dhqshgMjyBXsrjQUJpq02UaplxTyk5MLrJ/eP0MgmxWu4vX/l+s0tFiYRowQkNptTpUqFyK9RBk9penqewfnjhHt9/Du241BhQtrRT/Hoaj5RmC67b14XULh89VVltydHyBLp+bu6/dWtR76fC6eefNO/jKr97OH//4y8vep9VLkK/Y4dkLiw31BkAbgqKUo+0yNh1ad6XvFBHhjqs2860XZvJWSyhlhIVeeeUgI8EACw5+gNkU6yquJZY4XLSA8NxKPJdHsP7/WWgoTbXJN5NAKWVISxTJD0D+BiMnXcUWmSv0CiqHFiMJFlcSjgzBULcxJc7eVDa/bMTt3/KK0UwubdGJIYgVbv7LR4fXzctG+3jqpco8gqPjC1y3LVjWLI1yGezO7QWC8XmcnV1uaH4AtCEoiiUz4bSpLJpIcX4usi4/UAp37N9EKJrk6fO5y0hPT4c5O7vMGw5sIdjpLUvYrpjOUC2xkqSFKoeWswyBkSPIHxoq1kdQDfLNJLiwsMJSNFk0PwBGmADWnxRyzY7Ox2gVZgmfz4jNFTcELpewpXftDO9/eWqceDLNT916BUFTvsWpIciWl3DKwZ0DPDuxWLTQIB/RRIqTk0tcX6RAotpYxzxX2XAzJIpBG4KiDAcDRBNpx+GXs7PLpBXsLdMjAHjVniE8Lsk7n/bRY5cAeP2BrQQDPhZWSs8ROBE4qxVOxlVG4kk6bSeMgS4fkXiKWNYcg3oni2G9R+Cko9iiw+ump8OzLjRkhXmcGOZKhODCsSR/+c0X+bm/O4xLYL/DsuiR4KpnnE4rPvPkeW7eOcBVW3voC5RgCMoMDYEh1xJPpQtW4BTi+MVFkmlVdqy/XFZDQ+s9gmZIFIM2BEWx4qNOxees2atWM1k59HR4uWnnAF9/Prch+MqJKa7fHmRrX0fekEkxpkMxAl53RXX15dJRZEpZPJkmkVJ0eteGhmB9hVQ9k8VdeTyCE5NL5km1uCGA3LOLp0Oxol3F9nX0BUqbJbwQifPxx17gVb//NX7vS89z1ZYe/vEXb2Wnw3Gd9hne/zE2y7nLEX7qlUZuoCRDUGZoCODgFUZvzeEyw0PfNz3sV9TdEOQXnnuuCRLFoMtHi2LvLr4mj7ywnbGpEC6BXSXMw83FHfs38T8feZ7JxZXMFSAYYYjnLizyW3fvB6Cv08vCSiJnSWMhrFLFanULl0Kx0NBKRnBubbIYjISbPaFqDaXpzKH/X21WB9ivXfeJi/lnEORiqHt9OeHUUrRoV7GdkWDAUVlzKJrgk18b49NPnCMST/GGA1t43x17Sr4qHg4GmHpuklRa8Q/fPcdQt4+7zSbGUg1BuRcfg91+dg918dS5OaD0foJnJhYZ6etwlIepJj1+Dz63K6fMxHNNkCgG7REUJaP26DAxNzYT5orBLkdXdoW4w9TE+UZWGelXjhthIUv0KhjwEU+miSZKKyWcKjASsdZkksV5DMFyZl7xeo8gO08QMkMN9TBo+cZVnpjMP4MgF4ZHkBUactBVbGekr8ORIuiffeNF/uo/zvD6A1t49Fdu44GfPlhWaGSkr4NESvHsxAJfe36Kd9y0HZ/ZFNbpc+NxSVFDkEorIvHUmulkpXJwZz9Hzs2TLkOP6+j4fN3zA2AUgBhD7Nca/1A0wZmZZa7ThqD5Ger243GJY5Gv01PhsiuG7OzZ3M1oMLAuPPTo8Uvs3dzNblNKORMyKTFPkG9EZT0oNsDeKs8MZCWLIXdoqB5hIcidLF5cSTAxn38GQS4MyYHsZHFphtlpd/GzEwtcN9rHJ+59RUkyKdlYXumfPPYCCnjnzatKtCJG53cxQ2B9bl3+8i+SDl4xwEIkwZnZcEmvuxyOMT63Uvf8gEWuIfbWGNhrt2lD0PS4MxUTxa++Eqk0Z2eXK0oUW4gId+zfxH+OzWYSpHPLcb53dm6NrlDQdMvnl0vLE0yXGIqoJgGf8bXLFxqKZDyCtcliyOURFJagria5ksUnHXYU2xnq9jMfSWQawkrpKrYYDnawuJIoKIKnlOLkZMhx7qIQlmf8H6dnee1Vm9mWNQO8FENQieE+uLO8PMHRcSM/8PJtjTIE/nVzKKykd6MrhkAbAkeMBJ2JfJ27vEwyrariEYARHlqOpzhifukfPzlFWrHWEFhXyiV4BOFYkuV4qiGlo2BLFuc1BKtDaSzyJYuX6ugR+DwufG4XYdtMglIqhiyyG4zmS+gqthjN9BLk/17OhGLMLcfZP1y+J2BhCTACmSSxnd6At2hDWUZ5tILQ0K6hLga7fJnfhFOOji/gdgkva9DV91CO0NCzZs6ikfLTFtoQOGBrX4BLDro4x6yKoQp6COy88spBfB5XJjz0leOXGA0GuHZ09aSTqeEuoXLIySSsWtJRJEeQmVds8wj8HjedPjfzy+tzBPXoIbDo8ruJ2JLFTmYQZLMpq8FoqoSuYgsnJaQnLxlT8opJYzuhL+Al4HWzfSDA7XvXK6yW4hHkkqB2iohw4xX9ZsLYOUfHF9i3pafg7I1akkuB9FiTJIpBGwJHjPQZoaFiTWWW2NyVmyurGLLo9Hm4dfcgXz81zXIsybdOz/L6A1vWJEZXcwQlGIJQ47qKwUGy2DzRdmVV4fR3+phrYGgIWDfA3mlHsZ1smYlik+JyMeKgqex5M2zltFegECLCe19zJf/thw7klGYoyRBUWLJ8cGc/L12OMJNHzjubdFpxdHyhYfkBgIEuPyuJVOYiJxRNcKYJOoottCFwwHBfB/Fkuui4udPTYbb1B6p61XHHVZt4cWaZTz9xjngyvW7uQDBQeLB7LixDUEpMupqs9hHkrnRayZEshtwyE/VMFsPaKWXWDIJS8gOwGhqy5hKUI/expbcDlxSuZjs5ucRwX0cmfFgp/+V1e/POvXBkCKrU83EwI0DnzCs4M7tMKJqse/+AnexeAitR3KhQVTbaEDggn/piNqenq1MxZMcqI/3EV0/T3+nlpp1rB9Z0eF34PK6ScgSroaHGeARul+Dz5J9bvJwjWQzrZSZWh9LUMzS0OrfY6QyCbIask4LZaWp5BKWEl4xZwoVzV89fClUlLOQEyxAUKuu0xlSWKzFhce1IH36Py3HC2EoUN6J01GL1mBvf3+cmmidRDNoQOGKkb1WPPR+ptOLMTLgqFUN2dg51sWuoi0g8xZ1Xb8GTpWgqIgQDpekNTYdi+DwuegON6ycsNJwmV7IYDHlgu0cQiadIpVXdQ0NhM3TldAZBNt1+D37PaoPR1JLRVZxLjbYQI8GOvMnieDLN2HS4KmEhJ/QFvEWlqEPR6oSGfB4XL98e5IhDSeqj4/N0+z1cuam6v81SGOxaKzb43AUjUTzYBIli0IbAEU5mF0/MR4gl01VLFNuxZhTkc8tLlZmYNgfSNKKr2KLQcJpIPGl4DVlGr7/Tu8YjqKe8hEWXz53JETidQZCNiKwZWVlqV7HFcDCQN1k8Nh0mmVZ19QigsBR1tXIEYMhNHL+wWHTAEViKo311VRzNxip/tkJDxy4sNk1YCCo0BCLyOyJyQUSOmrc32R77sIiMicgpEbnLtv1uc9uYiHyokv3Xi8EuHz63q6BHsJoorv5Vx0/ecgXvOLidV+/LPcawVOE5o2a9MWEhC2s4TS6WY4byaLahCnb6WFxJkDLDD/UUnLOwJ4udziDIxVC3L1NXPl1iV7HFqNlUlquI4flLhrdydRVKR53Q60Bmwpg6567KCfngzn6SZhK4ENFEiucnQw1NFIMtR7AcZ6nJEsVQHY/g40qp683bIwAicgC4F7gGuBv4MxFxi4gb+BTwRuAA8E7zuU2NyyVs7esomCMYyzGnuFrs2dzNH7z9uryyFX0legSNGlFpp9AA+0g8uS4/AIZHoNTqyWapAR6BlSwuZQZBLuwewXSZch8jfR3EkumcA09OTi7h87hK9lbKxYneUCU6Q9ncuMNIGB95qXDC+NiFxiiOZtPp89Dpc3M5HOP4BbOjuMUMQS7uAR5SSsWUUmeBMeBm8zamlDqjlIoDD5nPbXqG+/LHY8HwCLb0+jM/iHoSDJQYGgrFGm4IAl5XwRxBLhE5S2bCCg9ZHkG9+wiWY0km5p3PIMjFoCk8l04rZsLleQRWEUOu8NDzl0JctaVnXU6pVjiZSRCqQII6m75OL/u2dBfNEzRDothisNvH5eV4U3UUW1TjW/IBEXlWRB4UEaukZRQYtz1nwtyWb3vTM1IgHgvGVLJa5AecEOz0Og4NRRMpQtFk3RUYswn4CuUIUusSxUBmaPlCxhDUbzqZRZffQ1rB980TTCUewdxynLlInERKsaUMw2wfKp+NIS1Rv++jY4+gikb74M4Bnj43nwkV5uLo+AKjwUDDKuTsDHQZ8uPPXlhkNBhomkQxODAEIvK4iBzLcbsH+HMMPdjrgUngj6u1MBG5T0SOiMiRmZnC83vrwda+/APDlVKM1aB01CnBTh/RRNrR5KZGDqSxEygSGsrVi9HfuVZXqRHJYuuK9vDZuZJmEGQz1O0nmVa8MGV0/5aTsxnuy91UNhOKMRuOsb9OiWJwaAiq6BGAkTAOxZKZzzAXjW4kszPUZQjPGR3F9Ts2TihqCJRSdyqlrs1x+6JSakoplVJKpYG/wgj9AFwAttveZpu5Ld/2XPt9QCl1UCl1cNOm9S3t9ebmnQMk04qPHjq+Ljk3uRhlOZ5qoCFwrgdfymzcWlI4R1BaaKiuHoFpoA6/NFfSDIJsrOShpVVUTmhooMuH3+NaJ4hY70QxOJOirmaOAOCmnYXzBLPhGBPzjVMczWaw28f5uUhTzCjOptKqoWHbn28Bjpn3DwH3iohfRHYBe4HvAYeBvSKyS0R8GAnlQ5WsoV7csX8z77n9Sv7xyfP8n2+fXfNYZipZowxBCd3FzeQRRPOEhpZjuZPF2TMJQtEkLlkvRVFLrGaoU1OhkmYQZLPJDAtkDEEZoQsRYTQY4EKWR3AyIy1Rv6tOJ1LU1Q4NbesPsLnHnzdPcNScSPbypjEE/owX+7IGqaDmo9Kj8ocicj2ggJeAXwJQSh0Xkc8CJ4Ak8H6lVApARD4APAq4gQeVUscrXEPd+M27ruKl2WXuf+QkVwx28foDxnCY06ZrundL43IEsBo7L0RG16bRhqBA+ehKPJXzSrvb78HjEuYjVmgoUbehNPY1AChVfn4AYMj8/K2mtHJHFQ7nUMZ9fjLEll5/pna9XjgyBFX0CESEW3YP8pXjU/zT985z703b13wXMoqjTXL1PWg7Hs2yJouKPAKl1LuUUi9TSl2nlHqzUmrS9tj9SqkrlVJXKaW+ZNv+iFJqn/nY/ZXsv964XMLH33E9Lxvt44MPfZ/jF43s/9h0mIEuX91/eBZWfNaJ8Nx0KIbXLZkwS6MolCNYjqdyXuWLCMFO35pkcW+dq7TsSexyK4Zg9aQwNh0uq6vYYqQvsK6s+WQdpSXsFJKiVkpVPUcA8F/ftJ/rtwf58Oef491/c3iNUTw6vsBVW3rKDt9VGyscOBoMNOxckQ/dWVwiAZ+bv/7pg/QFvPz83x5haina0EQxlCZFPbUUZVO3P6eCZD3p8LqJJtI5tWki8eSaecV2Brq8mWSxMYugvobAfiKrxCPo7/ThdknJcwiyGQ4GmAqtFjEY0hLVGUZTKoU8glgyTTKtqhoaAkOO+zO/cAsfu+caDp+d466Pf4vPHhknnVY8M77QFGWjFpbMRLN5A6ANQVls7u3g/7z7JkLRBD//d4d5YSrUsPwAlDacxhhR2fhSOusqLZZcW4UVT6ZJpBSdea6QgzbhOUOCur56SVaOoNQZBNm4XJK5Kqyky3s02IFShoEHODNrCOHVM1FsUcgQZCq8quwRgPFZvuuVO3n0V27j6pFefvPhZ3n7X3yHUCzZNIliWPUImklawkIbgjI5MNLLJ3/iFZy4uMRSNNlQQ9BlVmw4SRZfWoyytUHy03YCeaaUrWQE53KfMOx6Q/UeSgPQbSaxK/EGLKzwUCU17tkDap6frN4wmlIpZAiWqzCUphg7Bjt56Bdv5aM/ciCTe2mk9HQ2ezZ385ZXjPIj1400einr0IagAl67fwu//cOGQkYjrbwRO/dmkqiFuLQUzZw8Gkk+QxBJWBLUuT0CQ4raTBbH6juUBozOYp/HxXVVON6WR1HJyNCRrJGVJyeX8Lld7Bqqj7SEnT4zR5Ar3LcqOFfb4+VyCT/7ql18+YO38Yl7r29oyDYbv8fNx99xPTsGO4s/uc40Toe4RfjZV+3ijdcOZxRKG4VxNVY4NBSOJQlFkw1fK0CHzxpOs9YQWNPJ8iX4rGSxUqruQ2kAPG4XD7/nleyugqSxNamsnK5iC2tSmVVCevJSiL1buvHWSVrCTl/AS1pBOJ6kN8tAW6Ghrhwd47Vg51AXOxtgDDcq2iOoAs1wYu3v9BUNDVky2lubIUeQZ1xlJM9QGouBLi+JlGI5nmqIIQC4bluwKtUvmdBQBcej0+ch2OnNVMucnFxqSKIYbN3FOb6HlkfQU2OPQFMe2hC0CE5mEmQMQRMYrryhIStHUMAjAENWod5DaaqN1UtQ6chQq4R0NhxjJhRrSKIYCktRW9PJapkj0JSPNgQtQl/AV1Ri4pJZWTLcDIbAZ3z1skNDlkeQP1lsGILzlyNAfXWGqs3uoS68bmH7QGUx45FgBxcWVjh1qXGJYig8nCZcpelkmtqgDUGLYHgEhXMEl8yEYqOH0oBtgH0id44gf7LYONmcn7MMwcb1CF5/YAvf+dDrKlbGHAkGmFyM2qQlGuMRFBKes8Z7bmTD3cpoQ9AiBANeluMp4sn16qgWk4tR+ivoYq0m+XIElodQKFkMdkOwcU8sIlJRL4LFSDDA4kqCp87Ns6nH3zB544zUSZ7QkNsl+D36lNOM6KPSIjhRIJ1airK1CUpHYfVEn20Iloski7M9gnr3ETQjVqjvP07PNiwsBEU8AlNeopFzsjX50YagRegzr5QLlZBOLkabIj8AtmTxuhyB1VCWZyxnwItIa4SGqoU1oCYcS3J1g8JCUFiKOlRlwTlNddGGoEUIWsJzBSqHLi1GmyI/APYcwdpQViSexO0SfHnq4D1uF70dXsZbIDRULayRldC4RDEUlqION6jUV+MMbQhahFUp6nyiXykuL8ebxiOwYsW5ksWdPnfBEEJ/pzejUaQ9AqMhzdIQ3N+g0lGLfIZgOa49gmZGG4IWwRpOM5+ncsgaSNMMPQRgXD0GvO6cyeJ8PQQWVsK43kNpmhWP28WW3g68bmH3UGMlFfJJUYej1R1Ko6ku+si0CH1FksXWOMNm8Qgg9wD75Xju6WR2rISxTj6usq0/QH+nD1+Dq3L6At6cFyOhWJJtFfZLaGqHNgQtQm+HB3cBBVJLlKwZ5CUscg2nicRTeRPFFv2mNIMOC63yO2++BrVe663u9AW8vHR5ed32cDRZEwlqTXXQR6ZFsBJ1+WYSWHr1zRIaAujwunIYgiSd3mIegWUI9NfX4poK5idXk7zJYl011NToHEELEQzk1xuaXIzS7fc01VV0wLd+gL0jj8AMDWUrXGoaTy4p6lRaEYmndI6gianYEIjIL4vI8yJyXET+0Lb9wyIyJiKnROQu2/a7zW1jIvKhSvevWaWvM/9gkEuL0abyBiB/aKhYjiCoPYKmxS5FbWE1CWqPoHmp6MiIyB3APcDLlVIxEdlsbj8A3AtcA4wAj4vIPvNlnwJeD0wAh0XkkFLqRCXr0BgEA15mw7lDQ5OL0abKD4DRS2DJE1tEYsmiw8Z1aKh5sUtRWx6bFpxrfir1CN4L/L5SKgaglJo2t98DPKSUiimlzgJjwM3mbUwpdUYpFQceMp+rqQLBTl/BHEFTegTrqoZSRUtCrdBQM4W5NAa5pKjDdRhTqamMSg3BPuDVIvKkiHxTRG4yt48C47bnTZjb8m3XVIG+PDmCZCrNdCjWVKWjYOYIcvQRBIqVj3Zpj6BZySVFHdIeQdNT9MiIyOPA1hwPfcR8/QBwK3AT8FkR2V2NhYnIfcB9ADt27KjGW7Y8wU4voWiSZCqNxybRMBuOk0qr5vQIbIYgnkwTT6UdeAS6fLRZySU8l5lOpg1301L0yCil7sz3mIi8F/i8UkoB3xORNDAEXAC22566zdxGge3Z+30AeADg4MGDTVAh3fwEbT9CuxRxM/YQgJEjsPAnzWkAABETSURBVIeGVjKCc4W/lpt7/Pz8D+7i9Qc213R9mtLJ1di4miPQhrtZqTQ09H+BOwDMZLAPmAUOAfeKiF9EdgF7ge8Bh4G9IrJLRHwYCeVDFa5BY2JV02TrwTdjDwFYoaFV0blIwpxOVsQjcLmE3/7hA+zZ3FhdHc16cnkEyzpH0PRUemQeBB4UkWNAHHi36R0cF5HPAieAJPB+pVQKQEQ+ADwKuIEHlVLHK1yDxiSf8NyqvERzzCKwCHjdxFPpTCjLmk5WzBBompeuHFLUIcsQFMn9aBpHRUfGrPz5qTyP3Q/cn2P7I8AjlexXk5tgnpkElxaj+DyuTLVNs5CZUpZM0+12ZUJDxfoINM1LLilqKzTUVaRRUNM4dGdxC5FvJsGlJaOHoNkE2jp8a4fTWI1H2iPY2KwzBLEEAa97TQGDprnQR6aFKBQaarb8AKyfWxyxDIEuM9zQ9K4zBFqCutnRhqCF6OkwxjhmJ4svNdGISjuZcZUZQ2CFhrRHsJFZ7xGktPJok6MNQQvhdgm9HV4WbXrwSqlMaKjZCPjMKWWmAYiYyeJiEhOa5mZ9jiBBlzYETY02BC1GsNO7xiOYjySIJ9NNGRrqyPIIrByBThZvbNZ7BFqCutnRhqDFyJaitprJNlJoSHsEG5tsKeqQHlPZ9GhD0GL0dfpYsIWGLi1azWTN1UMAqyf8mC1Z7HZJZrC9ZmOSLUUdjunpZM2O/sW1GMHA2tDQJauruBlzBNmhoZgxuL7Zylw1pWGXogZdNbQR0IagxQh2rg0NXVqM4nYJm3r8BV7VGDKGIJ42/03pHoIWwC5FrZRiWecImh5tCFqMYKePpWiClBmfnVyMsrnHj9vVfFfZmYYyW7JYJ4o3PnYp6lgyTSKldNVQk6MNQYsRDHhRCkJRwytoxhGVFh2etQ1lKw7mFWuaH7vwnJag3hhoQ9BiZHcXN2sPAYDXLbhdskZiotOrTxgbHbsUtR5TuTHQhqDFyBiCleb3CERkzXCaiPYIWoJcHoE2BM2NNgQtRl/AnEkQiROKGj/EZuwhsOjINgQ6Wbzh6fK5cZtS1JkxlTo01NRoQ9BiBG1uudVDsKVJQ0NgyExEMxITSTp1snjDIyIEze7iZe0RbAi0IWgx7FLUVg9Bsw2ksWMPDS3HU1pwrkWwZCZ0aGhjoI9Oi2HFZ+cjcSYXjZNqM4eG7IZgJZ4ioD2ClsCSog7pMZUbAn10WgyP20WP38NCJIHL7NDd3Nt8zWQW1gD7RCpNPJXWHkGL0BfwshCJZ6qGevTg+qZGh4ZakL5O42rs0lKUwS4ffk/znlyNAfapjOCcHkrTGqyGhhK4XUKHV59qmpmKjo6I/LOIHDVvL4nIUdtjHxaRMRE5JSJ32bbfbW4bE5EPVbJ/TW4MmYl4U5eOWlihoYgeU9lS9GWSxSm6/R6tH9XkVDq8/h3WfRH5Y2DRvH8AuBe4BhgBHheRfeZTPwW8HpgADovIIaXUiUrWoVlLf6ePhZUE0USa0eDGMATL5lAabQhag76Al6VokqVoQieKNwBV8dfEMPc/DvyTueke4CGlVEwpdRYYA242b2NKqTNKqTjwkPlcTRXpC3hZjCS4tLjS9B5Bh8/NSjyd6S7W5aOtQV/ASyqtmFqKakOwAahW4O7VwJRS6rT59ygwbnt8wtyWb7umigQ7vUwtRZmPJJpWXsIi4DVyBKvTybRH0ApY1WsX5ld0xdAGoOgREpHHga05HvqIUuqL5v13suoNVAURuQ+4D2DHjh3VfOuWJxjwsWxeYTfjQBo763IE+uqxJbCkqC8uRLlisKvBq9EUo+ivTil1Z6HHRcQDvBW40bb5ArDd9vc2cxsFtmfv9wHgAYCDBw+qYuvUrGJ1F0Nz9xCAUTWUSqvMjFudI2gNLI8gnkprj2ADUI3Q0J3A80qpCdu2Q8C9IuIXkV3AXuB7wGFgr4jsEhEfRkL5UBXWoLFh/QiB5s8RmMNpLoeN8ZraELQG9u9gt877ND3VOEL3khUWUkodF5HPAieAJPB+pVQKQEQ+ADwKuIEHlVLHq7AGjY1gpy9zfyPkCAAuLxuGQA+maQ36bF6p9gian4qPkFLqZ/Jsvx+4P8f2R4BHKt2vJj9WaKinw9P0k6ECPsMpvRyOmX9rj6AVWOMRNPl3UKM7i1sSS3iu2fMDYPMIwnHcLsHv0V/JVsCSogY9nWwjoH91LYjlljd7xRDYcgTLcTp9bt2B2iKISMYr0B5B86MNQQsSNIfTbG1isTkLyyOYMw2BpnWwPFOdI2h+tCFoQXweF6/aM8gPXDnU6KUUxcoJXA7HdKK4xbB6CZo9T6XRMtQty2d+4dZGL8ERlkewHE/pRHGLYYWGerQhaHq0R6BpKFaOAHTpaKvRp0NDGwZtCDQNxe4FdPq1R9BK6GTxxkEbAk1DCdg8Ap0sbi1WQ0N6Olmzo021pqF0rDEE+uvYStyye4BnJoZ0aGgDoI+QpqG4XYLP7dLziluQV+/dxKv3bmr0MjQO0KEhTcOx5tkGtEeg0TQEbQg0DcdKGGuPQKNpDNoQaBqOlTDWfQQaTWPQhkDTcKyEse5A1WgagzYEmoZjeQK6fFSjaQzaEGgajhUa0uWjGk1j0IZA03AsQ6CTxRpNY9CGQNNwOqzQkM4RaDQNQRsCTcNZDQ1pj0CjaQQVGQIRuV5EnhCRoyJyRERuNreLiPypiIyJyLMicoPtNe8WkdPm7d2V/gc0Gx9tCDSaxlKpL/6HwP9QSn1JRN5k/v0a4I3AXvN2C/DnwC0iMgB8FDgIKOApETmklJqvcB2aDcxq1ZAODWk0jaDS0JACes37fcBF8/49wN8rgyeAoIgMA3cBjyml5syT/2PA3RWuQbPB6dAegUbTUCq9BPsV4FER+V8YRuUHzO2jwLjteRPmtnzbNW3Mm18+TKfPvUaJVKPR1I+ihkBEHge25njoI8DrgF9VSn1ORH4c+D/AndVYmIjcB9wHsGPHjmq8paZJ2bO5hz2bexq9DI2mbSlqCJRSeU/sIvL3wAfNP/8F+Gvz/gVgu+2p28xtFzByCPbt38iz3weABwAOHjyoiq1To9FoNOVRaY7gInC7ef+1wGnz/iHgp83qoVuBRaXUJPAo8AYR6ReRfuAN5jaNRqPRNIhKcwS/CHxCRDxAFDOUAzwCvAkYAyLAzwIopeZE5GPAYfN5v6uUmqtwDRqNRqOpgIoMgVLq28CNObYr4P15XvMg8GAl+9VoNBpN9dCdxRqNRtPmaEOg0Wg0bY42BBqNRtPmaEOg0Wg0bY4Yed3mRkRmgHMVvMUQMFul5VQTva7S0OsqDb2u0mjFdV2hlNpU7EkbwhBUiogcUUodbPQ6stHrKg29rtLQ6yqNdl6XDg1pNBpNm6MNgUaj0bQ57WIIHmj0AvKg11Uael2loddVGm27rrbIEWg0Go0mP+3iEWg0Go0mDy1jCETkx0TkuIikRSRvhl1E7haRU+Y85Q/Ztu8SkSfN7f8sIr4qrWtARB4zZzQ/ZqquZj/nDnPus3WLisiPmo/9rYictT12fb3WZT4vZdv3Idv2Rn5e14vId83j/ayIvMP2WNU+r3zfFdvjfvP/PmZ+Fjttj33Y3H5KRO4qdw1lruvXROSE+dl8VUSusD2W83jWcW0/IyIztjX8gu2xmswzd7Cmj9vW84KILNgeq9nnJSIPisi0iBzL87hIvWa/K6Va4gZcDVyFMd/gYJ7nuIEXgd2AD3gGOGA+9lngXvP+XwDvrdK6/hD4kHn/Q8AfFHn+ADAHdJp//y3w9hp8Xo7WBYTzbG/Y5wXsA/aa90eASSBYzc+r0HfF9pz3AX9h3r8X+Gfz/gHz+X5gl/k+7ip9Pk7WdYft+/Nea12Fjmcd1/YzwP/O8doB4Iz5b795v78ea8p6/i8DD9bp87oNuAE4lufxNwFfAgS4FXiyVp9Vy3gESqmTSqlTRZ52MzCmlDqjlIoDDwH3iIhgzFN42Hze3wE/WqWl3WO+n9P3fTvwJaVUpEr7z0ep68rQ6M9LKfWCUuq0ef8iMA0UbZopkZzflQJrfRh4nfnZ3AM8pJSKKaXOYsix31yvdSmlvm77/jyBMQCqHjj5zPJRq3nmpa7pncA/VWG/RVFKfQvjoi8fdZv93jKGwCH5ZiYPAgtKqWTW9mqwRRlDeQAuAVuKPP9e1n8R7zddw4+LiL/O6+oQkSMi8oQVrqKJPi8RuRnjSu9F2+ZqfF5O5mtnnmN+FosYn00tZ3OX+t4/j3FVaZHreFYLp2t7m3l8HhYRa5JhrT4zx+9rhtB2AV+zba7l51WMus1+r3QwTV2RAvOTlVJfrPd6LAqty/6HUkqJSN4yLdPav4y1U9s+jHFC9GGUkf0W8Lt1XNcVSqkLIrIb+JqIPIdxwiubKn9e/wC8WymVNjeX/Xm1GiLyU8BBVqcIQo7jqZR6Mfc71IR/Bf5JKRUTkV/C8KheW8f9F+Je4GGlVMq2rdGfV13YUIZAFZif7JB8s5QvY7hdHvPKztpe8bpEZEpEhpVSk+aJa7rAW/048AWlVML23tbVcUxE/gb4/+q5LqXUBfPfMyLyDeAVwOdo8OclIr3Av2NcBDxhe++yP68s8n1Xcj1nQowpfX0Y3yUnry0XR+8tIndiGNbblVIxa3ue41mtE1vRtSmlLtv+/GuMnJD12tdkvfYb9ViTjXvJGqhV48+rGBXPfndKu4WGDgN7xah48WEc+EPKyMB8HSM+D/BuoFoexiHz/Zy877r4pHkytOLyPwrkrDCoxbrEmC3tN+8PAa8CTjT68zKP3Rcw4qcPZz1Wrc8r53elwFrfDnzN/GwOAfeKUVW0C9gLfK/MdZS8LhF5BfCXwJuVUtO27TmPZ5XW5XRtw7Y/3wycNO/Xap65k+OIiOzHSLx+17at1p9XMeo3+73amfBG3YC3YMTKYsAU8Ki5fQR4xPa8NwEvYFj1j9i278b4sY4B/wL4q7SuQeCrwGngcWDA3H4Q+Gvb83ZiWHpX1uu/BjyHcUL7NNBdr3UBP2Du+xnz359vhs8L+CkgARy13a6v9ueV67uCEWZ6s3m/w/y/j5mfxW7baz9ivu4U8MYqf9eLretx8zdgfTaHih3POq7t94Dj5hq+Duy3vfbnzM9yDPjZeq3J/Pt3gN/Pel1NPy+Mi75J87s8gZHPeQ/wHvNxAT5lrvs5bNWQ1f6sdGexRqPRtDntFhrSaDQaTRbaEGg0Gk2bow2BRqPRtDnaEGg0Gk2bow2BRqPR1IhiwnJlvN+XRWRBRP4ta/suqUAEUhsCjUajqR1/S3U0kyz+CHhXju1/AHxcKbUHmMcoRXWMNgQajUZTI1QOYTkRudK8sn9KRP7DbGZz+n5fBUJZ71exCOSGkpjQaJoJEbkG+ASwA0PzaDNGt/Phhi5M0+w8gNE0dlpEbgH+jMr0lioWgdSGQKMpAxGxuop/DEMP/nngKW0ENIUQkW6MjuV/MS7kAWNuBSLyVnILJF5QSlV1uFE22hBoNOVxJ/B9pdRxyOgf/XFjl6TZALgwrt7XTc5TSn0e+HwZ71mRaKa1KI1GUzrXA98HEJERjElW/9nYJWmaHaXUEnD2/7V3t0YIBEEQhft5QsAjcSRADCjkJUBhCIEo0ARCBHgkghgwh1h+CseBQMz7Atha11uzVT3AInmuo5z+eObPJZAGgfSda15z2G3a/gPpDbBPazSdAGegS7JM0gHHtAK+T7e4BTikjSTn9/MeI6NNkjVwSvsz2A26p6Vz0nDAOO3VNUqrfJ4lufR9v/rrxaQvGASSVJyjIUkqziCQpOIMAkkqziCQpOIMAkkqziCQpOIMAkkqziCQpOJuFeTOHrUSc1IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lin = np.linspace(-1e-10, 1e-10)\n",
    "\n",
    "def plotting(alpha):\n",
    "    return nlml(np.array(m.x) - alpha*np.array(m.jac))\n",
    "\n",
    "out = [plotting(l) for l in lin]\n",
    "\n",
    "plt.plot(lin, out)\n",
    "plt.xlabel(r'$\\alpha$')\n",
    "#plt.savefig('Short_range.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3. ADAM**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def adams(grad, init, n_epochs=500, eta=10**-4, gamma=0.9,beta=0.99,epsilon=10**-8,noise_strength=0):\n",
    "    params=np.array(init)\n",
    "    param_traj=np.zeros([n_epochs+1,4])\n",
    "    param_traj[0,]=init\n",
    "    v=0;\n",
    "    grad_sq=0;\n",
    "    for j in range(n_epochs):\n",
    "        noise=noise_strength*np.random.randn(params.size)\n",
    "        g=np.array(grad(params))+noise\n",
    "        v=gamma*v+(1-gamma)*g\n",
    "        grad_sq=beta*grad_sq+(1-beta)*g*g\n",
    "        v_hat=v/(1-gamma)\n",
    "        grad_sq_hat=grad_sq/(1-beta)\n",
    "        params=params-eta*np.divide(v_hat,np.sqrt(grad_sq_hat+epsilon))\n",
    "        param_traj[j+1,]=params\n",
    "    return param_traj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running ADAM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "m_a = np.squeeze(adams(grad_nlml, np.random.rand(4)))\n",
    "t_adam = time.time() - t1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Comparison**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The time it took Nelder-Mead to converge:        208.37455773353577 seconds\n",
      "The time it took Conjugate Gradient to converge: 19.78248691558838 seconds\n"
     ]
    }
   ],
   "source": [
    "print('The time it took Nelder-Mead to converge:       ', t_Nelder, 'seconds')\n",
    "print('The time it took Conjugate Gradient to converge:', t_CG, 'seconds')\n",
    "#print('The time it took Adam to converge:              ', t_adam, 'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The inferred parameter with Nelder-Mead is:        1.99978331424706\n",
      "The inferred parameter with Conjugate Gradient is: 1.9998539316144641\n"
     ]
    }
   ],
   "source": [
    "print('The inferred parameter with Nelder-Mead is:       ', m_n.x[3])\n",
    "print('The inferred parameter with Conjugate Gradient is:', m.x[3])\n",
    "#print('The inferred parameter with Adam is:              ', m_a[-1,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}