jit-example-pln-jax.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4056e632-6640-450d-8e0f-f16b424a8cec",
   "metadata": {},
   "source": [
    "# FinistR : bootcamp R à Roscoff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52ac98d9-d91e-4019-897a-a2304e7d8fcf",
   "metadata": {},
   "source": [
    "JAX is a Python library, basically a wrapper around numpy for efficient scientific programming: automatic differentiation, parallelization, JIT, etc. Many numpy functions are rewritten in a low level API called LAX. It then uses the XLA compiler to opitmization computations, see [#numpy-lax-xla-jax-api-layering](https://jax.readthedocs.io/en/latest/notebooks/thinking_in_jax.html#numpy-lax-xla-jax-api-layering) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9fd3c6c-1e17-4a8b-9469-a805b8addf97",
   "metadata": {},
   "source": [
    "__This tutorial is a non exhaustive but dense introduction to JAX, especially JIT compilation with JAX. We try to point the reader to the main concepts with a lot of redirections to external resources and to JAX documentation.__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31d1bceb-0415-44f3-b75b-2d0f81e3f9a6",
   "metadata": {},
   "source": [
    "## JIT with JAX"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13797785-16b5-4390-bc6c-83147d8d2ca4",
   "metadata": {},
   "source": [
    "A practical definition of JIT compilation can be found in [JAX documentation](https://jax.readthedocs.io/en/latest/index.html):\n",
    "\n",
    "> When we jit-compile a function, we usually want to compile a version of the function that works for many different argument values, so that we can cache and reuse the compiled code. That way we don’t have to re-compile on each function evaluation.\n",
    "\n",
    "> For example, if we evaluate an `@jit` function on the array `jnp.array([1., 2., 3.], jnp.float32)`, we might want to compile code that we can reuse to evaluate the function on `jnp.array([4., 5., 6.], jnp.float32)` to save on compile time.\n",
    "\n",
    "> By default JAX executes operations one at a time, in sequence. Using a just-in-time (JIT) compilation decorator, sequences of operations can be optimized together and run at once.\n",
    "\n",
    "**In this tutorial we explore JAX JIT compilation on the same algorithm of the PLN model** as coded in [the PLN in pytorch tutorial](https://stateofther.github.io/finistR2023/torch_Python-PLN.html).\n",
    "Also, this tutorial is to be compared with [JIT compilation in pytorch](https://stateofther.github.io/finistR2023/jit-example-pln.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffd0991c-6ef8-4bb7-a974-b38c9b1df709",
   "metadata": {},
   "source": [
    "First, we make the necessary imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "35c14c1d-d801-447f-8375-78fffbd65fec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "#os.environ['CUDA_VISIBLE_DEVICES'] = '' # uncomment to force CPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a7a692fa-521d-4488-bb3c-7b5921da782a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import pyPLNmodels\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pyPLNmodels.models import PlnPCAcollection, Pln\n",
    "from pyPLNmodels.oaks import load_oaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "146a9584-e099-45be-acff-e6a0a4cbf518",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[gpu(id=0)]\n"
     ]
    }
   ],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "import optax\n",
    "jax.config.update(\"jax_enable_x64\", False)\n",
    "print(jax.devices())\n",
    "myfloat = np.float32"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "622d3141-5a6e-480c-8281-577c3bba2a2a",
   "metadata": {},
   "source": [
    "**Note**: We use a GPU since JIT compilation is particularly efficient on GPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "570dd0ec-5916-48db-b86a-3f03422c24e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "oaks = load_oaks()\n",
    "Y = np.asarray(oaks['counts']).astype(myfloat)\n",
    "Y = np.repeat(Y, 100, axis=0) # make data bigger to feel the speed up\n",
    "O = np.log(oaks['offsets']).astype(myfloat)\n",
    "O = np.repeat(O, 100, axis=0) # make data bigger to feel the speed up\n",
    "X = np.ones([Y.shape[0],1]).astype(myfloat)\n",
    "\n",
    "N_iter = 1000\n",
    "lr = 1e-4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c86dda3-2b04-4a78-a858-592d2dee7674",
   "metadata": {},
   "source": [
    "### JAX without JIT\n",
    "This section does not use any kind of jitting. This is our baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "58ef40bb-a938-4fa4-ac53-5c6b7d52784d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _log_stirling(integer):\n",
    "    integer_ = integer + (integer == 0)  # Replace 0 with 1 since 0! = 1!\n",
    "    return jnp.log(jnp.sqrt(2 * jnp.pi * integer_)) + integer_ * jnp.log(integer_ / jnp.exp(1))\n",
    "\n",
    "class PLN:\n",
    "    def __init__(self, Y, O, X): \n",
    "        self.Y = Y\n",
    "        self.O = O\n",
    "        self.X = X\n",
    "        self.n, self.p = Y.shape\n",
    "        self.d = X.shape[1]\n",
    "        ## Variational parameters\n",
    "        self.M = jnp.full(Y.shape, 0.0)\n",
    "        self.S = jnp.full(Y.shape, 1.0)\n",
    "        ## Model parameters\n",
    "        self.B = jnp.zeros((self.d, self.p))\n",
    "        self.Sigma = jnp.eye(self.p)\n",
    "        self.Omega = jnp.eye(self.p)\n",
    "\n",
    "    def get_Sigma(self, n, M, S) :\n",
    "        return 1/n * (M.T @ M + jnp.diag(jnp.sum(S**2, axis=0)))\n",
    "    \n",
    "    def get_ELBO(self, optim_params):\n",
    "      B, M, S = optim_params\n",
    "      S2 = jnp.square(S)\n",
    "      XB = self.X @ B\n",
    "      A = jnp.exp(self.O + M + XB + S2/2)\n",
    "\n",
    "      elbo = self.n/2 * jnp.log(jnp.linalg.det(self.Omega))\n",
    "      elbo += jnp.sum(- A + self.Y * (self.O + M + XB) + .5 * jnp.log(S2))\n",
    "      elbo -= .5 * jnp.trace(M.T @ M + jnp.diag(jnp.sum(S2, axis=0)) @ self.Omega)\n",
    "      elbo += .5 * self.n * self.p  - jnp.sum(_log_stirling(self.Y))\n",
    "      return -elbo\n",
    "\n",
    "    def fit(self, N_iter, lr, tol = 1e-8) :\n",
    "        ELBO = jnp.zeros(N_iter)\n",
    "        optimizer = optax.chain(\n",
    "            #adam(learning_rate=lr)\n",
    "            optax.scale_by_radam(),\n",
    "            optax.scale(-1.0),\n",
    "            optax.clip(0.1),\n",
    "        )\n",
    "        opt_state = optimizer.init((self.B, self.M, self.S))\n",
    "        \n",
    "        for i in range(N_iter):\n",
    "            loss_value, grads = jax.value_and_grad(\n",
    "                self.get_ELBO, 0\n",
    "            )((self.B, self.M, self.S))\n",
    "    \n",
    "            updates, opt_state = optimizer.update(grads, opt_state, (self.B, self.M, self.S))\n",
    "            optim_params = optax.apply_updates((self.B, self.M, self.S), updates)\n",
    "            self.B, self.M, self.S = optim_params\n",
    "    \n",
    "            ## update parameters with close form\n",
    "            self.Sigma = self.get_Sigma(self.n, self.M, self.S)\n",
    "            self.Omega = jnp.linalg.inv(self.Sigma)\n",
    "    \n",
    "            objective = loss_value\n",
    "            ELBO = ELBO.at[i].set(objective)\n",
    "        \n",
    "        return ELBO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aebf9f2-a55c-4423-8742-1d9d885f6744",
   "metadata": {},
   "source": [
    "Have a look a the way we update the ELBO vector in the previous functions. This is because JAX arrays are immutable [#in-place-updates](https://jax.readthedocs.io/en/latest/notebooks/Common_Gotchas_in_JAX.html#in-place-updates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "887fc6b2-e1ae-4318-9c06-64c85a01aac4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10min 53s, sys: 21min 15s, total: 32min 8s\n",
      "Wall time: 2min\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "pln = PLN(Y, O, X)\n",
    "\n",
    "with jax.default_device(jax.devices('cpu')[0]): # as opposed to the subsequent jitted versions of the code, running this code do\n",
    "    jaxELBO_no_jit = jax.block_until_ready(\n",
    "        pln.fit(N_iter, lr, tol=1e-8)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e705cd52-72d5-4adb-be8c-473b308d55b5",
   "metadata": {},
   "source": [
    "### JIT with JAX level 1: `jax.jit`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ead6c16-d951-48c1-943a-cd76b985ee6f",
   "metadata": {},
   "source": [
    "In this section, we will jit the functions used in the optimization process. Since it is not straightforward to JIT compilation class methods (see subsequent section), we will not use a class anymore. \n",
    "\n",
    "Thus, we first create an independant function to initialize variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5307021b-5159-4727-b3d7-9c1c53843c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_params(Y, O, X): \n",
    "    n, p = Y.shape\n",
    "    d = X.shape[1]\n",
    "    ## Variational parameters\n",
    "    M = jnp.full(Y.shape, 0.0)\n",
    "    S = jnp.full(Y.shape, 1.0)\n",
    "    ## Model parameters\n",
    "    B = jnp.zeros((d, p))\n",
    "    Sigma = jnp.eye(p)\n",
    "    Omega = jnp.eye(p)\n",
    "\n",
    "    return n, p, d, M, S, B, Sigma, Omega"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cfd17c7-4381-45cb-96ff-8df3d75104a8",
   "metadata": {},
   "source": [
    "Let's create the jitted functions. Note that `_log_stirling` will be automatically jitted when called the jitted `get_ELBO`. The actual _Just In Time_ compilation will actually happen at the first execution of the jitted function. Note that the call `jax.jit()` is equivalent to using the decorator `@jax.jit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d6f555ff-c562-4624-8801-41fd2397a2da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _log_stirling(integer):\n",
    "    integer_ = integer + (integer == 0)\n",
    "    return jnp.log(jnp.sqrt(2 * jnp.pi * integer_)) + integer_ * jnp.log(integer_ / jnp.exp(1))\n",
    "    \n",
    "def get_ELBO(optim_params, other_params): \n",
    "    B, M, S = optim_params['B'], optim_params['M'], optim_params['S']\n",
    "    X, O, n, Omega, Y, p = (other_params['X'], other_params['O'], \n",
    "        other_params['n'], other_params['Omega'], other_params['Y'],\n",
    "        other_params['p'])\n",
    "    S2 = jnp.square(S)\n",
    "    XB = X @ B\n",
    "    A = jnp.exp(O + M + XB + S2/2)\n",
    "\n",
    "    elbo = 0.\n",
    "    elbo = n / 2 * jnp.log(jnp.linalg.det(Omega))\n",
    "    elbo += jnp.sum(- A + Y * (O + M + XB) + .5 * jnp.log(S2))\n",
    "    elbo -= .5 * jnp.trace(M.T @ M + jnp.diag(jnp.sum(S2, axis = 0)) @ Omega)\n",
    "    elbo += .5 * n * p  - jnp.sum(_log_stirling(Y))\n",
    "    return -elbo\n",
    "\n",
    "jit_loss_and_grad = jax.jit(jax.value_and_grad(get_ELBO, 0)) # JIT !\n",
    "\n",
    "def get_Sigma(n, M, S) :\n",
    "    return 1/n * (M.T @ M + jnp.diag(jnp.sum(S**2, axis = 0)))\n",
    "\n",
    "jit_getSigma = jax.jit(get_Sigma) # JIT !\n",
    "\n",
    "jit_inv = jax.jit(jnp.linalg.inv) # JIT !"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5edaee77-7a3f-41d8-902f-d66d3b77f8f3",
   "metadata": {},
   "source": [
    "In the following, each function inside the for loop is jitted but the for loop\n",
    "is not jitted itself; it is inefficient to do so in JAX mostly due to the long compilation time it would induce, see [#jit-decorated-function-is-very-slow-to-compile](https://jax.readthedocs.io/en/latest/faq.html#jit-decorated-function-is-very-slow-to-compile)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "50ecdfe3-4035-42a7-b556-cc9b4d793ee8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaxfit1(optim_params, other_params, N_iter, lr, tol = 1e-8) :\n",
    "    ELBO = jnp.zeros(N_iter)\n",
    "    optimizer = optax.chain(\n",
    "        #adam(learning_rate=lr)\n",
    "        optax.scale_by_radam(),\n",
    "        optax.scale(-1.0),\n",
    "        optax.clip(0.1),\n",
    "    )\n",
    "    opt_state = optimizer.init(optim_params)\n",
    "    objective0 = jnp.inf\n",
    "\n",
    "    update = jax.jit(optimizer.update)\n",
    "    apply_updates = jax.jit(optax.apply_updates)\n",
    "\n",
    "    for i in range(N_iter):\n",
    "        loss_value, grads = jit_loss_and_grad(\n",
    "            optim_params,\n",
    "            other_params\n",
    "        )\n",
    "\n",
    "        updates, opt_state = update(grads, opt_state, optim_params)\n",
    "        optim_params = apply_updates(optim_params, updates)\n",
    "\n",
    "        ## update parameters with close form\n",
    "        other_params['Sigma'] = jit_getSigma(other_params['n'], optim_params['M'],\n",
    "        optim_params['S'])\n",
    "        other_params['Omega'] = jit_inv(other_params['Sigma'])\n",
    "\n",
    "        objective = loss_value\n",
    "        ELBO = ELBO.at[i].set(objective)\n",
    "\n",
    "    return ELBO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "350fb0ae-3d2f-4729-8d97-80994ae7cdee",
   "metadata": {},
   "source": [
    "Initialize the data for JAX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bf550788-0ffb-4a04-993d-142eb04dea66",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = jnp.asarray(Y)\n",
    "O = jnp.asarray(O)\n",
    "X = jnp.array(X)\n",
    "n, p, d, M, S, B, Sigma, Omega = init_params(Y, O, X)\n",
    "optim_params = {'B':B, 'M':M, 'S':S}\n",
    "other_params = {'X':X, 'O':O, 'n':n, 'Omega':Omega, 'Y':Y, 'p':p, 'Sigma':Sigma}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c915b5-e2ce-4a34-8a1a-af9db31c59e2",
   "metadata": {},
   "source": [
    "and run the learning process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "549f88cc-55ea-43e0-bb2c-1dcd69bb8f9d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.87 s, sys: 500 ms, total: 4.37 s\n",
      "Wall time: 4.13 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "jaxELBO1 = jax.block_until_ready(\n",
    "    jaxfit1(optim_params, other_params, N_iter, lr=lr, tol=1e-8)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "882cc6f9-9742-4ed2-949b-f4fd4e03c88d",
   "metadata": {},
   "source": [
    "__Note__: Because of the complex JAX internal mechanics, benchmarking with JAX should be done with caution\n",
    "[https://jax.readthedocs.io/en/latest/async_dispatch.html](https://jax.readthedocs.io/en/latest/async_dispatch.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2dd5f54-1bd4-406e-9b6d-fceda04a4986",
   "metadata": {},
   "source": [
    "### JIT with JAX level 2: `jax.lax.scan`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f11e6b9-9641-4f27-bed9-657212fbad12",
   "metadata": {},
   "source": [
    "In this section, we add resort to `jax.lax.scan`, which is the standard way to write `for` loops in JAX.\n",
    "`jax.lax.scan` automatically JIT compiles its content, it’s not necessary to add the calls to `jax.jit` here.\n",
    "\n",
    "For a complete tutorial on `jax.lax.scan` see [https://ericmjl.github.io/dl-workshop/02-jax-idioms/02-loopy-carry.html](https://ericmjl.github.io/dl-workshop/02-jax-idioms/02-loopy-carry.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2452c60f-2748-4743-ab56-4fc637d5115c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaxfit2(optim_params, other_params, N_iter, lr, tol = 1e-8) :\n",
    "    optimizer = optax.chain(\n",
    "        #adam(learning_rate=lr)\n",
    "        optax.scale_by_radam(),\n",
    "        optax.scale(-1.0),\n",
    "        optax.clip(0.1),\n",
    "    )\n",
    "    opt_state = optimizer.init(optim_params)\n",
    "    \n",
    "    def scan_fun(carry, _):\n",
    "        optim_params = carry['optim_params']\n",
    "        other_params = carry['other_params']\n",
    "        loss_value, grads = jax.value_and_grad(get_ELBO)(optim_params,\n",
    "        other_params)\n",
    "    \n",
    "        updates, opt_state = optimizer.update(grads, carry['opt_state'])\n",
    "        optim_params = optax.apply_updates(optim_params, updates)\n",
    "    \n",
    "        ## update parameters with close form\n",
    "        other_params['Sigma'] = get_Sigma(other_params['n'], optim_params['M'],\n",
    "        optim_params['S'])\n",
    "        other_params['Omega'] = jnp.linalg.inv(other_params['Sigma'])\n",
    "    \n",
    "        carry['optim_params'] = optim_params\n",
    "        carry['other_params'] = other_params\n",
    "        carry['opt_state'] = opt_state\n",
    "        return carry, loss_value\n",
    "\n",
    "    \n",
    "    carry, ELBO = jax.lax.scan(\n",
    "        scan_fun,\n",
    "        {\n",
    "            \"optim_params\":optim_params,\n",
    "            \"other_params\":other_params,\n",
    "            'opt_state':opt_state\n",
    "        },\n",
    "        jnp.arange(N_iter)\n",
    "    )\n",
    "\n",
    "    return ELBO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25c103d0-6d8b-4099-b9e3-57b3478f1630",
   "metadata": {},
   "source": [
    "Initialize the data for JAX:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "402d8997-3c76-4176-b496-531b4af8c567",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = jnp.asarray(Y)\n",
    "O = jnp.asarray(O)\n",
    "X = jnp.array(X)\n",
    "n, p, d, M, S, B, Sigma, Omega = init_params(Y, O, X)\n",
    "optim_params = {'B':B, 'M':M, 'S':S}\n",
    "other_params = {'X':X, 'O':O, 'n':n, 'Omega':Omega, 'Y':Y, 'p':p, 'Sigma':Sigma}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ccd9619-d192-45ca-8c45-b07859bf4f7a",
   "metadata": {},
   "source": [
    "**Note**: Let's us try the Ahead Of Time (AOT) compilation where the compilation is done beforehand in a particular explicit step and not at the first execution of the function (as in classical JIT compilation). By doing so we can get an estimate of the compilation time, which is, in our case, not a bottleneck."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6d046ba4-99ad-44a2-a74a-e56856b34e4b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 708 ms, sys: 40.1 ms, total: 748 ms\n",
      "Wall time: 502 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lowered = jax.jit(jaxfit2, static_argnums=(2,3,4)).lower(optim_params,\n",
    "other_params, N_iter, lr=lr, tol=1e-8)\n",
    "compiled = lowered.compile()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7c3c94-2291-4eba-ad32-aa5b81c23632",
   "metadata": {},
   "source": [
    "Let's run the complete fit (including compilation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "cc684410-91e3-4e91-b4ee-6bb784b8bd41",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.97 s, sys: 1.08 s, total: 3.05 s\n",
      "Wall time: 2.83 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "jaxELBO2 = jax.block_until_ready(\n",
    "    jaxfit2(optim_params, other_params, N_iter, lr, tol=1e-8)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc7ed2ca-ca9d-4b85-ae27-774f91921955",
   "metadata": {},
   "source": [
    "### JIT with JAX level 3: can we jit everything? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d929e407-98da-481e-af93-9ea6121e0595",
   "metadata": {},
   "source": [
    "**Now we will dive into some of the sharp bits of JAX**. Let's try to reuse the above jitted functions with a simple PLN `class` which acts like a container as we did in the tutorial [JIT compilation in pytorch](https://stateofther.github.io/finistR2023/jit-example-pln.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "32220f2c-30b3-4533-85fb-3b52c9b52dbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class PLN_container:\n",
    "    def __init__(self, Y, O, X): \n",
    "        self.Y = Y\n",
    "        self.O = O\n",
    "        self.X = X\n",
    "        self.n, self.p = Y.shape\n",
    "        self.n = self.n\n",
    "        self.p = self.p\n",
    "        self.d = X.shape[1]\n",
    "        ## Variational parameters\n",
    "        self.M = jnp.full(Y.shape, 0.0)\n",
    "        self.S = jnp.full(Y.shape, 1.0)\n",
    "        ## Model parameters\n",
    "        self.B = jnp.zeros((self.d, self.p))\n",
    "        self.Sigma = jnp.eye(self.p)\n",
    "        self.Omega = jnp.eye(self.p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8482744b-31f0-4ebc-9055-f595ea6808b9",
   "metadata": {},
   "source": [
    "Then we would like to adapt the `scan_fun` as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0fd7b8c0-3f28-485d-aab8-432b835fdf06",
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaxfit3(pln, N_iter, lr, tol = 1e-8) :\n",
    "    optimizer = optax.chain(\n",
    "        #adam(learning_rate=lr)\n",
    "        optax.scale_by_radam(),\n",
    "        optax.scale(-1.0),\n",
    "        optax.clip(0.1),\n",
    "    )\n",
    "    opt_state = optimizer.init({'B':pln.B, 'M':pln.M, 'S':pln.S})\n",
    "\n",
    "    def scan_fun(carry, _):\n",
    "        pln = carry['PLN']\n",
    "        loss_value, grads = jax.value_and_grad(get_ELBO)(\n",
    "            {'B':pln.B, 'M':pln.M, 'S':pln.S},\n",
    "            {'X':pln.X, 'O':pln.O, 'n':pln.n, 'Omega':pln.Omega,\n",
    "             'Y':pln.Y, 'p':pln.p, 'Sigma':pln.Sigma}\n",
    "        )\n",
    "    \n",
    "        updates, opt_state = optimizer.update(grads, carry['opt_state'])\n",
    "        updated_params = optax.apply_updates(\n",
    "            {'B':pln.B, 'M':pln.M, 'S':pln.S},\n",
    "            updates)\n",
    "        pln.B = updated_params['B']\n",
    "        pln.M = updated_params['M']\n",
    "        pln.S = updated_params['S']\n",
    "    \n",
    "        ## update parameters with close form\n",
    "        pln.Sigma = get_Sigma(pln.N, pln.M, pln.S)\n",
    "        pln.Omega = jnp.linalg.inv(pln.Sigma)\n",
    "    \n",
    "        carry['PLN'] = pln\n",
    "        carry['opt_state'] = opt_state\n",
    "        return carry, loss_value\n",
    "    \n",
    "    carry, ELBO = jax.lax.scan(\n",
    "        scan_fun,\n",
    "        {\n",
    "            \"PLN\":pln,\n",
    "            'opt_state':opt_state\n",
    "        },\n",
    "        jnp.arange(N_iter)\n",
    "    )\n",
    "\n",
    "    return ELBO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8638618-8f2b-4a75-b3ab-a4027f3621ce",
   "metadata": {},
   "source": [
    "And finally, we would like to start the optimization with:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9152d57d-95f0-470f-8f1c-d1dcdd50894c",
   "metadata": {},
   "source": [
    "```python\n",
    "pln_container = PLN_container(Y, O, X)\n",
    "jaxELBO3 = jax.block_until_ready(\n",
    "    jaxfit3(pln_container, N_iter, lr, tol=1e-8)\n",
    ")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b999685-6e19-4c59-90d3-52cc48f297c1",
   "metadata": {},
   "source": [
    "If we run the previous lines, we get an expected error: `class PLN is not a valid JAX type`. Recall that the `scan` loop automatically JIT compiles its content and we cannot JIT compile function whose arguments are custom classes without specific treatment that we will discover in the next subsection. __But what is a valid JAX type?__ From the [documentation](https://jax.readthedocs.io/en/latest/_autosummary/jax.jit.html) of `jax.jit` we read:\n",
    "\n",
    "> The arguments and return value of fun [the jitted function] should be arrays, scalars, or (nested) standard Python containers (tuple/list/dict) [pytrees] thereof [which themselves contain arrays or scalars].\n",
    "\n",
    "The variety of containers that JAX can handle are refered to with the term __[Pytrees](https://jax.readthedocs.io/en/latest/pytrees.html)__.\n",
    "\n",
    "We also read about __pure function__ and __static arguments__, more on that below !\n",
    "\n",
    "For the moment, as the class PLN acts as a simple container, a common work around can be to use a Python __native and jittable__ `namedtuple` class, or a `NamedTuple` from `typing` library which improves the class creation. See this introduction to `NamedTuple`: [https://www.geeksforgeeks.org/typing-namedtuple-improved-namedtuples/](https://www.geeksforgeeks.org/typing-namedtuple-improved-namedtuples/). NamedTuples are immutable class, so if we want to update one of its attributes we need to use the `_replace()` method which returns a new object!\n",
    "Note that we can use `jax.typing` to provide JAX related type hints [#jax.typing](https://jax.readthedocs.io/en/latest/jax.typing.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5dba1b03-12d6-4ade-8ecd-ae6c89920001",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import NamedTuple\n",
    "from jax.typing import ArrayLike\n",
    "\n",
    "class PLN(NamedTuple):\n",
    "    Y: ArrayLike\n",
    "    O: ArrayLike\n",
    "    X: ArrayLike\n",
    "    n: int\n",
    "    p: int\n",
    "    d: int\n",
    "    M: ArrayLike\n",
    "    S: ArrayLike\n",
    "    B: ArrayLike\n",
    "    Sigma: ArrayLike\n",
    "    Omega: ArrayLike\n",
    "    \n",
    "pln_namedtuple = PLN(\n",
    "    Y, O, X,\n",
    "    Y.shape[0], Y.shape[1], X.shape[1],\n",
    "    jnp.full(Y.shape, 0.0), jnp.full(Y.shape, 1.0),\n",
    "    jnp.zeros((X.shape[1], Y.shape[1])), jnp.eye(Y.shape[1]),\n",
    "    jnp.eye(Y.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6a385f8-4f5a-4f4a-a1da-fb1e0f79306d",
   "metadata": {},
   "source": [
    "Let's adapt the `scan_fun` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "41a092c8-6809-4b14-983d-b082d9a32c1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaxfit3(pln, N_iter, lr, tol = 1e-8) :\n",
    "    optimizer = optax.chain(\n",
    "        #adam(learning_rate=lr)\n",
    "        optax.scale_by_radam(),\n",
    "        optax.scale(-1.0),\n",
    "        optax.clip(0.1),\n",
    "    )\n",
    "    opt_state = optimizer.init({'B':pln.B, 'M':pln.M, 'S':pln.S})\n",
    "\n",
    "    def scan_fun(carry, _):\n",
    "        pln = carry['PLN']\n",
    "        loss_value, grads = jax.value_and_grad(get_ELBO)(\n",
    "            {'B':pln.B, 'M':pln.M, 'S':pln.S},\n",
    "            {'X':pln.X, 'O':pln.O, 'n':pln.n, 'Omega':pln.Omega,\n",
    "             'Y':pln.Y, 'p':pln.p, 'Sigma':pln.Sigma}\n",
    "        )\n",
    "    \n",
    "        updates, opt_state = optimizer.update(grads, carry['opt_state'])\n",
    "        updated_params = optax.apply_updates(\n",
    "            {'B':pln.B, 'M':pln.M, 'S':pln.S},\n",
    "            updates)\n",
    "        pln = pln._replace(B=updated_params['B'])\n",
    "        pln = pln._replace(M=updated_params['M'])\n",
    "        pln = pln._replace(S=updated_params['S'])\n",
    "    \n",
    "    \n",
    "        ## update parameters with close form\n",
    "        pln = pln._replace(Sigma=get_Sigma(pln.n, pln.M, pln.S))\n",
    "        pln = pln._replace(Omega=jnp.linalg.inv(pln.Sigma))\n",
    "    \n",
    "        carry['PLN'] = pln\n",
    "        carry['opt_state'] = opt_state\n",
    "        return carry, loss_value\n",
    " \n",
    "    carry, ELBO = jax.lax.scan(\n",
    "        scan_fun,\n",
    "        {\n",
    "            \"PLN\":pln,\n",
    "            'opt_state':opt_state\n",
    "        },\n",
    "        jnp.arange(N_iter)\n",
    "    )\n",
    "\n",
    "    return ELBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b0872a11-7d6f-4a03-b061-2664eeeb4e03",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.34 s, sys: 831 ms, total: 3.17 s\n",
      "Wall time: 2.91 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "jaxELBO3 = jax.block_until_ready(\n",
    "    jaxfit3(pln_namedtuple, N_iter, lr, tol=1e-8)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f20d9780-e348-4888-9aee-88401e14bee8",
   "metadata": {},
   "source": [
    "Ok, our code works but it is getting a little bit messy: we would like more than containers but real classes, with real methods! That's what we will see in the next section.\n",
    "\n",
    "First, let's try to summarize the **main points of JIT compiling functions with JAX**:\n",
    "\n",
    "- _Not all JAX code can be JIT compiled, as **it requires array shapes to be static & known at compile time**._ [#to-jit-or-not-to-jit](https://jax.readthedocs.io/en/latest/notebooks/thinking_in_jax.html#to-jit-or-not-to-jit)\n",
    "- _JIT and other JAX transforms work by tracing a function to determine its effect on inputs of a specific shape and type. **Variables that you don’t want _[or that cannot be traced because they are not of valid JAX type]_ to be traced can be marked as static**_ [#jit-mechanics-tracing-and-static-variables](https://jax.readthedocs.io/en/latest/notebooks/thinking_in_jax.html#jit-mechanics-tracing-and-static-variables)\n",
    "- Recall that JIT compilation in pytorch cannot handle dynamic shapes, control flows, etc. since it is lost when constructing the Intermediate Representation (IR) of the computational graph that is ued for jitting. The same behaviour happens in JAX when constructing the `jaxpr` which is the name of the IR in JAX, see [#why-can-t-we-just-jit-everything](https://jax.readthedocs.io/en/latest/jax-101/02-jitting.html#why-can-t-we-just-jit-everything). However, in JAX, **there exists way to not lose all the control flows**, to conserve conditional statements even inside the `jaxpr`, etc. Let's mention the `static_argnums` decorator, the special `jax.lax.cond` function, etc. To find more about control flows and jit, see [#python-control-flow-jit](https://jax.readthedocs.io/en/latest/notebooks/Common_Gotchas_in_JAX.html#python-control-flow-jit)\n",
    "- In fact, JAX is asked to JIT compile a function, it will _trace_ its arguments to create the jaxpr, more precisely [#jitting.html](https://jax.readthedocs.io/en/latest/jax-101/02-jitting.html):\n",
    "\n",
    "> When tracing, JAX wraps each argument by a tracer object. These tracers then record all JAX operations performed on them during the function call (which happens in regular Python). Then, JAX uses the tracer records to reconstruct the entire function. The output of that reconstruction is the jaxpr\n",
    "\n",
    "> The more specific information about the values we use in the trace, the more we can use standard Python control flow to express ourselves. However, being too specific means we can’t reuse the same traced function for other values. JAX solves this by tracing at different levels of abstraction for different purposes. For `jax.jit`, the default level is `ShapedArray`\n",
    "**We understand that the arguments of a jitted function should be traceable as `ShapedArray` class, or be a Pytree with either 1.Pytree nodes or 2.leaf nodes being either traceable as a `ShapedArray` object or leaf nodes treated as `static` arguments.** We saw above that scalars and arrays were the valid JAX types, it makes sense since scalars or arrays can be traced with ShapedArray objects.\n",
    "\n",
    "- **Note on static arguments:**\n",
    "\n",
    "> We can tell JAX to help itself to a less abstract tracer for a particular input by specifying static_argnums or static_argnames. The cost of this is that the resulting jaxpr is less flexible, so **JAX will have to re-compile the function for every new value of the specified static input**. It is only a good strategy if the function is guaranteed to get limited different values.\n",
    "\n",
    "- Last but not least, JAX requires the jitted functions to be pure, they cannot have any side effect [#pure-functions](https://jax.readthedocs.io/en/latest/notebooks/Common_Gotchas_in_JAX.html#pure-functions):\n",
    "\n",
    "> JAX transformation and compilation are designed to work only on Python functions that are **functionally pure: all the input data is passed through the function parameters, all the results are output through the function results**. A pure function will always return the same result if invoked with the same inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a99d2558-29fe-4b40-afb9-c3711cfd4be9",
   "metadata": {},
   "source": [
    "**JIT with JAX level 3.5: _dataclasses_ and _pytrees_**\n",
    "\n",
    "Python developpers also love to use `dataclasses` as simple data structures. However they are not supported by default in JAX. Many libraries offer ways to also add `dataclasses` to the Pytrees JAX can handle: [simple-pytree](https://github.com/cgarciae/simple-pytree), [tjax](https://github.com/NeilGirdhar/tjax), [jax_dataclasses](https://github.com/brentyi/jax_dataclasses), [chex dataclasses](https://github.com/deepmind/chex#dataclass-dataclasspy) etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b43b0530-261a-49ea-9980-10b8df2fcf6c",
   "metadata": {},
   "source": [
    "### JIT with JAX level 4: JIT on class methods\n",
    "In this section, we study the canonical way to JIT compile class methods, a introduction to this technique can be found here: [#how-to-use-jit-with-methods](https://jax.readthedocs.io/en/latest/faq.html#how-to-use-jit-with-methods). Note that the recipe we give here should be seen from the viewpoint of Pytrees, since a class is always somehow a container. We therefore want to [#extending-pytrees](https://jax.readthedocs.io/en/latest/pytrees.html). The main elements of making custom class jittable are:\n",
    "\n",
    "1. Decorate the class with `@register_pytree_node_class`\n",
    "2. Add the `tree_flatten` method\n",
    "3. Add the `tree_unflatten` method\n",
    "\n",
    "We here can draw a link with the elements we noted in the previous section: **we make our custom class a Pytree, so that it can be jitted provided it contains attributes being 1.Pytree nodes or 2.leaf nodes being either inherited `ShapedArray` object or leaf nodes treated as `static` arguments!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c7beb6be-fa55-4a42-8bc8-fb02454dc033",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax.tree_util import register_pytree_node_class\n",
    "from jax.typing import ArrayLike\n",
    "\n",
    "@register_pytree_node_class\n",
    "class PLN:\n",
    "    Y: ArrayLike\n",
    "    O: ArrayLike\n",
    "    X: ArrayLike\n",
    "    n: int\n",
    "    p: int\n",
    "    d: int\n",
    "    M: ArrayLike\n",
    "    S: ArrayLike\n",
    "    B: ArrayLike\n",
    "    Sigma: ArrayLike\n",
    "    Omega: ArrayLike\n",
    "    \n",
    "    def __init__(self, Y, O, X): \n",
    "        self.Y = Y\n",
    "        self.O = O\n",
    "        self.X = X\n",
    "        self.n, self.p = Y.shape\n",
    "        self.n = self.n\n",
    "        self.p = self.p\n",
    "        self.d = X.shape[1]\n",
    "        ## Variational parameters\n",
    "        self.M = jnp.full(Y.shape, 0.0)\n",
    "        self.S = jnp.full(Y.shape, 1.0)\n",
    "        ## Model parameters\n",
    "        self.B = jnp.zeros((self.d, self.p))\n",
    "        self.Sigma = jnp.eye(self.p)\n",
    "        self.Omega = jnp.eye(self.p)\n",
    "\n",
    "    def get_optim_params(self):\n",
    "        return {'B':self.B, 'M':self.M, 'S':self.S}\n",
    "        \n",
    "    def set_optim_params(self, op):\n",
    "        for k, v in op.items():\n",
    "            vars(self)[k] = v\n",
    "\n",
    "    def get_other_params(self):\n",
    "        return {\n",
    "            'X':self.X, 'O':self.O, 'n':self.n, 'Omega':self.Omega,\n",
    "            'Y':self.Y, 'p':self.p, 'Sigma':self.Sigma\n",
    "        }\n",
    "        \n",
    "    def set_other_params(self, op):\n",
    "        for k, v in op.items():\n",
    "            vars(self)[k] = v\n",
    "\n",
    "    def update_Sigma(self):\n",
    "        pln.Sigma = 1 / self.n * (self.M.T @ self.M + jnp.diag(jnp.sum(self.S ** 2, axis = 0)))\n",
    "\n",
    "    def update_Omega(self):\n",
    "        pln.Omega = jnp.linalg.inv(self.Sigma)\n",
    "\n",
    "    def tree_flatten(self):\n",
    "        children = (\n",
    "            self.B, self.M, self.S, self.Sigma, self.Omega,\n",
    "        )  # arrays / dynamic values\n",
    "        aux_data = {\n",
    "            'Y':self.Y, 'O':self.O, 'X':self.X,\n",
    "            'n':self.n, 'p':self.p, 'd':self.d\n",
    "        }  # static values\n",
    "        return (children, aux_data)\n",
    "\n",
    "    @classmethod\n",
    "    def tree_unflatten(cls, aux_data, children):\n",
    "        (B, M, S, Sigma, Omega) = children\n",
    "        obj = cls(\n",
    "            aux_data['Y'], aux_data['O'], aux_data['X']\n",
    "        )\n",
    "        obj.B = B\n",
    "        obj.M = M\n",
    "        obj.S = S\n",
    "        return obj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8782f509-3b4a-4f5a-8cda-ff6fa8dffa54",
   "metadata": {},
   "source": [
    "Let's adapt the `scan_fun` function, we are now free to use our new class methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b6672daf-f385-4899-82f6-778f5aa6134d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaxfit4(pln, N_iter, lr, tol = 1e-8) :\n",
    "    optimizer = optax.chain(\n",
    "        #adam(learning_rate=lr)\n",
    "        optax.scale_by_radam(),\n",
    "        optax.scale(-1.0),\n",
    "        optax.clip(0.1),\n",
    "    )\n",
    "    opt_state = optimizer.init(pln.get_optim_params())\n",
    "\n",
    "    def scan_fun(carry, k):\n",
    "        pln = carry['PLN']\n",
    "        loss_value, grads = jax.value_and_grad(get_ELBO)(\n",
    "            pln.get_optim_params(),\n",
    "            pln.get_other_params()\n",
    "        )\n",
    "    \n",
    "        updates, opt_state = optimizer.update(grads, carry['opt_state'])\n",
    "        optim_params = optax.apply_updates(\n",
    "            pln.get_optim_params(),\n",
    "            updates\n",
    "        )\n",
    "        pln.set_optim_params(optim_params)\n",
    "    \n",
    "        ## update parameters with close form\n",
    "        pln.update_Sigma()\n",
    "        pln.update_Omega()\n",
    "    \n",
    "        carry['PLN'] = pln\n",
    "        carry['opt_state'] = opt_state\n",
    "        return carry, loss_value\n",
    " \n",
    "    carry, ELBO = jax.lax.scan(\n",
    "        scan_fun,\n",
    "        {\n",
    "            \"PLN\":pln,\n",
    "            'opt_state':opt_state\n",
    "        },\n",
    "        jnp.arange(N_iter)\n",
    "    )\n",
    "\n",
    "    return ELBO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a8568fbf-173b-4d08-9bda-fcee5966abbe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.29 s, sys: 901 ms, total: 7.2 s\n",
      "Wall time: 4.34 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "pln = PLN(Y, O, X)\n",
    "jaxELBO4 = jax.block_until_ready(\n",
    "    jaxfit4(pln, N_iter, lr, tol=1e-8)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "949a28de-76ba-4721-ab7c-022e311aafed",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90cba7fa-862b-4268-b8ea-87fba86f8dbb",
   "metadata": {},
   "source": [
    "We check the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0be5ef9c-6830-4a08-ab44-01ca71773fba",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_54173/3050985271.py:1: RuntimeWarning: invalid value encountered in log\n",
      "  plt.plot(np.log(jaxELBO_no_jit), label='no jit')\n",
      "/tmp/ipykernel_54173/3050985271.py:2: RuntimeWarning: invalid value encountered in log\n",
      "  plt.plot(np.log(jaxELBO1), label='jit level 1')\n",
      "/tmp/ipykernel_54173/3050985271.py:3: RuntimeWarning: invalid value encountered in log\n",
      "  plt.plot(np.log(jaxELBO2), label='jit level 2')\n",
      "/tmp/ipykernel_54173/3050985271.py:4: RuntimeWarning: invalid value encountered in log\n",
      "  plt.plot(np.log(jaxELBO3), label='jit level 3')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABX8UlEQVR4nO3de3xT9f0/8NdJ2lza5tJ7KW1pBQYIIggIhQFVOijjt9GBXy+/boDiHNp+Admc4oYXHKu3zalzOHWCDhj98pUC47GfgxVbZBQcYIGKcrNAgaZcm/SapMn5/ZHm0NCWJmmaU8rr+XjkQZtzyecclbx8fy5HEEVRBBEREVEPppC7AURERESdYWAhIiKiHo+BhYiIiHo8BhYiIiLq8RhYiIiIqMdjYCEiIqIej4GFiIiIejwGFiIiIurxQuRuQKA4nU6cP38eOp0OgiDI3RwiIiLygiiKqK2tRWJiIhSKjusovSawnD9/HsnJyXI3g4iIiPxQWVmJpKSkDrf3msCi0+kAuC5Yr9fL3BoiIiLyhsViQXJysvQ93pFeE1jc3UB6vZ6BhYiI6CbT2XAODrolIiKiHo+BhYiIiHo8BhYiIiLq8Xwaw5Kfn4+NGzfim2++gVarxfjx4/HKK69g0KBBAIArV67g+eefx7Zt23DmzBnExsYiOzsbL730EgwGQ4fnFUURzz//PN5//33U1NRgwoQJWLlyJQYOHNi1qyMiol7B4XDAbrfL3Qzyg1KpREhISJeXHPEpsJSUlCA3NxdjxoxBc3Mznn32WUydOhVHjhxBeHg4zp8/j/Pnz+P111/H7bffjtOnT2PBggU4f/48/vd//7fD87766qt466238NFHHyEtLQ3Lli3DtGnTcOTIEWg0mi5dIBER3dzq6upw9uxZiKIod1PIT2FhYejTpw9UKpXf5xDELvwbcPHiRcTFxaGkpASTJk1qd58NGzbgxz/+Merr6xES0jYfiaKIxMRE/PznP8cvfvELAIDZbEZ8fDxWr16NBx980Ku2WCwWGAwGmM1mzhIiIuolHA4Hjh8/jrCwMMTGxnJh0JuMKIqw2Wy4ePEiHA4HBg4c2GZxOG+/v7s0rdlsNgMAoqKibriPXq9vN6wAQEVFBUwmEzIzM6X3DAYDxo4di9LS0g4Di9VqhdVqlX63WCz+XAIREfVgdrsdoigiNjYWWq1W7uaQH7RaLUJDQ3H69GnYbDa/e078HnTrdDqxePFiTJgwAcOGDWt3n0uXLuGll17CY4891uF5TCYTACA+Pt7j/fj4eGlbe/Lz82EwGKQXV7klIuq9WFm5ud1oyX2vz+Hvgbm5uSgvL8f69evb3W6xWDBjxgzcfvvteOGFF/z9mA4tXboUZrNZelVWVgb8M4iIiKhn8KtLKC8vD1u3bsXOnTvbXfe/trYWWVlZ0Ol0KCwsRGhoaIfnSkhIAABUV1ejT58+0vvV1dUYMWJEh8ep1Wqo1Wp/mk9EREQ3GZ8qLKIoIi8vD4WFhdixYwfS0tLa7GOxWDB16lSoVCps2bKl076qtLQ0JCQkoKioyOMce/fuRXp6ui/NIyIi6tUEQcCmTZsAAKdOnYIgCCgrK5O1TcHiU2DJzc3FmjVrsG7dOuh0OphMJphMJjQ2NgK4Flbq6+vxl7/8BRaLRdrH4XBI5xk8eDAKCwsBuG7+4sWL8Zvf/AZbtmzB4cOHMWfOHCQmJiI7OztwV0pERHSTq6qqwvTp0wEAycnJqKqqksaRFhcXQxAE1NTUyNjC7uNTl9DKlSsBABkZGR7vr1q1CvPmzcOBAwewd+9eAMCAAQM89qmoqEBqaioA4OjRo9IMIwD45S9/ifr6ejz22GOoqanBd7/7XXz66ac9Yg2WDx9dBtGhwqi56RiRkdn5AURERN3EPYwCcC3I1vr33s6nwNLZki0ZGRleLexz/T6CIGD58uVYvny5L80JCrvzLjSrDTB9ewLIkLs1RES3NlEU0Wh3dL5jN9CGKr2erZSRkYHhw4dDo9Hggw8+gEqlwoIFCzwmoZw5cwb//d//jaKiIigUCmRlZeHtt99uM2u2NUEQUFhYiOzsbJw6dQppaWn48ssvYTQacc899wAAIiMjAQBz587F6tWr/b7enqZL67DcCgS4/sNwcEloIiLZNdoduP25f8ry2UeWT0OYyvuvzY8++ghLlizB3r17UVpainnz5mHChAn43ve+B6fTiZkzZyIiIgIlJSVobm5Gbm4uHnjgARQXF/vctuTkZHzyySeYPXs2jh49Cr1e3+vWrWFg6YQgOgEAzTabzC0hIqKbyfDhw/H8888DAAYOHIg//vGPKCoqwve+9z0UFRXh8OHDqKiokNYR+/jjjzF06FD85z//wZgxY3z6LKVSKS3iGhcXB6PRGNBr6QkYWDojuiosTplKkEREdI02VIkjy6fJ9tm+GD58uMfvffr0wYULFwAAX3/9NZKTkz0WPb399tthNBrx9ddf+xxYbgUMLJ0Q4KqwsEuIiEh+giD41C0jp+vXIBMEAU6nU6bW3Py6vlZub+eusDTzXzIiIgqMIUOGoLKy0mOV9iNHjqCmpga33367X+d0Pwm59TIivQkDS6fcgaV3/gtARETBl5mZiTvuuAM5OTk4cOAAvvjiC8yZMweTJ0/G6NGj/Tpnv379IAgCtm7diosXL6Kuri7ArZYXA0unXJUVp6Pz6dpERETeEAQBmzdvRmRkJCZNmoTMzEzcdtttKCgo8Pucffv2xYsvvohnnnkG8fHxyMvLC2CL5XdzdATKyD2tWWxulrklRER0s2hvarJ7SX23lJQUbN682etzWq1WAEBERAQAIDU1tc26ZsuWLcOyZct8a+xNgoGlU6ywEBGRvCwWCzZu3AiFQoHBgwfL3RxZMLB0pmXQrcjAQkREMnn++eexbt06vPLKK0hKSpK7ObJgYOmUq8LCwEJERHJ544038MYbb8jdDFlx0G2nWmYHcVYzERGRbBhYOtXSJcTAQkREJBsGlk4IgqsriIGFiIhIPgwsnWKXEBERkdwYWDrVMuhWFGRuBxER0a2LgaUzQktpxcnAQkRE3ps3bx6ys7MBABkZGVi8eLHP5xAEoc2Cc90tNTUVf/jDH4L6md5gYOlUS2BhhYWIiHzw5ptvYvXq1QCAjRs34qWXXpK29dRQ4I2vvvoKs2fPRmpqKgRBCNp1MLB0RmBgISIi3xkMBhiNRgBAVFQUdDqdvA0KkIaGBtx22214+eWXkZCQELTPZWDpjBRYeKuIiMh7HXUJZWRk4PTp03jyySchCAIEwfv/Ia6srMT9998Po9GIqKgozJw5E6dOnQIAbNu2DRqNBjU1NR7HLFq0CPfee6/0+65duzBx4kRotVokJydj4cKFqK+v97oNY8aMwWuvvYYHH3wQarXa6+O6it/CnXFPa2aFhYhIfqII2OrleYmBWfF848aNSEpKwvLly1FVVYWqqiqvjrPb7Zg2bRp0Oh0+//xz/Pvf/0ZERASysrJgs9kwZcoUGI1GfPLJJ9IxDocDBQUFyMnJAQCcPHkSWVlZmD17Ng4dOoSCggLs2rXrpniyM5fm70xLYGGFhYioB7A3AL9NlOeznz0PqMK7fJqoqCgolUrodDqfulQKCgrgdDrxwQcfSFWZVatWwWg0ori4GFOnTsWDDz6IdevWYf78+QCAoqIi1NTUYPbs2QCA/Px85OTkSNWegQMH4q233sLkyZOxcuVKaDSaLl9fd2Fg6Yw7sLAYRUREMjp48CBOnDjRZixMU1MTTp48CQDIycnBuHHjcP78eSQmJmLt2rWYMWOGNJbm4MGDOHToENauXSsdL4oinE4nKioqMGTIkKBdj68YWDrDCgsRUc8RGuaqdMj12TKqq6vDqFGjPMKGW2xsLADX+JL+/ftj/fr1ePzxx1FYWCjNVHKf42c/+xkWLlzY5hwpKSnd1vZAYGDphKBgYCEi6jEEISDdMnJTqVRwOBw+HXPXXXehoKAAcXFx0Ov1He6Xk5ODtWvXIikpCQqFAjNmzPA4x5EjRzBgwAC/2y4Xfgt3RtHmByIioi5JTU3Fzp07ce7cOVy6dMmrY3JychATE4OZM2fi888/R0VFBYqLi7Fw4UKcPXvWY78DBw5gxYoVuO+++zxm8jz99NPYvXs38vLyUFZWhuPHj2Pz5s0+Dbq12WwoKytDWVkZbDYbzp07h7KyMpw4ccL7G+AHfgt3xl1hgVLWZhARUe+xfPlynDp1Cv3795e6czoTFhaGnTt3IiUlBbNmzcKQIUMwf/58NDU1eVRcBgwYgLvvvhuHDh2SZge5DR8+HCUlJTh27BgmTpyIkSNH4rnnnkNiovcDmc+fP4+RI0di5MiRqKqqwuuvv46RI0fi0Ucf9foc/mCXUCeuTY9ntiMiIu9ZrVZEREQAAIqLiz22jRs3DgcPHuz0HOJ1U6kTEhLw0UcfdXrc3r17O9w2ZswYbNu2rcPt7nVdOpKamtqmXcHAb+FOCFJhhbeKiIg619zcjCNHjqC0tBRDhw6Vuzm9Br+FOyFwDAsREfmgvLwco0ePxtChQ7FgwQK5m9NrsEuoE4LS3SfEMSxERNS5ESNGoKGhQe5m9DosG3RCULpukcjAQkREJBufAkt+fj7GjBkDnU6HuLg4ZGdn4+jRox77vPfee8jIyIBer4cgCG0ewtSeF154QXoAlPs1ePBgny6ku1yrsDDbERERycWnb+GSkhLk5uZiz5492L59O+x2O6ZOnerxlMeGhgZkZWXh2Wef9akhQ4cOlR4CVVVVhV27dvl0fHdRhLgrK6ywEBERycWnMSyffvqpx++rV69GXFwc9u/fj0mTJgGA9ECl66dwddqQkBCfHgIVLApWWIiIiGTXpW9hs9kMwPXkya46fvw4EhMTcdtttyEnJwdnzpy54f5WqxUWi8Xj1R0UoS2ZTmCFhYiISC5+Bxan04nFixdjwoQJGDZsWJcaMXbsWKxevRqffvopVq5ciYqKCkycOBG1tbUdHpOfnw+DwSC9kpOTu9SGjihCWgbdMrAQERHJxu/Akpubi/Lycqxfv77LjZg+fTr+67/+C8OHD8e0adPwj3/8AzU1Nfif//mfDo9ZunQpzGaz9KqsrOxyO9qjDFW1/MQuISIi8t68efOQnZ0NAMjIyJCGTPhCEARs2rQpoO3qTGpqKv7whz8E9TO94dc6LHl5edi6dSt27tyJpKSkQLcJRqMR3/nOd274ICW1Wu3xQKfuwgoLERH5480335SWsN+4cSNCQ0OlbampqVi8eLFfIUZu77//Pj7++GOUl5cDAEaNGoXf/va3uPvuu7v1c30qG4iiiLy8PBQWFmLHjh1IS0vrlkbV1dXh5MmT6NOnT7ec3xchLaFIZIWFiIh8YDAYYDQaAbjGeup0OnkbFCDFxcV46KGH8Nlnn6G0tBTJycmYOnUqzp07162f69O3cG5uLtasWYN169ZBp9PBZDLBZDKhsbFR2sdkMnk8Zvrw4cMoKyvDlStXpH2mTJmCP/7xj9Lvv/jFL1BSUoJTp05h9+7d+NGPfgSlUomHHnqoq9fXZaHqlkTMCgsREfmgoy6hjIwMnD59Gk8++aS09pi3Kisrcf/998NoNCIqKgozZ86UHla4bds2aDSaNuufLVq0CPfee6/0+65duzBx4kRotVokJydj4cKFHsuTdGbt2rV44oknMGLECAwePBgffPABnE4nioqKvD6HP3wKLCtXroTZbEZGRgb69OkjvQoKCqR93n33XYwcORI//elPAQCTJk3CyJEjsWXLFmmfkydP4tKlS9LvZ8+exUMPPYRBgwbh/vvvR3R0NPbs2eP1I7e7U4jKNYaFXUJERPITRREN9gZZXoF6QvHGjRuRlJSE5cuXS2uPecNut2PatGnQ6XT4/PPP8e9//xsRERHIysqCzWbDlClTYDQa8cknn0jHOBwOFBQUICcnB4Dr+zcrKwuzZ8/GoUOHUFBQgF27diEvL8/v62loaIDdbg/IjOEb8WkMizf/sF544QW88MILN9zn+kdXB2LgbncJUWsBMLAQEfUEjc2NGLturCyfvff/7kVYaFiXzxMVFQWlUgmdTufT+mMFBQVwOp344IMPpKrMqlWrYDQaUVxcjKlTp+LBBx/EunXrMH/+fABAUVERampqMHv2bACuGbY5OTlStWfgwIF46623MHnyZKxcuRIajcbn63n66aeRmJiIzMxMn4/1BR9+2IlQtbvCwjEsREQkn4MHD+LEiRNtxsI0NTXh5MmTAICcnByMGzcO58+fR2JiItauXYsZM2ZIY2kOHjyIQ4cOYe3atdLxoijC6XSioqICQ4YM8alNL7/8MtavX4/i4mK/wo4vGFg6odK60jQrLERE8tOGaLH3/+6V7bPlVFdXh1GjRnmEDTf3EIoxY8agf//+WL9+PR5//HEUFhZi9erVHuf42c9+hoULF7Y5R0pKik/tef311/Hyyy/jX//6F4YPH+7bxfiBgaUToRo1AAcgKNBssyNEFdrpMURE1D0EQQhIt4zcVCoVHA6HT8fcddddKCgoQFxcHPR6fYf75eTkYO3atUhKSoJCocCMGTM8znHkyBEMGDDA77YDwKuvvooVK1bgn//8J0aPHt2lc3mL/RydCAuPkH5uavR+FDUREVFHUlNTsXPnTpw7d85jEsqN5OTkICYmBjNnzsTnn3+OiooKFBcXY+HChTh79qzHfgcOHMCKFStw3333eaxZ9vTTT2P37t3Iy8tDWVkZjh8/js2bN/s06PaVV17BsmXL8OGHHyI1NVWaMVxXV+f9DfADA0snVNprJcDGuo4fFUBEROSt5cuX49SpU+jfv7/XM2LDwsKwc+dOpKSkYNasWRgyZAjmz5+PpqYmj4rLgAEDcPfdd+PQoUPS7CC34cOHo6SkBMeOHcPEiRMxcuRIPPfcc0hMTPS67StXroTNZsN9993nMWP49ddf9/oc/mCXUCe0EXoAJgBAY32DvI0hIqKbhtVqRUSEq0pfXFzssW3cuHE4ePBgp+e4fnZuQkICPvroo06P27u343E+Y8aMwbZt2zrcfv1MXl+3dxdWWDqhCbvWJWRtYIWFiIhurLm5GUeOHEFpaSmGDh0qd3N6DQaWTmjDIwDRCQCwNjTJ3BoiIurpysvLMXr0aAwdOhQLFiyQuzm9BruEOqEMCYEgOiAKClgb2SVEREQ3NmLECDQ08Psi0Fhh8YLQUmGxW1lhISIikgMDixcE0TVX3tbIwEJERCQHBhYvuAOL3W6TuSVERES3JgYWL7gDi6OJFRYiIiI5MLB4wT2Gpdlml7klREREtyYGFq+0VFjsDCxERERyYGDxgtQlZPftQVVERHTrmjdvHrKzswEAGRkZWLx4sc/nEAQBmzZtCmi7OpOamoo//OEPQf1MbzCweIVdQkRE5Js333wTq1evBgBs3LgRL730krStp4YCb2zcuBGjR4+G0WhEeHg4RowYgb/+9a/d/rlcOM4bLRUWZzMrLERE5B2DwSD9HBUVJWNLAisqKgq/+tWvMHjwYKhUKmzduhUPP/ww4uLiMG3atG77XFZYvCCAgYWIiHzTUZdQRkYGTp8+jSeffBKCIEAQBK/PWVlZifvvvx9GoxFRUVGYOXOm9DDCbdu2QaPRoKamxuOYRYsW4d5775V+37VrFyZOnAitVovk5GQsXLgQ9fX1XrchIyMDP/rRjzBkyBD0798fixYtwvDhw7Fr1y6vz+EPBhZvtMwScjY7ZW4IEdGtTRRFOBsaZHld/+Rkf23cuBFJSUlYvnw5qqqqUFVV5dVxdrsd06ZNg06nw+eff45///vfiIiIQFZWFmw2G6ZMmQKj0YhPPvlEOsbhcKCgoAA5OTkAgJMnTyIrKwuzZ8/GoUOHUFBQgF27diEvL8+vaxFFEUVFRTh69CgmTZrk1zm8xS4hr7gqK6KDgYWISE5iYyOO3jVKls8edGA/hLCwLp8nKioKSqUSOp0OCQkJXh9XUFAAp9OJDz74QKrKrFq1CkajEcXFxZg6dSoefPBBrFu3DvPnzwcAFBUVoaamBrNnzwYA5OfnIycnR6r2DBw4EG+99RYmT56MlStXQqPReNUWs9mMvn37wmq1QqlU4k9/+hO+973v+XAXfMfA4hVWWIiISF4HDx7EiRMnoNPpPN5vamrCyZMnAQA5OTkYN24czp8/j8TERKxduxYzZsyA0WiUznHo0CGsXbtWOl4URTidTlRUVGDIkCFetUWn06GsrAx1dXUoKirCkiVLcNtttyEjIyMg19oeBhavtFRYnIEpBxIRkX8ErRaDDuyX7bPlVFdXh1GjRnmEDbfY2FgAwJgxY9C/f3+sX78ejz/+OAoLC6WZSu5z/OxnP8PChQvbnCMlJcXrtigUCgwYMACA6+nUX3/9NfLz8xlY5Ca0VFhEBwMLEZGcBEEISLeM3FQqFRwO3yZy3HXXXSgoKEBcXBz0en2H++Xk5GDt2rVISkqCQqHAjBkzPM5x5MgRKWwEitPphNVqDeg5r8dBt15xV1hkbgYREfUKqamp2LlzJ86dO4dLly55dUxOTg5iYmIwc+ZMfP7556ioqEBxcTEWLlyIs2fPeux34MABrFixAvfddx/UarW07emnn8bu3buRl5eHsrIyHD9+HJs3b/Zp0G1+fj62b9+Ob7/9Fl9//TV+97vf4a9//St+/OMfe38D/MDA4pWWpMIKCxERBcDy5ctx6tQp9O/fX+rO6UxYWBh27tyJlJQUzJo1C0OGDMH8+fPR1NTkUXEZMGAA7r77bhw6dEiaHeQ2fPhwlJSU4NixY5g4cSJGjhyJ5557DomJiV63vb6+Hk888QSGDh2KCRMm4JNPPsGaNWvw6KOPen0OfwhioOZpycxiscBgMMBsNt+wVOaPD+e8hsawUQgLKcLDf1wR0HMTEVHHmpqaUFFRgbS0NK9nsPQUDz30EJRKJdasWSN3U2R3o3+O3n5/s8LilZZ+Rqf3i/sQEdGtqbm5GUeOHEFpaSmGDh0qd3N6DQYWbwgtRaheUYsiIqLuVF5ejtGjR2Po0KFYsGCB3M3pNThLyCvuheNYYSEiohsbMWIEGhoa5G5Gr8MKizeElkG3IgMLERGRHBhYvNKyDovI20VERCQHn76B8/PzMWbMGOh0OsTFxSE7OxtHjx712Oe9995DRkYG9Ho9BEFo89TIjrzzzjtITU2FRqPB2LFj8cUXX/jStO6lcAUWgRUWIiIiWfgUWEpKSpCbm4s9e/Zg+/btsNvtmDp1qsdjqRsaGpCVlYVnn33W6/MWFBRgyZIleP7553HgwAHceeedmDZtGi5cuOBL87qRa7StyMBCREQkiy6tw3Lx4kXExcWhpKSkzWOli4uLcc899+Dq1avSQ5c6MnbsWIwZMwZ//OMfAbiW+E1OTsZ///d/45lnnvGqLd25Dsu/Pt6Lo7vr0X90BLIevTug5yYioo7dzOuw0DWyr8NiNpsBuB6V7S+bzYb9+/cjMzPzWqMUCmRmZqK0tLTD46xWKywWi8eruxhj4gAAam1ggxARERF5x+/A4nQ6sXjxYkyYMAHDhg3zuwGXLl2Cw+FAfHy8x/vx8fEwmUwdHpefnw+DwSC9kpOT/W5DZwSFqyvIyaX5iYjIS/PmzUN2djYAICMjA4sXL/b5HIIgYNOmTQFtV2dSU1Pxhz/8Iaif6Q2/A0tubi7Ky8uxfv36QLbHa0uXLoXZbJZelZWV3fZZCqXrNvFpzURE5K0333wTq1evBgBs3LgRL730krStp4YCX61fvx6CIEjBrDv5tXBcXl4etm7dip07dyIpKalLDYiJiYFSqUR1dbXH+9XV1UhISOjwOLVa7fEEyu6kkCosfFwzERF5x2AwSD93ZehET3Xq1Cn84he/wMSJE4PyeT5VWERRRF5eHgoLC7Fjxw6kpaV1uQEqlQqjRo1CUVGR9J7T6URRURHS09O7fP5AUCjZJURERL7pqEsoIyMDp0+fxpNPPglBECAI3s9AraysxP333w+j0YioqCjMnDkTp06dAgBs27YNGo2mzXIiixYtwr333iv9vmvXLkycOBFarRbJyclYuHChx2xfbzgcDuTk5ODFF1/Ebbfd5tOx/vIpsOTm5mLNmjVYt24ddDodTCYTTCYTGhsbpX1MJhPKyspw4sQJAMDhw4dRVlaGK1euSPtMmTJFmhEEAEuWLMH777+Pjz76CF9//TUef/xx1NfX4+GHH+7q9QWEFFicDCxERHISRRF2q0OWVxcm1XrYuHEjkpKSsHz5clRVVaGqqsqr4+x2O6ZNmwadTofPP/8c//73vxEREYGsrCzYbDZMmTIFRqMRn3zyiXSMw+FAQUEBcnJyAAAnT55EVlYWZs+ejUOHDqGgoAC7du1CXl6eT9ewfPlyxMXFYf78+T4d1xU+dQmtXLkSgCsdtrZq1SrMmzcPAPDuu+/ixRdflLa5pzu33ufkyZO4dOmStM8DDzyAixcv4rnnnoPJZMKIESPw6aefthmIKxdWWIiIeoZmmxPvLSqR5bMfe3MyQtXKLp8nKioKSqUSOp3uhkMfrldQUACn04kPPvhAqsqsWrUKRqMRxcXFmDp1Kh588EGsW7dOChJFRUWoqanB7NmzAbgmrOTk5EjVnoEDB+Ktt97C5MmTsXLlSq+mju/atQt/+ctfUFZW5tuFd5FPgcWbdPnCCy/ghRdeuOE+7vJVa3l5eT4nvGBxD7plhYWIiORy8OBBnDhxAjqdzuP9pqYmnDx5EgCQk5ODcePG4fz580hMTMTatWsxY8YMaT20gwcP4tChQ1i7dq10vCiKcDqdqKiowJAhQ27YhtraWvzkJz/B+++/j5iYmMBeYCf4tGYvcNAtEVHPEKJS4LE3J8v22XKqq6vDqFGjPMKGW2xsLABgzJgx6N+/P9avX4/HH38chYWF0kwl9zl+9rOfYeHChW3OkZKS0mkbTp48iVOnTuEHP/iB9J7T6fpuDAkJwdGjR9G/f39fL80rDCxeYJcQEVHPIAhCQLpl5KZSqeBwOHw65q677kJBQQHi4uJuuCJsTk4O1q5di6SkJCgUCsyYMcPjHEeOHMGAAQP8avfgwYNx+PBhj/d+/etfo7a2Fm+++Wa3ronGxw97gQvHERFRIKWmpmLnzp04d+6cx5jOG8nJyUFMTAxmzpyJzz//HBUVFSguLsbChQtx9uxZj/0OHDiAFStW4L777vNYAuTpp5/G7t27kZeXh7KyMhw/fhybN2/2ekiGRqPBsGHDPF5GoxE6nQ7Dhg2DSqXy7Ub4gIHFC6ywEBFRIC1fvhynTp1C//79pe6czoSFhWHnzp1ISUnBrFmzMGTIEMyfPx9NTU0eFZcBAwbg7rvvxqFDh6TZQW7Dhw9HSUkJjh07hokTJ2LkyJF47rnnkJiYGNDr6w7sEvKCkoNuiYjIR1arFREREQBcDwRubdy4cTh48GCn57h+sktCQgI++uijTo/bu3dvh9vGjBmDbdu2dbi9vYkxN9J6jEx3YoXFCwIrLERE5KXm5mYcOXIEpaWlGDp0qNzN6TUYWLxwrUuIs4SIiOjGysvLMXr0aAwdOhQLFiyQuzm9BruEvOCe1iyyS4iIiDoxYsQINDQ0yN2MXocVFi9w0C0REZG8GFi8IK10y8BCREQkCwYWLyi4DgsRkawC9eBBkkcg/vkxsHiBT2smIpKHUula1dZms8ncEuoK95ie0NBQv8/BQbde4CwhIiJ5hISEICwsDBcvXkRoaCgUCv5/9s1EFEU0NDTgwoULMBqNUgD1BwOLFzjolohIHoIgoE+fPqioqMDp06flbg75yWg0IiEhoUvnYGDxQutBt6IoQhAEmVtERHTrUKlUGDhwILuFblKhoaFdqqy4MbB4wV1hAVxrsQhKBhYiomBSKBTQaDRyN4NkxM5AL7QOLOwWIiIiCj4GFi8wsBAREcmLgcUL7jEsAAMLERGRHBhYvKBQCEBLkcXBqc1ERERBx8DiJSWX5yciIpINA4uXuHgcERGRfBhYvMTF44iIiOTDwOIlBhYiIiL5MLB4ScExLERERLJhYPGSu8LiaOYYFiIiomBjYPGSMoQVFiIiIrkwsHiJs4SIiIjkw8DiJQ66JSIikg8Di5cUCgYWIiIiuTCweImzhIiIiOTDwOIlaZYQx7AQEREFnU+BJT8/H2PGjIFOp0NcXByys7Nx9OhRj32ampqQm5uL6OhoREREYPbs2aiurr7heefNmwdBEDxeWVlZvl9NN1KGsEuIiIhILj4FlpKSEuTm5mLPnj3Yvn077HY7pk6divr6emmfJ598En//+9+xYcMGlJSU4Pz585g1a1an587KykJVVZX0+tvf/ub71XSja11CrLAQEREFW4gvO3/66acev69evRpxcXHYv38/Jk2aBLPZjL/85S9Yt24d7r33XgDAqlWrMGTIEOzZswfjxo3r8NxqtRoJCQl+XEJwcJYQERGRfLo0hsVsNgMAoqKiAAD79++H3W5HZmamtM/gwYORkpKC0tLSG56ruLgYcXFxGDRoEB5//HFcvnz5hvtbrVZYLBaPV3diYCEiIpKP34HF6XRi8eLFmDBhAoYNGwYAMJlMUKlUMBqNHvvGx8fDZDJ1eK6srCx8/PHHKCoqwiuvvIKSkhJMnz4dDoejw2Py8/NhMBikV3Jysr+X4hV3lxCX5iciIgo+n7qEWsvNzUV5eTl27drV5UY8+OCD0s933HEHhg8fjv79+6O4uBhTpkxp95ilS5diyZIl0u8Wi6VbQwsrLERERPLxq8KSl5eHrVu34rPPPkNSUpL0fkJCAmw2G2pqajz2r66u9ml8ym233YaYmBicOHGiw33UajX0er3HqzspGViIiIhk41NgEUUReXl5KCwsxI4dO5CWluaxfdSoUQgNDUVRUZH03tGjR3HmzBmkp6d7/Tlnz57F5cuX0adPH1+a1604S4iIiEg+PgWW3NxcrFmzBuvWrYNOp4PJZILJZEJjYyMAwGAwYP78+ViyZAk+++wz7N+/Hw8//DDS09M9ZggNHjwYhYWFAIC6ujo89dRT2LNnD06dOoWioiLMnDkTAwYMwLRp0wJ4qV3DLiEiIiL5+DSGZeXKlQCAjIwMj/dXrVqFefPmAQDeeOMNKBQKzJ49G1arFdOmTcOf/vQnj/2PHj0qzTBSKpU4dOgQPvroI9TU1CAxMRFTp07FSy+9BLVa7edlBd61lW4ZWIiIiILNp8Aiip1/WWs0Grzzzjt45513vDqPVqvFP//5T1+aIQt2CREREcmHzxLykoJL8xMREcmGgcVLnCVEREQkHwYWL7FLiIiISD4MLF7iLCEiIiL5MLB4SZol1MzAQkREFGwMLF5ilxAREZF8GFi8JHUJOVlhISIiCjYGFi9xlhAREZF8GFi8xC4hIiIi+TCweEnqEuKgWyIioqBjYPESnyVEREQkHwYWL7FLiIiISD4MLF7is4SIiIjkw8DiJc4SIiIikg8Di5fYJURERCQfBhYv8VlCRERE8mFg8RJnCREREcmHgcVLCgW7hIiIiOTCwOIlzhIiIiKSDwOLlziGhYiISD4MLF5SumcJNbNLiIiIKNgYWLzkrrCIIiA6WWUhIiIKJgYWL7kDC8BuISIiomBjYPGSe+E4AHBwphAREVFQMbB4yT1LCGCFhYiIKNgYWLykUDCwEBERyYWBxUuCIEihhYvHERERBRcDiw+4FgsREZE8GFh8wMBCREQkDwYWHyhCXLeLs4SIiIiCi4HFB6ywEBERyYOBxQdSYGlmYCEiIgomnwJLfn4+xowZA51Oh7i4OGRnZ+Po0aMe+zQ1NSE3NxfR0dGIiIjA7NmzUV1dfcPziqKI5557Dn369IFWq0VmZiaOHz/u+9V0M/ficZwlREREFFw+BZaSkhLk5uZiz5492L59O+x2O6ZOnYr6+nppnyeffBJ///vfsWHDBpSUlOD8+fOYNWvWDc/76quv4q233sK7776LvXv3Ijw8HNOmTUNTU5N/V9VNlOwSIiIikoUgiqLf374XL15EXFwcSkpKMGnSJJjNZsTGxmLdunW47777AADffPMNhgwZgtLSUowbN67NOURRRGJiIn7+85/jF7/4BQDAbDYjPj4eq1evxoMPPuhVWywWCwwGA8xmM/R6vb+XdEPrX9qLy+fq8cOFI5B8e1S3fAYREdGtxNvv7y6NYTGbzQCAqCjXl/f+/ftht9uRmZkp7TN48GCkpKSgtLS03XNUVFTAZDJ5HGMwGDB27NgOjwEAq9UKi8Xi8epu7i4hzhIiIiIKLr8Di9PpxOLFizFhwgQMGzYMAGAymaBSqWA0Gj32jY+Ph8lkavc87vfj4+O9PgZwjacxGAzSKzk52d9L8RpnCREREcnD78CSm5uL8vJyrF+/PpDt8drSpUthNpulV2VlZbd/JgMLERGRPPwKLHl5edi6dSs+++wzJCUlSe8nJCTAZrOhpqbGY//q6mokJCS0ey73+9fPJLrRMQCgVquh1+s9Xt2Ns4SIiIjk4VNgEUUReXl5KCwsxI4dO5CWluaxfdSoUQgNDUVRUZH03tGjR3HmzBmkp6e3e860tDQkJCR4HGOxWLB3794Oj5ELZwkRERHJw6fAkpubizVr1mDdunXQ6XQwmUwwmUxobGwE4BosO3/+fCxZsgSfffYZ9u/fj4cffhjp6ekeM4QGDx6MwsJCAK6nIC9evBi/+c1vsGXLFhw+fBhz5sxBYmIisrOzA3elAcAuISIiInmE+LLzypUrAQAZGRke769atQrz5s0DALzxxhtQKBSYPXs2rFYrpk2bhj/96U8e+x89elSaYQQAv/zlL1FfX4/HHnsMNTU1+O53v4tPP/0UGo3Gj0vqPu5nCbFLiIiIKLi6tA5LTxKMdVj++UE5Tuy7gO/+10DcOaX7ZyURERH1dkFZh+VWwy4hIiIieTCw+ECaJeRklxAREVEwMbD4gBUWIiIieTCw+ECpYGAhIiKSAwOLD9yzhBzN7BIiIiIKJgYWH7BLiIiISB4MLD5gYCEiIpIHA4sP+CwhIiIieTCw+IAVFiIiInkwsPiAgYWIiEgeDCw+ULZ0CTnYJURERBRUDCw+YIWFiIhIHgwsPmBgISIikgcDiw84S4iIiEgeDCw+YIWFiIhIHgwsPnAHFi7NT0REFFwMLD5Qhri7hFhhISIiCiYGFh+wS4iIiEgeDCw+UCpZYSEiIpIDA4sPFCHuCgvHsBAREQUTA4sP3NOaOeiWiIgouBhYfKAM4RgWIiIiOTCw+ECqsDCwEBERBRUDiw+kCgu7hIiIiIKKgcUHrLAQERHJg4HFB9I6LKywEBERBRUDiw9ar3QriqyyEBERBQsDiw/cFRYAcDoZWIiIiIKFgcUH7goLADibGViIiIiChYHFB+6VbgEuHkdERBRMDCw+UChadQlxphAREVHQMLD4QBCEVk9sZoWFiIgoWHwOLDt37sQPfvADJCYmQhAEbNq0yWN7dXU15s2bh8TERISFhSErKwvHjx+/4TlXr14NQRA8XhqNxtemBYUixP08IVZYiIiIgsXnwFJfX48777wT77zzTpttoigiOzsb3377LTZv3owvv/wS/fr1Q2ZmJurr6294Xr1ej6qqKul1+vRpX5sWFEpWWIiIiIIuxNcDpk+fjunTp7e77fjx49izZw/Ky8sxdOhQAMDKlSuRkJCAv/3tb3j00Uc7PK8gCEhISPC1OUGnaLUWCxEREQVHQMewWK1WAPDozlEoFFCr1di1a9cNj62rq0O/fv2QnJyMmTNn4quvvur0sywWi8crGNwVFs4SIiIiCp6ABpbBgwcjJSUFS5cuxdWrV2Gz2fDKK6/g7NmzqKqq6vC4QYMG4cMPP8TmzZuxZs0aOJ1OjB8/HmfPnu3wmPz8fBgMBumVnJwcyEvp0LVBt6ywEBERBUtAA0toaCg2btyIY8eOISoqCmFhYfjss88wffp0KBQdf1R6ejrmzJmDESNGYPLkydi4cSNiY2Px5z//ucNjli5dCrPZLL0qKysDeSkdUkqDbllhISIiChafx7B0ZtSoUSgrK4PZbIbNZkNsbCzGjh2L0aNHe32O0NBQjBw5EidOnOhwH7VaDbVaHYgm+0QZ6goszXYGFiIiomDptnVYDAYDYmNjcfz4cezbtw8zZ870+liHw4HDhw+jT58+3dU8v4WEKgEADhsDCxERUbD4XGGpq6vzqHxUVFSgrKwMUVFRSElJwYYNGxAbG4uUlBQcPnwYixYtQnZ2NqZOnSodM2fOHPTt2xf5+fkAgOXLl2PcuHEYMGAAampq8Nprr+H06dM3nFUklxCVu8LikLklREREtw6fA8u+fftwzz33SL8vWbIEADB37lysXr0aVVVVWLJkCaqrq9GnTx/MmTMHy5Yt8zjHmTNnPMa0XL16FT/96U9hMpkQGRmJUaNGYffu3bj99tv9va5uE+LuEmKFhYiIKGgEURR7xXQXi8UCg8EAs9kMvV7fbZ/z6XvlOHngAiY+MBDD7wnOzCQiIqLeytvvbz5LyEdSlxArLEREREHDwOKjEJVr0C1nCREREQUPA4uPro1h4aBbIiKiYGFg8VEI12EhIiIKOgYWH7nHsDhYYSEiIgoaBhYfKUM5hoWIiCjYGFh8xHVYiIiIgo+BxUfuZwnx4YdERETBw8DiI6VSAAA4nb1ivT0iIqKbAgOLjwR3YHGwwkJERBQsDCw+cj8DyelghYWIiChYGFh8pJAqLAwsREREwcLA4iMGFiIiouBjYPERAwsREVHwMbD4SKFsGcPCWUJERERBw8DiIwVnCREREQUdA4uPFAp2CREREQUbA4uP3BUWkV1CREREQcPA4iP3GBYHKyxERERBw8Dio2tdQhzDQkREFCwMLD6SuoRYYSEiIgoaBhYfcR0WIiKi4GNg8REDCxERUfAxsPio9cJxosjQQkREFAwMLD5yV1gATm0mIiIKFgYWH7UOLOwWIiIiCg4GFh+5pzUDfJ4QERFRsDCw+IgVFiIiouBjYPGRoGBgISIiCjYGFh8JgsAHIBIREQUZA4sfrq3FwuX5iYiIgoGBxQ9cPI6IiCi4fA4sO3fuxA9+8AMkJiZCEARs2rTJY3t1dTXmzZuHxMREhIWFISsrC8ePH+/0vBs2bMDgwYOh0Whwxx134B//+IevTQsaRYj7ic2ssBAREQWDz4Glvr4ed955J955550220RRRHZ2Nr799lts3rwZX375Jfr164fMzEzU19d3eM7du3fjoYcewvz58/Hll18iOzsb2dnZKC8v97V5QaF0BxY7AwsREVEwCGIX1pcXBAGFhYXIzs4GABw7dgyDBg1CeXk5hg4dCgBwOp1ISEjAb3/7Wzz66KPtnueBBx5AfX09tm7dKr03btw4jBgxAu+++65XbbFYLDAYDDCbzdDr9f5eklf++uvdsFxqwqynRqFPf0O3fhYREVFv5u33d0DHsFitVgCARqO59gEKBdRqNXbt2tXhcaWlpcjMzPR4b9q0aSgtLb3hZ1ksFo9XsEgVlmZWWIiIiIIhoIFl8ODBSElJwdKlS3H16lXYbDa88sorOHv2LKqqqjo8zmQyIT4+3uO9+Ph4mEymDo/Jz8+HwWCQXsnJyQG7js4oQxlYiIiIgimggSU0NBQbN27EsWPHEBUVhbCwMHz22WeYPn06FIrATkhaunQpzGaz9KqsrAzo+W+EY1iIiIiCKyTQJxw1ahTKyspgNpths9kQGxuLsWPHYvTo0R0ek5CQgOrqao/3qqurkZCQ0OExarUaarU6YO32BbuEiIiIgqvb1mExGAyIjY3F8ePHsW/fPsycObPDfdPT01FUVOTx3vbt25Gent5dzesSd5eQk4GFiIgoKHyusNTV1eHEiRPS7xUVFSgrK0NUVBRSUlKwYcMGxMbGIiUlBYcPH8aiRYuQnZ2NqVOnSsfMmTMHffv2RX5+PgBg0aJFmDx5Mn73u99hxowZWL9+Pfbt24f33nsvAJcYeNcqLFw4joiIKBh8Diz79u3DPffcI/2+ZMkSAMDcuXOxevVqVFVVYcmSJaiurkafPn0wZ84cLFu2zOMcZ86c8RjTMn78eKxbtw6//vWv8eyzz2LgwIHYtGkThg0b5u91dStliGulW3YJERERBUeX1mHpSYK5Dsv2D7/CsS+qMeG+ARiRmdKtn0VERNSbybIOy62Cg26JiIiCi4HFD5zWTEREFFwMLH5ghYWIiCi4GFj8oAxtGXRr7xXDf4iIiHo8BhY/KFhhISIiCioGFj+E8FlCREREQcXA4geOYSEiIgouBhY/MLAQEREFFwOLH7g0PxERUXAxsPiBS/MTEREFFwOLHxRcOI6IiCioGFj8wDEsREREwcXA4gclpzUTEREFFQOLH/gsISIiouBiYPEDu4SIiIiCK0TuBtyMrs0S6t5pzXarA45mJ0SnCKdTRKhaiVCVEoJC6NbPJSIi6mkYWPzQHRUWp8MJ07dmnD9eg/PHa3Cxsg5Ndfa2OwqASq2EVqdCmEGFcIPa9TKqoY/RQB+jhT5GA3VYaMDaRkREJDcGFj8EetDthdMW7Pj4G1w+V9f5ziJga3LA1tQI88XGDndTh4VAF62BIUYLXYwW+mgN9LEtf0ZrpWsgIiK6GTCw+CFQFRa7zYEvtnyLg0WVEEVXyEgeEoXEgUbEp+lhjAuDUqWAQhAAAWi2O2FrbIa9yYEGiw31ZisazDbU11hRd7UJlstNsFxqRGOtHdaGZlgb6nCpsm0IEgQgIkoDY5wWhrgwGOPCYIjTwhgXBl2MBkolwwwREfUsDCx+cAcWZ7MIURQhCL6PKblqqsfWdw7B0lIlGTgmHhPvHwitTtXhMaEq1xgWGABjfFiH+9mtDlguNUoBpvZSE8yXGlF7uRGWS02wWx2ovdyE2stNqPz6qsexgkKALtoVZlxBJkwKNrpoDRQcP0NERDJgYPFD6+4UZ7MIZahvX+L1Zis2vfElGsw2RESqMfn/DkLqHTEBa1+oWonovhGI7hvRZpsoimistaPmQgPMFxpQc6HR489mmxOWi42wXGzEma+ueByrUArQx2hbVWa0MCaEIzI+DGEGlV/BjYiIyBsMLH4ICbkWWJqbnT6PB9mz6SQazDZEJoThRz+/64ZVlUATBAFhehXC9CokDjB6bBNFEfU1tpYA0wDzhUbXnxcbYb7QCEezEzXVDaipbgBw2ePYUI0SkfFhMCaEITI+HJEJrp+NsWEcL0NERF3GwOIHRYhrTAlEoNnmgFrr/W2sqW7AN6UmAMC9c4YENax0RhAERESqERGpRt9BkR7bRKeIuhqrZ5CpbsDV6gZYLjbC3uTAhdO1uHC69rpzAroYrSvAxIchMj4MkQnhMMaHQasLZVWGiIi8wsDiB0EQEBKiQLPd6fNqt1/vrgIA9BsWjYTbDN3RvG4hKAToojTQRWmQPNhzm8PuhPliI65W16OmugFXTa5XjaketiaH1MV0+rBnVUYdFnItyLSEmMiEMOhjtRz4S0REHhhY/KRUuQJLs837wOJ0OPHNHldgGTK+T3c1LeiUoQpEJYYjKjHc431RFNFgsaHG5KrEuP6sx1VTA2qvNMHa0AzTtxaYvrV4HKdQCNDHuqoyUX3CEdknHFF9wmFMCHMNOiYiolsOA4ufQlVKWOub0Wx3eH3MmSNX0GC2QRMeitThgRtk21MJgiAtbHd9F1OzzYGaC424arpWlXH9WY9m27WxMhUHL7U6IaCP1iAywR1iwqSffemWIyKimw//lveTeyBpsw9dQt+Uuqor3xkbL02NvlWFqJSISYpATJLnTCbXwF8rrlY14IqpHlerXBWZK+fr0VRvh+VSEyyXmnC63LN7Kdyovq4iE4bIPuHQRvScMUJEROQ/BhY/hYS6uiYcXnYJNdbZpGpBb+oOCjTXwF8NIiI1SL49ymNbY60NV031uFLVgCtVLWGmqh71LYvn1ddYcfYbz3VltLrQNhWZqD7hnIZNRHSTYWDxU4jKXWHxrkvo2BfVcDpExCRHICZJ151N67W0OhW0OhUSB3p2L1kb7K4qTFWrikxVPWovN6Gx1o7GWtfzmVpTaUPajJGJTAiDLkrDh0sSEfVADCx+CvGxS8jdHTRkfGK3telWpQ4LRcJthjazruxWB2qqXeGldZgxX2iArbEZ1RUWVFd4DvgNUSlaKjItYaalIqOP1XKVXyIiGTGw+EnZ0iXUbOu8wnKxshaXKuugCBHwnbvju7tp1CJUrURsig6xKZ4VLYfdiZoLbSsyNdWulX4vnqnFxTOe68koQxQwxodJY2PclRlDHKdgExEFAwOLn6QuIS/GsBz7ohoAkDY8Bprw0G5tF3VOGapo99EFTocTlktN1yoypnpcrWrA1ap6NNuduHyurs0TtRUKAYY4rWfXUp9wGOO10jgnIiLqOp8Dy86dO/Haa69h//79qKqqQmFhIbKzs6XtdXV1eOaZZ7Bp0yZcvnwZaWlpWLhwIRYsWNDhOVevXo2HH37Y4z21Wo2mpiZfmxc03nYJiU4Rx//jCizfGZPQ7e0i/ymUriqKMT4Mt42Ild4XnSJqrzR5dC1daQkydqtDWigPX16UjhEEtKwl4wox7spMZEI4QtUMMkREvvI5sNTX1+POO+/EI488glmzZrXZvmTJEuzYsQNr1qxBamoqtm3bhieeeAKJiYn44Q9/2OF59Xo9jh49Kv3e02dwuAOLo5NBt+dP1KC+xgqVNgQpw6JuuC/1TILC9dBHfYzW4yGVoiii7qrVo1vpakuosTY0w3zB9QymU4cueZxPF61pM/06KiEcKq4lQ0TUIZ//hpw+fTqmT5/e4fbdu3dj7ty5yMjIAAA89thj+POf/4wvvvjihoFFEAQkJNw8FQilyj2G5cYVlmMt1ZX+I2PZRdDLCMK1xxWkDI2W3nev8Nu6EuPuYmqstaP2chNqL7e/lsz1Y2Si+oSzG5GICN0whmX8+PHYsmULHnnkESQmJqK4uBjHjh3DG2+8ccPj6urq0K9fPzidTtx111347W9/i6FDh3a4v9VqhdVqlX63WCwd7tsdvOkScjicOLn/AgBgIAfb3jJar/CbNPi6tWTqbK5F8VpVY65fS6by6+vWktGrPNaQcQcZPjySiG4lAQ8sb7/9Nh577DEkJSUhJCQECoUC77//PiZNmtThMYMGDcKHH36I4cOHw2w24/XXX8f48ePx1VdfISkpqd1j8vPz8eKLLwa6+V7zJrBUf2uBtaEZmohQ9P1OZIf70a1DG6GCdqAKiQONHu+3Xkvm2qJ4rmcuNVpsOGex4dzRGo9j1OEh1wJMqzATbuSieETU+3RLYNmzZw+2bNmCfv36YefOncjNzUViYiIyMzPbPSY9PR3p6enS7+PHj8eQIUPw5z//GS+99FK7xyxduhRLliyRfrdYLEhOTg7sxdxAiMq90m3HY1gqv7kCAEgeHMk1POiGOlpLxtbULK0l4+5iulJVD8ulRljrm1F1woyqE2aPY1QapWuAb0uQca8pw0XxiOhmFtDA0tjYiGeffRaFhYWYMWMGAGD48OEoKyvD66+/3mFguV5oaChGjhyJEydOdLiPWq2GWq0OSLv9cW2l244rLGdbSvtJQzjYlvyj0oQgrp8ecf30Hu832xy4Wt3QZsBvzYVG2JocN1wUz1WJCeOieER0UwloYLHb7bDb7VAoPBfSUiqVcDq9f0igw+HA4cOH8f3vfz+QzQso9wDajgJLs92BC6dcXxhJg9gdRIEVolIiNlmH2OTrFsVrdi2Kd/04GS6KR0Q3O58DS11dnUflo6KiAmVlZYiKikJKSgomT56Mp556ClqtFv369UNJSQk+/vhj/P73v5eOmTNnDvr27Yv8/HwAwPLlyzFu3DgMGDAANTU1eO2113D69Gk8+uijAbjE7iE9rbmDLqGrVQ1wOkVowkOhi9YEs2l0C1OGKBCdGIHoxLaL4pkvNraZfn3V1AAHF8UjopuAz4Fl3759uOeee6Tf3eNI5s6di9WrV2P9+vVYunQpcnJycOXKFfTr1w8rVqzwWDjuzJkzHlWYq1ev4qc//SlMJhMiIyMxatQo7N69G7fffntXrq1bdbbS7aWzrr/8o5PCOQCSZKdQtjwjKSHcY1E8p1NE7eUmjxlL7iDT2aJ41wb8clE8Iup+giiKotyNCASLxQKDwQCz2Qy9Xt/5AV10uvwytv7xIGKSI/DAr+5us33X/x7HwX9VYvg9SZj4wHe6vT1EgdR6UTyP1X1NrkXxOuJeFM+YEIbI+LCWkBQGTQSnYBNR+7z9/ubSmn66ttJt+xWWy1KFJaLd7UQ9WWeL4rWeeu3Nonjq8BBExrvCizGhJcjEh0Efo4GC42SIyAsMLH5SdtIldKWqHgAQlRgetDYRdbfWi+IlX78oXq0NV03XKjE1pgZcrW5A7eUmWOubYfrWDNO3nlOwFUoBhpZnLrmCTBgi410/q/moAiJqhX8j+OnaLKG2g26dDicaLDYAgD5aG9R2EclFq1NBq1MhcaDnrDi7zQHzhQZpPEyNqR5XqxtQY2pAs915bZzMdcIMKldFJj68pXvJVZ3RRXI9GaJbEQOLn2600m2DxQ6IrofmaSP4HBi6tYWqlIhJ0iEmyXMKtugUUXu1SarEuP50DfhtMNuk1/Ur/IaEKqQxMsaWMTLup2yHqjjol6i3YmDxU6jG9Rej3eqAKIoeAwobLK5nHIXpQvl/gkQdEBQC9NFa6KO1HuNkAMDa2OwRYGpMrm4m88VGNNuduFRZh0uVdW3OqYvStBknY0wIQ5iejysgutkxsPhJpWm5daIrtEi/A6ivcQWWcKN8K/ES3czU2hDEp+kRn+Y5Y8DpcMJyqcm1yq97nIzp2uyl2itNqL3ShDNHrngcp9IoYYwPgyHOXY3Rwtjyc+v/domo5+J/qX4KUSkgCIDYXmAxu8avhBkYWIgCSaFUSN0/acNjpPdFUURTnV0KLzXVDS2hpgG1l1yPK7hwuhYXTte2OWeYQSWFF2OcqyJjjNNCH6OFMoQzmIh6CgYWPwmCgFC1ErYmB+xNDqDVM+vqzS0VFoNKptYR3VoEQWg16Nfosa3Z7oDlYhNqqhtQc6HB48/GWrs0Vub88RrPcyoE6KM1LQHGHWi0MMaHIdyoZhcTUZAxsHSBShsCW5MDtibPhbQaWGEh6jFCQpWISgxvd4mBpno7zBca2wQZ97OXzBcbYb7YiNPwXFcmRKW4VpGRgozrsQXqMA60J+oODCxd4F6G3NbkObWZFRaim4MmPBSatNA2Y2VEUUR9jc0jwLh/tlxqQrOt44G/Wl3otSDTEmoM8VoYYvkMJqKuYGDpgtCWcSv2Dios4aywEN2UBEFARKQaEZHqNk9bdzicqG0Z+OsOMuaWn+vNNjTW2tFYa0bVSfN1JwX00RoY4sJgiHUN+jXEuYIMx8sQdY6BpQtUmk4qLJwlRNTrKFsN/L2eranZ1cV0XfdSTXUDbE0OWC41wXKpCZXXHScIrucweYSZWC0MHPxLJGFg6QJVOxUWp1NEo8U9hoVdQkS3EpUmBLEpOsSmXLdIniiisdaOmup61FxohPlCI8wXG1w/X2xEs9WLMBOrdU3LZpihWxQDSxeEtlNhaay1QRRdf8lodQwsROTqYgrTqxCmb/voAvcDJc0XGjoPM19fve68nmGmdVcTwwz1NgwsXSBVWKzXAot70TitXgUFV7klok60fqBkx2HGNZPJfLHxWrBhmKFbDANLF0gVlsZrXUIccEtEgeIZZowe2zoKM+aLjai50HmYiYjSwBinhSH2WogxxGqhj9XymUzUIzGwdIE06LZ1hYVTmokoCLwNM1L3Uqufm60O1F5uQu3ltmEGcI2/M8RqYYhxBRh3kDHEaqEJD+WieSQLBpYuaG/QrbQsP2cIEZFMfA0zlpYF8iyXGmFtaJZW/606YW5zbpVGKYUX95Rsd6CJiNSwK5y6DQNLF7Q36FaqsOhZYSGinudGYQZoWf334rUQY7507ef6GitsTY4OF81TKAXoY66FGKkyE6OFPkaDEHY1URcwsHSBSt22wiKNYWGFhYhuQprwUGjCQxGfqm+zrdnmGhfTOsS4KzOWy41wNovSujPtCTeqPULM9V1NRDfCwNIFodp2Kiwts4T4HCEi6m1CVB0/l8npFFFfY5UG/loutQo0F11PzK6vsaK+xtrmQZMAoA4L8ehecnc36aM1iIhUQ6HkrKZbHQNLF6i1rttnbWhdYeGgWyK69SgUAnRRGuiiNEga7LlNFMU2XU3uQGO52Ih6sw3WhmZcPFOLi2dq2z13RJTaFWBiXSFG6nqK0UIdHsKBwLcABpYu0ES4SphNdXaIoghRBBosnNZMRNSaIAjQRqigjVAhIc3QZrvd5nB1K7WqyJgvNbqmZbd0NbmnaOObtrOaQjVKqRrj7m5yr0Gji9bwoZO9BANLF2gjXFUUR7MTdqsDdquj1Sq37I8lIvJGqEqJ6MQIRCdGtNkmOkXUm60tgaVlvExLkHFXZ+xNDlw+W4fLZ9sOBAZcFW93RUYXo2kZBOwaCBxuUEPgzKabAgNLF4SoFFCGKuCwO9FUZ0dTvR1Ayyq37G8lIuoyQSEgIlKDiEhNu7Oamm0O1F5pgvliI2ovuwYE17oHBl9qhL3JgXqzDfVmW9snaANQhAjQR7vCi+vPlp9bup/cXf8kP/6T6AJXmTMUdVetaKyzo+5KEwBAF6WRuWVERLeGEJUSkQnhiExoOxBYFEVY65ul8OKazdQES8sYmror1k5nNrkHA0shJuZauNFFa/iIgyBiYOkijTuw1NpQ6w4s0QwsRERyEwQBmohQaCLan6btdDhRd9XqCjGtu5tafm6std9wMDAEIMKohi5aA130tRCjj9ZAF61FRJQaSlbbA4aBpYvCjWpcqqxDfY0VtZdZYSEiulkolAqpaoJBkW2225qaUSuFmaY2VZpmuyvw1F21trsqsCC4viP0LYOA3WHGHWw4Xds3DCxdpIt0hZPaK02oueAqKRpitXI2iYiIAkClCUF03whE921nMLAoorHWDsvlRum5TJbLTai97Ao3tZeb4Gi+FmhwvO35BYUgVWj0Ma6qjL4l2OiiNYgwMtC0xsDSRRFRrunLdVesuFrlCiyRfdr2pRIRUe8hCALC9CqE6dufqi06RTTU2lqFmcaWQHPtd2eziNorTai90oTz7QQa9/ozuuhrYcbd3aSL1iDcqL6lnt3EwNJF7u6fS+fqpDEsUe0M/iIioluHoLj2zKaE2zoINBbbtaqMO8hcaqnYXGmC09Fq/RnUtDmHQikgIkojVWX00a2rNFqEG1S9aso2A0sXuUuF7vn/umiNtKAcERFRewSFgHCjGuFGNfr0bz/Q1JttrcJMq1BzuQl1l1sCTctCe+1RhAjQRbaEmRitayVid5dT1M1XofE5sOzcuROvvfYa9u/fj6qqKhQWFiI7O1vaXldXh2eeeQabNm3C5cuXkZaWhoULF2LBggU3PO+GDRuwbNkynDp1CgMHDsQrr7yC73//+z5fULBFJoRJa7EAaPdfPCIiIl+41p9RIyJSjT4D2m53P7uptlWYcQcbV4XGNWXb/TwnoO0KwVKXU9S1EOPuftJFaXrcLCefA0t9fT3uvPNOPPLII5g1a1ab7UuWLMGOHTuwZs0apKamYtu2bXjiiSeQmJiIH/7wh+2ec/fu3XjooYeQn5+P//N//g/WrVuH7OxsHDhwAMOGDfP9qoJIoVQgaXAkTh++DADof1eczC0iIqLervWzm9DOgnpOhxP1ZpvUxWS51CiNl6m93ORag8bZusupHe5p260CzZ1TkqHVyfOsPEEURdHvgwWhTYVl2LBheOCBB7Bs2TLpvVGjRmH69On4zW9+0+55HnjgAdTX12Pr1q3Se+PGjcOIESPw7rvvetUWi8UCg8EAs9kMvb7tfPvudLGyFts++AqJ3zFi8kODbqoSGxER3XqcThENZuu12U0tQab1n+6eg9bmvTIh4M/K8/b7O+BjWMaPH48tW7bgkUceQWJiIoqLi3Hs2DG88cYbHR5TWlqKJUuWeLw3bdo0bNq0qcNjrFYrrFar9LvFYuly2/0Vm6xDzovjZPt8IiIiXyhaPfKgvS4n97Rt94ym2iuuqkyYTNUVoBsCy9tvv43HHnsMSUlJCAkJgUKhwPvvv49JkyZ1eIzJZEJ8fLzHe/Hx8TCZTB0ek5+fjxdffDFg7SYiIiKX1tO249OC22vRkYCPpnn77bexZ88ebNmyBfv378fvfvc75Obm4l//+ldAP2fp0qUwm83Sq7KyMqDnJyIiop4joBWWxsZGPPvssygsLMSMGTMAAMOHD0dZWRlef/11ZGZmtntcQkICqqurPd6rrq5GQkJCh5+lVquhVge2H42IiIh6poBWWOx2O+x2OxQKz9MqlUo4nW0H77ilp6ejqKjI473t27cjPT09kM0jIiKim5TPFZa6ujqcOHFC+r2iogJlZWWIiopCSkoKJk+ejKeeegparRb9+vVDSUkJPv74Y/z+97+XjpkzZw769u2L/Px8AMCiRYswefJk/O53v8OMGTOwfv167Nu3D++9914ALpGIiIhudj5Pay4uLsY999zT5v25c+di9erVMJlMWLp0KbZt24YrV66gX79+eOyxx/Dkk09CEFzTfTMyMpCamorVq1dLx2/YsAG//vWvpYXjXn31VZ8WjpNzWjMRERH5x9vv7y6tw9KTMLAQERHdfLz9/u45a+4SERERdYCBhYiIiHo8BhYiIiLq8RhYiIiIqMdjYCEiIqIej4GFiIiIejwGFiIiIurxAv60Zrm4l5OxWCwyt4SIiIi85f7e7mxZuF4TWGprawEAycnJMreEiIiIfFVbWwuDwdDh9l6z0q3T6cT58+eh0+mkRwAEgsViQXJyMiorK7mCbjfifQ4e3uvg4H0ODt7n4Omuey2KImpra5GYmNjm4cmt9ZoKi0KhQFJSUredX6/X8z+GIOB9Dh7e6+DgfQ4O3ufg6Y57faPKihsH3RIREVGPx8BCREREPR4DSyfUajWef/55qNVquZvSq/E+Bw/vdXDwPgcH73PwyH2ve82gWyIiIuq9WGEhIiKiHo+BhYiIiHo8BhYiIiLq8RhYiIiIqMdjYOnEO++8g9TUVGg0GowdOxZffPGF3E26aeTn52PMmDHQ6XSIi4tDdnY2jh496rFPU1MTcnNzER0djYiICMyePRvV1dUe+5w5cwYzZsxAWFgY4uLi8NRTT6G5uTmYl3JTefnllyEIAhYvXiy9x/scOOfOncOPf/xjREdHQ6vV4o477sC+ffuk7aIo4rnnnkOfPn2g1WqRmZmJ48ePe5zjypUryMnJgV6vh9FoxPz581FXVxfsS+mxHA4Hli1bhrS0NGi1WvTv3x8vvfSSx7NmeJ/9s3PnTvzgBz9AYmIiBEHApk2bPLYH6r4eOnQIEydOhEajQXJyMl599dWuN16kDq1fv15UqVTihx9+KH711VfiT3/6U9FoNIrV1dVyN+2mMG3aNHHVqlVieXm5WFZWJn7/+98XU1JSxLq6OmmfBQsWiMnJyWJRUZG4b98+cdy4ceL48eOl7c3NzeKwYcPEzMxM8csvvxT/8Y9/iDExMeLSpUvluKQe74svvhBTU1PF4cOHi4sWLZLe530OjCtXroj9+vUT582bJ+7du1f89ttvxX/+85/iiRMnpH1efvll0WAwiJs2bRIPHjwo/vCHPxTT0tLExsZGaZ+srCzxzjvvFPfs2SN+/vnn4oABA8SHHnpIjkvqkVasWCFGR0eLW7duFSsqKsQNGzaIERER4ptvvintw/vsn3/84x/ir371K3Hjxo0iALGwsNBjeyDuq9lsFuPj48WcnByxvLxc/Nvf/iZqtVrxz3/+c5fazsByA3fffbeYm5sr/e5wOMTExEQxPz9fxlbdvC5cuCACEEtKSkRRFMWamhoxNDRU3LBhg7TP119/LQIQS0tLRVF0/celUChEk8kk7bNy5UpRr9eLVqs1uBfQw9XW1ooDBw4Ut2/fLk6ePFkKLLzPgfP000+L3/3udzvc7nQ6xYSEBPG1116T3qupqRHVarX4t7/9TRRFUTxy5IgIQPzPf/4j7fP//t//EwVBEM+dO9d9jb+JzJgxQ3zkkUc83ps1a5aYk5MjiiLvc6BcH1gCdV//9Kc/iZGRkR5/dzz99NPioEGDutRedgl1wGazYf/+/cjMzJTeUygUyMzMRGlpqYwtu3mZzWYAQFRUFABg//79sNvtHvd48ODBSElJke5xaWkp7rjjDsTHx0v7TJs2DRaLBV999VUQW9/z5ebmYsaMGR73E+B9DqQtW7Zg9OjR+K//+i/ExcVh5MiReP/996XtFRUVMJlMHvfaYDBg7NixHvfaaDRi9OjR0j6ZmZlQKBTYu3dv8C6mBxs/fjyKiopw7NgxAMDBgwexa9cuTJ8+HQDvc3cJ1H0tLS3FpEmToFKppH2mTZuGo0eP4urVq363r9c8/DDQLl26BIfD4fEXOADEx8fjm2++kalVNy+n04nFixdjwoQJGDZsGADAZDJBpVLBaDR67BsfHw+TySTt094/A/c2clm/fj0OHDiA//znP2228T4HzrfffouVK1diyZIlePbZZ/Gf//wHCxcuhEqlwty5c6V71d69bH2v4+LiPLaHhIQgKiqK97rFM888A4vFgsGDB0OpVMLhcGDFihXIyckBAN7nbhKo+2oymZCWltbmHO5tkZGRfrWPgYWCIjc3F+Xl5di1a5fcTel1KisrsWjRImzfvh0ajUbu5vRqTqcTo0ePxm9/+1sAwMiRI1FeXo53330Xc+fOlbl1vcf//M//YO3atVi3bh2GDh2KsrIyLF68GImJibzPtzB2CXUgJiYGSqWyzUyK6upqJCQkyNSqm1NeXh62bt2Kzz77DElJSdL7CQkJsNlsqKmp8di/9T1OSEho95+Bexu5unwuXLiAu+66CyEhIQgJCUFJSQneeusthISEID4+nvc5QPr06YPbb7/d470hQ4bgzJkzAK7dqxv9vZGQkIALFy54bG9ubsaVK1d4r1s89dRTeOaZZ/Dggw/ijjvuwE9+8hM8+eSTyM/PB8D73F0CdV+76+8TBpYOqFQqjBo1CkVFRdJ7TqcTRUVFSE9Pl7FlNw9RFJGXl4fCwkLs2LGjTYlw1KhRCA0N9bjHR48exZkzZ6R7nJ6ejsOHD3v8B7J9+3bo9fo2Xxy3qilTpuDw4cMoKyuTXqNHj0ZOTo70M+9zYEyYMKHN1Pxjx46hX79+AIC0tDQkJCR43GuLxYK9e/d63Ouamhrs379f2mfHjh1wOp0YO3ZsEK6i52toaIBC4fn1pFQq4XQ6AfA+d5dA3df09HTs3LkTdrtd2mf79u0YNGiQ391BADit+UbWr18vqtVqcfXq1eKRI0fExx57TDQajR4zKahjjz/+uGgwGMTi4mKxqqpKejU0NEj7LFiwQExJSRF37Ngh7tu3T0xPTxfT09Ol7e7ptlOnThXLysrETz/9VIyNjeV02060niUkirzPgfLFF1+IISEh4ooVK8Tjx4+La9euFcPCwsQ1a9ZI+7z88sui0WgUN2/eLB46dEicOXNmu9NCR44cKe7du1fctWuXOHDgwFt+um1rc+fOFfv27StNa964caMYExMj/vKXv5T24X32T21trfjll1+KX375pQhA/P3vfy9++eWX4unTp0VRDMx9rampEePj48Wf/OQnYnl5ubh+/XoxLCyM05q729tvvy2mpKSIKpVKvPvuu8U9e/bI3aSbBoB2X6tWrZL2aWxsFJ944gkxMjJSDAsLE3/0ox+JVVVVHuc5deqUOH36dFGr1YoxMTHiz3/+c9Futwf5am4u1wcW3ufA+fvf/y4OGzZMVKvV4uDBg8X33nvPY7vT6RSXLVsmxsfHi2q1WpwyZYp49OhRj30uX74sPvTQQ2JERISo1+vFhx9+WKytrQ3mZfRoFotFXLRokZiSkiJqNBrxtttuE3/1q195TJPlffbPZ5991u7fy3PnzhVFMXD39eDBg+J3v/tdUa1Wi3379hVffvnlLrddEMVWSwcSERER9UAcw0JEREQ9HgMLERER9XgMLERERNTjMbAQERFRj8fAQkRERD0eAwsRERH1eAwsRERE1OMxsBAREVGPx8BCREREPR4DCxEREfV4DCxERETU4zGwEBERUY/3/wEv7UP4GHLJyQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log(jaxELBO_no_jit), label='no jit')\n",
    "plt.plot(np.log(jaxELBO1), label='jit level 1')\n",
    "plt.plot(np.log(jaxELBO2), label='jit level 2')\n",
    "plt.plot(np.log(jaxELBO3), label='jit level 3')\n",
    "plt.plot(np.log(jaxELBO4), label='jit level 4')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2538caeb-51f7-4709-9fe2-ccce8823d6a9",
   "metadata": {},
   "source": [
    "- **Jitting is more involved in JAX than with pytorch**, the best performances require some specific JAX knowledge that we exposed in each subsection of this document.\n",
    "- When jitting with JAX we fall back on quite the same difficulties as when jitting with pytorch trace.\n",
    "- However, JAX offers clear ways to overcome the difficulties that pytorch tracing does not offer (scan, static_argnums, pytrees, ...)\n",
    "- In this particular example on **GPU**, **jitted JAX can perform twice as fast as jitted pytorch**. But in this particular example on **CPU**, JAX is slower than pytorch; is it link with the points raised in [#benchmarking-jax-code](https://jax.readthedocs.io/en/latest/faq.html#benchmarking-jax-code) ? However, this also recalls what we can often read online: JAX is particularly suited for GPU optimization. We need to test JAX on other examples involving neural networks, mini-batches, etc. in order to understand JAX potential on CPU optimization too.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}