{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f3755af6-a293-4b3f-b032-2226af185804",
   "metadata": {},
   "source": [
    "# Cross-validation analysis (advanced)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b038303d-a7c4-4520-b472-66c5815fe0bc",
   "metadata": {},
   "source": [
    "In the [previous example](tutorial_1.ipynb), we used a single length scale parameter. However, this should be changed depending on whta your data is; if the intrinsic spatial resolution of the data is high, a smaller length scale is needed to reproduce smaller structures. To do this, the cross-validation analysis is useful."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5762592-768a-44ea-9014-5cda945e0a16",
   "metadata": {},
   "source": [
    "## 1. Import modules and data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c45f31f-3a06-4be1-b2e7-752cc0f20d94",
   "metadata": {},
   "source": [
    "As we introduced in the previous section, we first import modules and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0f078f1a-fafc-4b2c-8fe5-799e67a2fb72",
   "metadata": {},
   "outputs": [],
   "source": [
    "import frap\n",
    "import numpyro\n",
    "import pickle\n",
    "import dsharp_opac\n",
    "import jax.numpy as jnp\n",
    "\n",
    "numpyro.set_host_device_count(10)\n",
    "\n",
    "bands = [ 'band3', 'band6' ]\n",
    "data = {\n",
    "    'band3' : 'HD169142_band3.pkl',\n",
    "    'band6' : 'HD169142_band6.pkl'\n",
    "}\n",
    "\n",
    "obs = {}\n",
    "\n",
    "for band in bands:\n",
    "    with open(data[band], 'rb') as f:\n",
    "        obs[band] = pickle.load(f)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed630a8f-972e-48a4-8bdc-b007e5167013",
   "metadata": {},
   "source": [
    "In this example, we try a five-fold cross validation. We randomly split the dataset to 5 subsets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b24def8b-ac15-4a98-82ab-53eb549e1020",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_folds = 5\n",
    "train_obs, test_obs = frap.cv.split_data(n_folds, bands, obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54a44627-8ecf-47e7-ba19-63ba271a18d9",
   "metadata": {},
   "source": [
    "You can see the splitted data like this. Here, 5 sets are plotted with offsets. The blue dots are the training data sets and oranges are the test datasets. The model will be fitted to the blue ones first, and the chi-squared of the best-fit model will be measured against the orange data. This procedure will be iterated for the 5 sets, and the mean chi-squared and its standared deviation will be calcurated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "104acf14-4f09-4550-a589-b35a7426720f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGvCAYAAACjACQgAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAl9tJREFUeJztnQeYFFX29k/35DxMgpkBiSIqMAgiEgVFETCCCsKa1tVdRXcVXcVPJAiumFFAMWNiSWL4KyKKIIiYQMKKoGQY0szA5Mh0f8+5PdVTXV3VVdVdPfH9Pc+I3V1dfetW1b1vnXuCzel0OgkAAAAAoJFjr+8GAAAAAABYAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSRBKTQCHw0FHjhyhuLg4stls9d0cAAAAABiA8/8WFRVRRkYG2e2B21mahKhhQdOmTZv6bgYAAAAA/ODQoUPUunVrCpQmIWrYQiN1Snx8fH03BwAAAAAGKCwsFEYJaR4PlCYhaqQlJxY0EDUAAABA48Iq1xE4CgMAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJqvqJk3bx61a9eOIiMjqU+fPvTTTz9pbrt8+XI6//zzKTExkWJiYqhHjx703nvveWxz6623ihh1+d/ll1/uT9MAAAAA0EwxnXxv8eLFNHHiRJo/f74QNLNnz6Zhw4bRrl27KC0tzWv7pKQkevTRR6lLly4UHh5On332Gd12221iW/6eBIuYt99+2/06IiKC6h1HNVHOeqKyo0RR6USpA4nsIfXdKgAAAACoYHNyNSkTsJDp3bs3zZ07111MklMc33vvvTRp0iRD++jZsyeNHDmSZsyY4bbU5Ofn08cff0z+pllOSEiggoIC6zIKH1pOtOlfRKWHa9+Lbk3U60WiNqOs+Q0AAACgGVNo8fxtavmpsrKSNm3aREOHDq3dgd0uXm/cuFH3+6yfVq9eLaw6gwYN8vhs7dq1wnpz1lln0V133UV5eXma+6moqBAdIf+zFBY066/zFDRMabbrff4cAAAAAA0KU6ImNzeXqqurqWXLlh7v8+tjx45pfo8VWGxsrFh+YgvNnDlz6NJLL/VYenr33XeF4Hnqqafo22+/peHDh4vfUuPJJ58Uyk76s7RCNy85sYWG1AxYNe9tus+1HQAAAAAaDHVS0JKrb27ZsoWKi4uFcGGfnA4dOtDgwYPF52PHjnVv261bN+revTt17NhRWG8uueQSr/098sgjYh/KKp+WwD40SguNB06i0kOu7Vq62g8AAACARiZqUlJSKCQkhI4fP+7xPr9u1aqV5vd4iapTp07i/zn66ffffxfWFknUKGHBw7+1e/duVVHDTsRBcyRmp2ArtwMAAABAw1t+4uWjXr16CWuLBDsK8+u+ffsa3g9/h/1itDh8+LDwqUlPT6c6h6OcrNwOAAAAAA1z+YmXfW655RaRe+aCCy4QId0lJSUiTJu5+eabKTMzU1hiGP6Xt+XlJBYyK1asEHlqXnnlFfE5L0lNnz6dRo8eLaw9e/bsoYceekhYduQh33UGh21zlBM7Bav61dhcn/N2AAAAAGi8ombMmDGUk5NDU6ZMEc7BvJy0cuVKt/PwwYMHxXKTBAueu+++W1hfoqKiRL6a999/X+yH4eWsbdu20TvvvCPCujMyMuiyyy4T4d71kquG89Bw2DZHObGA8RA2/JqIes1GvhoAAACgseepaYjUXZ6aNi5Bgzw1AAAAQIObv+sk+qlRwsIl82pkFAYAAAAaCRA1ektRCNsGAAAAGgWo0g0AAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkAUQMAAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkAUQMAAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkAUQMAAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkQWt8NaAxUO5y0fXcF5RVUU3JCCHXrFEEhdlt9NwsAAAAAMiBqdFj3aynNW3qKcvKr3e+lJobQhOtb0KDzovW+DgAAAIA6AstPOoJm2uu5HoKG4df8Pn8OAAAAgIYBRI2PJSe20Phi3rJTYjsAAAAA1D8QNRqwD43SQqMk51S12A4AAAAA9Q9EjQbsFGzldgAAAAAILhA1GnCUk5XbAQAAACC4QNRowGHbHOXki9QWrvBuAAAAANQ/EDUacB4aDtv2xYTrWiBfDQAAANBAgKjxAeehmXZHipfFhi00/D7y1AAAAAANByTf04GFS/+sKGQUBgAAABo4EDUGl6J6dI4M/tkAAAAAgN9g+QkAAAAATQKIGgAAAAA0CSBqAAAAANAkgKgBAAAAQPMVNfPmzaN27dpRZGQk9enTh3766SfNbZcvX07nn38+JSYmUkxMDPXo0YPee+89j22cTidNmTKF0tPTKSoqioYOHUp//vmnP00DAAAAQDPFtKhZvHgxTZw4kaZOnUqbN2+mrKwsGjZsGJ04cUJ1+6SkJHr00Udp48aNtG3bNrrtttvE35dffune5umnn6aXXnqJ5s+fTz/++KMQP7zP8vJyqk+4AveWP8pp9c8l4l9U5AYAAAAaLjYnm0lMwJaZ3r1709y5c8Vrh8NBbdq0oXvvvZcmTZpkaB89e/akkSNH0owZM4SVJiMjgx544AF68MEHxecFBQXUsmVLWrBgAY0dO1Z3f4WFhZSQkCC+Fx8fT1aw7tdSmrf0lEelbk7Cx1mGkXQPAAAACByr529TlprKykratGmTWB5y78BuF6/ZEqMHC5jVq1fTrl27aNCgQeK9ffv20bFjxzz2yQfI4klrnxUVFaIj5H9WwoJm2uu5HoKG4df8Pn8OAAAAgIaFKVGTm5tL1dXVwooih1+zMNGCFVhsbCyFh4cLC82cOXPo0ksvFZ9J3zOzzyeffFIIH+mPLUVWwUtMbKHxxbxlp7AUBQAAADTH6Ke4uDjasmUL/fzzz/TEE08In5y1a9f6vb9HHnlECCXp79ChQ5a1dfvuCi8LjZKcU9ViOwAAAAA00jIJKSkpFBISQsePH/d4n1+3atVK83u8RNWpUyfx/xz99Pvvvwtry+DBg93f431w9JN8n7ytGhEREeIvGOQVVFu6HQAAAAAaoKWGl4969eol/GIk2FGYX/ft29fwfvg77BfDtG/fXggb+T7ZR4ajoMzs0yqSE0Is3Q4AAAAADbSgJS8d3XLLLSL3zAUXXECzZ8+mkpISEabN3HzzzZSZmSksMQz/y9t27NhRCJkVK1aIPDWvvPKK+Nxms9F9991HM2fOpDPPPFOInMcee0xERF1zzTVU13TrFCGinHwtQaW2CBHbAQAAAKARi5oxY8ZQTk6OSJbHjry8RLRy5Uq3o+/BgwfFcpMEC567776bDh8+LBLrdenShd5//32xH4mHHnpIbHfnnXdSfn4+DRgwQOyTk/vVR0VuDtvmKCctJlzXQmwHAAAAgEacp6YhUmd5alqECEGDPDUAAABAw5u/TVtqmgssXPpnRYkoJ3YKZh8aXnKChQYAAABomEDU+IAFTI/Odb8EBgAAAADzoEo3AAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEqNJtgmqHk7bvrqC8gmpKTgihbp0iRCVvAAAAANQ/EDUGWfdrKc1beopy8qvd76UmhtCE61vQoPOig3V+AAAAAGAQLD8ZYN3mYpr2eg7l5J/2eJ8FzrTXc4XgAQAAAED9AlGjQ/WB5TTvvd+JyElE6ktN85adEktTAAAAAKg/IGp8cWg5bf/8WcqpaOmzq3JOVQtfGwAAAADUHxA1WjiqiTb9i/KqUg11JDsPAwAAAKD+gKjRImc9UelhSg47YagjORoKAAAAAPUHRI0WZUfFP90SfqLU8CNsutHY0EmpLVzh3QAAAACoPyBqtIhKF/+E2Bw0ocP0mjeVwsb1esJ1LZCvBgAAAKhnIGq0SB1IFN1aRDwNSllJ07rcRanhxzw3icihaX9LRp4aAAAAoAGA5Hta2EOIer1ItP46t7Dpn7yKthdcQHlVLYWvTbeRD1JI2951esIAAAAAoA5EjS/ajCIauExEQbHTMC9F9Uj8gSi6DVGv2a7PAQAAANAggKjRg4VL5tWuaCh2HmZfG16aYksOAAAAABoMEDVGYAHTcnDQTwYAAAAA/AeOwgAAAABoEkDUAAAAAKBJAFEDAAAAgCYBRA0AAAAAmq+omTdvHrVr144iIyOpT58+9NNPP2lu+/rrr9PAgQOpRYsW4m/o0KFe2996661ks9k8/i6//HJ/mgYAAACAZoppUbN48WKaOHEiTZ06lTZv3kxZWVk0bNgwOnFCvfDj2rVr6cYbb6Q1a9bQxo0bqU2bNnTZZZdRdna2x3YsYo4ePer+++9//+v/UQEAAACg2WFzOp1OM19gy0zv3r1p7ty54rXD4RBC5d5776VJkybpfr+6ulpYbPj7N998s9tSk5+fTx9//LFfB1FYWEgJCQlUUFBA8fHxfu1Dta0OJ23fXUF5BdWiCjcXrQyx2yzbPwAAANCcKbR4/jaVp6ayspI2bdpEjzzyiPs9u90ulpTYCmOE0tJSqqqqoqSkJC+LTlpamhA8F198Mc2cOZOSk5NV91FRUSH+5J1iNet+LaV5S09RTn61+73UxBCacH0L1HoCAAAAGvvyU25urrC0tGzZ0uN9fn3smGexRy0efvhhysjIEEJIvvT07rvv0urVq+mpp56ib7/9loYPHy5+S40nn3xSKDvpjy1FVguaaa/neggahl/z+/w5AAAAAJpxRuFZs2bRokWLhFWGnYwlxo4d6/7/bt26Uffu3aljx45iu0suucRrP2wpYr8euaXGKmHDS05sofHFvGWnqH9WFJaiAAAAgMZqqUlJSaGQkBA6fvy4x/v8ulWrVj6/++yzzwpRs2rVKiFafNGhQwfxW7t371b9PCIiQqy9yf+sgn1olBYaJTmnqsV2AAAAAGikoiY8PJx69eollokk2FGYX/ft21fze08//TTNmDGDVq5cSeeff77u7xw+fJjy8vIoPT2d6hp2CrZyOwAAAAA00JBuXvbh3DPvvPMO/f7773TXXXdRSUkJ3XbbbeJzjmiSOxKzj8xjjz1Gb731lshtw743/FdcXCw+53///e9/0w8//ED79+8XAunqq6+mTp06iVDxuoajnKzcDgAAAAAN1KdmzJgxlJOTQ1OmTBHipEePHsICIzkPHzx4UERESbzyyisiauq6667z2A/nuZk2bZpYztq2bZsQSRzWzU7EnMeGLTu8zFTXcNg2Rzn5WoJKbeEK7wYAAABAI85T0xzi3KXoJy2m3ZGCsG4AAACggc3fqP2kwqDzooVwYYuN0kIDQQMAAAA0TOo0pLuxCRsO20ZGYQAAAKBxAFHjAy6J0KNzbT4dAAAAADRcsPwEAAAAgCYBRA0AAAAAmgRYfvIDVO8GAAAAGh4QNSZB9W4AAACgYYLlJxOgejcAAADQcIGosbh6N28HAAAAgLoHosYgqN4NAAAANGwgagyC6t0AAABAwwaixiCo3g0AAAA0bCBqTFbv9gWqdwMAAAD1B0SNiZIJE65v4XObCde1ENsBAAAAoO6BqDEBqncDAAAADRck3zMJqncDAAAADROImgDKIwzuFY3lJgAAAKCBAFFjEJRHAAAAABo28KkxwNrNJTTt9VzKya/2eJ9f8/sseAAAAABQv0DU6LB2UwnNfDPP5zYojwAAAADUPxA1PmALzONv5pFeOaecU9XC1wYAAAAIKo5qouNrifb/1/UvvwZu4FMTQAFLf8ooAAAAAH5xaDnRpn8RlR6ufS+6NVGvF4najEKnwlITWAFLf8ooAAAAAH4JmvXXeQoapjTb9T5/DmCpscLyUp/lEaQw85z801RQ5KCEWDultggV7UF2YwAAaALwEhNbaEjNF4LfsxFtuo8o82oie/N+wMbykwWWF63yCMq8NlYLDbUwc4m4aBv16x5FPc+KhMgBAIDGTM56bwuNB06i0kOu7VoOpuYMRI1OAUtfS1CsTybfniyyDNd1XhsOM3/8De2orKJSJ335Q6n4Y2KjiPpnRVOPzhFUVOKkhDg7pSbCogMAAA2esqPWbteEgajRKWDJeWi0mPzXZBrcM0ZV0Mi/Z6dq6pbwEyWHnqCPlqQROUfQoJ5xgYWZv+U7zFxJcRl5iByJuGg7jb44lsZfnoDlKgAAaIhEpVu7XRPG5nQ6dQKWGz6FhYWUkJBABQUFFB8fb+m+VS0uLULEkpOaxYWXnMZNPuLefmDyFzShw3RKi6hV0LlV6dRi8BwKaTvar/b4Elr+Eh9tp4njkyyxIhlaH85ZT47SI7Q3L4UOOC6k5MRw+AEBAIDWmPlpO5dTsKpfjc0VBXXVvkbnU2P1/A1RYwAzvjFb/iinibNPuAXNtC53iYvQJtvc4bSJ17aBy0yF4XE7xk8+ROnOjZQcdoLyqtJoe8EF5CDrLuIpf1O3PgUzJPFERTrN2zuVtpSNoFFD4ugvw2E1AgAAr7GTo5wEcmFTM7mYnE8aChA1ddApgbD65xJ64u08seS0sHd/Sg0/6iFoJJxkI5tJZb3vx0UUs2Oih9VHEgTr84Zb0n7hJ8TLar2sFzbVBz4k+4bra45eIfKIaNrOV8RxxMfYaeK4OrIaAQBAo85T04ao1+xGKWgYiJoGJGrULDgfrCygBZ8VUlbCRnqh21j9nVyyxpi3+qHl5GSV7lSx+sgEgVVMuyPFUlGxbnMRnbPtLEoOUxd5fBy5Fa1o3C8b3JanoFuNVJbEhKMdr0unDmx0ZtxgR9sBABoAQRirqutx7LB6/oajsIW+NuyXUljqEP/Py0OWeavLchQoBYHd5hSCYEL76bQh7zLLlqK4nlX/rChLLmwWNNtWPEODOmgfKx9HWuRR4VDNS2r874ZPTlBiSSfq0f8yywWGPL9PTM7H1LN4EkVW83p1DeFJRGf9i+jcR1V/2z0I5FdSW/sP1CE5l+zRGfUjhhzVtPX7r2jthj20Pz/FvSSJsH4AmqAI4e/qPAhX64gU+efZOVX02XcllCubyxpzAAl8avzAiLOupZYaru+xeojuru7fvoi2FvQlq3j+vjTq0Tky4CWnU2vupZRwY6GGSw/fThelrvBYYiuiDNqe8BSVpFxjSc4duSDV8nuSqLQn0S9xL1JR0tUiTD4h1kmOY+tp584DFO/YQyNb/ZfSIo65ty8PzaSjbZ+lM3qPCf5gwIPlb09Q1W+zKcxxSndJkkVOXfksKZNC8m9z/8XF2ERKgbhYOxUVOyghxkntQ3+sV1Go1VYr0h64933qNBXw8fq5T0tEdLAc9BX73Xe6DxWU2Mwdq0bbzukQTjv2VvptQdCa3C21TPgSLb7KGnCivCBYXN7/ooCWrykS17BaOhFf+c2URIXb6IZLgztmYPmpDjrFF8roJi0kn5qU8GPCChGIt/rvaxfQ2Udu023bzJ0v0feFV1N5JVnC6CGxNOH6pMB9aBRLZr6QYvHUltim75xHBaeTqF1CDp11djsKSx9kelCW5/fR83uSt2fazvniX2UkG3+u1tZZ+16lzJ7Xm3rSUU6ASgHgFgaxdorL/ZjOL/oXhTtOeu1Hb0kyKsJGYy6Ns/QpTD7pnj66Tog+udVIDbXIQAprQdT6anK0vIT2nkyzZOL1Ei1Sv0bbaMufFfT91jIqKtMOAnVbvDqHCQHWLimH9p9M9WgboxQd27OT6MkV59CJfO996uWskrd5885y2rC1jM6LXuHVX+UhmbQ5dhaVJF/pIQ6rkwfQ9r2nxaTdIs5Oiac+poxD/6bI09le4vf30yNp5IBYykgJ9RBekqBwZyuXBGjSCbJX5hJFpBIV7yHa/RpRmfd+pWtP6r8eZ0a4Hwzah9T2Y2nhCeqSN4nCq7z38V3ecOESK6XF0Lv39SwQKQl26toxgjbtLPeY9BOiHfT3gf+jy7oXeohFXxnbmYM/L6b0Aw969Kv7waZVGIXU+BDKcXkUOqnKnuRx/yofiHSFlySmOCKqIod+PxZPH3wbRz8cP999z7nTidQElnQ5/xJavNozrYcRgunn2CBEzbx58+iZZ56hY8eOUVZWFs2ZM4cuuOAC1W1ff/11evfdd+l///ufeN2rVy/6z3/+47E9N2Hq1Kli2/z8fOrfvz+98sordOaZZzY4USOPbtJDsgKIm1ImbCQ3WUf/pbph3ayqP1ryf4asPlvP+IK69hvmNRgW+xiwg+Vbo+dDo4QroTvJTnZyqDtWO4kcZKcQm2t5Tz7wbSsfSf8a20LXuVnK7yNVXTdqTePti6oSKT6sQNOio+UfFBsdph8qr7F85Pu6+of4f622qPkoKYkMJxpzabx4CmPkA6hyQnMLAZkgkCwvW/6ooO+3lVGPKO9JV8tqpGchU35/c+kIurRPDGUkh3qIO6XYk7+WW9WM9Ksv1ASY1LYfC13H1ideexstf7drLoqhjJQwt2DgiX719gSav74bFZTWdozPSEpyUuHpREoIq1VPuZXpNGeP63e1v+uKm1l25Hb6/uSl7v6RJsKU8BOUW+mKsOyfvMpbgHom6TckqtX60deDDO+D0erXLaUjqF9WlEgqytfhhq2lIieXcjL/raAXnZuwSTVqVK1NLBa/qJxJC7Ze4iF+JCLCiPq1+IImn6l1ToiKTidQXFi+R98YOWZ+IDqdfg1t31OpvSSU/ZG3BUjRN4zZ67GufS0bhKhZvHgx3XzzzTR//nzq06cPzZ49m5YuXUq7du2itLQ0r+3Hjx8vREq/fv0oMjKSnnrqKfroo4/ot99+o8zMTLENv/fkk0/SO++8Q+3bt6fHHnuMtm/fTjt27BDfaUiiRopuCmgwLE+nefum0s7qK3Sf1tgqxE9/vqw+QiRFtyabitVH7el/y64KWvVjqWq2AyWck2fhjAxDT8nSb323tZT2bP7K2PKbbGA2as3RGjxvGBpH/xjVQlPQPP6m53m7OOUTmtzln4Z/T2mVMbMcOOV274gy7i9+0kvZ8wDF0RFDA48R65JWG7QICSEKsRFVnq59j/dtZmTwPel6TnBmjsHVBptpR3hfIsTsgG7k2Bijx2+2rbqRlDrWzbs7zDDU1/yb3+RcRRenfurRlgIW9KEuwWT0+lcT1Vr9qHVf8T6KTvNvu5ZWjfarWn9WO9UfiPw9b0bOidnxzOjDyMUtv6BHJTGl8rlTcd9aGVhiZj5oNKKGhUzv3r1p7ty54rXD4aA2bdrQvffeS5MmTdL9fnV1NbVo0UJ8n8UR/3xGRgY98MAD9OCDD4pt+OBatmxJCxYsoLFjx9abqFEz//Fro5YaCeVTg/KJUUv9quW8UVp9/M15ozbJB+Jbo1ynNSMYWOR9mzeCrs98kwIdBK69KIYGnhfjYarlJaeZb9ZaaEz7PfkJLwd+k3u1+H9uydUXxdCg82KEFWTRqkLK3ryUJrX/u6kB1Wyb5W0IFnoDvPIcmT0GIwO9vwLLqmMjm5NSw48Zjuwz01a2wvhznfI+WJC0CK/1t/KFfCZQig7le0aRRLVZMa5sl9F+NSqc5BYuFmxmz1uwxw6thxGj/ej0cc7M3k/B8LUM5vxtN7NxZWUlbdq0iYYOHVq7A7tdvN64caOhfZSWllJVVRUlJbl8Nfbt2yeWseT75ANk8aS1z4qKCtER8j+r4UmarSQsKtgyw//y6/yiarEebga+cPgC5cmF/1VeSBxpxALKV6VwHoR5gMutbOWxDV+c29PfN52jgK0Go4cYK9XAS1hqcJtZeL28JFcskXULXS5udr7xWLwZYe7ex8TNxeZvf5BHTTEffVviPld8DiUfGpXuFeKSn9iClVNb3gf8Ex/XtO3yfx6mdz87RbdnTtWMaOPtOaKN+1KO4ag6lTYEC+57firWGmSV58jsMSi/73NbqhZP6Wb7NdBjY2dxo8dvtq0pYbWO6GbgfRgVNIzUBmVbxEOTnw/m0rnW60cj7dLrV1/9qdW/vGRn9rzJjytYaO3faD/afJwzM/eTmfmgoWAqpDs3N1dYWtiKIodf79y509A+Hn74YWGZkUQMCxppH8p9Sp8p4aWq6dP54q3b6Ca2QrB1g5c5lnxdZNnv5ZyqFhYgpfpVVgpnYcNh20qrz7MD/Kv3wSHbH67RPw7ehi0fcmuSZJnpEvKZy9TbzdN0/vLex8S/Wktm/LRwsiqduo/4N3V12mnmmy6Boe1YbW4Q4HOlF6HG4pJN0MIC5sOZmQWRGWur9CTE50YLaXDSQj7wyJ/YjIoUFmo5Fek+22AVRgd4aTt/hZaR3/G3XwP5TaMo92W0rQnh5uq8BYK/4kUL6VwHUwQohZNRjN7Tam03eg2bHTv09h/M69EoavNBQ8KUpSZQZs2aRYsWLRI+NUZ8ZbR45JFHhKlK+jt06JBlbWTrA0/WvlizqVQkhjNrsfGF3CqjrBTuy+qT3KI2+sIsavvXQm5NYusHCwYWNCwI2BQqh4XJ1C4TxPq8rWaS9ziGmiWzFkNeEoU92Wo05Y6W9N7x6arbG8HfiVKygLEZWg1pGYDN+EbaJW3PPlO+TLtmhYDSuuSrLZIfil4brMJo30vb1R6D9b/jb78G8ptGUe7LaBsKKpN1z3lDg9vKy8qSqPanH41aUIMtnNTarncfupb+WqiOZ3xcWsem7DcjbfGXvAD2NXfpSdXVhUYnalJSUigkJISOHz/u8T6/btXKc1lEybPPPitEzapVq6h79+7u96XvmdlnRESEWHuT/1mFK3KoWteykhgbQgtnZtCwC61Rq0qrDMM+IXddpz7ZSnBhTX+dtqRK5GasSSKC6M08Q6bzi1M+FY6KyiUzttAoI79Y9d/30N9of6eFVGLztDyxk5/Dz0HAqLAZ/eNmevvA/V7ihi0uHM79/O5ZmgOU9/b6PhtmhYDSuuRL/BVWt7A8w7QvjAzw8nMkPwYjE5eZc+xvvwZ6bCcqWhk+frNtyK1qpXnOfU2Q9YWasNfrR+UxuLZz+QQZ7Vd/J2l/BIav+1A6/ud3PynGDuX4V3jaNeZqfc/Xw4jRBwJnAMLJCLn5DpE9v9GLmvDwcBGSvXr1avd77CjMr/v21TblPv300zRjxgxauXIlnX/++R6fcbQTixf5PtlH5scff/S5z2ChZjHR2o7XFr/8wXzMv5pHuZq1hZd4XlmWr/kdK8Lr+Pucj8YI6351ORfzDWXU14Dzyoz7eYNwfGOnVf53R/edqqHsLLLa9xlLcWMPivD0Fw/NFdvP2DmXK2X5NQgY5ZrBCXTe6P9QzLgcj99mfx8WB1o+TRVhmfRrxCP03P45HttbLQTkaLWFEwVSt+kUO+4EXXvTLfTIrUlCdHNemmBiZIBXniPXMczXtJB5OHiaOMeB9GsgxzZv7zRTx2+2rVrnXGuC9BetidDXJGlE2PvqR7X9SvvQepjwRzipWlMq+frzb2zx5ecoHT//Kce/0T9uUhU7Rh6IavuR26y+jVP2fjDHTC4HxHNUQ8OvkO5bbrmFXn31VZFrhkO6lyxZInxq2A+GI5o4VJv9XqRw7SlTptDChQtFaLdEbGys+JO2YSuOPKR727Zt9RLSbTQPzbP/TKWn3j1pKCujHmriRC9rsVqIcF3k3pEwGt0kj74RxTJvN17PSa+UgRQaH4hFQquAp2aW2ZrEYcqMrrz9eysK6N0Vxp3WfUa0GYjScSckS8ylwf07Ula/SzVLOnCW0SVfFVFZZfAe69VCactCM+nXmFnujMxqGYW7lj1HrfLmka3SO5GgP+fY336VJ4kTSfm2lblzlPhKzSDty8g2gbY1IcZB15y9hfJzst25d3zlkPGFMr+MNBO42qKfC0d+jJ8dv5GOlLUT1pJOPYZQcYXdo//kx6vWRy/ve0w8AGnlkbm343SPrORa/aqdG0w7lw5j9ryZiW61+nvScWqd8xM1bQ/0uOoqxLveQ7oZDseWku/16NGDXnrpJRGtxAwePJjatWsnwrEZ/v8DBw547YOT7U2bNs0j+d5rr70mku8NGDCAXn75ZercuXOdd4qRjMF8Ih++KYkefCnH0D55+7tGJdIrH+Z77Jff5+UjpaAx2gar8gXw7934aDblFtTmcdDDaEijXq4WU9Rk0Pztt3305lfRtC2AZGoSVhfN/KZmec4oaoNTEWXS9oRZVJR8tWZGYSEIZNlNjeYRYpPxolVFVB4kcZMY46S/D9hOlyoys+qiyI6640g8LVzvmR3VMoEl9atOOQSvbLKyBHlfb0+guWvOpdJKz7bFRDjo+u5b6NSJWtHB7ef7dUivaFrzS6nXfe1LDG0pG+HKZnxWpMe5Vmub47gr0eCB/CSa1HkipYQf18xrZePaZiFRRGW1ydvKQjLpQMRoalfxocfDg9Rv7qzFSSdo794jtGJTFO3Lr52MleOZ3oOB1I+vru9KBaV2r/EyMS6kNp1Gh1AKyftOlFH4alu8V3JCZX8qRRDZQoictf1eHtqajrZ9hjJ73eAqw6CTCZsTVfJMWVHl/Xss1K+5KJZWbCihPBNjqBy2qK7bXGbqoUMSRSPPK6CdfxwWvle5Va4ABa2MwnrCSWuuCmaId4MQNQ0NqztFz0rClpWq007DSfgkSwzf5Fv/KBeZLxnOgpnVOdJrQjJqObEyX8A7n+fTO58btzLolYGQ50JgZ2Y18RYI8nIH/mDWamSG+ctPmYqOkwae1nE5NHRQJ02Li1VI4mYxW24qfN/+butF5wjNjML82oo6SWrt9EgcqZNBWC2jsJpVzSp83c96NYeUxyQXTMoSDP7UOrJnL6duR/8i3mMrSy01++K8Vlq1hwwWYLSqfpI/+9EqfeEW+jUiyH0Myf2I8r43fExqtboYrZIJ3F4jNQF9PaAyfF9++E0xFdUURhafJ4ZQ147h9MvvFZ7vy0SkWj2n1JrPGa1aT7zvkQNiKDM1zKPvzRzLo7cl0yW9/R9HIWrqoFMYXxcJX0RGhcetV8TTzSMStfepUgPGaNbiQC+mQDIlGykDsa3le+RsMzpoZezNFGYL5vKdGuxQ/eKiU1RQov/kFhlBNGaoq1xBXVbEVZuUu3aKCKiAIGhAqBZTbEPUa7bpvFbA/zGJ6yYV+hgHlO4H/hbh9PV5tR+FVY0+6MJSEwTqMqOw/CLRWyJKSbTTf2dmGlK+8gu7Piw1/vjVaK7t1uHAKS/NsHxNcb0WZtNqGw8km3eVe/kYcFtGDamp5QLhAIKBQasLoKDOGxxU4ushuSFSzW4Jk4941J9qsj41DY26rP0kx6hQMesjw9uPn3yI0p0bNddCra7BYbT6uK/lkzsuK6VzzmlfbwOnL8tNQxAQVpnrAQCNj8Z4/68z8TDuLxA1ddApVi5T+WV5ObScKjb+kyJk5eyVxfiCUS3V3zXhYPqnmMUfMysAAAD/57hAgKipg06xWoGb8pFp9SXR+utqghBVEjodepV6DxsfNJMlX8DPLzzpcw24rv1TAAAANE0rU6HF87ep2k9A+wT78m1Ryxasuh2fT3bsUwia2iy9NprYZQbZsv4atFPBYolrQom8Jl/7jo6pS/8UAAAA9UOI3WZpZe5gAlFjEqMRTGo1lvR8arrF/+QZqaBAhGaWHnI5/rUcTMG8gG8ZmSiicVjcLF9T5OHgyuGTowbH1Xm0DgAAAOALiBqLqnfz+1p+LlKNJV/+KqKGU4V6VXIvOJKhDpCLm8bm4AYAAKD5UadVuhszRqp3yytZK2Gxw6JHWRWbrR63jIwXSz4i5NIIRrez2PTIOXH4XwgaAAAADRGIGourd/N2WrCw4crenJAvLtrV9ZwFkxMccTj1uuxe5Ixq7U5e543NlQOGQ6YBAAAA4AFETRCqd/uCkzBxdVN5umv3EtYbp+ip36YIP2HvSrM1rzmpHZJnAQAAAF5A1BjEcASTj+2MLGGtOjJMtZx9eWimq14L0psDAAAAqsBR2CCGI5hqCp/5u4TFcIK9DXmXeVRXPWrvSx9c1ybAmtQAAABA0wWWGoNIEUy+EBFMPqKCjC5hMVwSYWtBX/om92rx74lTLlEEAAAAAHUgakygFcHEFhojZQuMLmFZIYoAAACA5gaWn/zMuOtP3hYjS1jBFEUAAABAUwaWGjM4qomOr6WQg4uoR8IPdEmvSFN5W4wsYfnrrwMAAAA0d2CpMcqh5a66TPIyBtGtiXq9aCoiiS09nKeGw7rNoOevAwAAADR3YKkxQPWBD8m5/jpyKusylWa7Kmqz4DFBZmqYqe1HD4lD0UgAAABAB4gaHdZtLqJTa+4lcqrl+a0pibDpPtfSVJB8Y0QJBQAAAAD4BKJGp4DlR0tXUEr4UbJprvzIKmcbRHIYNgJ8aQAAAMiTuG75o5xW/1wi/tWqN9hcgU+NTvbfbmEnLK+cbaRqtwR8aQAAAEgP2jwvySNo+QGZ5xO9lCLNBVhqdLL/cjbfYFTO1sp5Yzb3DQAAgOYhaPhBWJkSRNQNfD1XfA5gqdFNdLe94AI6UZFOKeHHyG7zNvOxp42No6D8qJwtz3mTc+o0FRQ7KCHOTqmJoYZz3whfnpz15Cg9QnvzUmjf6T5UUGIzvx8AAAANEiN1A+ctOyXmk+Y+3mP5SceZl8sVzNs7laZ1uUtUzpYLG34tfG00KmfzhaiXpI9fc66bQMPM2eTWiYjiK9JFe7l+FBMXbafRF8fS+MsTmv3FDgAAjREjdQNzTlWL7fyeT5oIEDUGsv+yQODK2RM6TKe0iFrfmZNV6dRiyEsUopKnJthrnxxmbt9wfY2tqBa2KLEA4/Zyu4tKHSInzqJVhXRRz2jqeVYYtQ/5kTok55I9OsNlYVIRZAAAABoGRkvk5KGUDkSNUWdetcrZ114/gga1jdNc+1QirX1OuT2ZBveKMXVRe+x/cxGds+1eSg5zekVlsSWJLUgT2k8X7WVLE1NeSVT654eUVe0pzMpDM+lo22fpjN5jYMkBAIAGiNE0IMkopQNHYTPOvFLl7O3Vo+jaG66kQT3j/Fr7nPlWHq3dXELBCjNnYZMWeVQIMImByV8IC05quGeUVnjVEWq3exzNfvqNOnM0k0ISv/qpmJatLqSvfixGaCIAAASQBgTpP1xg+cniApZG1j45rcDjb+SR/Q6bqaUoFgNzl5yk7gbDzNmixNipWiydiaUqDcvOTWlTadzrQ2ny39JocE//rUh6qC3LScRF22jUkDj6y3D4/wAAgJk0IEj/4QIh3QaQnHkv6R2jW8DSzJome6ubSZz0wcoCyi1wGA4zl7Zjiw0vORmx7Mx8M4/WbvLPimRk2eyjJf9H3UKXU1bCRiG25BSVOumdzwtp9MPZCE8EAAADaUCQ/sMTWGoCQC26KTunyvD3zXirs4VDKoKpF2bOlpfcilZiO7nFRg/eTliR3swju92cFcmIYzP7AQ3qVrv8dUIRqSVRWOJw+R79LTmoViP5OfQ7pL4hURPeLxJBct4kOIED0KxXDoxiJFK3sQBRY+EySny0nQpLHab2Y8Syo/TT0Q0zZyvQvqluJ2Gzlh2rcx5IkVrs2Ew+IrWUzHgjj/aPqKKbRli/HMV9+v4XBbR8TZGwEKkthfXrHkU9z4qk1BahdE6HcNqxt7L2pu8QSiF537nzAx1wXEjJieH1MhjwsRz8eTGlH3iQIk9nu993hifRseQJtCPmwXprW1MePAEwck1bfc0bSQNSbaJN/CD+2XcllCubyxpzKhCb0+ls9IUjCgsLKSEhgQoKCig+Pj7ov6cV3eQPz9+XpnuBslPtxNne1hZ2/lWGmZ8oTxeCRi4SeJlnYe/+upadWX8+R0lhuULcsJXn2fvSA8554IrUOouSw9SXv6TfHvfLBrcIUxIVwVajWoER6KDA5+/5hSeFRcgo3HbpTuF+v7fjdOGsLXGiohV9fuxGKrR3pMH9O1JWv0tNh8obHfzk1qXNu8rJcXA5TWr/D1WfKaagKpGWH7ktoLaZxd3G/NNUUOQQInHLnxX0/bYyDxEZG2UT4rlH5wgqKnHWi6VMra3cxoRYp0h/0C4ph/afTPVLuFplCbRiP8rjTIi1W3I/ee03xkntQ/XTRsiv9xZxLk+IvMJqj7Z5PUz42W/K7xsWGjrWT73UHb4+t8ri4nEs8TwYraNvv99L+/NTxDjO46pem7SICrfRDZcG18/R6vkbosaPC2jc5COGLgg9+EJbODND92LhC/DDNUWqn7FgkYeZd+pxMRVX2OjbTaVULlsJk6KfeDrxtuw4qfB0IiWE5XssDR1r9yx1HzLO7+PjCK9Pln5GL3Qbq7vt/dsXicgyIyREO+jvA/9Hl3UvNJ1rh9vETtr+IvWjUkCw4JG/LqIMyu34nKFQeS2rkXScQ7sVuCbV6t5UdWwD7fp9P+0vSHUvL7Jg5ag2LZ+pQNpmFKXQUooXs4hjH/A/Oistj05WpVFIq0HU7cxovyYe9+el2UQVOeQIT6G9J9NU+1MS1qoPDDXLpVtKR1C/LJfIlkJopcnYLYji7HQk5zR9/l2x8IPTc4rXElZxsXbasqucNmwto+Iyp9d93y4xl7p0aUv2lgNdmcTVhBi3UWWic7clyuY6ns5hQoxI391X3YcKim0UF2MTglPepiIWVrF2OpLreYxq/VYekkmbY2dRSfKV7v1/vT2BXl3flQpKfbt18tmWX0Xutmo84OhZINISicadt5X27T3o0RexUUT9s6KpZ5fIWrGY/ZE7ualaCoz1W0o9xhLlWNzl/Eto8WrviFJpu9axOXS4uPa6U1pHjFhb2M/yw2+KRT4yX9fs+pqH3OsuiaVlq4vJLPExdpo4LikoZXsahKiZN28ePfPMM3Ts2DHKysqiOXPm0AUXuAZYJb/99htNmTKFNm3aRAcOHKAXXniB7rvvPo9tpk2bRtOnc3ROLWeddRbt3LmzwVlqtKwm/sDLVRPH+75QzFiFbhkZT7eMTFS94Bm1i76gqgXFh7qWtuQTn5Qt2TZwGZFKckFf8G+/90U+vbeiiIYkf0KTu/xT9zszd75E3+Rerbud2jFQRApR75eJzuBkhNqwAzSH0/tb1FayePkSEEqx+Hne36hNz9Ga1hEtq5HacVY77RRic3gMWJ8fG0u3tX3B1HFIS5RP73uZunfNpC5pJymvMo32n76Q4mJDXRaTmqdtvQmOJ1HHsfW0c+cB1QnTH7QG5zezp5H9jFHU48wI9+Qal/sx9SyeRJHV2V6TaFHy1RSX94nX5776kycARk24Sv2mtVxqBmliy4zNodT01vTxjvOooNQWUP98k3MVXZz6qeL+do0HygcWpS+bkQnRSLu0+039wcnM/uW4RV1CDp11djsKaTWQNu+qog1bS6m4TLt9Ro/x0oyVLuunzTO5qTRbvp89kd7Zf48hEWy2nyPDbdS3WyRt31OpuSTEIvf5D0663R189z1Zcs0ywahHWO+iZvHixXTzzTfT/PnzqU+fPjR79mxaunQp7dq1i9LSvH03fv75Z1qyZAn16tWL7r//fnr44YdVRc2yZcvo66+/dr8XGhpKKSkpDU7UcLn3J972/0nfzIVixiqkZfXxMlvHOqn66Hra9Os+OlmZQg93nkip4cdUJ2l3Xaur9hmyhEgWhyVfFVFZpeuy4ignqyw1Wjeum7P/TXTe09qRV0tXuJ+i/Jl8jR6LGiW2DIrs9xKFtB2tazUyag2SJgs9gaUG78tB2pO6VQO0WcwMznrbLs6+k8ZkvqZ5vWj1J0++8aH5fi+XGjlGf/vN17Uhofa+r760YkLUE/xG2hHM/jNzjEYfXnjJed7eaeL/jezbKuERFkpUdbr2tV57HbJrlpFbk8yOgxxptXCG/upCoxI1LGR69+5Nc+fOFa8dDge1adOG7r33Xpo0aZLP77Zr104IGjVR8/HHH9OWLVv8OYagiRo18x+/tspSo3ehmLEKmVXQbLX4ZJmxpSG6ZA1Ry8H6YubrIiqr8LycjPrz6E0SuoNmzb/bW75DjjY3eKydH/hxEcXvekDhA2N+8r04xZjVSQ1pwtwS+Qh1O+9c+np7PD23qitVO0P8tga59svLiX41SbONEnUxQMsxOzj73pZbZic7OUyJPqP9aWa5VE4g/WZEOBg9Vqkv//LLOnq/9yBDfe7r/vRX8JsVif70n5nrittg9Fj4WuHdGRHBrn4eqPkAGYhYNtretw5MpCta/TfghxAjfqD1OX+bylNTWVkplpGGDh1auwO7XbzeuHFjQA35888/KSMjgzp06EDjx4+ngwcPam5bUVEhOkL+ZzW8LMBWEhYVbJnhf/l1flG1bmZHs0ih3f7mvBk9JNa0SZDLNNw5TMNOq+D33/erf+Copq3fraTZs+bRlu++pIqK096b1ERq2WpuXI/PVCK1tNDLtWOr+et67Db6+L23adRDh+nxN3LohaffoPZ7xwtHZbXIKx4kjWI0ikwNSdCdV/Ekhf7wF7q85Cr64Pz+Xr+vd5ze+3X9a4W7v1x0qiVp5E+5/EYoVfpM5ihtp8xDpIeZfEr625KwQpm1YhkViEbTJHjsWycJpl6/6d4DJo5V6surMt413OdW94eZ/QfSf2auKzPHwtcK/yYvqent+/5Oj1BaxDHDbTCD0fbedsbzXhnl/RkHeemrIWNK1OTm5lJ1dTW1bNnS431+zf41/sLWnwULFtDKlSvplVdeoX379tHAgQOpqEjdOfbJJ58Uyk76Y0uRlUh+LMplH37NOVyGnG+9s5SagDFax4Md3Pzh7LPbGdrutS+jvJPhHVpOFcvaUtbB4fRAu3vFkwI/DandHFJB0NzKVh7vc0FQx4ClNOSav+i2wfhA46CpXe6mHlEraN3mYrq55VTLJl8pP5BSnBlF2Qa1AcXfycEqYSP8qGx1MxEqMZNPKZB+sgJ/BK7ZyVVJMI45M0L74dHMbwci+I3sP5D+M3NdmT0WoyL48pbLTLXBDGbaa7NgHOSglboqqdNoMwoPHz6crr/+eurevTsNGzaMVqxYQfn5+cIXR41HHnlEmKqkv0OHDlnWFiO1m9ZsKhWJ4ay02KgJmKDX++AokejWijrftfDkzSHiPJnLsx9z3hnn+usovCrbsOpnYTPu5w3CbM9Owfzvju47hY8JW414+Ywdp60aNPkm7Z7wg6WTry+rkz+oDSj+Tg5vH5wozOB1gVUToRIz+ZQCnUT9FYjye8IsZidXJcE45uyKMwxtp/fbgQp+I8fmb/+ZzdMV6LEEgj/nWK+9DhV/pkAfQuYuPWkqG36DFTXsuBsSEkLHjx/3eJ9ft2rl+RQeCImJidS5c2favXu36ucRERFi7U3+ZxVGajfxclFirMsxd9iFgVtttIQJ+4TcdZ3viSqgeh/s/NvrRfG/ektDfMxb/yin9z7Po1Nr7hWjvlnVLxUEXZt3NV19/RUeBUF5+ezDpzPp1iviRV4atRv3VGWSocPidvFNel7CRlODYNaZ4SKqxxdaVid/rSTKAcXsgCpNsgsP3UOjf9xMbx+430vcWJ2JyqqJ0PzgXCsojAzkHN1kdmLi7Qsquf9sAS2XquFPEkw5esds5jxLffnpkZsN97m/gp/bZYVI9Lf/zFxXymMx2qda26k5SGtt569YNrLEb/VDSG6+Q0TXNnpREx4eLqKYVq9e7X6PHYX5dd++5p3mtCguLqY9e/ZQeno61TVG/Vh4O15b/PKHwM1wWsKETXyvLKsNgQxKvY82o2h7+vtekzQ7rSmd7ibPz6HN331tukK4ksl/VS9/wH1w84hE+vS51iI8XS4w+MZ9cc9MnwOkEqNjPA+CU25Pphfub0XLn24tHOEeuTVJCFbOXyGHj1tudXrx0Fz6NeIREd0UCNKAojc5+JpkY2JCyd59CsWMyxHO3Y6+H9DR9Me82uaa7NXbYWQCsmIiVPumGf8r/W1ttDT7Dp8Tk1Z/Pr9nlqpwVbsnzGB2cvX63MC1oTyvaudT3penKdwSnzdfgr+wJqw80P3723/++PVJx5KjOBazItgM/opleXtzVa7ZBQfvD8pDCJftaYjLUH6FdN9yyy306quvitw0HNLNy0ScU4Z9azjcOzMzU/i9SM7FO3bsEP8/YsQI4QTMf7GxsdSpUyfx/oMPPkhXXnkltW3blo4cOUJTp04VkVD8vdTU1Dr1njYacfTsP1PpqXdPGk7Cx5P0ig0lnpklW4QIQaMmTPTy0/AkzMs2VsDH/ODso4ZC/YxGAKnlnTGbwEmekGzzTldCtxvTnqAxma8acop8YPtCerjzAz4jr9ivh5fB5FYjXxFwmhlO2Wn6+6/o0OYPaWTyGzVJDskwymgaI3lqykIy6VdOapZ6je/ssjVt+/rb3SLZV0LoSZraZYJmIkaGlyTVym8oo5/U91G7nVRuQsotI2XClfpRSjXA+W+2/OHKOMz+UEYyZWv1k3xb1bxGWv0Zmkm/xsyioqSrvTIKr96eQPPXdzOcS0aLi1t+QY920u83LWKjiW7L+oaGh0/2yL3Dx/xNrkqeGp5wbYo8NTX9s6l4BJXWRCsazU5uhNEXRVKnyJ88chf1T14V0P6ljN5GrzurQsHZ4jy+zRy69YwXvO5p+W8yasf32fEb6a9tnzcUmfT+oX9RoNgVCQDlCToDjUANVoh3vYd0MxzOLSXf69GjB7300kvC2ZcZPHiwCN1mx19m//791L59e699XHTRRbR27Vrx/2PHjqV169ZRXl6eEDEDBgygJ554gjp27FjnnWIkNwyfyIdvSqIHX8oxdeIZo2nwjbTBqnwB/Hs3PpqtmvlUiT95ZyIjiMYMjQ841bYkMOyHllC347fxs5HqdvKblAdTzUGQtUj/pR65YwKF27ju0/fp3FMPewxwWuG2agPK5NuSKCkhlPLyK6mt/Yfa7LDVvaltyM+66ed9tU1KyNgj6nPNSYYxMgGpTRKGhZaPNop7JL+S2tg20tHDh2jV1gT68cT5mgOuNJCnhJ+g3EpvQS593j7xBI3oVUYdOmS4Mwqb6U+lyFYmeouJtNFlfaKFUFBmVJYeYDg1vlqdLql/t5SN8BCB8uy9Hhl0a7IkS7XH9p3u484ozNfMF9/sdgsKRjnRJbcIp/emp7uEpaK8gZGMwq4sx57Hr3xI0yqfIM92fDjH4ZV1mfdz16hESowL8XqY4H3F5HgnXFRen/KyJuIasNWE68syMXPJkK4XDqXte0/rZsIeccaXdFfbaRTjPOJ9zspH0KjBcdS2lZ0+WfaFKUHBOKNa03thW2jx16Xu/F5muLBrJP3wv3LNzwf3jKKwYx+JZIJ6YlDq+1c+zDf8wB5oiHeDEDUNDas7Rc9Kwss+VaedhpPwyZeJjNQcMWotsjJfwDuf59M7n+uHxpvJOxMREUpjLo0LTlG0A8uINlwvblL5nrWStHmlbg9tTZF9XzSdLdkonOxv1acrKdJxjDKi9tNt7ic97QElmKnIfQmHsKpjAWUUZmuGv0LLTHvVSgjIJ3w1608w60n5upd173OFKLG6IKqRMSzQ68yqQo1+7Ucp6mquT62aUUZrSPmsscUPUTrnTK2u0mUZX9LDHf4u/l+eB8o9ctVkbVfLAs9wsEjXjuH0y+8Vnu/LRKRqjakWtZ9rFr2NbkP705+mvfYrPfrFTCb7R29Lpkt6+79qAFFTB53C6F0kRoUHO76yn4jmPmWFxsxmLQ70YvI3U7KeCfipfa9SRs/rg1/h9dByr9osWiZt+dP8yKGdqEf/y+qkqKMRy4j0pBfMonGg+aE3hoHgoCrSVOpIUXQbol6zvR6s/C3CaUgcOnTqpPnxoAtLTRCoy4zC8icxvSWilEQ7/XdmpiHlK39yqg9LjdmaVmrWj5zKdCo66zlq22ds3U3Osie2r7bF05w1Xam0Utv/nUPx1ZyU68oy4rGcZPHTOQDBsqYACzAhKBoC1eyWMPmIR/2pJutT09Coy9pPcowKFbM+MlwTaOabvgsvWl2Dw5/q40qntGuvH6HpcFtXaJpw8YQKAAANbhkToqYOOsVqE68ZywtXazaylhmMaqlm1lHl1JUviBnwhAoAAA1/GdPq+Ts04D00c/ikckSDLxOv0dw37Jz2xie+ExrxbjnPSzAEBO+TxdLzC08KcaUHH+FNI+LpphENzxeE22Nl0TUAAGiODDIwxzUkIGossgL4mkCN1nBib3u95R9ekuJQx2BfwFpVt+U8ZmGuHAAAAA2TkEb0kAhRYxKjEUxqNZz0fGo4fNAIRi0/gVzAt4xMFNE4LG6WrylSzbnRkJabAAAAAIgaC3xOWKzw+1p+LiwSWPT48ldhkcC+KUYwavmxUtw0FtMjAACA5gtEjYXVu7mSNS/dqE34kr+KZOVRZrbMyrqU+N2UxBDdEDq/q3I3A9MjAACA5gtEjcXVu3k7LQEg+at4ZXY8yJXHWtPOFk9TZeXA4FXlBgAAAJowpqp0N2fMVO/2BWeWbL9nnEeqasZZmk3nHB5PWVGfq36Pl6aCEcYNAAAANBUgagxi1I/F53acUZJTZXvU/yB3TRB+d0L76WJpSkl4GAkrDwAAAADUgagxiBTB5AtdfxdOkS2v/aE8GTYnpUUeFb42SnLzHWJpCwAAAADqQNQYRIpgCsjfhWt+GKBnwgZVa02wQ7kBAACAxgxEjQmkCCalxYYtNIb8XbiImQFuOmMOLezdXxSNrI9QbgAAAKAxgoKWdVlXiH1qPm1HVMpOwr7riDqcNlGGYNrOV2h93nDLC1gCAAAA9Y3VtZ9gqTEpZrg45dpNpeL14F7RInzbsNDgMvO9XhRyhkWLz01tno7DCOUGAAAAfIM8NUEsj6BKm1G0v+NCitkxkdIifPvYSI7Dz43bRVnntTf+GwAAAEAzBJYaA6zdXCJKHCiT70nlEVjwmGGv/Uoa9/MGevfgvYa2z2rtO5MxAAAAACBqdFm7qYRmvpmnWx6Bl6aMwn44DgqhXwv6G/uCQQdjAAAAoDkDS40P2ALz+Jt5pKdXpPIIZnPebC+4gE5UpGv61zjZVTi6DVGq79IJAAAAmpdv5+qfS8S/Zh6omwPwqQmggKW/OWTkVbvn7Z1K07rcJYQN+9DIBY2QOr1muxyMAQAANGss8+1swsBSE0ABy0ByyEg5b3ZWXyHCtnMrW3l8botuTTRwmXAsBgAA0LxhQWOlb2dTBZYaCywvuuURNJCqdm/ffQttz7+R2tp/oA7JuWSPznAtORmw0Eg5c3LyT1NBkYMSYpzUPvRH0/sBAADQeFcO2Lezf1ZUs89lBlFjgeVFK4eMkSR9/Jpz3RDx32WGf5P3/cHKAvrwm2IqKnWI9zgD8YQO0z1CxctDM+lo22fpjN5jmv3FDgAATXXlQPLt7CHmk+YLRI2OM6+vC4n1yeTbk1XXMtXWPuOibTRqSBz9ZXhCQAKD9/38ByepsEbMSIKGfXOUmYrDq45Qu93jaNaqArK3HU09OkdQUYmTEuLslJoYajwbMgAAgAa9cpCH+oAQNUacebWY/NdkGtwzRnPtU0lRqZPe+byQPlpbTBPHJfnl2KW2b5FxuMN0l3uxQp+w8zE7Id+eOZ3G/XAZfflDqfs7XA28XWIuDe7fkbL6XYplKgAAaMQrB8moDwhRY8SZ18vbvEWIWHJSEyVG1j4LSxxCmBgqgqnY99wlJ73eZ3HiKzuxlJmYt9ta0Nd7meogUcWRTIro+1KdOCZ7+QHF2im1BaxGAADg78qBv76dTQ0sPxl25jVWwNJM1JRZxy72ocktqF1ykkgOO2Ho+7yd1jJVWNURcq6/jhz9l1JI29EULNSW5WA1AgCAwFYOUB/QBUSNAWqdea1d0zTj2MViYMFnheq/WZVm6PdOVaXQw50f8LlMdWrNP2lH98toUM84Cka5icff8MzOXN9WIwAAaKorB80RiJoAUItuys6pMrUPIyJIb0lLykycEn7MI4GfBIuV3IpWwjajt0yVEn6EPlq6gsh2pWU3Cbf/vRUF9N4KT1FW31YjqW3CunbqNBUUOxqlA7VyOS8uxuZyBseyHgDNeuXAKEYidRsLEDUWLqPER9s9IpKscuzSW9LiOlJamYn5NV+a8/ZNpaSwXGNtCjthWc4DEam18KTwIzLr3Hxyzb30tn0AnXdWjPU+N45q2vr9V7R2wx7an58ihCH3ozxSrV/3KHe0GL9mR29JMMTF2qmo2FH7fs3ruhYSatehkrhoO42+OJbGXx5Y1F1AOKqJctYTlWYTVeQQRaQSRWc2yDxKTWmAb3RI10nZUVfNuwZ4fZi9Xqy+ngytHDh896O8Tfwg/tl3JZTrEanbAMYMP4Go8QOt6CazgsaoYxdfeJLfCQsOXm5STsLr84aLzMTKPDVsoWFBw59nJWw01C7ef05B4DkPtPrJqHNzavhROrp9La38oa+1N9qh5VSx8Z+UdTqbstoQURsSli4WhtxPDAsVjhSTosWMEpB/kGLid4Sn0N6TaXSguje1DfmZ2iXl0P6Tqe7XFUXZ9NGqaMoruICHOs3dch4jXrpcvKqIbrg08JQCZo7FUXqEjh/YSS2Ov0mR1Ue8NnOGJ9Gx5An0v6gHqKDEVqeWMjUL15Y/KmjD1jIqLnN6CdyeZ0VqClatiStQS6Dy+24rnB/9pOmg3yGUQvK+81tIePVjjdDXaqNaO7jPEk99TBmH/k2Rp7Pd2zrDkuhY6gTaEf0gJSeG6x6v8jyc0yGcduyt9Hpttq15+ZWqyVH1yhaofZ6SGEJXDIih9ORQn9eEUTEk365FnF21H4kz1Pd6USzpG3kQKimtol/Xf0k5v+bQ0Is6NaroWJvT6Wz01bAKCwspISGBCgoKKD4+Pqi/xRfQuMlHTJVQ0GLK39RDwpXs+3ERxeyY6CEClJMw0/vsCBraO5Kqj6+n/23fT9klqR7ihyfchb376y5Tjftlg/jO6CGxNOH6JL/7afzkQ5Tu3KgqxC5O+YQmd/mn7n5m7nyJvsm92uO9yDCii3pF+5xkNNt14EOyb7ieyOlpIZIsWiwM5X1qBrXkhxWh+v5B3FcHf15M6Qce9ByIpM+ddgqxOTRfq10LvoiKsNGYS+MMiUOjk7XcchV38hPqWTyJIqtlk5OTvCxycgqqEun53bPcxyCERLcoSk4MIbuNn0wjqGunCNcEJU0wSSfIXpnrJf6U7++r7kMFxbZaQRDjFPfIrt/30/4Cz3vECHFRNuqX5bLi8bH/tqeCftlZQWUVnkKoV5dI+t+eClXn/tgoov5ZfA2HUfuQH92CVWqr0TYmRDvpHwO30yXdCjy+L1+C5PPG52r52tpEnfJr9t6O0yklXJawMySTNsfOoqKUa1znVLJOqlglj+Seps+/K1Y9RrnA79KlLdlbDqTNf1TR99vKxH60lqLVrhPp+thSOkL0vfLe5+vx/S8KaPmaIo998+Utr/fIu/Y14Um5xG4cFk+LVhWK5KY9oj73uq9L7Rn0Lc2kVVujNR80r78klpauLvbxa9oPQuu2VujWeNJKwDpNpR9FcWQi+i3jffrnkgFBGccayvztl6iZN28ePfPMM3Ts2DHKysqiOXPm0AUX8NOiN7/99htNmTKFNm3aRAcOHKAXXniB7rvvvoD2WZ+ihquiTpxtLNpID0OFyA4tF/4lepMw72vhzAz3JKV2wcsvej7pastUykndbNi5a2fVdPTb6RRzcA7Fh+WrTr5sNXqh21jdXd2/fZEIQ/dFQoyd/jW2BQ3upS0QRX+sOEkjcrtRcthR1YFTKerMoDWYiH7lgVXFP0g6R9mbl9Kk9n/XHNCVgkD52qwgkwbRzNgc6npuG/GUmBSeS1VhreiA40LXpBhtoy1/VnhNQPx+x9ZhtPtQlYc1Q68f9JBGoWk753sdg9xKmRG1n65otZDSIo557UMp9tSuO7UB26wotAqttnyTcxVdnPqpbhsDPRaf12wQBL5W26SHLbbMal0zruvD5tUmFpgdWofRrgOVVF7pV1N1j0NVJNRcr/L3/LmO1Pop35FOL+zS3g8/DLPQVyZg1etHB49vla1o3M/a45veOGYLQj3Cehc1ixcvpptvvpnmz59Pffr0odmzZ9PSpUtp165dlJbmHYXz888/05IlS6hXr150//3308MPP+wlaszusz5FDZd7f+JtzwieQNEUDmzC/7QdUelh1e/JJ+Epd7TUzJujNF/H5n5MZ+U+7PF0dqI83b1MpVwiWzijViz5hNv72xNEO54hqvZ+QpEPlhvyLjNlNTLCDUPj6B+jWmj69bQP3WCZkJJjZDDJq0qn37r9ToN7xbvFzKJVRVRZeVp3QDeC0f5SG0TlBDLBG5mcfMFP1LkV6R7HoNZeLauPr/d5QlycfSeNyXwtKJO4WbQnD5c1QbTYRxsDFSSGJsBgCHyVthl9wOHzmKO4PoKJXh9pPVxM3zmPCk4naVpwjPWTk5Zl307fn7rUsCUx0AdFu97xks1VaPmqfZYuRVk9f5uu0v3888/THXfcQbfddhudc845QohER0fTW2+9pbp97969hQVm7NixFBERYck+65NgZGxkp1ye6LwQ/hXqgkaeVO+5cbs0rSmSU9mlfWLpukviadiFcdT/iptoR/dd4uLm5R3+lwcKtUFQCjv3Bbedl8gql6QRbZ+qKmik9vJRTmjPDsIkJk9bzU0sR+7cbGbwWvJ1Ec1ZnCesaVJ/chg5+/Wwo7KZfD5mkPyDtCZyl3/QEfpk2Rf018eP0JUTDwkfl/JKp+53jSJPsKiFNIjyoKUFi0zehrc1S6DHwrpZfgxa7dXav6/3+Wq4PvN1Tcd06brkgT3Y+HaSr22zVhtDqdKnk72RYzFyzepdT+aPTb1tRu833p8/bfIXvT5SPz4nPdblHiEueHmd/2WhoLyf9PqJ37u+9Zua31cj0PGtm97x8tkrPeSalxowpkRNZWWlWEYaOnRo7Q7sdvF640ZjTqhW7LOiokKoO/lfXWd2tBJN4cCOewbIau07g7EanIemU8/LhL8Kq3Zf4oEdJ9Vg4fDuinx6+j9viPpSYdXe2Y59DZaSc3NeZSuPbfgJ0d+n5o++LRHLg9c9dIAWvPkRbfjkLfEEw4OI0Xw+RrfzZzDZf+y0h5ncrIDyty2+BlGrJnirjoX3Y7S9RuHj4qUpqyfxupgslW28KuPdgAVJfQp8ZduCdb8Fij+/w6JUuQSq9qBg5gHA6INGoONbssHj/f33/dRkRE1ubi5VV1dTy5YtPd7n1+wL4w/+7PPJJ58U5irpr00bDmGp28yOVqOar4YjEYxgdDsFHLJthA/XFIklHDn8evRD2fTuZ6fo9syppicfvoHiY+x0yaib6PcsY1Yjo/DN/+q5/ejWqFH06Fm1T0sJoSfF8orSMiTB7/MyHJt7zRDIYGJ2QPe3LWYGUX8neKuOhfdjlQXLLHUxYQb6G5kRBwP+nYYg8JV5ttSM1Va0yV+s+h21BwUz14DRB43afrT5Nb4ZPd7Xvozymg8aEqaXnxoCjzzyiFh/k/4OHTpUL5kdrbTYqC5rpQ4kZ1Rrt+e6Nzai6Dau8MIgW53kS2TrNhfRR0v+j86P/oiuyXjbr8mnd6/29OFTmaIv2Wp07Q1X0i+l1+pajfTQWrLgp52pXSYIJ0wrl7wCHUzMDuha6A1Y/kykZr+j1w9mjqGunsbrY8IM9DeyK84I+HcCnQD9+U2t7aQ8WyIyycd94G+b/CXQ69nXg4LZa8DIg4a8Hx1+jG/61wRRQWUi2clBLy/LVXeZaGyiJiUlhUJCQuj48eMe7/PrVq08lxCCuU/2zWGHIvlfXcOT8V3XJVqyL618NRzW9/yuKcIhwPtCq3nda7bfTltmrE5iiezPUqreNp3O+98Z7jXjezrMMPWbfBs4o9vQsCuHezgfc39++HQm3XpFvMhH4w9G1vMvTvlUOPJxFICcvKpW9HH1W7SlbITp3w1kMKn9rk1T2CgHeuVr5W9cOziGhl0YTRFhgU2kZr/jqx/0whGUx2C1uOD9c3SUVh/X5YSpN3k4ddr46ZGbAxYkgU6AVosl11L0fCo8rT6mBtImfwnketZCEuv+CiY9sb++Zkk/148lfb3j5eE6ITyfnus2jl7qdKFIQ9Hok++Fh4eLKKbVq1fTNddcI95zOBzi9T333ONXA4Kxz7qAVeory2rDlQNBrRBZbeK6YVRY4p1Urzw0kyL7upIpBQKLCc5H8+GaYl0LyJk/TaKQkHyKCyBloxiaNIQY98HNIxJF/hQpOdfmneW0dlMpVRioPmG0WjlHJnBYozyZ4bXXj6Bre8bRVRolB/QzCo+n73OiqXPOw5SqkfxQC63EiRIOslMIOTRfS7+x03EFTbujNkXAg39xRVkt/qpIt5SGWuSLPxO81rFw1Mpnx2+kYns7uv7cnyi97HOiqlOa/WSmvYaiotiRPPsOEf3EwkYr67Z8whQOxkF4GPWdAbwml4pGZA238TSFG8ogrjf5G0nYae2x+W4b/x5HRY5rM5euy3jLIx1EIG0KBF/3pus6Mrc/Saz76icj39dr84a8y3wma/V3LJLg+zJ1zziijPAGV6PPr5DuW265hV599VWRR4bDrzlke+fOncIPhkOzMzMzhd+L5Ai8Y8cO8f8jRoyg8ePHi7/Y2Fjq1KmToX02pJBuf/LVsCXmrlGJ9MqH+YYKkakl+FNmFD5q70sfzGhjSdZVvWNxLen8Q/y/rzwSuktQEclEF7xm+ibQyrmjxJ+EfuzXM3FckmV1rtZuKhRRTmYHE/k5bhN3nK66oJw6dMjwmVHY/dpxoc9sq1JisiO/LqFJ7f/hlaNIjlXhzcrEa6HpgzzbKMs4vDcvhfad7uOVmE0rp5LWtaaZp0aWrkAtRLwsNJN+jZlFRclXe9TMkmefZXGtljTOLJHh5HYUV83lUp5O3+R656kpC8mkXzkhXvLVIusxt6VH1ArV7xud/KMjbFRa4XSfq5TwE5Rbae6aVePai6KpV+xKrwSMyrbxA0GHzDDac7iSimWxCNz//7ohnga32Syuj6+2xdP89d2ooNQ6Byt3GZQzI0Q+pnW/lnkkT1QSEVpNPVr8TDHO4+K+Zv+8uzvM8Oh7vv54acZMiLxeigW971tNaosQGtEvht77/BR1T/iBpnW5m+JC84Ma4l3veWqYuXPnuhPl9ejRg1566SWRX4YZPHgwtWvXjhYsWCBe79+/n9q3b++1j4suuojWrl1raJ8NUdSYyVcj5aExmvbaqGB6/r60gMoYSHC7rns4mwoU9ZmsyD0iCI0lOvvfROc+GtDFL/Xful9L6ONvSwLK07CtsC/dPDye/jLC+pIB85efEuHlZuFMv5xrJ1hlDPQyF/uaFLWsFrHRRP27R1OPsyI9LFf+ljuQp9A/VVhNSQWf0IWVj1Cs84hHxttTLf9KLc/o7HdG4fahP3qlvDfTNsmKx7mfOF2APOvxf78s9MpsKz3ASMUI3WUTYp2aGYV9tVGZvl/+/VNFRB+sLPQpwLg9701PVy8hoMwSrZJReMuuctqwtdRDjHg9pKkIV7VSGEbGRXnfs8BUlrKQ4IeUXl0iaPueSo9aRuw7OHJADGWmhmnWZ1LLSMzHPmqw655kPMo7KM5RdVkO2b8fUxv+bPBBIa0F0aTLf6PsX5fTiOQ3akQ8Gf7+hV0j6Yf/lWue68E9o+jH/5VTWaX29cDWes5wLblB8EN1hvM7Q+MpXbKGqOVgatSipqHRkC017CPCSypmCrYZFUyP3pZMl/TWL7NgBE7JzVFO/goFNarsSRR27r8CFjNqcP6Zx9/w7CMzZSAm357mMwNxwO3bVELPvn9SPAnrERlBNGZofN3UZFJONrkthJhVyygsFyfKmjl1WrizCRY5DDa+6q75nSm8AR2jZh0rE0IpKMd0aDnRpn955Bcrokx69vfHNK1n8odetYcO9kH8lp6gF74b4mGplotI1RpULWo/17J2a60W8P7WffSmIcs39VtI1O5G8heImjroFKtqQKUk2um/MzNrbwqVC15eaKy+LDW+ftPoko7E3L2PUYU9jYYOCn4RNLUbWa8MxPOHXqXew8ZbttwUyLIZP1WOGtI4K+GCho/eRAfqToSr1nLSOhcaIt6KauDVJkTb1u9WUtZBA0vQsNQ0DVFj+mmIBQ3XcPIqp1ZzQclqavCFd+Oj2apF4vwqXxCASDOTwjy3Kp2+SNlO40ck1csTm2SWPi/a29egPLQ1HW37DJ3Re0ydC4hAKzUD0FgtRqARnwtHNTk/aUdUlu2xnOaZUqSJ+NQ0NOpL1Bh+GtKp4aS8OKRaRbxWH0zzsRGRprekw7jetakWbaxrlL4GZv0mAAAAkOJhnHH6fBj3F4iaOugUyxX48bVEq4fo7+iSNbTuyAU+rT9WR+uo+arMfDPPI5+H1pJOoFFNAAAAGjiH1Nwm2rhSc1gw5ls9fweQcQQoi0ZqYrCGEztustXHFxFhNsPlDfxhcM8YIcgffzNPN3dBpT2JwoPkCAwAAKAB0GYUUebVjcZZH6KmLtZGDdZm4rBHX47HDH/Ov22Vg7AaHBVkt9s8ltXkCZ0498jg/h2D7ggMAACgAWAPCShsuy6BqLHChyYxRJQb0FwSYlXLPjOlHKan7XDF4bRE+f4Vv7QYPhZ3Tg2PsMmrG76DGwAAgGYJRI0F0U4scPh9TeddVrkcti0crkQSdNUaTsll4f4Xv6yPZTUAAACgAdEoq3TX15KTnr+LvJK16roke4pHZ3q8zTWc9nVcSNWZ1woLSIpO1Wyt4pcAAABAcweWGoO4lmGq9StZ+/J3qXG42vr9V7R2wx7an5/irrOS+skRGnJ+NFX6SGWtVfwSAAAAABA1hjHqx6K3HWeWnLbwHCLiv1pYMPmqFxTsUG4AAACgsYPlJ4MY9WPxtZ2RJSwtwsMoqKHcAAAAQGMHosYg7MfCUU6B+LsYWcLSIjffIb4PAAAAAHUgagzCfiwcth2Iv4va0hSXIeD6Slw4kv/l12a+DwAAAAAXcBQ2AfuzcNi2v5VvlUtTXH5AmaX3REU6zds7VbVMfV2FcgMAAACNEYgaE7BPDDvs/u3qBL+qLUtLWCyIpHpKymR8XDiS3+eyBHJhg1BuAAAAwDcQNRZkEjYaYi0tYT3++nFhoWFBY1N8lQtGOpw2mtB+uihLwOHeDEK5AQAAAN/Ap8ZEJmGlk6+USZg/NwovUT03bpdYclIKGvdJsTkpLfKoqLPE3DIyHqHcAAAAgA4QNcHOJKxCVmtjYd3JYSfEv63TwgzvGwAAAGiuQNRYmEnYMAardudVpYl/4SAMAAAA6AOfmjrKJKxWtdtZmk02lard7FOTW9FKlFCAgzAAAAAJXhXgh2iec/iB12igSnMBoqYOMglrVe22rb+OXYU9hA0LGr485+2bKpyE4SAMAABAL2AFJXRcYPmpDjIJ+6rabVNU7WYLDYdz73RcIXLi4EIFAABgZcBKUwaWGoNh2HzRaOG3NaWmajflrCdH6RHam5dCBxwX0rUDwg2bFN2myPxKamv/gdol5dD+k6liP8mJxvcDAACgcQes9M+KavbjPURNMDMJO6qFYKGyoy7nYPal4aUnOfy65WBhMuvUnqiTiYv8g5UF9OE3xdQj6nOPzMS8j/iazMRTy0ZQv+5R1LNLpKlEgQAAABpfwEqPzpHUnIGoMQgLF1bBhh20Di0n2vQvotLD7rfKQzPpaNtn6YzeYwISFmxmfP6Dk1RY6jCUmfjLH4bTVz8Uibw37RJy6Kyz21FY+iBYcgAAoLkGrDRRIGpMwELEkApmQbP+Oi+hEV51hNrtHkfPrymi3sPG++UvI62rMlz80khmYhs56O4OM2prTJUQnfhNsuSMpNEXx9L4yxNgwQEAgOYSsNJEgaOw1fCSE1toVEK1WWjwuzelTRWlEsw6dvGS09wlJ92v2fJiJDPx1C53U2p4bdFMuSWHl64WfFZIox/KrjNHMz6OLX+U01c/FdOy1YX01Y/F4rWZBIYAANBcCFrAShMElhqrYR8a2ZKTrxII85YNMOXYxT40uQUOr4zDRvBlySk6HU9JYbn00ZI0IucIGtQzjoKB3A+oqLT2OCTiom00akgc/WU4rEYAAFAnAStNDIgaq2GnYAOwINmaa9yxi60obFFRyzish54l5/lu49zv5W5Pp+rkORTSdjRZidwPSIuiUie983khLfm6iMZcGoclMQAACDRgpZkBUWN1ZseiPw19VxIkRhy7tML5OOPwiYp0sZTEAkUJr+aYFe5JocfItuF6IvsyV8i5BazdVEKPv5lnePuyCqcQcMu/KaaJ45OCe7PWRKhJIfX7qvtQQbGNEmLtlNqicUSLSddhTv5pKihyUFyMjYpKnI3qGABoqliZAdh0wIpRjETqNhIgaizM7DjyjC9p4hnTREZgLeQlEIw6dmmF83HGYXb2Zd8Y3q9c2LgyE5v3UXH5/diINt3nyqET4IW9dnMJzXzLuKCRU1xaRR8t+T9KyS+lc85pb+mNxgPNwZ8XU/qBBynydLYrpF4WCr8+b7h7SYxD4nucGSEsSZJgiIu1U1GxQ3wu3q95XddCQu06VBIXba8/Z3ClaDzdhwpKbPXWX4YG9OR+RHnfN6gB3tLU+BZPYFa2ra5KANRVm31lAPZXnBgJWKnWOz75NcAP4rtfJyqzPlK3PrA5nc5G751ZWFhICQkJVFBQQPHx8UH/PXkEkgRHIi3s3V845Got90g9PW3nfDFpstlw4YwM3Ytm9c8l9MTb2sKAw7rleWqYE+Xp9PK+x0TUk5YlR5dL1ogcOoEImsff8E/QqB0ThScRnfUvonMfDWgQ5vP385cf0MQ2f/eKHJPKVHAovCRs/MEfIeG2uJw6TQUsmHQE1JHcKvr42xLDbYoKt9ENlwbXZ0luNYrJ+Zh6Fk+iyOps9+cnFKJRIjaKqH9WNPXoHOGyMsXZLc+rpBzoz+kQTtmblriFrYST7CJiUKI8JJM2x86ikpRrNAWY1iSiPKdmj0vLD03pf6Z2bDv2VorXLeJc8SD8/zG53ufEGdWa9mc8Q7vpCkNtlP/W4RNV9Pl3xcLXj8dAkTYiMZe6dGlLoRppI5SWRRa23OZtuyvoo7VF4hqXSIh20N8H/o8u615I9ugMXcHpdR46hFJI3nfu7auTB9D7XxbT8jWev6Plz6fWVvk1oCZapHv/jFZhPse/yHAblVc6tccMHfGpdm0lx4fQ1t3l9PFavl6c6qUUVNKNKJHGwVn75lPGeTcEdcywev6GqDEJX0jjJh/xejLOSthIL3Qbq/v9tw/cT+8duk/8/5S/JdPgnjG633nn83zha+ILHlD+2ncb3dC/mL7enkCvru9KBaV2dx4bvrzNCpvfM96mswffSv4uObGFxp+AJnnuHTWBWGlPop9jZ9Pv1VeQ3cZPLRGU1TnS0E0nLEdvnPApQCVr2rhfNghrWCDER9t1l9D4mnr/iwKvgdZqpEmndVwOdT23neak44HWwKqwwByo7k2nj22gnTsP0P78FEoIPUlTu0wISDRye/uk/ULDehRQ/z4dKSS1v5jQlEuFSrHnfi0ThVt2ldOGraVUXKZ/nfHDh682SxPt0G4FtD8vif63p4L27D1Mh4tShQWWrxne5ppztlD+icO0v6D2fffEX5MvKqTVQM1j2PJnBX27uZTKK7X7KCrCRme1DaM9h6s8rh1uv/JxVet4tc6Jx3GeTBVWts1/VNH328q8rlPVB6saAbulJgEoWzv5mNS+r4baPpWCk6IyiTrdSY7YTrRqW7x73JO+f2/H6ZQii/zMqUinuSqiWiImwkETBv+PLuuWT3v3HqEVm6JoX36a+/y5t4u00bkdwumnHRWa7ec+dfpzf8bm0Ojex6ldyQIP6wkXQeaagewWYMRCq8ZLN3xHXY/8xVDL5ONgbEwYTRwXHFeABiFq5s2bR8888wwdO3aMsrKyaM6cOXTBBa7lFDWWLl1Kjz32GO3fv5/OPPNMeuqpp2jEiBHuz2+99VZ65513PL4zbNgwWrlyZYOz1HDo8cTZ3lFHF6d8QpO7/FP3+zN3vkTf5F5tuBCZmlVIDd7Xwpm1Vh+5ild7OjPC/dsX0bU3XGn6Qg7EQuOPxYuJj7H7vOm4P95bUUDvrSik7gYFKB//1oK+ZAVTbk+mwb0UAtZRTVu//4q+XrfbY0IMBr4mnV+KR9AFZ4fRjedtpQ7Jue4JLCbv/7yf6sOTKD/mYooq/J4iq4+436922inE5vB4bSeH36JRrb1sC+B9Ktvvj0XNyHWm1mbJ+ulhQZTBbfom5yq6OPVTr77Wet/fY7DyeJXnxNf1IrVVmoT7JX1F12W8Kd6zyuppVHDyUCA/HKmNjBkBp2kd1jh2q1EXcJ7HJtwCiOi3zPfpn4sHmNq/naopK+EHmtLlbooPy/fpIuFrHAxGPcJ6FzWLFy+mm2++mebPn099+vSh2bNnC9Gya9cuSkvzjsb5/vvvadCgQfTkk0/SFVdcQQsXLhSiZvPmzdS1a1e3qDl+/Di9/fbb7u9FRERQixYtGpyo0VoKygpgotS6ULSsQmroXmyyJ+uvtsZQr/x/UFLYcQ0H49oBLrlFuKElMqVw8NfeYLQf2QKUW5HuNTGqiQcRebXwJBWWOPwWoIHCA+ujf02mi7ltfC5+e4KqfptNYY5TQR849Z7QF2ff6TXZFlQlUnxovrvtSrwmF8VrI2iJRn8tKMG4zpRIo6XWsfJ1KX2kbLtEMJY79TAzPvF51xMEjJYACNTqaVZwKn+Pv1J4OkEch1FRrWcddp0/W1DOk95vK9tReDqRpu98mbYVXOjRp5LITA47IQJRpIckX2LNCB4P4gZdJsxg9fxt2lH4+eefpzvuuINuu+028ZrFzeeff05vvfUWTZo0yWv7F198kS6//HL697//LV7PmDGDvvrqK5o7d674rlzEtGrViho6Wo69bG5XPq36chA2UojMSL0P5tYr4vXVs6zG1LD2RNUHTosoJ3UHY6J5+6aKG8JIPRFp+YRDsTlyKRCM5t7hrpLy/cgnxhlv5dGBo1X0lxGuNeB1m4voo6Ur6PyoE5QX6rrRjYbCG91Ot6082MT/RN9/coJs2w9T/8jXKdxxisIU28lLW1g1cOplneZBckzmq17f8yVo1N43O/lonWtf7fWVa2lD3mWmrFxmcjwp8XWsfF2qCTxpOcjKYzCD0eNNCTtGd7R/ymeW8okdJ1F8WIGhJQx5Xi6jVk8pqag/SNd0Qli+4Tb5uuYk+H0WrFafJyO/rWwHHxun4ZA/BGlZ1r7JuYrGZL5mciHMk8TwHPEgKITSqQsafH0pU6KmsrKSNm3aRI888oj7PbvdTkOHDqWNGzeqfoffnzhxotfS0scff+zx3tq1a4Wlh60zF198Mc2cOZOSk5NV91lRUSH+5EqvrjM7ysUGX1CS/4CvpzdJKCjREg5G63hkpiqnR304D822ve9Tq/0PetwILLy4nfJJdcPWMtWLWHJiXPyVf2KGb9DuncJo14Eqt9+AWSGhHKx5QHtnRSEt/rqQ/tZzNQ1yPEqDunne6C/vfUwnFF5bgFpiVubTqjKABWOC05sgjIqWYKB2rs1OaP5Mmlq/bQSjE4+Z9/09BjMYPd6E8Dyf/c9tTQjPN22ZMyMiAxGcjNF2Sb9j9JoLxnkKRMBJD0FsaVUTLinhx9wPLP7cz3yOecn3ng4z3O/xuHks+1mizrW5zRq1qMnNzaXq6mpq2bKlx/v8eufOnarfYb8bte35fQm25IwaNYrat29Pe/bsof/3//4fDR8+XAiikBDvgZ2XsqZPZ3Vb/5kdjShtdm57fOdcn0/fagIm2PU+HJmjaNzSvqomSzkfrikSYk6yBkliZtGqIg/vfbM8xktFPWPE/t5dkU/vf1Hkzr1j1PSsNVj3jvuCrg25iyjE+0ZnASoNBHqWqkDQKjbq+0nf2oEz0AkiGPgSjf621+z39HI8KfFnea0hnSv9nFauc1JQqf4gqcRsX5gRkVZZSI3+jtl+t/I8BbIv10MQ0fWZr2ta1vxFWi6V+7AxfP2kHv0L0aFIy/KYNcnaT2PHjqWrrrqKunXrRtdccw199tln9PPPPwvrjRpsKeL1N+nv0KFD9ZLZkS02evWXGF6SKjidZFqYsJBICWK9D/4eR7/w5MlrpvyvFKHBa/BscuR/+TUvkbH4YP8UrhPFyfH8FTRsopdHfrFQvO2KFkLkSLl32JThy9uLB2EOW1ebGHWXXIRfzac0fec8yq30XPLkgd2K5R+zZuVgDZx1NUEYRU80+ttes9+TrjNbTZvkKK875efBIpjnytfxys9JbpW1LgC+7lM9ARZov2tFXirbZLbfrTxPge6Lx1KeXwIV3E7lNV8jDbSEkpPzmLFvYGMXNSkpKcJywk69cvi1lj8Mv29me6ZDhw7it3bv3q36OfvfsEOR/K+uYWFz13WJhicfX9tpCRNe9qnUEQ6B1PuQrE5K6wI76bFTITvT8r/8uov9M3p/RYGwUPkqdWCEyX9VD2Xn91gs7qy+QggLdojzZ2I0WuiThea4nzcI50h2hlsZ/X+09dzfKfrM0SJvSiAYEbvBGuyuHRxDwy6Mpogw6yYIfwdK9jMzIxrNttefSVOC28BtUQpbaUCXt5kdM422S0uMa70fyDFYcbzyc6LX/2ZSNPhr9TQjOLV+t6Cqha6Ak9pUe8z6+7X6PNXn/Sknp6KVSDfC4+DcvY/5FEqc1NVWekhEbjb65afw8HDq1asXrV69WlhUGIfDIV7fc889qt/p27ev+Py++1y5WRh2FOb3tTh8+DDl5eVReno6NVTYavHKsnzKcAbudKomTPRCufVCmI3C3x89JJY+XFOsuVzidmD9kV/5b8HgQ5xcs+Tkqz2uTJu30C/5N9I5pc9S8rG5wrHWl9+PHDNCkwc2tlDJo6Yu7aNSekAl94kyId6vO8tp1Y+lAVlaeNDmXBr+DJzKGjAP/sW1TPjm5uk0qf3fNbNOu3IYqbeFMSLMtCKTZuycK8SjlJDN3mEgdWtro341uWOUOUt8ZcnW+o1Algr5GmL/JbEEG36C8irT6LeCXnRuwiZKCT9BuZW1S7K8jKzWLjX/ObXoMOlzj8SuFi53mj5elSVnI1nKOTIuLrRAd3lD7z7VaycLLbWw/hDZkojWNfH87ifFa+X31drkecxOzXshGOfJzPVuNU6NSCq20Bvh/1b9SQUxgxpczSm/QrpvueUWevXVV0VuGg7pXrJkifCpYV8ZDvfOzMwUfi9SSPdFF11Es2bNopEjR9KiRYvoP//5jzuku7i4WPjHjB49Wlhv2KfmoYceoqKiItq+fbuwyjS0jMLyfDVS+KHWWjXnFrBFt6Z1Z/xG85YVGipEZiSUOyXRTv+dmWlJaB0fy4OzjwY9IZ1qrhYj1ORzWbthj0jsppfPxUz46r7T/S1NKvUNJx18M8+vkGEpbJSzeNrPGOXOrquVUVh6bSQDrLwkhAQ/dX6T64qOUCZnlCawsuoYig71zlqsHHCVkX9loZn0a4zvLLxaWZS3/FFBjgPL6fbW0zwmJOVvcPv9mTSVien4N0cNjqMbh8W7s/DKs/KysN280yXAekSt8Bkey4J0tUo+mrKQTDoQMZraVXzomWHZ4DHwQ8yoIbE09rJ4+u+XhZqJGo0mfIuOtFF0pJ1yPTLh2qhja1ciP7XjlNrKqCf0dP26o/N9tCr7Eo9EeGrw73XI5N+r9EiKyFl7770+kZISQikvv5JOH13nTuooCU4WZBlR++mKVgspLeKYVxu/OzlcnGMpzFkpUOX9uWhVocjY3CPqc81zW0SZNGf3FPr62OWWJduTM3PkOupW8DDFUW3upxPlrWjF8bF0bcYCigvNVxdb7oSEPNOozD81bykfXLRSCfCc9P+Gb6esg8MNjZ9H7AMCDvGu9zw1DIdjS8n3evToQS+99JLIWcMMHjyY2rVrRwsWLHBvz3lsJk+e7E6+9/TTT7uT75WVlQmrz6+//kr5+fmUkZFBl112mQj9VjoYNyRRI89Xo5W1V1w4NiLbQFdxSKP1RrQS/Cl5/r40S0LruF0zn1xGU9vfEJSEdFZZleTWk82/l9OXNVYRJXpCk89LXlU6rUzdRuOGJ1me/nv+8lO07Ot8n21QgzMlZ7efF7x6K/JcRdviae7arlRS4co6rTWB8VP9+DZzaHTm2+4wb+nzz47fSEX29nRWl3YU1qo/tQ35WSTvE+nsragn9Gcp7dvyDe3dc0AkJ5QmtHaJOeI3fWXjVbOqSSnu5SUEzNTccd+/+ZXU1v4DtUvKERmFC0oclByRS23anSGOe/ve0+5tvPpDo4CqWcGqlb5fOjbOnszWVy14mVer9pDqcdYkZOS6XdymTs7PqN2Rf5PNI+NtG6Jes90OpFrWTuUxGRkXzZZmUJ5jvXOuPOYOSSfIXplLFJFKFJ0pzl812X0eT35xtV9JR925X8ih+vA2SETXaolIIjr7QaLfn615z+nxOb860fI+Sshb4vVAw/c3L/WPHBAjImjd/cKWsE/bkbM0W1UoKR9wA52HGoSoaWjUp6VGQmtiKDn3eWrfp/aJ3cgNrFfrSeLR25Lpkt5+WD5U+H3tAjr7iCv3kFUJ6TiF+5hL44JWSFGZVE+OntB09F8qwtqDBZeJ+PGLhfRQ+797tUFp5WAxE35u4DWtzMLX4tY/yoVlxOaspv4Zv3hkFJYmMDEB1dTQcZdFcFyoX2LBwnbWRZHDpoZqMUUN63Bjq+zcUK8JtT7nhzq1MUorcarRGmUeIlKtnlO07HNlWRO9+/fQcnKuv85l7VKJDpVbeAKdhyBq6qBTjKC2RKTM6HjU3pc+mNHGUMVW+UVd15YawfG1RKuHWGKpCbaYUU38x7lyFE7VakKzPLQ1RfZ11U8JNr6WfVxWjg40uH9Hyup3ab1XgQZNk4Y6+Te3Puegj4AFpp6ItFhk7vtxEcXsmKhqwZUvWcFS00REjRFnXrkCN7Mt3xQ3PpotzKtaWJ6umm8IEyZHLTFzw9DgVoE2W81Ybp6uNwFh9ikJANDkaGwCs9rhpPGTD1G6c6NmHjMr5iFYauqgU6w28Rpx/JVfHL6WVSSCUVjMjMlRDn9204h4uqmmNEF9onQ81XOiBQAAoI6Zh3F/gaipg06xWoGbWU5iIVMXodxaVB/4kE6tuZdSwn2bHC2JagIAANCsfbMK67ugJfBeuwxJHejTt8VoDSe2LrzxCReK0yYizCaiFoIFO8/u6H6ZKALpq3RCXQgsAAAA9csgd+6wxrF0BlFjFlUv89ZEvbSdT43WZuLlEr2q3Px5sKukDuoZR2S70qXOVQSZlNejPnxnAAAA1C0hdluDrswtB6LGrKBZf513iqXSbNf7NflojFT2VsLmPPb/MIJRy49V6lyZC6Mhq3QAAADNF4gaM0tObKFRzRnJ79mIuMhX5tVe0TXKyt5q8PokL+cEsyp3U1bnAAAAQIOo0t0oYB8a+ZKTF06i0kOu7XQqeyuXcm4ZGS+sIsGuyg0AAAA0ZWCpMQo7BQe4nbSkI8+nwmnc3/m8kFZsKKEh50cHtSo3AAAA0JSBqDEKZ2i0YDvOLLngs0Kv99nfZsnXRZrfQ6QRAAAA4BssPxmFU05zlJNURMwLm6vWBm/nI58NRxT5Q3gYBTWUGwAAAGjsQNQY7qkQV9i2QClsal5z8TAfKfhdkUT+RS7l5jvE9wEAAACgDkSNGThcm8O2uRS9HLbgaIRzWxmKXReh3AAAAEBjBT41ZmHhwmHbflRDDTQUu65CuQEAAIDGCESNP7CAaTnY9NeMJOHTAqHcAAAAgG+w/FSHSEn4/AGh3AAAAIBvIGrqGM5Vc+sV5iqRjh4Sh6KRAAAAgA4QNfVAZmqYqe0Ryg0AAADoA1FTD5hx+IUvDQAAAGAMOArXA2YchuFLAwAAQJ7ElXOWcYoPfkDm+QSlc2qBqKkHjFTtRlkEAAAActb9Wiqy0ssfiPkBmecT9tcEWH6qN/Sqdn/4VCYuUgAAAG5Bww/CSgs/v+b3+XMAS029mgmlqt2BmBKl38jJP00FRQ6Ki7ZRUamTEmKd1D7kR+qQnEv26AzDCQIBAAA0LIzUDZy37JSYT5r7UhSWn+rZTMgXYI/OkeZ/wFFNW7//itZu2EP781Noe8EF5CCXaBmY/AVN6DCd0iKOujd3hiXRsdQJtCP6QUpODMc6LAAANBKM1A3MOVUttvNrPmlCQNQE0UyoRDITTrk9mQb3ivH/Bw4tp4qN/6Ss09mU1YZLNxCdqEineXunio+ndbmLZYzHV2xVJyn9yAyKrppDz++eRVPLRtCoIXH0l+EJzV7ZAwBAQ8Zo3b881AdESHd9mAlnvpVHazeX+PcDh5aTc/11FF6V7fF2SvgxIWYmdpokBI1NwwIZH5ovtusRtYLe+byQrnrgML27Il+0GwAAQONNA5KM+oAQNfVhJmT98PgbeeYduxzV5PzlX6qixW5jUeKkhLB8TUHDuD5z0oT208lO1VRW4aQFnxXS6Iey68zRjAXUlj/KafVPhbT7l1Xk2LeQ6PhacXwAAADU04D4AjnNXCD5nsWYMf+xY5cZCwn70NjKDpOWZvElZpTbpUUepW4JP7nfKyx1iKUxvy1IBuBjZavQqIey6aP33qFuv3WhTn8MI/vG8USrh1D5sra078dFsBoBAIDJuoHIaeYCoiYQ2LLAFob9/3VbGrJzqgx/XXLsMgJbUf5v1Z9kJclhJ7zem/lmHq3dZL2w4fazNYitQj2iPhdLYKnhtY7MTHjVEWq3+0b6Yu6dtPW7lbDcAACAThoQttDw+4HkqamWrOc/l4h/G7M7AhyF/eXQcqJN/yIqPex+qyI0k/Zsf4yIhltq2ZH8dDKcaYb2ydejkai+vKo09aWxN/PIbrdZlieHRRLvk+ElL47M0l5CI7oi5Q2ig29QxZFMiuj7ElGbURQs3CHxp05TQbGD4mLtVFTsoIRYO6W2CNWMEmtIWT29wvpjbFRU4tQ9BgBA3d2fVowVRtOAVJv4TXmkLo/PbMFvl5hLg/t3pKx+lza6VCAQNf4KmvXXeUUYhVUdERaIaTtfofV5wy1z7JL8dPLoAhHlxE7BkgCQ43DaqOh0IsWHniKnU3s5irfLrWglwsCDnfOAl7PYMVqCbxh5qLkvuD/ZKXpby/fI2Wa0pZMz3/Tvf1FAy9cUibw+SqSbu33iCRrRq4w6dUgnqswjR3gKrdqeSK+u70oFpbWGzrgoG/XLiqKeZ0XWqZBQSx2gPIb6HKC8RGON4JLEozuvUpydUhPrR4ApJ4BzOoTTjr2VDUKweliFc9YTlR0likr3K++UUvxaJXqt2G+w2hbM/jQiIHyl9vA3R5mvNCDVGuOaMp2I1OYNW0vpwzXF6qlADlKdPFhajc3p5OmvcVNYWEgJCQlUUFBA8fHxwf0xvhk+bedhoVETDON+2eDOG6MFmw0XzsjQvZDZJPjE23nuC4+FE580ubDh3+W9sKBiOAqKnYbV2idtpye8nr8vLaCcByxo2CFazsUpn9DkLv80vA95f0ZHhYqBoMdZkYasKVrwQPP8wpNUWOJQ/Vwtz48SKYReqw/jou00+uJYGn+58ZB5owJAen0kt4o+/rbE8DEUUQbldnyOzug9xvKJQs3itWVXuRg0i8uM7yc2iqvSRwd8jo1MmCxe/vtloXsCkERgSvgJyq1Mc+d+Soh20N8H/o8u614oEllWJw+g7XtPe05G5BCTpKP0CO3NS6EDjgtd+aA6hFJI3ndEpdlEFTlEEalE0ZnGJ1Eeb357gmjXi0SVJ91vl4dk0ubYWVSUfLXrGlGKQ+l3y45SdUQrWri5Oy1bU0ZFpd7XfEqCnUYOiKWMlFBx7qR9KMWd9Frej0dyT9Pn3xVTboH3ft1Cv3MYtQ+tTQTq7r/8Smpr/4Eqi7Jp6fextD77fK8xk4+rX/co6nFmhMe1r3ldKMSK8lzxMWRvWkLpBx6kyNPZunm89MSWmmiR7v0zWoV5jX9yIsNtVF7p1Bwz9KwtynvuWN5pWvlDCZWWa0/pvEzFKNsszStKC7qYL2xEjv5LKaTtaGoM87dfombevHn0zDPP0LFjxygrK4vmzJlDF1yg/dS/dOlSeuyxx2j//v105pln0lNPPUUjRoxwf85NmDp1Kr3++uuUn59P/fv3p1deeUVs2+BEDfvOrB6iu9n92xfR1oK+Hu9Jgyb7svDSz7XXj6BBPeN09/XO5/ki/NrXhHWiPJ1+inqSwtqNcmcU7lr6HLXKm0c22WDI283bpz0Zyxk9JJYmXJ9E/i45sYVGuTSblbCRXug21vT+1PpTIiHGTv8a28JQ7p91m4voo6Ur3OdAnrTQ182txHXX2HTFYXy0nSaOT/K5lKdnNTKLzwGKiGbtm0/2M0ZRzy6RmtYRvSUt9wR36jRt3lVO32/jCTN4z0cxUTa6tE8MZSSHeok7zdcxNtryR4WhtqneUxXp9E3OVXRx6qce7+dWptOcPVNpQ95l4n4e3PIrujTtY4qmPK/vDk37lFIUvmMeoiTlGk+LVU0m8HZJOZRzcBclH5tL4Y5Tph5O+Fju7Tjd43f1RLgavH/x8KQh9oyg1q/cf1+f8O7XwqpEWnbkr7Tw0D2G9y8XTl3Ln6NWOfNETi7luZKOW+/+LqhKFHm8tpSNoI6tw2jP4Sqva4f748KWv9CgLvn0xZZ4zf4IoWrqKhvvjfYbi52+3SJp+55KytVI3urLQusL6TqTH0v3hB9oWpe7KS5UPXKWr7WTVem0o/tOQ/NVoxM1ixcvpptvvpnmz59Pffr0odmzZwvRsmvXLkpL8/bR+P7772nQoEH05JNP0hVXXEELFy4Uombz5s3UtWtXsQ2/5s/feecdat++vRBA27dvpx07dlBkZGTDEjXsFPz9ON3NZu58ib7Jvdr30390a6JeL/o07Wkl8lMKpKP2vvTBjDbeT7M1Ty78BPnVtniau7YrlVQY9w/3xwFNzUIjb/fC3v01l9CM9qcaNwyNo3+MUo8Q4En6wI+LKH7XA5qDvdQ2dmA2EklmxiqnlnCR2/TBygJa/FWRCK23Ar1jUGtzVISNzj87gs7tEEEtYokcx9fTzp0HvDJVU4ATXENFWwS6JnbGWxw6qfB0oqo1lJGPqmrnQU0Uq2YCN7CMLD+XeoLWzNK4VpuMCiSz/aoUFkbbyb9jxDLNIlTv/tZ7WDHaH4H0m95YH91mEG38n/GAFC2MWKSVD5bX3nCl5TUJ613UsJDp3bs3zZ07V7x2OBzUpk0buvfee2nSJE785smYMWOopKSEPvvsM/d7F154IfXo0UMII/75jIwMeuCBB+jBBx8Un/PBtWzZkhYsWEBjx45t9JYa7aeDmhcDl6kKG570xk0+YkiNGxUfvM/3VhTQeysKFR5BJsWSH/uW0FpC89dSI2fyX5Po4vNjvYThz19+QBPb/N3nYM+TlJVWJHlfnqpKoyuvG06De8W7xcyiVUUe5mcrMGoJU2uz3kBs5UDdUNATgVqiQho19ax5vj6Xi5L+yasMWQh9nUt/BK0vAhFI/var9Bkzbed8v4WT2nHP+vM5er6b/gMp/35ORbpXPxntDyuFZbDuuYEGLdLKB8vt1aMMuUyYwer521RId2VlJW3atImGDh1auwO7XbzeuHGj6nf4ffn2zLBhw9zb79u3TyxjybfhA2TxpLXPiooK0RHyvzqD18LZwqKRLYYvXF7ikZxwfUX7uB2NN92nGr5sJJEfc+sV8YbVM1+Mt16RSI/dnqx6ofNAxJMi+73wvy91upAO/rzY5z55kuYlMld2Yt+ChuGbkW/s3MpWuu1V9qcW3M88oX//f2/TlvVfuPtTWI1eP043tZyqGXHF7eVkhClhx8iq0HhlXz7XbSyds60LvfXym3TlxEMitN1qQaPVFiPbSYOcMsxeylR9Z7v/+Pycv98YkRzXtQZ2X+/rTQZ6n/O1x/mishJ+8DFGGD+Xesci/Z48P5Vm23SiFKV7hrezsl+lz/T2r9dGteM+L0F9PlH7fWU/Ge2PUKoMqN/M3JP+3nN2g/2mhB9yzaQhqS9MiZrc3Fyqrq4WVhQ5/JqFiRr8vq/tpX/N7JOXqlj4SH9sKaoz2LmPl4wEnleEpMTZZ0VS+Ho3t7h9Sw+5nNv8TOSXmRpm7hiIxFLI6CFxhm6gdnvGuSK+NJLpsZhhnx8zSygsbMb/soEWlC2nz3L+Jp6OlP43av2phoeAOOuf1OPQCKpckkqblz1CT7x5wvBgnxCu7dRnJjTeV1/elnAH9Y4LngBQC9PX287IgH195uuWDNQNDaMiMJiwGPc9Rhg7l/4KWjUCFUiB9isbAvQEmP7Y6onZRwj5MRjtj6sy3rVEWAYqKn1htt+UD5Ybtprw/K8HGmXyvUceeUSYqqS/Q4cO1W0DeKmIl4w4ikFGUVUCLTh4n1i7NX1zs7d+Hdf74Egio7lj5NYkScxcOfGwsDj46w/y6O1pdOvt19Lwe1+jfR0WCmc0OWwy1jPVagkIdq7sWTmLll3Qk/q1+MpQewoqk4Vpl29iI6hZkYI5GBmB2+LrGNTabGTADrE5LLEANDSMisBg4q8hnx8E5OfSH0GrRaACyap+9dUOs8KJl+hc94ax7eXHYPS3MiMOGtpOb39WWt0C6Te1B8sP1xQFJUFrvYialJQUCgkJoePHj3u8z69btVJfSuD3fW0v/WtmnxEREWLtTf5X57CwuWo/HU1/THjtMwnh+XRb2xeE1UAyDRq+uTlXggKOSklLdD3JcSg0/6ucDAOp9yHVE9E1FcusSfLMwP4un/BT2JS/JdPgnjHuJbEOfW+kHd13Cf8AXrvlf3lN25egMWJG5QKe12W+aahduVWtxFq1iPrQOTThE6BiRQrmYGQEbot0DEpho2X5sspa0RCsHlaLwEASXuh9VxKYvxrwFVPft83jXPojaLUIVCBZ1a++2mG0jSxi+Li3Flxo6P5W6yejv5VdcYah7fT2Z6XVLRDBqfVgGVBR5oYkasLDw6lXr160evVq93vsKMyv+/ZVvzH5ffn2zFdffeXenqOdWLzIt2EfmR9//FFznw2G7E+o1dGZIhROa81Turmd2hWbiKLbuHx1FIRkf0Tv9Ojn4eMiF0yB1vuQ6okYvTF2/L5PRGJxnahAmPzXWkEjh8MF2bt+++lR4qlKz5nRiBlV+qzaadd8QuNBLKcyndLPHkz9rhhPGxPfFdmhfcGOhHyzf3dyeJ0NRkZJ6XadOIYSmzHLl1VP1Q3B6mGtCHT9q5wA+bUR0Svfh9rnksDcVnChKQshU3i6hde59EfQahGoQPKnX83s30gbpd+QH7fLn2++CArQ+l21fjLaH58eudkSYWml1U2J/rGw1TqRHti+UPPB0u+izHWAXyHdt9xyC7366qsiNw2HdC9ZsoR27twp/GA43DszM1P4vUgh3RdddBHNmjWLRo4cSYsWLaL//Oc/XiHd/Lk8pHvbtm0NM6RbkYTPWapeYNIrsuHsu2q2cxqLftLIWizddM8fepV6DxtvSXgd11nKOqjvSf/+wQm0qWCg32G8rL0m364uaLRypGzeqZ3AzWwiPzEZkc07aaFacikpiZeUNC082Z1ReO/JNHdyNXlCsk2/l9Pxnav9jj4KFLbasch1XxOOalEE9et1u+lwUarmedMLs+c+4n6zkUO1/IaRqBp3ErXOEV4JBbf8aSyXTDDhB4V7OkynVFmUSU5lBq0+caVXPpWCqhaurN2Ka0kOT1zf5HrnYtHKF6UVEaiMEuLJ+MPs2+iDQ/dq9rVWHiuj+ank+9FL9GkkOkmtLdw3w1suVg3DdtZMtkb37yuSkifn5/d4h4fzNT+uzVy6PvMtj4dSX/1ktD+s6Dcj96SZSDbj/Sbs8jTtd2MRWkYTyDbokG6Gw7ml5Hscmv3SSy+JaCVm8ODB1K5dOxGOLcF5bCZPnuxOvvf000+rJt977bXXRPK9AQMG0Msvv0ydO3c21J56ETUmQruP2AfQwjt/oJBf7/PMRMwWml6zvQWNTtZiYfWJbk22q/ZZk/beUU3OT1igZRsKsfY3pFAtV4sRlJkzs3Or6JNvS0wn8luafTtdlLLCY4AtD21NkX195woyw9pNhSLKKZDBiPPGDDpPXQCoJZnTKzFgJLmf3kC8OPtOGpP5mu5A7RYvJjMCK88xZ0v+6sdSKvGRHdUs3DZOqLb7UKWHSOYkfaMGx9FfhsW6s/B6ZKPNr6TTR9d55O7hBxWlCDpVmUQbCkdRXuxI+nhHD1FGQ15uo1e7Itq0P4725avn99ESAJ8dv5GK7B1EqYuuFw4VbVImRVSKQ2VqBuXvxcfYKTzM5pHcTeofecI5tUR+RZRJz/7+mOb9z/d5YlxIbRbeGKdq7qOEGAdNufA16uF41SNhHo+L1ee9QNsrRtYep+zad2WqLqPiMqd29uzTibQ0uzaRn4iokl1KfPyjhsTS+Mtc55zzeK3aFu9V/kQiMtz1/QvitQXjd3nD3Y+gZoXlsAujad3mMiqTLekbFUfctvJKMs1A1dxpbWhf+tN0+zv9DO8n0MzzDULUNDTqRdSYSMI36NrbXU/PRuuNGBRMdMkaopaDyQr2/biI2u0eZyh3jCt5ln42XfkAMnGc76y6ftWUeuNEzdPMUUMFPFlg8oAqJY8bObQT9eh/meX1kH779r90zuHxpp/UWMxwAsG/DDdeWsEMcuHAWYDX/Vrm4eStNsiVhWbSrzGudPxxeZ9Qz+JJFFldm16+LCSTfo2dRSWp11heu0mr2KiZjMIiC7JC9PlbYNDrex1Cxf18aP9ByqtMo5BWg6jbmdE+f0PvmOQZhfefTK0tt2CyjfKaPlp5rbRqD6kdp1zs8bi1bmuFd00jpaVQr/+kY/KjDpNX+YIYp3opBpP1vJTnR37tMKL9NeUdlOcov7jaI+monrBUWjsYzl/14TfF7nIWQlR24nQT3uJop+MK0d/ifKokaJUY3DOKNu0sV32gSWtB9MjwHZTV+pS777mspdH8aMyjtyXTJb3NP6xKQNTUQacYwqDw2HrGF5Q14HL3a0MDqkHBRP0WErW7kayA60ut+/R9mtBhGqVF6OdrYWGTq5KgSjlJj7k0zlT9IzPUJtW7U7w2k3nVX6uRUaoPfEjl3/+TYpxHdJ/UIiOIxgyND5qY0Wyjw0lb/ygXpQQcTqcQn5xRWD45eE0wFhYCBMFFtZiijvAwQ0OqVN+Q+5zvK606c2qJU7VEpVddMb3CmS1qz7UvsaZ2zrQy2asBS00QqE+fGuFzoZIBQW2JyFfFVo9Bph4sNVv+KKeJs0/QeQkb6DkDWTd9+YYEW8zI4ZuVkwNm7ptA4Q6ZCVvDOhIMq5EmNT4tazfsUS054DaB10E/geYJhEfD6HNeLgumwNT63ZAAxhWt+n0N3acmNOA9NFekJHzCmVcECco+dE2kwl9GJmjUlC9f5Py+h1pPHUjOqNZEZdk1bltKXIJJLWLKX6Tw7hahxtS5WhRPsJdP1ODfad9nLFHv61WrGbOFhq0jW8tH0q1X1LGAsIcIK13XfuaekgCwCr6+AvF3ANb0OY/tWst9DfVcD2ZLts0V5aRFING3wQKixookfJv+pXAAbu3hAMwKmlW6L+YtOyUuelHOfmsF/bxriqhV5PCKsKi5gGSCyQqk8O6PlpgLEeS1Ym7RTSPi6aYR9Whx4L7oNoXo3EfdBTwlU+21A8JpSj0KCEwsAIDGOA4M7hlD9jtsQbcyWQlETaCwcMm82qefgZEaTlJNDV57dVl0hlFhyStejpvloZmWRuvIEReocwTlbkunpDDfVbQlnxpeUnksyP4ppuB+bzlYJGDq1J6oU323BwAAGjGD6sDKZCUQNRZOpFoYreHEyxNvfFLgfs0+IFxywatq9nVt/MhMYAxOgFedPIdsG64XC1/qOXhcC2zvnZhOU+5o2SDVOgAAgOZnZYKoqQOM1mZifwulRYcdS5WOuKyYg3mBiSR0dpVltRoqw1rT0bbP0H03jmmwah0AAEDzA6KmDpCccH0tQfEaJTuQGsGo5ceyZTUpq25EqijiGZk6kNojjBcAAEADA6KmDpCccH3F/d81KlFk4QxmVW6rl9UAAACARlvQEvgP+51w2DZbbNR45cN8yi+qDmpVbgAAAKApA0tNHQsbztyqFvfPS1NrPn6f3uj6OMXSEc06Sw0xLwAAAADQEIClpg7hfDWvLPOuSisvXiZPq89wYUR+f+QZX3ql0wYAAABALRA1dYhWvhpeYuJ8NKK4gsIII3LF2IgmdplBg7Kw7AQAAABoAVFTh2hFLXEeGk6wp1WQkUsl2EoPuSKRAAAAAKAKRE0dohW1JK+f5BPOWAwAAAAAVSBq6iFfjRLOFGwILsEAAAAAAFUgauohX40Srp/EUU4Op/r6k1NU5W5jaVVuAAAAoKkBUVPHcPTSrVfEe5VC4LBtljRKYeN+bXFVbgAAAKCpAVFTD2Smhnm9x3lopu18hXIrW3m8n1vRiranvx+UqtwAAABAUwLJ9xqQw7BaVW5emnp2AHxpAAAAAD0gahpYgUtlVW6URQAAAACMgeWnBuQwrAbKIgAAAADGgKhpoAUu2UKDsggAAACAcbD8VM/Cpn9WlKt8wqnTVFDsoIQ4O6UmhoolKhSuBAAAAIwDUVPPsHDp0TmyvpsBAAAANHqw/AQAAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkAUQMAAACAJgFEDQAAAACaBBA1AAAAAGgSQNQAAAAAoEkAUQMAAACAJkGTyCjsdDrFv4WFhfXdFAAAAAAYRJq3pXk8UJqEqCkqKhL/tmnTpr6bAgAAAAA/5vGEhAQKFJvTKnlUjzgcDjpy5AjFxcWRzWazVEGyUDp06BDFx8dbtl+Afm+o4JpHvzc3cM3Xf9/v2LGDzjrrLLLbA/eIaRKWGu6I1q1bB23/LGggauoe9Hv9gb5Hvzc3cM3XH5mZmZYIGgaOwgAAAABoEkDUAAAAAKBJAFHjg4iICJo6dar4F9Qd6Pf6A32Pfm9u4JpvWn3fJByFAQAAAABgqQEAAABAkwCiBgAAAABNAogaAAAAADQJIGoAAAAA0CRo9qJm3rx51K5dO4qMjKQ+ffrQTz/95LPDli5dSl26dBHbd+vWjVasWFFnJ6u59vvrr79OAwcOpBYtWoi/oUOH6p4nYE3fy1m0aJHI2H3NNdege+ug3/Pz82nChAmUnp4uokM6d+6M8aaO+n727Nkiw21UVJTIeHv//fdTeXm5vz/fLFm3bh1deeWVlJGRIcaNjz/+WPc7a9eupZ49e4rrvVOnTrRgwQLzP+xsxixatMgZHh7ufOutt5y//fab84477nAmJiY6jx8/rrr9hg0bnCEhIc6nn37auWPHDufkyZOdYWFhzu3bt9d525tTv48bN845b94856+//ur8/fffnbfeeqszISHBefjw4Tpve3Pre4l9+/Y5MzMznQMHDnReffXVddbe5trvFRUVzvPPP985YsQI53fffSf6f+3atc4tW7bUedubW99/8MEHzoiICPEv9/uXX37pTE9Pd95///113vbGzIoVK5yPPvqoc/ny5Rxh7fzoo498br93715ndHS0c+LEiWJ+nTNnjphvV65caep3m7WoueCCC5wTJkxwv66urnZmZGQ4n3zySdXtb7jhBufIkSM93uvTp4/z73//e9Db2pz7Xcnp06edcXFxznfeeSeIrWya+NP33N/9+vVzvvHGG85bbrkFoqYO+v2VV15xdujQwVlZWenPz4EA+p63vfjiiz3e44m2f//+6Fc/MSJqHnroIee5557r8d6YMWOcw4YNM/VbzXb5qbKykjZt2iSWMiS49gS/3rhxo+p3+H359sywYcM0twfW9LuS0tJSqqqqoqSkJHRxHfT9448/TmlpaXT77bejv+uo3z/99FPq27evWH5q2bIlde3alf7zn/9QdXU1zkGQ+75fv37iO9IS1d69e8Wy34gRI9D3QcSq+bVJFLT0h9zcXDFA8IAhh1/v3LlT9TvHjh1T3Z7fB8HrdyUPP/ywWKdV3gDA+r7/7rvv6M0336QtW7age+uw33ki/eabb2j8+PFiQt29ezfdfffdQsxzBlYQvL4fN26c+N6AAQN4JYNOnz5N//jHP+j//b//h24PIlrzK1fyLisrE/5NRmi2lhrQOJk1a5ZwWP3oo4+E0x8IHkVFRXTTTTcJR+2UlBR0dR3icDiEdey1116jXr160ZgxY+jRRx+l+fPn4zwEGXZWZavYyy+/TJs3b6bly5fT559/TjNmzEDfNwKaraWGB+mQkBA6fvy4x/v8ulWrVqrf4ffNbA+s6XeJZ599Voiar7/+mrp3747uDXLf79mzh/bv3y8iGOSTLRMaGkq7du2ijh074jxY3O8MRzyFhYWJ70mcffbZ4mmWl1TCw8PR70Hq+8cee0yI+b/97W/iNUe5lpSU0J133imEJS9fAevRml/j4+MNW2mYZnt2eFDgJ6DVq1d7DNj8mtey1eD35dszX331leb2wJp+Z55++mnxpLRy5Uo6//zz0bV10PecumD79u1i6Un6u+qqq2jIkCHi/znUFVjf70z//v3FkpMkIpk//vhDiB0ImuBd85LPnlK4SOISpRKDh2Xzq7OZh/px6N6CBQtECNmdd94pQv2OHTsmPr/pppuckyZN8gjpDg0NdT777LMitHjq1KkI6a6Dfp81a5YIyVy2bJnz6NGj7r+ioqLAL4Jmhtm+V4Lop7rp94MHD4oIv3vuuce5a9cu52effeZMS0tzzpw5088WNF/M9j2P69z3//3vf0WY8apVq5wdO3YU0a/AODw+cxoO/mOp8fzzz4v/P3DggPic+5z7XhnS/e9//1vMr5zGAyHdfsCx8GeccYaYNDn074cffnB/dtFFF4lBXM6SJUucnTt3Fttz+Nnnn3/uz882e8z0e9u2bcVNofzjwQcE/5qXA1FTd/3+/fffi5QRPCFzePcTTzwhwutBcPu+qqrKOW3aNCFkIiMjnW3atHHefffdzlOnTqHrTbBmzRrVcVvqa/6X+175nR49eojzxNf822+/7TSLjf9jrREJAAAAAKDuabY+NQAAAABoWkDUAAAAAKBJAFEDAAAAgCYBRA0AAAAAmgQQNQAAAABoEkDUAAAAAKBJAFEDAAAAgCYBRA0AAAAANFm3bp2oAZeRkUE2m40+/vhjMsO0adPE95R/MTExZDUQNQAAAADQhAt6ZmVl0bx588gfHnzwQTp69KjH3znnnEPXX389WQ1EDQAAAAA0GT58OM2cOZOuvfZa1c8rKiqEcMnMzBTWlz59+tDatWvdn8fGxooq3NIfV9/esWMH3X777WQ1EDUAAAAA8Jt77rmHNm7cSIsWLaJt27YJC8zll19Of/75p+r2b7zxBnXu3JkGDhxIVgNRAwAAAAC/OHjwIL399tu0dOlSIVI6duworDYDBgwQ7yspLy+nDz74IChWGiY0KHsFAAAAQJNn+/btVF1dLSwvyiWp5ORkr+0/+ugjKioqoltuuSUo7YGoAQAAAIBfFBcXU0hICG3atEn8K4d9adSWnq644gpq2bIlBQOIGgAAAAD4xXnnnScsNSdOnND1kdm3bx+tWbOGPv30UwoWEDUAAAAA8GmN2b17t4c42bJlCyUlJYllp/Hjx9PNN99Mzz33nBA5OTk5tHr1aurevTuNHDnS/b233nqL0tPTRTRVsLA5nU5n0PYOAAAAgEbN2rVraciQIV7vs1/MggULqKqqSoR8v/vuu5SdnU0pKSl04YUX0vTp06lbt25iW4fDQW3bthXi54knnghaWyFqAAAAANAkQEg3AAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAIAmAUQNAAAAAJoEEDUAAAAAaBJA1AAAAACgSQBRAwAAAABqCvx/bxWnfaf8gX0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "band = 'band6'\n",
    "\n",
    "offset = 0\n",
    "\n",
    "for fold_id in range(n_folds):\n",
    "    \n",
    "    plt.scatter( train_obs[fold_id][band]['q'][0], train_obs[fold_id][band]['V'][0] + offset, color='royalblue'  )\n",
    "    plt.scatter( test_obs[fold_id][band]['q'][0], test_obs[fold_id][band]['V'][0]+ offset, color='orange'  )\n",
    "    \n",
    "    offset += 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a357187-b450-46fa-aeaa-94948e762681",
   "metadata": {},
   "source": [
    "We construct the model in the same way as before. However, we make a function this time as we are going to iteretively run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "16922921-423b-4fa6-98b5-1c39fe819f8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dust_opacity = jnp.load(dsharp_opac.get_datafile('default_opacities_smooth.npz'))\n",
    "\n",
    "\n",
    "\n",
    "def set_model( obs, lengthscale ):\n",
    "    \n",
    "    disk = frap.model( incl= 6.28,\n",
    "                         r_out = 1.3,\n",
    "                         N_GP = 260, # smaller\n",
    "                         flux_uncert = True )\n",
    "    \n",
    "    \n",
    "    D = 114.9 # pc\n",
    "    \n",
    "    disk.set_parameter('log10_T',\n",
    "                       free = False,\n",
    "                       profile = lambda r: jnp.log10(77.0*(r*D/10)**(-0.5)) )\n",
    "    \n",
    "    \n",
    "    disk.set_parameter('q',\n",
    "                       free = False,\n",
    "                       profile = lambda r: 3.5 )\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    disk.set_parameter('log10_Sigma_d', free = True,\n",
    "                                        bounds = ( -6.0, 3.0 ),\n",
    "                                        lengthscale = lengthscale)\n",
    "    \n",
    "    disk.set_parameter('log10_a_max', free = True,\n",
    "                                        bounds =  ( -2.0, 2.0 ),\n",
    "                                        lengthscale = lengthscale)\n",
    "    \n",
    "    \n",
    "    f_s = { 'band3': 0.025, 'band6': 0.05 }\n",
    "    f_mean = { 'band3': 1.0, 'band6': 1.0 }\n",
    "    \n",
    "    for band in bands:\n",
    "        disk.set_visibility( band = band,\n",
    "                                q = obs[band]['q'],\n",
    "                                V = obs[band]['V'],\n",
    "                                s = obs[band]['s'],\n",
    "                                nu = obs[band]['nu'],\n",
    "                                Nch = obs[band]['Nch'],\n",
    "                                f_s = f_s[band],\n",
    "                                f_mean = f_mean[band] )\n",
    "    \n",
    "    \n",
    "    disk.set_opacity( opac_dict = dust_opacity )\n",
    "\n",
    "    return disk\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b802787e-32b1-4104-aadb-de85a9d50c0f",
   "metadata": {},
   "source": [
    "This is the core part of this experiment. We vary the length scale parameter from 0.01 to 0.2 arcsec with 6 steps. We fit the model to the training datasets using the Stochastic Variational Inference (SVI). Note that this aims to just get the MAP estimate rather than the posterior probability distribution. The best-fit model (contained in `map_estimate`) is then passed to ``calc_chi2()`` and calcurated the chi-squared between the test data and trained model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b69a9c5e-ba7f-4c1a-b32a-bbf845e4d6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be332c80-f0d3-4785-88ee-042f80cbc8cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████| 10000/10000 [00:02<00:00, 3636.10it/s, init loss: 2267982.4833, avg. loss [9501-10000]: -946.7230]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3575.05it/s, init loss: 2780281.4840, avg. loss [9501-10000]: -1501.2962]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3589.67it/s, init loss: 3204897.3779, avg. loss [9501-10000]: -421.0165]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3449.84it/s, init loss: 3134429.5749, avg. loss [9501-10000]: -227.6692]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3595.64it/s, init loss: 3166669.0292, avg. loss [9501-10000]: -155.5769]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3476.86it/s, init loss: 2303864.9237, avg. loss [9501-10000]: 2584.8933]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3650.19it/s, init loss: 2825341.3065, avg. loss [9501-10000]: 1046.0936]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3640.46it/s, init loss: 3257051.8628, avg. loss [9501-10000]: 4434.7356]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3516.19it/s, init loss: 3189747.6306, avg. loss [9501-10000]: 3704.0627]\n",
      "100%|███████████████| 10000/10000 [00:03<00:00, 2979.97it/s, init loss: 3211385.7718, avg. loss [9501-10000]: 2661.0913]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3701.40it/s, init loss: 2374480.8694, avg. loss [9501-10000]: 17719.4199]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3535.45it/s, init loss: 2905487.1398, avg. loss [9501-10000]: 9934.0114]\n",
      "100%|███████████████| 10000/10000 [00:02<00:00, 3426.57it/s, init loss: 3353864.3141, avg. loss [9501-10000]: 8755.2184]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3691.94it/s, init loss: 3289188.5347, avg. loss [9501-10000]: 21970.7591]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3588.52it/s, init loss: 3290168.2765, avg. loss [9501-10000]: 19431.2547]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3674.87it/s, init loss: 2446200.3768, avg. loss [9501-10000]: 17891.6687]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3651.14it/s, init loss: 3003678.9450, avg. loss [9501-10000]: 24149.0978]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3551.20it/s, init loss: 3454946.6006, avg. loss [9501-10000]: 19812.0084]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3607.67it/s, init loss: 3391084.1780, avg. loss [9501-10000]: 18368.9768]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3695.26it/s, init loss: 3355895.4526, avg. loss [9501-10000]: 18125.9190]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3615.70it/s, init loss: 2507411.4472, avg. loss [9501-10000]: 41976.0602]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3641.97it/s, init loss: 3035349.6109, avg. loss [9501-10000]: 5057216.7438]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3623.31it/s, init loss: 3482559.4627, avg. loss [9501-10000]: 7948036.6270]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3573.53it/s, init loss: 3415938.8750, avg. loss [9501-10000]: 1820248.4722]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3646.44it/s, init loss: 3358206.0644, avg. loss [9501-10000]: 33879.2780]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3675.28it/s, init loss: 2482746.2140, avg. loss [9501-10000]: 2404236.5582]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3634.88it/s, init loss: 2951094.8391, avg. loss [9501-10000]: 5066346.5318]\n",
      "100%|████████████| 10000/10000 [00:02<00:00, 3623.11it/s, init loss: 3417574.6154, avg. loss [9501-10000]: 2789714.0355]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3615.19it/s, init loss: 3361061.3063, avg. loss [9501-10000]: 81569.4412]\n",
      "100%|██████████████| 10000/10000 [00:02<00:00, 3648.19it/s, init loss: 3364994.4475, avg. loss [9501-10000]: 74550.8066]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "lengthscales = np.logspace( np.log10(0.01), np.log10(0.2), 6 )\n",
    "chi2_arr_dict = {}\n",
    "\n",
    "\n",
    "for il, lengthscale in enumerate(lengthscales):\n",
    "\n",
    "    chi2_arr_dict[il] = np.array([])\n",
    "\n",
    "\n",
    "    for fold_id in range(n_folds):\n",
    "    \n",
    "    \n",
    "        disk_train = set_model( train_obs[fold_id], lengthscale )\n",
    "\n",
    "        disk_test = set_model( test_obs[fold_id], lengthscale )\n",
    "\n",
    "        inference = frap.inference( disk_train )\n",
    "\n",
    "        map_estimate = inference.SVI_MAP( num_iterations=10000, adam_lr=0.1, seed = 1)\n",
    "        \n",
    "\n",
    "        chi2 = disk_test.calc_chi2( map_estimate )\n",
    "              \n",
    "        chi2_arr_dict[il] = np.append( chi2_arr_dict[il], chi2 )\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af1caaec-416c-4e4b-8fe0-e53ddbf63caf",
   "metadata": {},
   "source": [
    "Here we plot the results. We can see that the chi-squared satuates at a length scale of ~0.03 arcsec, although there are some outliers. In principle, any length scales smaller than this value can reasonably fit the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3db6d229-dd92-4e9c-8f17-bbec7a093bde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGnCAYAAABsCPeYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJtxJREFUeJzt3X9w1PW97/H3biihIT9IDEkgBrFO6iFQQSFJ4WprlAJ6xKqnyj+3DXivds5knHrSzhT8A2TGgp12OGlLepA718FTO1OoDpQz46RW1KKFkgClVREbvIzDAPklJEvCMSnJ3nl/QiK7SRaS/e7u5/v9Ph8zO+n38/24+80m3bz4/Hh/A+FwOCwAAACWCKb6AgAAAK5GOAEAAFYhnAAAAKsQTgAAgFUIJwAAwCqEEwAAYBXCCQAAsArhBAAAWIVwAgAArEI4AQAAViGcAAAAf4eTjz76SBYsWDD8+OIXvyh79uxJ9mUAAABLBVJ547/u7m6ZPXu2fPLJJzJ16tRUXQYAALDIpFS++N69e+Xee+8dVzAZGBiQs2fPSlZWlgQCgYReHwAAcIaOhVy8eFFmzpwpweA1Jm7C4/THP/4x/MADD4RnzJihIy7h3bt3j+izdevW8E033RROT08PV1RUhA8dOjTqc33zm98Mv/rqq+N6/dOnT5vX5cF7wO8AvwP8DvA7wO+AuO490L/j1zLukZOenh6ZP3++PP744/LII4+MOL9z506pra2Vbdu2SWVlpdTV1cny5cvNWpOCgoLhfqFQSA4cOCC/+c1vYr5eb2+veVydvNTp06clOzt7vJcPAABSQP/ul5SUmJmPhK450WmV3bt3y0MPPTTcpoGkvLxctm7dOjwNoxfz1FNPydq1a4f7/epXv5Lf//738vLLL8d8jWeffVY2btw4or2rq4twAgCAi8JJTk7Odf39dnS3Tl9fnxw5ckSWLl36+QsEg+b44MGDEX137dolq1atuuZzrlu3znwjQw8dMQEAAN7l6ILYjo4O6e/vl8LCwoh2PT5x4sTwsYaMxsZGefXVV6/5nOnp6eYBAAD8ISW7dXRYp7W1NRUvDQAALOfotE5+fr6kpaWNCB56XFRU5ORLAQAAj3I0nEyePFkWLlwo+/btG27TBbF6vHjx4rieu76+XsrKysxiWwAA4F2TJlLV9eTJk8PHp06dkmPHjkleXp7MmjXLbCOurq6WRYsWSUVFhdlKrNuP16xZE9eF1tTUmMfQal8AAOBN4w4nhw8flqqqquFjDSNKA8mOHTvMDpz29nZZv369tLS0mPvnNDQ0jFgkCwAAYN29dRK9TxoAAPi8zgkAAICrb/w33gWx+tA6KgAA/xgYCMu55k7pCfXK1Ox0mVE6TYJBbvyaCAMD/XLmww+ku/OCZE7LleI5cyUYTJNkY1oHAGCtj//SJu/sbJaezs/vsTZ1WrrctapUbrn98/u1IX7Nhw7Imzu2S/f5juG2zLx8uWf1k1JauSTu52daBwDgiWDS8ML7EcFE6bG263k4F0z2btkUEUyUHmu7nk8m1pwAAKycytERk1je3dVs+iH+qRwdMYnlrZe2m37JQjgBAFjHrDGJGjGJ1n2h1/RDfMwak6gRk2gXP+0w/ZLFNeGECrEA4B+6+NXJfhibLn51sp+vwolWhz1+/Lg0NTWl+lIAAAmmu3Kc7Iex6a4cJ/v5KpwAAPxDtwvrrpxYMnMHtxUjPrpdWHflxJJ1Q77plyyEEwCAdbSOiW4XjuXOx0qpd+IArWOi24Vjqap+Mqn1TggnAAAraR2TFd+dN2IERUdMtJ06J87ROiYP1j4zYgRFR0y03Yk6J+NBETYAgNWoEOuNCrHjKcLmmvL1AAD/TvEU35q8xZh+FgymScnc21J9Ge6Z1mErMQAA/sC0DgAASDjurQMAAFzLNdM6AADAHwgnAADAKoQTAABgFcIJAACwCuEEAABYxTXhhDonAAD4A3VOAABAwlHnBAAAuJZrpnUAAIA/EE4AAIBVCCcAAMAqhBMAAGAVwgkAALAK4QQAAFiFcAIAAKzimnBChVgAAPyBCrEAACDhqBALAABcyzXTOgAAwB8IJwAAwCqEEwAAYBXCCQAAsArhBAAAWIVwAgAArEI4AQAAViGcAAAAqxBOAACAVVwTTri3DgAA/sC9dQAAQMJxbx0AAOBarpnWAQAA/kA4AQAAViGcAAAAqxBOAACAVQgnAADAKoQTAABgFcIJAACwCuEEAABYhXACAACsQjgBAABWIZwAAACrEE4AAIBVCCcAAMAqk1J9AQAAwA4DA/1y5sMPpLvzgmROy5XiOXMlGExL+nUQTgAAgDQfOiBv7tgu3ec7ht+NzLx8uWf1k1JauSSp75BrpnXq6+ulrKxMysvLU30pAAB4Lpjs3bIpIpgoPdZ2PZ9MgXA4HBYXCYVCkpOTI11dXZKdnZ3qywEAwPVTOf+n5n+NCCZXy7ohX/731v8b1xTPeP5+u2bkBAAAOM+sMYkRTNTFTztMv2QhnAAA4GPdnRcc7ecEwgkAAD6WOS3X0X5OIJwAAOBjxXPmml05seiaE+2XLIQTAAB8LBhMM9uFY6mqfjKp9U4IJwAA+Fxp5RJ5sPaZESMoOmKi7cmuc0IRNgAAIBpAbimvpEIsAACwRzCYJiVzb0v1ZTCtAwAA7MKaEwAAYBXWnAAArDYwEJZzzZ3SE+qVqdnpMqN0mgSDgVRfFhKIcAIAsNbHf2mTd3Y2S09n73Db1GnpcteqUrnl9oKUXhsSh2kdAIC1waThhfcjgonSY23X8/AmwgkAwMqpHB0xieXdXc2mH7yHcAIAsI5ZYxI1YhKt+0Kv6QfvIZwAAKyji1+d7Ad3IZwAAKyju3Kc7Ad3IZwAAKyj24V1V04smbmD24rhPYQTAIB1tI6JbheO5c7HSql34lGEEwCAlbSOyYrvzhsxgqIjJtpOnRPvoggbAMBaGkBunj+dCrE+QzgBAFg/xVN8a26qLwNen9Y5deqUVFVVSVlZmXzlK1+Rnp6eVFwGAACwUEpGTlavXi3PPfec3HXXXXL+/HlJT2crGAAASFE4+eCDD+QLX/iCCSYqLy8v2ZcAAAC8NK2zf/9+WblypcycOVMCgYDs2bNnRJ/6+nqZPXu2TJkyRSorK6WxsXH4XHNzs2RmZprnuOOOO2TTpk3xfxcAAMC/4UTXh8yfP98EkNHs3LlTamtrZcOGDXL06FHTd/ny5dLWNnj3yMuXL8s777wjv/zlL+XgwYPyhz/8wTzG0tvbK6FQKOIBAAC8a9zh5L777jPrRR5++OFRz2/ZskWeeOIJWbNmjVnwum3bNsnIyJAXX3zRnC8uLpZFixZJSUmJWWty//33y7Fjx8Z8vc2bN0tOTs7wQ/87AADgXY7u1unr65MjR47I0qVLP3+BYNAc6yiJKi8vN6MoFy5ckIGBATNNNGfOnDGfc926ddLV1TX8OH36tJOXDAAAvLwgtqOjQ/r7+6WwsDCiXY9PnDgx+IKTJpl1Jl/72tckHA7LsmXL5IEHHhjzOXV0hd08AAD4R0q2EuvUkD4AAAASOq2Tn58vaWlp0traGtGux0VFRU6+FAAA8ChHw8nkyZNl4cKFsm/fvuE2XVeix4sXL47ruXV3kC6w1TUrAADAu8Y9rdPd3S0nT56MKEWvu220mNqsWbPMNuLq6mqzI6eiokLq6urM9mPdvROPmpoa89CtxLprBwAAeNO4w8nhw4fNfXGGaBhRGkh27Nghq1atkvb2dlm/fr20tLTIggULpKGhYcQiWQAAgNEEwrplxkWGRk50W3F2dnaqLwcAADj89zsldyUGAABwfThhQSwAAP7AtA4AAEg4pnUAAIBruWZaBwAA+APhBAAAWIVwAgAArOKacMJuHQAA/IHdOgAAIOHYrQMAAFzLNdM6AADAHwgnAADAKoQTAABgFdeEE3brAADgD+zWAQAACcduHQAA4FqumdYBAAD+QDgBAABWIZwAAACrEE4AAIBVCCcAAMAqrgkn1DkBAMAfqHMCAAASjjonAADAtVwzrQMAAPyBcAIAAKxCOAEAAFYhnAAAAKsQTgAAgFUIJwAAwCqEEwAAYBXXhBMqxAIA4A9UiAUAAAlHhVgAAOBarpnWAQAA/kA4AQAAViGcAAAAqxBOAACAVQgnAADAKpNSfQEAAMAOAwP9cubDD6S784JkTsuV4jlzJRhMS/p1EE4AAIA0Hzogb+7YLt3nO4bfjcy8fLln9ZNSWrkkqe8Q0zoAAPhc86EDsnfLpohgovRY2/V8MhFOAADw+VTOmzu2x+zz1kvbTb9kcU044d46AAA4z6wxiRoxiXbx0w7TL1lcE05qamrk+PHj0tTUlOpLAQDAM7o7Lzjaz1fhBAAAOE935TjZzwmEEwAAfKx4zlyzKyeWrBvyTb9kIZwAAOBjwWCa2S4cS1X1k0mtd0I4AQDA50orl8iDtc+MGEHRERNtT3adE4qwAQAA0QByS3klFWIBAIA9gsE0KZl7W6ovg2kdAABgF9acAAAAqxBOAACAVQgnAADAKoQTAABgFcIJAACwCuEEAABYhXACAACsQjgBAABWIZwAAACruCac1NfXS1lZmZSXl6f6UgAAQAIFwuFwWFwkFApJTk6OdHV1SXZ2dqovBwAAOPz32zUjJwAAwB8IJwAAwCqEEwAAYBXCCQAAsArhBAAAWIVwAgAArEI4AQAAViGcAAAAqxBOAACAVQgnAADAKoQTAABgFcIJAACwCuEEAABYhXACAACsQjgBAABWIZwAAACrEE4AAIBVCCcAAMAqk1J9AQAAxDIwEJZzzZ3SE+qVqdnpMqN0mgSDAd40DyOcAACs9fFf2uSdnc3S09k73DZ1WrrctapUbrm9IKXXhsRhWgcAYG0waXjh/YhgovRY2/U8vIlwAgCwcipHR0xieXdXs+kH70nJtM7s2bMlOztbgsGg5ObmyltvvZWKywAAWMqsMYkaMYnWfaHX9Cu+NTdp1wWPrzk5cOCAZGZmpurlASBu4f5+uXT4iFxub5dJ06dLxqKFEkhL4511gC5+dbIf3IUFsQAwAaHXX5fWTZvlckvL5x+oRUVS+Mw6yV62jPc0Trorx8l+8Piak/3798vKlStl5syZEggEZM+ePSP61NfXm6mbKVOmSGVlpTQ2Nkac1//u61//upSXl8uvf/3r+L4DAEhBMDnzvacjgom63Npq2vU84qPbhXVXTiyZuYPbiuE94w4nPT09Mn/+fBNARrNz506pra2VDRs2yNGjR03f5cuXS1vb56uq3333XTly5Ijs3btXNm3aJH/729/i+y4AIIlTOTpiIuFRFmJeadPz2g8Tp3VMdLtwLHc+Vkq9E48adzi577775LnnnpOHH3541PNbtmyRJ554QtasWSNlZWWybds2ycjIkBdffHG4T3Fxsfk6Y8YMuf/++02IGUtvb6+EQqGIBwCkilljEjViEiEcNue1H+KjdUxWfHfeiBEUHTHRduqceJeja076+vrMiMi6deuG23RHztKlS+XgwYPDIy8DAwOSlZUl3d3d8uabb8pjjz025nNu3rxZNm7c6ORlAsCE6eJXJ/shNg0gN8+fToVYn3E0nHR0dEh/f78UFhZGtOvxiRMnzP9ubW0dHnXRvjrKomtPxqJBR6eJhujISUlJiZOXDQDXTXflONkP1zfFw3Zhf0n6bp0vfelL8te//vW6+6enp5sHANhAtwvrrhxd/DrqupNAQCYVFpp+ACyoEJufny9paWlmdORqelxUVOTkSwFASmgdE90uPHgQdfO5K8d6nnongCXhZPLkybJw4ULZt2/fcJuuL9HjxYsXx/XcujtIF9jGmgICgGTQOibFP6szIyRX02Ntp84JEJ9AODzauOTYdBHryZMnzf++/fbbze6cqqoqycvLk1mzZpmtxNXV1fLCCy9IRUWF1NXVya5du8yak+i1KBOha05ycnKkq6vLlMAHgFShQiyQmL/f415zcvjwYRNGhgwtVtVAsmPHDlm1apW0t7fL+vXrpaWlRRYsWCANDQ2OBBMAsIlO3UytrEj1ZQCeM+6Rk1Rj5AQAAPcZz99vR9ecAAAAxMs14YQFsQAA+APTOgAAIOGY1gEAAK7lmmkdAADgD4QTAABgFcIJAACwimvCCbt1AADwB3brAACAhGO3DgAAcC3XTOsAAAB/IJwAAACrEE4AAIBVCCcAAMAqrgknbCUGAMAf2EoMAAASjq3EAADAtVwzrQMAAPyBcAIAAKxCOAEAAFYhnAAAAKsQTgAAgFVcE06ocwIAgD9Q5wQAACQcdU4AAIBruWZaBwAA+MOkVF8AAACww8BAv5z58APp7rwgmdNypXjOXAkG05J+HYQTAAAgzYcOyJs7tkv3+Y7hdyMzL1/uWf2klFYuSeo7xLQOAAA+13zogOzdsikimCg91nY9n0yEEwAAfD6V8+aO7TH7vPXSdtMvWQgnAAD42BldYxI1YhLt4qcdpl+yuCacUIQNAADn6eJXJ/v5KpzU1NTI8ePHpampKdWXAgCAZ2ROy3W0n6/CCQAAcJ5uF9ZdObFk3ZBv+iUL4QQAAB8LBtPMduFYqqqfTGq9E8IJAAA+V1q5RB6sfWbECIqOmGh7suucUIQNAACIBpBbyiupEAsAAOwRDKZJydzbUn0ZTOsAAAC7sOYEAABYhXACAACsQjgBAABWIZwAAACrEE4AAIBVJrnpxn/66O9P3i2bAbfqH+iXo21Hpf1Su0zPmC53FNwhaUms7ggA8QiEw+GwuEgoFJKcnBzp6uqS7OzsVF8OYJ03PnlDNh/aLG3/3TbcVvDFAllXuU6W3rQ0pdcGwL9C4/j7zbQO4LFg8m9v/1tEMFF6rO16HnCbgYGwnPnogvy9qcV81WN4m2umdQBceyrn2QPPxuyz8cBGqSqpYooHrvHxX9rknZ3N0tPZO9w2dVq63LWqVG65vSCl14bEYeQE8IjDrYelq68rZp/Ovk7TD3BLMGl44f2IYKL0WNv1PLyJcAJ4ROO5Rkf7AamkUzc6YhLLu7uameLxKKZ1AK8IONwP1xTu75dLh4/I5fZ2mTR9umQsWiiBNHZFOeFcc+eIEZNo3Rd6Tb/iW3P5bfUYwgngEeWF5bJdtl9XP8Qv9Prr0rpps1xuaRlum1RUJIXPrJPsZct4i+PUE+p1tB/chWkdwCPKi8olJz0nZp9p6dNMP8QfTM587+mIYKIut7aadj2P+EzNTne0H9yFcAJ4hBZZe3Zx7N06GxZvYKeOA1M5OmIio5WIutKm57UfJm5G6TSzKyeWzNx00w/eQzgBPESLrP373f8uhRmFEe16rO0UYYufWWMSNWISIRw257UfJi4YDJjtwrHc+Vip6QfnDAz0y+kP/iYf/umP5qsepwJrTgCP0QCitUwoX58YuvjVyX4Ym9YxWfHdeSPqnOiIiQYT6pw4q/nQAXlzx3bpPt8x3JaZly/3rH5SSiuXSDIRTgCPTvGwtiQxdFeOk/0QmwaQm+dPH9y9E+o1a0x0KocRE+eDyd4tm0a0a1DR9gdrn0lqQCGcAMA46HZh3ZWji19HXXcSCMikwkLTD87QIMJ24cTRqRsdMYnlrZe2yy3llRJM0g1EWXMCAOOgdUx0u/DgQdR6hyvHep56J3CLMx9+EDGVM5qLn3aYfslCOAGAcdI6JsU/qzMjJFfTY22nzgncpLvzgqP9nMC0DgBMgAaQrHvvpUIsXC9zWq6j/XwVTurr682jn9oBwHXdoZjdOomnUzdTKyv4jYSrFc+Za3blxJrayboh3/RLlkA4PNqKLnuFQiHJycmRrq4uyc7OTvXlANZ545M35PnG56X1UmtEnZO1FWupcwJgXLt1hjixW2c8f79ZcwJ4LJjUvl0bEUxU26U2067nATfeofjMRxfk700t5qsew1kaPDSA6AhK9IhJsrcRK0ZOAA9N5Sx/dfmIYDIkIAEzgtLwLw2UsIdrfPyXthFF2LSsvVaPpQhbYrYVm907nRfMGhOdynFq+zAjJ4AP6RqTsYKJCktYWi61mH6AW4JJwwvvRwQTpcfarufhLA0iJXNvkzn/4+vma7Lqmoy4jpS8KgDHtV9qd7QfkEo6daMjJrG8u6uZKR6PIpwAHjE9Y7qj/YBUMuXqo0ZMonVf6DX94D2EE8Aj7ii4w6wp0bUlo9H2oowi0w+wnd5Hx8l+cBfCCeChm/3pdmEVHVCGjn9Y8UMWw8IV9AZ/TvaDuxBOAA9ZetNS2XL3FinIKIho1xEVbdfzgBvonYd1V04smbmDdyiG97imQiyA66MBpKqkigqxcP2diHW7sO7KGcudj5WafvAe6pwAAFxV50RHTDSYUOfEXcZT54SREwCAtTSA3Dx/+uDunVCvWWOiUzmMmHgb4QQAYDUNIsW3Ju+OuEg9FsQCAACrEE4AAIBVCCcAAMAqhBMAAGAVwgkAALAKu3UAYILC/f1y6fARudzeLpOmT5eMRQslkJaaW8wDXkI4AYAJCL3+urRu2iyXW1o+/0AtKpLCZ9ZJ9rJlvKdAHJjWuaJ/ICwHP/5UfnfsjPmqxwAwVjA5872nI4KJutzaatr1PICJY+RERBrePycb/+u4nOv6bPiNmZEzRTasLJMV82bE8fYC8OJUjo6YSHiUf8BoWyBgzmfdey9TPMAE+X7kRIPJv758NCKYqJauz0y7ngeAIWaNSdSISYRw2JzXfgAmxtfhRKdudMRktAmcoTY9zxQP3KZ/oF+aWprktf/3mvmqx3CGLn51sh8Ai6Z1Ll26JHPmzJFHH31UfvrTn6bkGhpPnR8xYhIdUPS89lt8yw1JvTZgot745A15vvF5ab3UOtxWmFEoayvWytKblvLGxkl35TjZD4BFIyc/+tGP5Ktf/aqkUtvFzxztB9gQTGrfro0IJqrtUptp1/OIj24X1l05urZkVIGAOa/9ALgonDQ3N8uJEyfkvvvuk1QqyJriaD8glXTqRkdMwqNMVA61/bjxx0zxxEnrmOh24cGDqIBy5VjPU+8ESGI42b9/v6xcuVJmzpwpgUBA9uzZM6JPfX29zJ49W6ZMmSKVlZXS2NgYcf4HP/iBbN68WVKt4uY8sytnjH//mHY9r/3gEF37cOodkfdeGfzKWgjHHG07OmLEJDqgtFxqMf0QH61jUvyzOplUWBjRrsfaTp0TIMlrTnp6emT+/Pny+OOPyyOPPDLi/M6dO6W2tla2bdtmgkldXZ0sX75cPvroIykoKJDf/e538uUvf9k8Dhw4cM3X6+3tNY8hoVBInJIWDJjtwrorR4PI1f/eHAosel77wQHH94o0/FAkdPbztuyZIit+LFL2IG9xnNovtTvaD7FpANHtwlSIBZwXCIdH26x/nf9xICC7d++Whx56aLhNA0l5ebls3brVHA8MDEhJSYk89dRTsnbtWlm3bp28/PLLkpaWJt3d3fKPf/xDvv/978v69etHfY1nn31WNm7cOKK9q6tLsrOzxQnUOUlSMNn1nagIqK4Ev8f+k4ASJ92V8/jvH79mvxeXvyjlReXxvhwAjIsOLuTk5FzX329Hw0lfX59kZGTIK6+8EhFYqqurpbOz04yaXG3Hjh3y/vvvx9ytM9rIiYYdJ8OJ0u3CuitHF7/qGhOdymHExCE6dVM3L3LEJEJgcATl6fdEgtyXJJ41J8tfXW4Wv4627iQgAbNrp+FfGiSN9xmAxeHE0QWxHR0d0t/fL4VR87B63BKraFEM6enp5pu4+pEIGkR0u/A3FxSbrwQTB31yIEYwUWGR0JnBfpgwDRy6XXgoiFxt6PiHFT8kmACwXkrL169evTqVL49k6W51th/GpHVMtty9ZdQ6JxpMqHMCwHfhJD8/36wlaW2N/COjx0VaFwD+lFnobD/EpAGkqqTK7MrRxa/TM6bLHQV3MGICwDUcndaZPHmyLFy4UPbt2zfcpgti9Xjx4sVxPbduTy4rKzOLbeEyNy0ZXFMSa9N2dvFgPwCA74175ER32Jw8eXL4+NSpU3Ls2DHJy8uTWbNmmW3EugB20aJFUlFRYbYS6/bjNWvWxPVm19TUmMfQghq4iC6+1O3CZrfOGJu2VzzPYliHUL4egNuNe7fO22+/LVVVVSPaNZDo7hul24h/8pOfmEWwCxYskJ///Odmi3GyV/vCDXVOigeDCXVOHC1fH71bZ2hBrK5HYd2Jc8L9/dQ5AWzbSpwKhBMPbCvWXTm6+FXXmOhUDttaHd1KPFaVWLYSOyv0+uvSummzXL5qJ6LeU0dL11MhFrBoKzFwTRpEbr5L5CvfGvxKMHEM5euTG0zOfO/piGCiLre2mnY9D2DiXBNOWBALxEb5+uRN5eiIiYw26HylTc9rPwAeDye6GPb48ePS1NSU6ksBrKRbhp3sh9GZe+nEKioZDpvz2g+Ax8MJgNi0lokWW4uuDjtE24syikw/TNzl9nZH+wEYiXACeATl65Nj0vTpjvYDMBLhBPBg+fqCjIKIdh1RYRuxMzIWLTS7ciQwRlHBQMCc134AXHhvnfEuiNWH3lgQwNgoX59YgbQ0s11Yd+WYgHL1wtgrgUXPaz8AE0OdEyQXdU7gpTonP9pktg8Poc4J4EydE9eMnMCrFWJnDpa2p0IsXCh6M7HLaloC1mLNCZIXTPTeOlcHExU6N9iu5wGXFWHrj7oDe39bG0XYAAcQTpCcqRwdMRnx70z5vK1h7WA/wHIUYQMSj3CCxNN76USPmEQIi4TODPYDLEcRNiDxCCdIPL3Jn5P9gBSiCBuQeK4JJ9xbx8X07sNO9gNSiCJsQOK5Jpxwbx0Xu2nJ4K6cMcqqm/bs4sF+gOUowgYknmvCCVwsmDa4XTiWFc8P9gNcUoRt8CAqcFOEDXAE4QTJoXVMljwlEoj6ldNjbafOCVwke9kyKf5ZnUwqjJyK1GNt1/MAJo4ibEgOrWNy4BcjtxNr0Sptv7GcgAJX0QCSde+9g7t32tvNWhSd8qFsPRA/wgksqHMSGKxz8k//zNQOXEWDyNTKilRfBuA5TOsg8ahzAgAYB8IJEo86JwAAL4YT6py4GHVOAADjEAi77Daa47nlMixac1I3b/Amf6OuO9E6JzNFnn6PNScA4FHj+fvtmpETeKXOSXQhtivH1DkBAFxBOEFyaB2Tx/5TJHtGZLuOmGg7dU4AAFewlRjJowFEtwvr7h1dJKtrUbRkPZVhAQBXYeQEAABYhZGTK/oHwtJ46ry0XfxMCrKmSMXNeZIWHOtGdZhwlVgtxhY6Gzmto+tRmNYBAFxBOBGRhvfPycb/Oi7nuj4bel9kRs4U2bCyTFbMi1ojgYkHk13fGblbR3fwaDvrTgAAV/h+WkeDyb++fDQimCg91nY9j0SXr9eEuHawHwDA94J+n8rREZOxCr1ou57XfogD5esBAF4MJ4moEKtrTKJHTKLpee2HOFC+HgDgxXBSU1Mjx48fl6amJsees6Xrvx3thzFQvh4A4MVwkgjne/oc7YcxaC0T3ZUzojrs1eXriwf7AQB8z9fhJC8z3dF+GAPl6wEA4+DrcFKUPcXRfriO8vVZRZHtWTPYRgwAiODrcKKF1rSeSSx6XvvBIYGoqR3q3AEAovg6nGgFWC20Fouep1Ksg0XYrq4Oe3URNj0PAIDfwwmShCJsAIBx8HU4GSrCNhadcaAImwMowgYAGAdfh5NrFWHTurAUYXMARdgAAOPg63CidyB2sh/GQBE2AMA4+Dqc5E9Nd7QfxkARNgCAF8NJIu6tc93bWNnuGh+KsAEAvBhOEnFvnY7uXkf74TqKsGXPiGzXsvbarucBABCRSX5+FwqypjjaD9egAeSf/nlw944uktW1KDrloyMrAABc4etwMlQhtqXrM7MzZ7TZnCIqxDpLg8jNdzn8pAAAL3HNtE6iK8RGLysZOqZCLAAAyeXrcKJWzJsh//E/7zAjJFfTY23X8wAAIHl8Pa0zRAPIN8qKTFE2rWmia0x0yod76gAAkHyEkys0iCy+5YYU/AgAAMDVfD+tAwAA7EI4AQAAViGcAAAAqxBOAACAVQgnAADAKoQTAABgFcIJAACwCuEEAABYhXACAACs4roKseHw4P2DQ6FQqi8FAABcp6G/20N/xz0RTurr682jr6/PHJeUlKT6kgAAwDhdvHhRcnJyYvYJhK8nwlhkYGBAzp49K1lZWRIIBEacLy8vl6ampgmnOg09p0+fluzsbAeuFk6L5+frVm76nm251mReR6Jfy+nnd+L5+Jz1vvIE/F5r3NBgMnPmTAkGg94YORmi39CNN9445vm0tLS4g4X+94QTOznx83UbN33PtlxrMq8j0a/l9PM78Xx8znpfWoJ+r681YuLZBbE1NTWpvgQkkB9/vm76nm251mReR6Jfy+nnd+L5bPk5I3FS/TN23bROIum0jqa6rq4uK/71BwBew+csfDlyEo/09HTZsGGD+QoA4HMWqcHICQAAsAojJwAAwCqEEwAAYBXCCQAAsArhBAAAWIVwAgAArEI4mQAtb3/33XdLWVmZ3HbbbfLb3/7W+Z8MAEAefvhhyc3NlW9961u8Gz7CVuIJOHfunLS2tsqCBQukpaVFFi5cKH//+99l6tSpzv+EAMDH3n77bXM/lpdeekleeeWVVF8OkoSRkwmYMWOGCSaqqKhI8vPz5fz5807/bADA93SUWm/0Cn/xZDjZv3+/rFy50tz5UO9cvGfPnhF96uvrZfbs2TJlyhSprKyUxsbGCb3WkSNHpL+/39zNGAD8JJmftfAXT4aTnp4emT9/vvk/xWh27twptbW1plT90aNHTd/ly5dLW1vbcB8dGZk3b96Ix9mzZ4f76GjJd77zHdm+fXtSvi8A8ONnLfzH82tONM3v3r1bHnrooeE2Te/l5eWydetWczwwMGBGPp566ilZu3btdT1vb2+vfOMb35AnnnhCvv3tbyfs+gHAz5+1Q+tO9DlYc+Ifnhw5iaWvr89MxSxdunS4LRgMmuODBw9e13Nonlu9erXcc889BBMASNBnLfzLd+Gko6PDrBEpLCyMaNdj3XlzPf70pz+Z4UqdX9UhSX289957CbpiAPDnZ63SMPPoo4/Ka6+9JjfeeCPBxicmpfoC3OjOO+80w5MAgMR64403eIt9yHcjJ7rtNy0tzdQpuZoe67ZgAACftUgt34WTyZMnm6Jp+/btG27TURA9Xrx4cUqvDQC8gs9axMOT0zrd3d1y8uTJ4eNTp07JsWPHJC8vT2bNmmW2tlVXV8uiRYukoqJC6urqzJa4NWvWpPS6AcBN+KxFwoQ96K233tLt0SMe1dXVw31+8YtfhGfNmhWePHlyuKKiIvznP/85pdcMAG7DZy0SxfN1TgAAgLv4bs0JAACwG+EEAABYhXACAACsQjgBAABWIZwAAACrEE4AAIBVCCcAAMAqhBMAAGAVwgkAALAK4QQAAFiFcAIAAKxCOAEAAGKT/w82b7Nw8l5wsQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for il, lengthscale in enumerate(lengthscales):\n",
    "    plt.scatter( [lengthscale]*n_folds,  chi2_arr_dict[il])\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')"
   ]
  }
 ],
 "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.13.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}