{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "049ed0ad",
   "metadata": {},
   "source": [
    "# Optimization Tutorial\n",
    "\n",
    "Trey V. Wenger (c) June 2025\n",
    "\n",
    "Here we demonstrate how to optimize the number of cloud components in a `bayes_spec` model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f1c46dd4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:06.989214Z",
     "start_time": "2024-08-05T19:59:04.772064Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pymc version: 5.22.0\n",
      "bayes_spec version: 1.7.9-staging+2.gb86248c.dirty\n"
     ]
    }
   ],
   "source": [
    "# General imports    \n",
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import pymc\n",
    "print(\"pymc version:\", pymc.__version__)\n",
    "\n",
    "import bayes_spec\n",
    "print(\"bayes_spec version:\", bayes_spec.__version__)\n",
    "\n",
    "# Notebook configuration\n",
    "pd.options.display.max_rows = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70d808d8",
   "metadata": {},
   "source": [
    "## Model Definition and Simulated Data\n",
    "\n",
    "Like in the basic tutorial for `GaussNoiseModel`, our model is a Gaussian line profile with the spectral noise as a free parameter. We generate synthetic data as in the basic tutorial, and the data key must be `\"observation\"` for this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "da706708",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:10.588998Z",
     "start_time": "2024-08-05T19:59:06.991731Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHGCAYAAABke8+ZAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAmsBJREFUeJztnXl8VNXZx3+TZBKykSEBVBYlE3d9VbJUizvJoEXrRkLUulUlAa19W60JobZKrSWhWFtfFxKor7XVShKXVkXNROpaW0kGsFRRyYAVwQWSScCQZJLM+0feczn3zl3n3lnzfD8fPkzuds6cufec332e5zzHFggEAiAIgiAIgiBkSYp2BQiCIAiCIGIZEksEQRAEQRAqkFgiCIIgCIJQgcQSQRAEQRCECiSWCIIgCIIgVCCxRBAEQRAEoQKJJYIgCIIgCBVILBEEQRAEQaiQEu0KJAKjo6PYvXs3srOzYbPZol0dgiAIgiB0EAgEsH//fkybNg1JScr2IxJLFrB7927MnDkz2tUgCIIgCCIEPvvsM8yYMUNxf9yKJZ/Ph0WLFqGyshLl5eWKx61cuRL79u1DXl4eurq64HK5VI8PhezsbABjjT1x4kRLrx1v+P1+tLW1Yd68ebDb7dGuTkJDbR0ZqJ0jA7VzZKB2FtPX14eZM2cK47gScSeWKioqkJubCwBobW1FZWWl4rHV1dUoKChAQ0ODsM3lcqG7uxtVVVWW1Ym53iZOnEhiye9HRkYGJk6cSA9imKG2jgzUzpGB2jkyUDvLoxVCE3diqaWlBQDg9XrR1NSkeJzH40FTUxOk6wQ3NDSgqKjIUrFEEARBEETikrCz4RobG1FYWBi0nW1rbW2NdJUIgiAIgohDElYstbe3w+l0yu5zOBxwu90RrhFBEARBEPFIwoolr9crxDZJyc3NRUdHR4RrRBAEQRBEPBJ3MUt68Pl8qvsdDofmMWoMDg5icHBQ+Luvrw/AWOCc3+8P+bqJAPv+470dIgG1dWSgdo4M1M6RgdpZjN52SEixFG5WrFiB5cuXB21va2tDRkZGFGoUe5CbM3JQW0cGaufIQO0cGaidx+jv79d1XEKKJYfDobrfjFUJAOrq6nD77bcLf7M8DfPmzaPUAX4/3G43XC4XTUsNM9TWkYHaOTJQO0cGamcxzDOkRUKKJS26u7sVg7/1kJaWhrS0tKDtdrudbr7/h9oiclBbRwZq58hA7RwZqJ3H0NsGCRvg7XA40N3dLbvP5/OhuLg4wjUiCIIgCCIeSVixtHDhQni9XsX9LpcrgrUhCIIgCCJeSVixVFFRAY/HExSf1N7eDgAoKyuLQq0IgiAIgog34lYsMRGk5GorKytDeXk5VqxYIdre0NCAlpYWzSBwgiAIgiAIIA4DvGtra+H1euHxeIS/3W43cnNz0djYKDq2paUFK1euRG1tLfLy8tDV1YXq6mqUl5dHo+oEQagwPDyMhoYGlJWV4fTTT492dQiCIATiTiw1NDQYOr6mpiZMNSEIwkrWrFmDu+66C3fddVfQAtgEQRDRJG7dcARBJBYffPCB5jGBQADPP/88duzYEYEaEQRBjEFiiSCImCA5OVnzmI0bN+Lyyy/HzTffHIEaEQRBjEFiiSCImCApSbs7+vzzzwEAX3/9dbirQxAEIUBiiSCImEDOsrR161Y88sgjGB0dBXBoHafh4eGI1o0giPFN3AV4EwSRmMiJpcLCQvj9fthsNixZsgTffPMNAGBkZCTS1SMIYhxDliWCIGICObHk9/sBjKUBASCIJbIsEQQRSUgsEQQRE6jFLH366acADrnhyLJEEEQkIbFEEERMILUssTgl4JBYIssSQRDRgGKWCIKICXixNDIygr6+PtHfX3/9NQV4EwQRFciyRBBETMCLpaGhoaB1H3fs2EEB3gRBRAUSSwRBxAR8zJKcWPL7/eSGIwgiKpBYIggiJuDF0uDgoKxYIjccQRDRgMQSQRAxAe9ak7MsDQ0NkRuOIIioQGKJIIiYQEsskWWJIIhoQWKJIIiYgBdAWjFLZFkiCCKSkFgiCCIm4AWQXMwS74YbHR0V5WEiCIIIJySWCIKICfRYlpgbDiDrEkEQkYPEEkEQMYGemCVmWZIeTxAEEU5ILBEEERNILUu9vb2i/UNDQyLLEgV5EwQRKUgsEQQRE0jFks/nAwDk5OQAGFtE1+/3C8eQZYkgiEhBYokgiJhAGuDNLEtTpkwBAEE8MciyRBBEpCCxRBBETKBkWSKxRBBEtCGxRBBETMBblg4ePIi+vj4AymKJ3HAEQUQKEksEQcQEvKVo3759wufJkycDIMsSQRDRg8QSQRAxAW8p+uqrrwAAEyZMQHZ2NgAEzY4jyxJBEJGCxBJBEDEBbyliYiknJwd2ux0A0NPTo3g8QRBEOCGxRBBETCBnWXI4HIJYIjccQRDRgsQSQRAxAS9+vv76awBjlqXU1FQA8m44csURBBEJSCwRBBETGLUsPfnkk5g4cSJefPHFiNWRIIjxSUq0K0AQBAEoxywxy5KU+vp6AMB3v/tdBAKB8FeQIIhxC1mWCIKICXjL0oEDBwCILUsEQRDRgsQSQRAxgVzANj8bjiAIIlqQWCIIIiaQC9bOzc1VdMMRBEFECopZIggiJpCzLE2fPh2jo6NRqA1BEMQhyLJEEERMIGdZmj59OrnhCIKIOpZZlvr6+tDd3Q2fzwen04mJEydadWlLWLlypbDelM/nQ1FREaqqqqJcK4IgGEqWpe7u7ijUhiAI4hAhi6W+vj40NTXB7Xajvb0dAETTd202GxwOB4qLi1FRUYGbb77ZfG1DpLq6GrW1tXA6ncK2pqYmVFdXo7GxMWr1IgjiEEqWpW3btkWhNgRBEIcw7IbbuXMnFi5ciPz8fLz33nsoLy9HR0cHuru7MTo6Kvzr7u5Ge3u7sP/oo49GZWUldu7cGYavoUx7ezscDodIKAFAVVUVOjo6ghLdEQQRHaSWpaysLGRnZ1OAN0EQUceQZWnNmjVobGzEsmXL0NzcrHpsTk4OZs+ejdmzZ2PRokUAgNbWVlRVVWHevHn4yU9+EnqtDeDxeOD1emX3OZ1OeL1eFBYWRqQuBEEoIxVLQ0NDAEAxSwRBRB3dlqWlS5eit7cXHR0duOKKK0IqrLy8HG1tbcjJyUFdXV1I1zCK0+lEa2srmpqagvZ5PB4SSgQRI0jdcEpiKSWFJvESBBFZdImlTZs2obKy0jJr0KJFi7Bw4UJs3rzZkuupUV5eDqfTierqarhcLvh8Pvh8PlRUVKClpSXs5RMEoQ+pZSknJwcAgtxwsTZ5hCCIxEfXK9rs2bMtLzgc11Sis7MTFRUVaG9vx6RJk1BYWIjXXnsNDocjpOsNDg5icHBQ+Luvrw8A4Pf74ff7rahy3MK+/3hvh0iQaG3NLEuNjY148skn8fOf/xx+vx82m0103MSJE4NmyIWzDRKtnWMVaufIQO0sRm87jAt7tsPhQEVFBRwOB1pbW+HxeLBo0aKQLUsrVqzA8uXLg7a3tbUhIyPDbHUTArfbHe0qjBsSpa33798v/H/77bfjwIEDWL9+vWLMIc/69evDXb2EaedYh9o5MlA7j9Hf36/rOFvAwHLda9euDSkFwJIlS/Doo48aPs8qXC4XqqurUV5eDq/Xi+rqarS3t8PpdKKzs9OwhUnOsjRz5kzs3bt33LsI/H4/3G43XC4XBeaGmURr66OOOgp79uzBe++9h9NOO03Y/sEHH4j+Puuss/D222+LzmXxTeEg0do5VqF2jgzUzmL6+vowefJk9Pb2qo7fhixLjY2NIYkllocpGlRXV6OiogLl5eUAxgK+3W63kGeptrbWcK6ltLQ0pKWlBW232+108/0/1BaRI1HamsUspaeni76P1FrLYpl4IvH9E6WdYx1q58hA7TyG3jYwJJY6Ozvx3HPP4fLLL9d1/M6dO1FRUaHLjB4umpqaIGc8q6qqgs/nw7p166JQK4IgpLCYJelsNz7AOy0tDZmZmRGtF0EQhOGklDU1NboSS65atQoFBQXo7OwMpV4RgVmbCIKIPsyylJycLNrOv/nNmjWL3oYJgog4hsRSTU0NPvnkE6xevVpRMO3cuRMlJSWoqalBfn4+WlpakJ+fb0VdQ6KwsFDRDdje3o7KysoI14ggCDmULEu8ODrqqKOCxBRBEES4MSSW6uvrhf9Xr14tTJlnrF27VrAm1dTUYPv27ViwYEFU8xm1tLSguro6yBXo8XjgdrtRU1MTpZoRBMGjZFni3XDTp0+npJQEQUSckHud+vp6LF26FMuWLUN3dzcqKirg8XgEaxKfRymSOZWksIDuhoYG0faCggJKSkkQMYQey9L06dPx5Zdfyp5LFieCIMKFqVe0+vp6LF68GGvWrEEgEEBNTY1gfYolnE6n4RlvBEFEFj0xS9OnT8e+ffuCzh0aGkJ6enp4K0gQxLjFcIC3lNWrV2P27Nlob2+PSaFEEETsMzo6KnyWWpb4v5XccOHMs0QQBGFILCktftvR0YHm5mbD5xEEQQDideGkliV+uZNvfetbsu42EksEQYQTQ2JJLV9SQ0ODoiiKZp4lgiBiHxavBARbloCxWbb//ve/cdhhh5FliSCIiGMoZqmlpQWbNm2SzaALAD6fD62trUHLh3g8npArSBBE4qNmWQLGUgao7SexRBBEODEc4L19+3bDhUhXDScIguDRsizxyO3n12okCIKwGkNuuMLCQoyOjhr+F83UAQRBxD5aliUesiwRBBFpDIklp9MZUiHFxcUhnUcQxPiAiSWbzYakJPVuiWKWCIKINIbE0po1a0IqRJoQkiAIgoe54fQkliSxRBBEpDEUs6QU2B2u8wiCSHzuuusubNq0CYB2vBJAbjiCICKPLsvSjh07LM+VVFdXh82bN1t6TYIg4ovR0VHcd999WL9+PQCyLBEEEZvoEkv5+flYuHAhLrjgAnz66aemCty5cycuuOACuFwunHbaaaauRRBEfCOdxRaqZYlmwxEEEU50u+Fmz56NdevWYeHChbDZbKitrcXcuXN1F7RhwwasXr0avb29WL16NfLz80OqMEEQiYNU5JBliSCIWMRQzJLD4UBbWxtee+01rF69GuXl5SgoKEBZWRny8vLgcDiQm5uL7u5u+Hw+7Nu3D+3t7fB4PCgsLMTSpUuxYMGCcH0XgiDijFAsSySWCIKINIaTUgJAaWkpSktLAQDPPPMMNm7ciPfeew8+nw9erxcOhwNOpxO5ubmoqqpCWVkZWZIIgggiFMsSBXgTBBFpQhJLPAsWLCBrEUEQISEVSxMmTNA8JzU1NWgbiSWCIMKJoTxLBEEQViIVS3rSjJxwwglB20gsEQQRTkgsEQQRNaRiaeLEiZrnnHrqqZrXIQiCsBISSwRBRA2pyMnOztY8h9xwBEFEGhJLBEFEjVAsSwCCJoyQWCIIIpyQWCIIImpIRY5eseR2u3HppZeipKRE9joEQRBWQmKJIIioEaplqaCgAM8//zzOPfdcACSWCIIILySWCIKIGqHELPGw+CUK8CYIIpyQWCIIImqEalliMLFEliWCIMIJiSWCIKIGiSWCIOIBEksEQUQNs2IpLS0NAIklgiDCiyViafPmzbjggguQl5eHtWvXCtsXL16MDRs2WFEEQRAJiFUxSySWCIIIJ6bF0qZNmzB37lzk5OSgvr5etG/16tXo6enB5s2bzRZDEEQCQm44giDiAdNiqb6+Hp2dnWhubsaiRYuC9i9YsADt7e1miyEIIgGxSiw9//zzuOeeezA6OmpZ3QiCIBimxVJ+fn5QNl2CIAg9WCWWAGD58uVobW21pF4EQRA8psXS5MmTRX8HAoGgY/bt22e2GIIgEhCrYpYYH3zwgek6EQRBSDEtlrZv344tW7YIf9tsNtH+VatWmS2CIIgERSqWsrKyDJ3PZsMxDhw4YLpOBEEQUlLMXqC+vh5OpxMulwslJSXo6upCbm4uvF4vGhsb4XA4sHHjRivqShBEgsGLpeXLlyM5OdnQ+VLLEoklgiDCgWmx5HA40NHRgerqatTU1AAAGhsbAQA1NTVBM+QIgiAYTCytXLkSd955p+HzSSwRBBEJTIslAHA6nXC73ejt7UVHRwdyc3Mxe/ZsKy5NEEQCw8SS1J2mF6lY2r9/v+k6EQRBSDEtlp555hk0Nzdj3bp1yMnJQWlpqRX1CgtNTU3o6uoS/i4oKEBVVVUUa0QQ4xurxVJfX5/pOhEEQUgxLZYaGxvh9XrR19dneNpvpPD5fKioqEBFRQUaGhoAAB6PBxUVFSguLkZhYWGUa0gQ4xOWTDJUsSQ9b+/evabrRBAEIcX0bDiXy4Xt27erCqVoz4hbtGgRnE6nyIrU3d2N7u5uOByO6FWMIMY5VluWvvrqK9N1IgiCkGLaslRWVoZVq1ahqqpKUTBFczacx+NBa2uryP0GjNW7p6cnSrUiCAKwXizt3bsXo6OjSEqiNcIJgrAO02KpubkZPp8P+fn5cDqdyM3NFVlrfD5fVJc7WbFiBRwOB5xOZ9TqQBCEPFaLpdHRUezbtw9TpkwxXTeCIAiGJTFLubm5KCoqAjCWwTuWLDZerxdOpxM+nw9NTU0AxjKKU3A3QUQfq8USMJY+gImlnp4eTJo0KfQKEgRBwAKx5HQ60dHRoXrMwoULzRYTMh6PB2VlZWhqahLyQAFARUUFOjs7hZxQRhgcHBQl02MzcPx+P/x+v/lKxzHs+4/3dogEidDWAwMDAIDk5OSQvod0xQB2Tb/fj+effx4LFy5ETU0NfvnLX4Zcx0Ro53iA2jkyUDuL0dsOtoDcYm4GeO211zTTBWzatCkqeZd8Pp/wVtnV1SVyxXm9XhQUFMDtdqOsrMzQde+55x4sX748aPtTTz2FjIwMc5UmiHHEokWL8PXXX2PlypU49thjDZ8/ODiIyspK0baHH34Y06dPx8033yzMjnv++eetqC5BEAlGf38/rr76avT29qpOVDMtlvTw7LPP4oorrgh3MbLYbDY4nc6gAG+2r7y8HC0tLYauKWdZmjlzJvbu3Ruz6RMihd/vh9vthsvlgt1uj3Z1EppEaOuZM2fiyy+/xMaNG3HqqacaPn9kZATp6emibZs2bcJJJ52Ek08+GR9//DGAQykKQiER2jkeoHaODNTOYvr6+jB58mRNsWRJBm8tVqxYETWxBEAxPYDD4YDX6zV8vbS0NNkYC7vdTjff/0NtETniua2ZC3vy5MkhfQe5c5KSkmC325GZmal6XChlxWs7S/n888+Rnp6O3NzcaFcliERq51iG2nkMvW1gWiyVlJQo7vP5fOju7jZbhClYcLccStsJggg/Q0NDOHjwIAAgJyfHsusODw8DALnEFfD5fJgxYwaAsQk5BEFoY1osdXV1obi4GA6HQ/SW0t3dDY/Hg6KioqhO2y8rK0Nzc7Pi/uLi4gjWhiAIRm9vr/DZSvf1yMgIABJLSjDXJDDWVsnJyVGsDUHEB5bMhmtra1Pc/8wzz0R16m5tbS2ampqEFAIMj8cDAKiuro5W1QhiXMMsu9nZ2ZYO2GRZUod3O/T39yM7OzuKtSGI+MB0mts1a9ao7l+wYIEgTKKB0+lETU1NkChatGgRampqaF04gogSzLJk9ZJDcpYlcjcdgs9u/s0330SxJgQRP5i2LEUjJYBRGhoa0NTUhIqKCuTm5qK7uxt1dXUoLy+PdtUIYtzCLEtWxisB8palgYGBoFlz4xU+rwyJJYLQR0Rmw0VzbThGVVUVZewmiBgi3JalCRMmCNv2799PYun/4dMokFgiCH2EfTac1+tFQ0OD2WIIgkgwwm1Z4l1vfX19mDp1qqXlxCsklgjCOJbMhnM6ncIiulJWr16tmeGbIIjxR7gtS0w0AWOWJWIMEksEYZyIrA1HEAQhJVyWJSaW2P/AoeSXBESrD5BYIgh9hH02HEEQhBxWWZb+/Oc/4zvf+Q5OPPFEAIcsSmRZkocsSwRhHNNiSWk23M6dO+ltjiAIWT7//HM0NjYCMG9ZuvLKK7F+/XpMnjwZgLxlicTSIXixdODAgSjWhCDiB9Niqa6uLmhbb28vurq64Ha7sWrVKmzYsMFsMQRBJBB33HGH4A6yKmaJJbaUsyzRi9shyLJEEMaxJMBbSk5Ojiioe9WqVZg7d67ZogiCSBDWrVsnfE5NTbXkmkwskWVJHRJLBGEc05Ylm82meYzb7TZbDEEQCcRZZ50lfD7ppJMsuWZKyti7H1mW1KEAb4IwjmHL0tKlS+H1eoXgzI6ODlxwwQWKx3d0dFAySIIgRIyOjgIYW3ZILVebEdQsSwMDA5aUkQiQZYkgjGNYLNXX1wMAWltbUVVVBZvNprjuksPhQH19PRYtWmSulgRBJBRsyY1LLrnEsmuqWZZ44TQe6e7uRmdnJ0pLS0ksEUQIhByzVF5ejsLCQixduhTNzc1W1okgiASHDdh2u92ya6pZlsa7WPrWt76Frq4uNDU1kVgiiBAwFbPkdDpRWVlpVV0IghgnMMuSVcHdQLBlicTSIdhEnNbWVhJLBBECpmfDLViwQHHfr3/9a9hsNlRVVWHixIlmiyIIIkGIhGWJ3HDBJCUliQK8Kc8SQejD9Gw4Ne6880785Cc/wYoVK8JZDEEQcQZZlqJDcnIyWZYIIgRMW5YYzz77LLxeL/bt2yfa7vP54PV6rSqGIIgEgCxL0SEpKYnEEkGEgCViqbi4GB6PR/ibZeT1+XxwuVxYvXq1FcUQBJEghNOyJBfgzVIVjHeklqX+/v4o1oYg4gfTbrilS5eirKwMPT09GB0dRUtLC7q7u9Hd3Y3R0VEhvQBBEAQjnJYlSh2gjFQs8Z8JglDGkpil+vp6YTFMp9OJzZs3C/sWLFiA1tZWK4ohCCJBiLRlicTSGNIAbxJLBKEP02KJrfTNcDqdonWfCIIgpJBlKTpILUtMtBIEoY5psbR3714AwIYNG7Bz507k5OSgs7MTn376qXAMrQ1HEARjdHRUEC9WWpYoKaU25IYjiNAwLZaqq6uxePFilJWVobq6GgBQVVWFwsJCLFmyBBdccIEQ8E0QBMFbM6y0LNFyJ/Lwwe1kWSKI0DA9Gy4/Px+rV69GRUUFnE4ngLGlULq7u7F06VLk5eWhs7PTdEUJgkgM+AGaLEvhhxdH0tQBw8PDGB0dRVJSWFPuEUTcY9kTUlpaivz8fOHvqqoqdHd345NPPqHs3QRBCPCDNVmW9BEIBETfxwgDAwPC5+TkZFGAN0DWJYLQg2mx9Mwzz9D6cARB6IYfnJnAsYJEtixdeOGFcDqdIuGjxUcffYTi4mI8/fTTwjabzRYUp0RiiSC0Md1TNTY2wuv1oq+vjyxIBEFowgZnu91uaQ42NctSvCelbGtrAwC8/fbbKCsr03XODTfcgM7OTlEYhN/vDxJLFORNENqYtiy5XC5s375dVSitWrXKbDEEQSQIbHC2Ml4JSGzLUihIl54C5MUSWZYIQhvTYqmsrAyrVq1CX1+f4jEbN240WwxBEAkCb1myEmZZeuyxx3DLLbckVMwSIxAImDqfLEsEERqm3XDNzc3w+XzIz8+H0+lEbm6uKFWAz+dDe3u72WIIgkgQwm1ZAoBHH30UWVlZwt+JIpbMuhOHh4cpwJsgQsCSmKXc3FwUFRUBGHvz6enpMV0xgiASk3BblhiJYFl67bXX0NzcLPxtxLIkdyxZlggiNEyLJafTiY6ODtVjFi5caLYYgiAShEhYloDEiFmSBnOHww1HliWC0MZ0zFJDQ4PmMXV1dWaLIQgiQYiUZYkXAfEqlqRYKZbYTESyLBGENqbFUmlpqeYxs2fPNlsMQRAJQqQsSzyJIpbM4vf7BRHJYrrIskQQ2liSwXvz5s244IILkJeXh7Vr1wrbFy9ejA0bNlhRBEEQCUK4LEtqYine8ywxzFqW9u/fL3zOzs4GQJYlgtCDabG0adMmzJ07Fzk5OaivrxftW716NXp6erB582azxRAEkSCEy7Kklg08USxLcmLJyHfz+XzCZzZrmSxLBKGNabFUX1+Pzs5ONDc3Y9GiRUH7FyxYQKkDCIIQiIZlKdHEEmvDpqYmOBwOvP3227rOZ2JpwoQJmDBhAgB5y9LWrVs1J+4QxHjCtFjKz88XLaAbL7BUBwRBRBayLIXO6OgoPvjgA2RnZ+Ouu+5CdXU1Dhw4gKuuukrX+UwsZWVlCe0vFUv79u3Df/3Xf6GkpAQHDhywtP4EEa+YFkuTJ08W/S1nJpZLux9Nqqur4fF4ol0NghiXkGUpdEZGRrB06VIMDg7ivvvuUz1WLb4pMzNTaH+pG46PO1VbmYEgxhOmxdL27duxZcsW4W/pwpixti6cx+Mh8zJBRBGyLIXO8PCwKNmmHH6/H7W1tejq6lI8JjMzU9Gy9Ic//EH4nCjtRhBmMZ2Usr6+Hk6nEy6XCyUlJejq6kJubi68Xi8aGxvhcDhiam24devWobKykixLBBElyLIUOsPDw5rf5Te/+Q1WrlypeoyaZWnPnj3C50RpN4Iwi2nLksPhQEdHB7q7u1FTU4PGxkaUl5ejpqYGCxYsiCmhtHLlSkqQSRBRhixLoTMyMqJpWfrTn/6keR2lmKVAICCKU0qUdiMIs5i2LAFjS5643W74fD50dnYiNzc35hJRejweOJ1O0SK/BEFEHsqzFDpabji/34+tW7dqXod3w/GWpaGhoYRYU48grMYSscRwOBwoKSmx8pKWsW7dOl1Ls+hhcHBQtHI3C4Lks+OOV9j3H+/tEAnita0HBgYAjFmCIlX3kZGRkMuKpXYeHBxUrIff71eNU+JJT09HUtKYY+HgwYPCNaWLoA8MDETse8dSOycy1M5i9LaDZWKprq4OTU1NwtRUh8OBZcuW4Y477rCqiJCx2v22YsUKLF++PGh7W1sbMjIyLCsnnnG73dGuwrgh3tr63//+N4Cx2Jj169dbfl05Dh48aLqscLVzIBAImhijxObNm7F3796g7ez7ffrpp7qu093dLVjb/vWvfwlt89VXX4mOe/3117Fjxw5d17SKeLuf4xVq5zH6+/t1HWeJWCouLobH40F5eTmcTicAoLOzE3feeSfcbjdeeeUVK4oJiXC43+rq6nD77bcLf/f19WHmzJmYN28eJk6caFk58Yjf74fb7YbL5bLczUKIide2/sc//gEAOProozF//nzLrqv2jCcnJ4dcVjjb+bbbbsMrr7yCjo4O5OTkaB5//PHH47333gvanp6ejvnz5+ue6XvCCSegv78ff/vb3+B0OoW2+eCDD0THzZkzB6eeeqqua5olXu/neIPaWYze9BimxdLSpUvhdDrx2muvBT3sPp8PVVVVWLt2LW6++WazRYWEle43RlpaGtLS0oK22+12uvn+H2qLyBFPbd3T0yMsi5SWlmZpvVlGajlGRkZMlxWOdl6/fj0+++wzbNu2DWeddZbm8YFAQDH+ym6363YpZGdnC9cZHR0VvhcfXgAASUlJEb+34ul+jmeoncfQ2wamxZLX60Vzc7PsPofDgebmZixZsiQqYqm1tRUejwfV1dWi7ezti21vaGigwG+CiAArVqwQPls9Gy4eUwcwcXPw4EFdxyvNhmNuPL3XyczMFNwP/Gw4acbuWG03gog0psUSc7uZPSYclJeXo7y8PGg7y+Dd2NgYhVoRxPjl3XffFT7z+XysQC51gM1mQyAQiNlBnwkVvXETWrPhWPC8FpmZmejt7RXVASCxRBBKmM6zpAc5sbRhw4ZIFE0QRAyRl5cnfJ46daql15azLDHrVawO+kyo6LUIaYklI5YluaSUJJYIQh7TYsnlconWEpKyYcMGTJo0KWh7NK063d3dUSubIMYzbLbVscceizvvvNPSa8tZlphYitU8S1ZblvSKJaWklCSWCEIe0264xsZGvPbaa2hsbERubq5oX3d3N7xeL4qLi0VB1t3d3VFZbqSpqQlutxutra0AgKKiIhQXF5M7jiAiBBNLjz32WNAi3GZRsyyNjo4amqKvl3vuuQdTpkzBrbfeavjcQCBguWXJiBtOzrL0zTffBJVHEIQFYqm9vR1OpxOTJk0KWuV60qRJKCoqCtquthp2OKmqqkJVVVVUyiYI4pBYstoFB8hblvhZq6Ojo6pB4EbZtm2bkG/tlltuMSzEeCESacuS0kK6ZFkiCHksCfDWm9uDZ+HChWaLJggijjh48CD2798PIDxiSc2yBIwN/FaKJRYgDYxZZLKysgydz4sUs5Yltk3tOunp6cJ+peVOSCwRhDymY5ZCzWFEC9oSxPiCWZVSU1PDkrxVy7IU6sB/4MABWWs4L1r0Jrbj4cWSXsvSyMiI7Pdg+ZHUxBIvHLOysgQ3HFmWCEIb02KptLRU85hnn302aFusLbRLEER4YWLpsMMOszx2CNBnWTLKli1bkJubKxvXyKxk0s960RJLckHpSpYlJpbUYpb4Mo488kiyLBGEASKSOoACqAmCCGe8EiAfC2lWLN19990AILtkE29NUrIseb1exTXrtNxwcvUdHh6WzdKtx7LEnzdx4kRZy5I0wJvEEkGMYcnacKtWrcK6deuERXSleL1eK4ohCCKOYWJpypQpESvTrFjat2+f4j49YqmgoAAAsHfvXlGOKSA0y5Lf7xedxxgeHsbo6Kiu2CcmkljbvPPOOxgcHERaWhpZlghCAUvWhmtqakJxcTGKioqC9u/bt4/yGhEEIVgtsrOzw3L9I444Atdffz2efPJJwVUVTbHEW7r+85//qIolvZal/v5+xdnEg4ODusTSUUcdBeCQaDpw4ADuuOMOPPTQQySWCEIBS9aG0xJDNPONIAjmBgrn4p2PP/44SktLcd111wEAMjIyhH2hJKY0I5a0FrXVsizJCRWpmOEZHBzUlWeJiSV+ZYX//Oc/AMgNRxBKmI5ZKikp0Twm1BlzBEEkDpEQS4B4Btxpp50mBH5H2rLEiyG5gHYty5KcuNMSS0YsS7Nnz8Zdd90lKouJNiYySSwRxBgRCfDesWNHJIohCCKGiYZYmjNnjimxpHYOn2dJbjYcb1nSEktWWZaYWLr88suD9t94441ITU3FsmXLhG3HHHMMgENpEFg9mKuUxBJBjGFaLFVVVWHVqlXYuXOn4jE0G44giEiJJV5QnHHGGSKx9NZbb2HdunWmy3jggQfw2GOPCX9rueHk4oz4/XpjltRSFPBuuBtuuAFtbW244oorAIyFQqxduxY9PT1C0DlwKDcViSWCUMd0zFJOTg727t2LgoICOBwO5ObmwuFwCPt9Ph/NhiMIImJi6bDDDhM+T548WSSWzjnnHABAYWGhYFVRghc40oSXt99+u+hvLbEklxsplNlwei1LGRkZKCsrw9lnn40NGzbgvPPOg81mE8VwAYe+FxNFJJYIQh7TYmnx4sVobm5GaWmpKGCQQbPhCIIAIieWSktLsWbNGhQXFwM4lKySd5sppTnh4Y/XmsGnFbMkF+wdymw4LcsSu056ejoAYMKECZg/f77iOaxthoeHMTIyIlimSCwRhBjTYqm7u5tmwxEEoUmkxJLNZsPNN98s/M0Ewa5du4RtfFyTEiwvFH8NJbQsS1piSW/MktqyKrxYmjBhgmp9GbwbjhdsJJYIQozpmCWXy6V5DM2GIwgiUmJJSlLSWDf32WefBdVFDX4mnFbaAbNuOL2z4ZRyLAHimCVmWdKCd8Pxgo0tCkxiiSDGMC2W9JizaTYcQRDREktyliW5LNhSeLGjJZa0ZsNpWZb8fn+QoDIqVOTccFrwbjgmltLT04XfSE7kEcR4hGbDEQQRVnbs2IEvv/wypsSSHsuSEbHExzfJna8lloBg65IZsRSKG47PsWQm3QJBJCKmY5aqqqrg8/lQW1tLs+EIghCxb98+YeLHNddcAyB6YsmoG05JLPGfZ8+ejU2bNsnGHBkJ8AbG4pb4QHKjGcf3798vuOnMuOFILBFEMKbFktvtRnFxMRYsWIDc3Nyg/TQbjiDGLx999JHwOd4sS7wLihcu/PaHH34Yc+bMweDgYND5RmKWgOAgb6NChXcF6glgB6yxLLW0tGDWrFm6VnMgiHjFtFhyOp1oa2tTPYZmw0WG3t5e2Gw2TJw4MdpVIQgA4szVTBzEu2WJFz6ZmZkA5GOgIu2G48USv4CwGnIxS0bEUmdnp9C/qwWfE0S8Yzpmac2aNZrH0Gy48OP3++FwOJCTk0NBmUTMwIslJgaiJZa03GJS+GN40SAnlrQsS3rdcDxqbri8vLygbUwsJSUlaaY6YJh1w23btk1XOQQR75gWS7Nnz9Y8Jj8/32wxhAZff/218Fktyy9BRBI2bR+Ivlji0TMbTo8bjoklv98fJG74MvS44YxYlqZPny58Zu3JxJJeqxIg74bLzMzULZaMxlURRLxiyUK6mzdvxgUXXIC8vDysXbtW2L548WJs2LDBiiIIDfhOS27RTiJxCQQCsgHGsQB/L7I6RivPEo9Vbjh++ZBjjz0Wf/rTn2TPD8WypCZUpk2bJnyePHkygNDEklk3HIklYrxgWixt2rQJc+fORU5ODurr60X7Vq9ejZ6eHmzevNlsMYQGfLwAdWDjixtuuAGZmZn45JNPol2VIGJBLMlZlowGeAcCAeEZY+empKSIAqm7urpw7bXXypaxe/duHH300Vi+fLmwTcuypPQcZ2VlieISmUvOjGUpVDcc9TXEeMG0WKqvr0dnZyeam5uxaNGioP0LFixAe3u72WIIDZTefInE589//jMA4MEHH4xyTYLhRXy8iSXpMSMjI/jd736HGTNmABgTGmrfhT//l7/8Jbq6unDPPfcoXl+vZWnq1Kmics1YlszOhiOxRIwXTIul/Px8ikmKAfhOjXKjELECL9wTQSz96Ec/Ev5OSUlBUlKS4vfhLUdyM8VCtSxNnTpVEDnAIbHEllwhNxxBWI9pscQeVIZcp8CvsUSEB6WZO8T4IRanbvP3YiyJJaMB3kDwc8UEi1JOIy1BFmrMkpJYCocb7sEHH8QLL7ygeH4s3nMEEQ5Mi6Xt27djy5Ytwt/S4OJVq1aZLYLQAd+xk1giYgU5yxI/0EcCKy1LPFaLJb2z4aZOnSpKVRBONxwAXHLJJYrn83Uk4UQkMqZ7rfr6ejidTrhcLpSUlKCrqwu5ubnwer1obGyEw+HAxo0bragroUKsWJZ27tyJ9evX4/zzz4+4BYGIPeREfCxYlmJRLOnNs3TYYYeJMqOz5aWsnA0nFW5KSJeB0ZvfiSDiDdNiyeFwoKOjA9XV1aipqQFwaOHcmpqaoBlyRHiIFctSYWEhDhw4gKysLNx///1RqwcRG8jdi/EilqRuOOnfVomliRMnoq+vz5Blibfms/LNWJYCgYCQny0jI0OXmxIQi6WRkRESS0TCYkgsrV27Fl1dXfD5fMjJyYHNZsOKFSvgdDrhdrvR29uLjo4O5Obm6kpWSVhHrFiWWIdL+bUIQH5mZryIpVAtS0w0aAkOtj8nJwd9fX2GYpYGBgaEv1n57HqhiCVgbLkkYCwppTSx7fDwsKz7VCqWCCJRMSSW6uvrUVRUhKamJuTk5ATtz8nJQWlpqWWVI4IJBALYunUrjj76aGFl8Y0bN+Lhhx8WjglH6oBAIICvv/4aU6dOtfza0nIoqWbiEAuWJSuSUgL6xVJ/fz+ys7N1W5YcDgc+++wzQ7Ph5MQSw0j78kKyp6cHAJCdnY29e/eKjjt48CCys7ODzufjlEgsEYmM4QDvNWvWyAolIjI888wzOOWUUzB37lxh27e+9S08+eSTwt/h6LTuvvtuHHbYYXjssccsvzbj448/xpQpU6Luuq2ursbFF19M06ItIFYtS1bOhpNacr755hsA+t1wrD/dv38//vd//xfvv/++bHmMqVOnioSVVCyFalliYmnixIlBbaYUw0SWJWK8YEgsOZ3OkFa055dAIczBFi7+xz/+oXhMODqte++9FwBwyy23WH5txo9//GPs27cPdXV1YStDD01NTXjppZco87wFxIJlySo3nFRgKVmW3nzzTfzoRz8SrdcoB3N7MWttc3MzbrzxRtx0000AjLvhGGbFUnZ2dlCbKS2nIw3wJohExZAbjs26MEpnZyduvvnmkM4lxOhxUYXzDS+cHaLeoNJIIbeSPGGMWLUshRLgLY3jURJLlZWVuurFXF1HHXWUaHtHRwf27duHRx99NOicU089FZMnT7ZMLPFtw8SbnGVJj1giyxKRyBiyLIUaS+L1ekM6jwhGLv5CSiyIpVByrugZwMINX29aNsY88WBZ2rNnD370ox9h8eLF+Oqrr2SPAQ651xhas+HUCAQCglg68sgjRfumT5+OhQsXBi0T1dDQgM7OTiQlJeH6668HAMyZM8eUWEpKSgrq141YlmjlAGK8YMiy1NnZiQsuuMCQhcnn80V1bTifz4cVK1bA5/PB6/Wiu7sbdXV1KC8vj1qdzBBty1I4rx0OsfTNN9/g73//O8477zxdgzS9KVtLPFiWHnnkEfzud78DABx//PHCkibS+1GvZUmLQCCA3t5e4fpSsTQ4OCg7mzQlJUX4LkuXLkVxcTHmzJmDrVu3io4zIpbYdfnvKmdZ8vl8+Mtf/oLzzz9fFIoRK7NwCSLcGM6z5Ha7DRcSrdlNPp8PtbW1aGhoEASex+NBUVERysvL0dLSEpV6mUGvWBocHMRf//pXzJ07V1iVPNYJh1iqqKjAyy+/jGXLluG+++7TPJ4XS7Fg6Yp34sGy9OWXXwqf2fpqQLDQY3mMGKGKpZGRESGeKSsrC7m5uaL9ehJC2u12fOc735Et36hYSk5OFrVHZmZmUJv94Ac/wEcffYTvfOc7WL9+vbCdxBIxXjAc4N3T02PoX3d3d9QW2l2xYoVIKAFjSRMbGhrQ2toaVYtXqOhxww0PD+NnP/sZFi5ciIsvvjgCtbKGcIiTl19+GQCwevVqXcfzYonccOaJBcvSpEmThM9MSPDxcSywGRALlXC54YaHhwWxNGXKFCEFCIOPR9KDWbHEB3lnZWUhKSkpSCyxjOHseWLESjJcggg3hsSSw+FATk6OoX8OhwOFhYXhqr8qra2tKCoqCtpeVlYGAHFjWbr66qvx7W9/G8PDw7pjlpqamgCoz5qLNcJpyZEOSErEiljq7+/H66+/HveCLRYsS06nU/jMpunz95rP5xM+q4klqRuOfY9QxBKLV5oyZQoyMjJE+5VEh1IcoJViibnY9Gbi5tuIZsMRiUxEArylZuZI4XQ60d3dHbSdWZrk9sUif/7zn/GPf/wD7777rm43HJvZEi1iLcA73sRSZWUlzj//fCxfvjxqdbCCWLAs8ZZtObGkZFnSOxvOqDjRsiwZxUqxxBJPhiKWyLJEJDKGxFKos9r0ukCsxu12izpChsfjAQCUlJREukqmGBwcjHqAdzgJpziJN7H04osvAgD+53/+R/W4cK30fsMNN1iS7kN6LyYnJ0c8hjEUy9Jf//pXvPDCC6LrhEssTZgwwdD5UqyIWWKYsSzFa79DEHowFODd09OD5557Dpdffnm46hMRGhsb4XA4UFVVFdL5g4ODohw8LCjU7/dbbh3hB+/+/n7R4KhUljQ/kNV10nO90dFRw+Xyx1td5wkTJui6Jt92AwMDMRHkrVQHv98fNHvPivru2bMHf/jDHwCMTVcPJREtQ3ov2u32iLfpzJkzhc/MjT00NCTUg3+h+uabb7B//35ceumlQdfhg7/Ztfx+v2GRcPDgQSGoPDc3F5mZmcI+m82mKICVfl+pa14asK2FNGbJ7/erinD+2lY+L+zcWHjmEhlqZzF628GQWMrJycGCBQtQUVEBl8sFp9MpWnYjHmhvb0d7eztaWlpCTrK5YsUKWfdIW1tbUPyBWfgf8t133xXlgeFnpfC89957or+VjgsVPdc7cOAA/vKXv2B4eFi3VYefbWR1nfv7+3Vdk69DZ2enKaFgBX6/X7Xe/KD26aefWtJufObpl19+WXZNML188MEHor9tNpvlv60Rdu3aBWDsO65fvx6BQEBkWfrPf/6jGMv48ccfi/5m1zBqcf/LX/6CP/3pTwDGElO++eabWLlyJVJSUvCzn/0sKJCc8eGHH8q2nVSQer1eQ23MB7uz50Qtez1/7Z07dwqf33jjDXz66ae6y1UilBnXhHGoncdQyiEmxXCeJUZvb29cJpusqKhAY2OjqTxLdXV1uP3224W/+/r6MHPmTMybN8/ywZXvOP/rv/4Ln3zyifD3/PnzZd8ATzvtNOFzVlYW5s+fb2md9FwvKysLP/zhD7Fnzx709PSI3p6V4N+Qra7zjBkzdF2TX0D05JNPtrweRrHb7Yp18Pv9eOmll4S/jzrqKEvq29XVJXw+55xzcMQRR8ge9+c//xkjIyO45pprFK8lHXTT09Oj2qZ5eXnYs2cPMjIyMH/+fPT19Ymsc9nZ2aLnh0f6cjVz5kzMnz8fr732mqE6vPPOO8J9tmTJEhQVFQltsnLlSkWxdPzxx8u2ndSydcoppxhq4+zsbEEgH3300Zrn8vvXrVsnfJ4zZ45i2+nB7/fD7XbD5XJFPK5tPEHtLEZqMVbCcJ4lRk5ODmbPnh3q6VGhoqICdXV1IbvfGGlpabIzYOx2u+U3Hx/fMTo6KoolsNvtmi6AiRMnWl4nPdez2WzYs2cPgLE34m9/+9ua5/AxQlbUmR8EMzIydF2Tb9/R0VFT9QgEAqbjcwKBgGod+O+YlJRkSbvxv4NS+X19fUIW6fLyct3Wp3A8I3r4wx/+gN///vdYvHgxrr76agwPD8NutwcJk4GBAcXOU/oGmpqaGtJ3Ydbha6+9FmeccYZon5oVNjk5WbY8u92O5ORkoS9IT083VC/eDedwOGC321Vn+PHX5vsfq+6/aN0j4w1q5zH0toGhAO94pra2FiUlJaipqYl2VQzBD1x+vz8oPkFOLPFm+Wi5kXiLl94pxVb70PnBTW8QLV9XM2vV3XbbbZg5c2bYZ1yGY7o2PyNMaX083gKnFggv3Retzvm6667DG2+8gRkzZgCAbLwSMCaW9u3bJ3sNpaSURoPs2X150kknBe0LNdibFzdWzIbjt6lBAd7EeGFciKWmpibk5eUFCSWWiyiW4QebgYGBIEuFXAfFD2Rm4k0AYPfu3bLBrlrwgzg/mKilNLBaLPGDm94BxCqx9NBDD+Hzzz/HmjVrQr6GHqIllvi2VRskpfui/SbLymf3Gh+vBIx9d6NiyehvwMSS3D0ZDbHEW1ONpg6gpJTEeCHhxVJ7ezt8Pp+sRUnaUcYifGd08OBBXWKJDwI3O7V88eLF+Otf/2r4PL7ebDB56aWXMGnSJPz2t7/VPMcK+Kneejtyq8SS3PVCQev3M3t9ufP1iCX+2YkHy5K0fKlliYkDNbGklDrA6G/A2ldO1ISac8kqyxKLLVQTS0ozckksEYlMyDFL8YDX60V1dTXKyspQW1sL4FAnz/bFOvxg09/fr0ss8bOZzHZg/GwXI/DlssFk06ZNCAQC2Lhxo+HrvfrqqwgEArjwwgt1n8NbAvQKMavFUrjyIDH4+hodtD/44AOceeaZqK2txdKlS4XtesQS715MBMvStGnT8Nlnn8mKJTadX2m5k1DdcLFiWeLFEpvNqyaWhoaGhPJILBHjhYQWSy6XC16vV9Hd1tDQEOEaGUcqlviYpUAgoCmWzFprjAzA/KAhZ1li4kNpto8S/f39gkjav38/srKyZI97+eWXMTQ0JLgNx5tYMjpY/eQnP4HP50NdXZ1hscTH+ii17e9//3v85je/EW2LFbHEflsmlg4//HBBLEnjzCZMmICDBw8GuZCZyDjhhBNky7ruuuvw4osvBl1PTSypWZbU7iWr3HB6xNLBgwdlxRItd0IkMgktlvgp0PGKmmVpZGREdoDk34zNiiWtATgQCKC9vR0nnngiDj/8cNlyP/roI5xyyikhiyVe9AwMDAhiaceOHXj66aexZMkSZGZmClOa9+7di7y8vJDccPxx8SCW+PoaFUtKg6pRsaRUrlwG8GiLJfad2SDP7pHJkycDkHfDpaWl4eDBg0Gz4ZhYWrJkCfbt24eUlBTcfffdAMasUY8//ji2b9+OY489VnReqJalo446SnHf4YcfLuSBCrcb7uDBg0IaBYpZIsYLpmOWNm/eLPxjbNiwAUuWLEFlZSWeffZZs0WMa9TE0vDwsKwY4qc+h9uy9Prrr2PevHlwOp2isvjPt9xyC6ZMmSIrln7/+98LA4wSSt/h7LPPxrJly7BkyRKRsGFv8mYtS0pCwQhmxZKRmCWjg5XSwGzUDWfkHjO7tIdZpG44di8ysTQ8PCxk12Yo1Zlf7uQXv/iFKEFvamoqbDabrDhkbarXsjRr1izcc889qrnheOuWUUFq1A3Hi0ZywxHjBdNiafXq1aitrRUSVL722mtwuVzo7u5GVVUV3nvvPRJMIdLX14fPP/9c+FtOLMl1UFaKJa0OkMUfDQ0NYWBgQPU8Jmj4zvbmm2/GL37xC9WMwbwQ4q/L2oa53xjsOyuJpUAggOXLlwtrf+3evVsQJfEc4G10sFLKpWOVZUmOaC2qzVASS1OmTBGO4Z85QFssyR3HhJDaFHy9lqW5c+fi7rvvDkobwnPiiSeqXlcNXhjptSw98cQTWLx4sej+ILFEJDKm3XAFBQWihXJra2vhdDqFzK6lpaVYtWqV2WLGJS6XS7R0SX9/PyZNmiT8rSSWQrGoKKE12E+fPl34/OGHH6qWyzpWNkDxQkAtHxHfIctdl1/niz+Gd8Px573yyiu45557AAAPP/wwbr31Vvz0pz/FL3/5y3EVs2RGLIVqWYoVsTQ6OoqRkRHhXszLyxOOkYolJfGhJpZYOWpWHr2WJT0pNXjLUrgDvA8ePCgkJOUhsUQkMpanDvB4PKaWEiEOIe30pKkD9FiWzHZgcmLJ6/XizTffBCDuyDs6OkR1k8IGJvY/f66aqOAtVkuWLEFra6tov3QBY3a8kmjcvXu38PnWW28FANx3330AxleAtxk3XLxbloCx+4bdi1lZWYJ4lP5mSu0knZkqZ1lSW+ZHTkjJlWVULBklFLEkB4klIpExLZZ4S8drr70Gm80Gl8slOibcA0aiIu1M+/v7g2JqeGsOg++0wmFZKigowLnnnoutW7eKhMxtt90mfJbr4NnMIzZA8QOx2oDAl/HCCy+goqJCtH94eFhWLCm5CNTcGYkklvbt24d7771XcXFTJRHAt7eZ2XByRFss8d95YGBAuBczMzMVZ6Lx50ydOlX4LP1t5WakqS2sbaVlibfwKs0WVSIUN5wcZlzOu3fvxmeffYY///nPtMArEZOYFkvbt28XPjc0NMDhcIgCHXfu3CkSVIR+pJ1pf3+/aECsrKzExRdfrHqNcLrhNm/erDiYym3nxVIgEBAdw3fA7I197dq1KCoqwo4dOzTryQsbdi0l0ThexNL111+Pn//85zj77LNlz+UHd35AtiLPktL3jrZYSktLE56r/fv3GxZL/Iw06bPBH8esNUlJSYrWJb0xS3rEks1mw4YNG/D444/juOOO0zyeh/+tmLjjnxHpAsJKq7SHalkaGhrC9OnTUVBQgHXr1uGiiy4K6ToEEU5Mi6Xq6moUFxfjmGOOQXt7u7C8w2uvvYbFixejoKAgLjJlxyJalqV33nlH8xpmxJJSagK+frwVQgtmjRgZGcHQ0JBoIOZdh2ygWbRoETweD3784x8HXUs6GMsN9vx3//vf/44LL7wQ+/fvTyixxF9f+lu1t7cDAD777DPZc/mBmR8A9YglrfQUSiI72mIJOLReYl9fn2GxNGvWLOGzmruO36dk6bHSsgQA559/vmwskRb89Zmw43+/u+66S3S81W446TIyBBGLmBZL+fn56OjoQHNzM3p6enDFFVcI+yoqKtDW1obZs2ebLWZcomVZ0kOoYqm3txdHHnkkvvjiC8VjUlJSQhJLwJh1iT+XT/gnFTN8kk2GdBDXEkvAWBbwBx54QNXFYDZ1QCAQEAUIx3LqAD5WxYhY2rt3r+ZyJ0p1iTWxxCYBqIkl/jmcNm2a8FlNLPG/ixGxJGdZsnoZICn8b8zqxJe5ZMkSvP322zjvvPMAWG9ZksZ+EYfYu3cv/vrXv4b9HiC0sSwppVQQlZaWWnXpcYucZSlSYunpp58WBULLkZKSInS0GRkZip0ogx9g+/v7RZ00v29wcBB/+tOfhL/lrBTSbMpyYkmurXp7e1U7Z7NJKRcvXizKGB9NsWTkXCNiSRonJ9fO8SKW+ABvPZYldi4QfF/y4jNUsWTGshQq/G/Mng2+35gwYQLOPPNMIb0CbwXmCVUsUeZvZU4//XR4vV6sWrUKd9xxR7SrM66hpJQxjJxlyWjHEmoHpjexHbMO5eTkaB7Ld8rffPONolgCgGuvvVb1WlKxJBezJCcUk5KSVEWQWTecdGmdaMYsaZXNH29ELH3wwQeiv+XaWUmkx0L8ohk3HC98pO3Li/BQ3XDS+CAgsmKJcdJJJ+GUU05BWVmZYOll7aMUVhFqX0NWE2VY/sKWlpYo14SgpJQxjBWWpZGRkZAGbD2Zlv1+vyCW+DduPWiJJS2kQaxys+Hk2ircYkmKXNv/+c9/xsyZM0WpFkLFzHInoYqlRLQsZWZmKt7z/HOoJpZ4QrUssUziPOEWS3L3eXJyMjZt2oS2tjZhG4tnUrI49/b24tFHH8WePXsMlU9iSRtKyxB9TIulgoICvPrqq0KsEp+UsrS0FPX19YKQIowh7UxHR0dFiRb1EsqDppSwkGdwcFAYTPVYlnjMiiUp4bYshdpZyQ2oV199NXbt2oWFCxfqPkcJq9xw/BI0VliWlNrLqKgOB6wOPp9P+K6ZmZmKyRz5GDde+KhZefl9RmbD8ckxGdGwLAFjzwpvLZszZw4ACAmHpfz4xz/GLbfcgrKyMkPly90/lG5GDIml6ENJKWMYOVeYWqZrJUJ5c9OTBZhf4iQUscQHeJsVS3wchdVi6bXXXkNOTg4ef/xxw/VS6/SV6mGkYwxFLD311FN47rnnQrYsSWfXGbEsxUIwLxNLO3fuFLZlZmYGPW+pqak44YQTFMWSXsuSnsBxRjQsS3onMlx22WXIyMjQfGGTimktjLhxxysU1xV9KCllDCPXmUZKLOk5hxdLZt1w0pXejcLPtFML8FYSS/wyGAx2XFlZGb755ht8//vfN1wvfs25F198UTS7UEk4GPm9jM6G27t3L773ve/hiiuuEIlVI2KJ1Y+JCL2WJaV8T5GG3asPPfQQgLHfIT09XfS8paeno6enB++//74ocDsUsaRkpZV7vuViuiI5G06NrKwsUQ49q1Bawog4BIml6GN6NhwlpQwfcpYlXhToJZTOVs/bLC+W5AJT1ejv7xfFiMilBzCCnFiS+942m012cFATS2ZgouHxxx/HTTfdJLLA6RFLRmfDjYyM4KabbsK3vvUt2XN5UcRb44yIJfadUlNTcfDgQVlhJG37J598EldeeaXqd4kUUmGfnp4Om80mEi92u1126Y/s7Gzhs9pvw7eJkliSe77l0lrcdNNNiuVYgZEUGUYtyHpQWkdSbamY8Qa54aIPJaWMYeTePENpy1AeND1iycqYpS+//NJYBSXwFrdQArzlxFIoeZaksDq8+OKLAMSz+JTEkt4184BgsfTiiy/iD3/4A2699VbNc0MVS2xwY2JXj2XJ4XCoJgONJFKxxL47L154a5JZN5xc4HhKSooul+Trr78um5TVSoy8FOiZ+AGMvUQ//vjjuvoesixpQ2Ip+pi2LLGklJs2bYLT6RQNmhUVFUHreBH60Tt9X4tIWJZCEUv89zNrZubFUigxS6wuanmWtCykctdldTASq2PGDSdNqSCF/135Y/kAb/57qFmWmMVET8ySWiLQSKPkMpZalhi8yMvKysJxxx2Hjz76CJWVlYplaLnh9MQEAsC5556r6zgzGHkpUIq/knLMMccAGLuXqqqqVI+Vu39ILIkhN1z0oaSUMYzeDlUOu92OQCCA4eHhiIilUGKW9L6l6kGvGy6UmCUGS8qnBC84GKwOclYVPW44rd9OKpaUcv0wlMQSb1nij5FrK94Np1RH6QDIW2qijdK9ygsk/jPfpllZWdi0aRP27NkDp9OpWAbf9mbEUiS49NJL8cwzz+C0007TPNboM/vuu+9qiiWyLGlDYin6WGIX37x5My644ALk5eVh7dq1wvbFixdjw4YNVhQxLjFjWcrIyBAGKKUBV82NoFcsWeWGM4sVAd6jo6Oi2KmhoSHRsfyK83LIzRLq6enBgQMHZIWRHrGklSfLaIC3llgKBAKi66i54ZgI4Ou7a9cu/OY3vwmKrYtnyxLfBizTt5pQArTdcLEkltauXYsHH3wQr7zyiuaxei1LDD3fk8SSNuSGiz6mxdKmTZswd+5c5OTkoL6+XrRv9erV6OnpEWX3JvRjpkPNzMwMmq3k9/uFAX3Hjh2YOnUqfvGLX8ierzdmyYwbzkqxxM8yU7MsAcpJ+C6++OKghUh37dolfNb6jnILgra0tCA7O1v2zVBPzBKgbl1SsyzJwX93ObEkLduoG660tBR33HFHkDUhlixLLHCbcfXVVwMQP298ffkYLr2WFTNuOKvc73pxOBy47bbbcNhhh2kea9SyRGLJGkgsRR/TYqm+vh6dnZ1obm7GokWLgvYvWLBAWP2cMIbRTpN3hfKWpe9///vYs2cPSkpKkJOTgxNPPBFOpxN79+7F3XffLZzDd1rhjlmSrg1nFj5njpplaWRkRFEAvPzyy0Hbd+zYIXzWMoWr5Z+RS/mgN3WA2m9hVCxpWZZCEUt8fT/++GMAwPvvvy86J5YsS/yAzwKRAWU3HO+i1Bt7ZsaydN111wEIDm2IBYxalvQktyWxpI1RN1wgEMATTzxhOOcVoYxpsZSfn4/8/Hwr6kJIMGJZmjt3ruhtPjMzUxBLb731Fm666SZs2bIFo6OjQctVAMC9994Lh8MhDHKRCPAOVSzJBdby8UKsTkrJ7tQCsaXwIkzr7U5NLMm1p16xZKVlia8HX1/WJtJ6DgwM4O233xbNwpS64fS89caSZem0007DnXfeidWrV6OgoEAQRkpuON6ypAV7BpcvXy5sk3P7qT3bv/vd79DY2Ij169frLjdSkGUpOiiJpW3btmHx4sX49NNPRdtfeOEFXH/99TjppJMiUb1xgekeTJpxVi6+wmzCwfEK32FPmDBBlERQSnJysujtnRdLAPD222+rlvXzn/8cAHDXXXfhr3/9a0RiltS+jxpyWY551NxwRsUSn61aK9haTSypDbj79+/H6Oio0IbStg+XZYmHd9Xy7Nq1C2effTaOOeYYwWqkZllSIpYsSzabDStXrgzaridmSYtHHnkEt912m2iQKi8vR1NTE958803ZsqRkZmZqBkVHi0hZlqy0OicCSi8kc+bMQU9PD9577z14PB5hu9SyS5jHtGVp+/bt2LJli/C3tLNetWqV2SLGLXyHKrdmFI9ULPFuOEA+nkYNozFLRjtRM5YlacyJFDU3nJJYUhI6e/fuFT5rWVDU2ph35TBsNhtGRkYwceJEOBwOoT3CaVlSemNnv7dSWZ988omoHODQ/RlvliUllPIsGbEsJScn4+STTxb9DmlpaXjjjTfwq1/9StgWSwHeRjBqWeLFUiAQgMfjCZo1SpYlbZSeMTaRYtOmTaLtvDWTVtCwBtM9WH19PZxOJ1wuF0pKStDV1YXc3Fx4vV40NjbC4XBg48aNVtR13MF33pMnT8bnn3+ueKyWZUkNfpBnnaEesTQwMCAcZ7QTNSOWtGK5QrEsyYkZQGwVNWNZUhJL/EC8Z88ezJo1K2IxSzxKliUp/Gw5taSUUmLJsqSEFW44NfjnMV7FkpnZcM888wwqKipwzjnn4I033hC2k1jSxmjMEi+W9u/fHxMLWMc7pi1LDocDHR0d6O7uRk1NDRobG1FeXo6amhosWLCAhJIJjFqW+M6Ynw2nBR+XoxToK4ecyNKLGbGkJQJDsSwpoWZZGh0dFXVioYgl/prsWmYsS1oo/a5su9rv/vHHH8su45GIliWlAG8zJIJYMvqc8+L9d7/7HQCI3JFAZMRSR0cH/vu//ztuV5MwKpb4e83sUlLEGJb0YE6nE263G729vejo6EBubm5MzuSIN/gOWyt7tJYbTg1+xhebgq80aJ522mlCKgh+uQw9sQk8VoqliRMniuqi5lIaHh42VC5vWeJFwcjICAoLC5Geno53330XNpstJLHEt/PIyAi8Xi/OOuss0XFGLEtawkVpENJjWTruuOPw05/+VPib/eZffPEFzj33XNx0002KsXXxbFm68847ce2115pejYC/byOdHsAqjIol/n5UWgQ8Ehm8S0pKAIzNAGWzH+MJo6kD+Pbbu3cvCgoKrK7SuMPS172cnBzZzN3PPvssrrjiCiuLGhfwnbfWopJm3HByYknJmvH222/jt7/9Le666y5BoCQlJRm2HPT394cc4C0ta9q0aSKxNDIygtHRUUNuOCWU3HAffvihEER58OBBZGRkhCSW+GsODQ3hRz/6UdBxRsSSlktMrxvObrfLHnvfffcJn5lYevTRRwGMWQzy8vJkf9d4sCwp5VmqrKzEN998E5SDyyiJYFky6obj70clsRRJN9w///nPsFw33JgRS7FuWers7MTHH3+Mq666KtpVUSUiK1s2NjZGopiEg3/71OqkzFiWeDfcl19+idHRUWGgvPHGG0XHZmZm4sQTTwRwyLI0YcIEQ2ufAcqWJT2dsfStXM7q5vf7ZTsYv98fslhSektmbj81d41ckKXUsjQ4OChKrskw4oYLVSxJ3XBaQfSA/ICvJOrjwbKk5IYDgMMPP9z0d0gEsRQOy1IkxZJVLtVIY9QNJ7UsxTLFxcW4+uqr8dZbb0W7KqpYIpZWrVqFkpIS2X/HHHMMJaUMEb5D1Rq8zFiWvvzyS+Hz8PAwuru7hUHzuOOOU6wXE0tGXXCA2LLECyQ9Ykn6veTSFgwNDVliWeJn7vDX++qrr4TP7HsYtZRJxdLAwIDsgBxJyxL7X49YkvvdlQZTo2I6Gii54ayCjzuMVzecGcuSlguYR+8zum7dOkybNg3vvPOOruP7+/tRXl6OM888M66yYmuJJenzFU9iifHvf/872lVQxbRtfOnSpWhqakJxcbHsekn79u1TfKMg1OE71HCKJWnQ4xdffCFyx0hhgyRzOxntQKXl/uxnP8OyZctw6qmnYvfu3Zrn6RFLSpYlqVjiY7C04K/Hz0xklqVQ3IpSy5Lcb2bEshRqzFIoliU5saTUqcfDwKRmWbKCefPmCZ87Ozstv34kMGNZUsJMnqUrr7wSAHDzzTfLJtuVsnfvXjzzzDMAxp7hI488Ulc50SaR3XDxgmmx5PV6NcXQwoULzRYzLjFjWcrIyJBd6V4OftkLAOjq6lIVS1IXQmFhoa5ypDAX1+zZs/HJJ59g4sSJOOWUUzTPs8qytH79epx55pm6E2ry1+NFXahiSc6yJCeWohGzlJqaiuTkZNVOWu7eULq+0UE2GijFLFlFVlYWJk2ahJ6eHs3ZrbGK0Rcjdv/w8Xy8FWTPnj2y7hejbrhQrNvRyD/U19eH7Oxsw5bWWHLDNTU14Te/+Q1efvnlcbV6h2k3nMvl0jymoaHBbDHjEn4w0hPgLU0doCf/DXAoluDYY48FAGzYsMGQWJo7d65mGXxcEetwWblpaWk4+uijMXXqVF0Cz6xlib21zpo1S7Ndefjr8WLJjBuO/40GBwdl3XB6LUt8rJkSemOW7Ha7alxNSkqKrKCQG+QefPBBXYu0Rptwu+EAYMuWLbj22mvR1NQUluuHG6Oil927e/bsEbYFAgFcfvnl2Lp1K6ZPn46nn3466DyjYmn69OmGjgcib+388MMPkZeXhwsvvFBXahYjYk7NDWe1Zam6uhofffQRfvjDH1p63VhPnmlaLOnJW8HPtiL0w3feRgO8MzMzdT2QwKGHiU2NbmtrMySW5GZASpkyZYqobsAhccG/FeoRS9I6hWpZSk1N1W19Y+cy5NxwRpMXWm1ZApQHGdYR6bUspaSkqAoG6f2mVNe5c+fitttuU7xOLBFuNxwAzJw5E0888QROPfXUsFw/3ITqhpNOXHj++efhcrkUB0g9Yom3iE+bNs1Qvfi6RYotW7ZgeHgYbW1tssvtSJE+20bCWfj2C1dQu9rs30TEtFiqqqrCqlWrRDOqpMTCbLiVK1eitrYWK1euRHV1NVpbW6NdJU2knbdaB56SkhLkhtMrlpg77PLLL0dycjK2bdsm/J5qMUvAmPA5+eSTNct44IEHAABLliwJcinyQtBKy5IesWSz2XS7XJQsS1a54QYHB2XN81aIJVZ3vTFLWveb1JKpVNd4mAXHiIRlKd6R+83V2ordd3LPhtzMT4YescSPOaH8XpESS0wQ8v0Rv44bT29vL7773e/iqaeeCqpfXl6eapvx8O2ndxwwilz7GXUXxhOmHfNVVVXw+Xyora2Fw+FAbm4uHA6HsN/n88Hr9ZotxhTV1dUoKCgQuQNdLhe6u7tjdsFKQNx522w2pKamKt74cpYlo6bso446CjNnzsTOnTsFa6CWZWnWrFm6BM68efOwd+9e5Obm4vXXXxft48WSnsFVb8ySVuoAJvpSUlJ0uSzDLZYGBgZk3wL1uuEA5cDY4eFhpKSk6J4NZ5UbLp7EktLacIQ66enpuu8rvRgVS6GkGjBap1C45pprsHnzZvzzn/8U1VHpOV2xYgVefPFFvPjii7jsssuC9v/tb3/TlY+ILytcaRik7dfZ2YnS0lLce++9IVmTY90NZ7pHcLvdKC4uxoIFC5Cbmxu0P9qz4TweD5qamoJ+iIaGBhQVFcW0WOI770AggNTU1KBFKBlm3HCMnJwcTJo0CTt37hQGbS2xpDWb5KijjsKvf/1rpKSkCEGt0jgh3rSv17J0xx134P7778dtt91myLLECxr2Pex2uy6hw64XCARES71YFbM0MDAg+/uy3/Gpp55CXl4eLrjgAmGfVBAqdcLsOCMxS1a44eJJdJBlKTQyMjJESWF52H0XDrHEh3doHZ+enh7kJo+EZenJJ58EADzxxBOivk2pvnx8kVz9lNZ4U4tZipRYuvnmm9Hb24sf/vCHusVSPFmiTPdkTqcTbW1tqsdEczZcY2Oj7Gwttq21tRXl5eWRrpYu+M6biSUlkpKSQnbDAWOdid1uDxK8WmLpqKOOEj7LTcO/8cYbg5aJsEIs1dfX48orr8Rpp52Gt99+O+gYpQBv3nLDvodRN5xUkPT39+Pvf/+7KPeSHuTccHJxAMPDw9i9eze+973vARj7bqzO0pcANcsSO1fP/lDdcNIONJ4sSySWQkMtnlLvAs1S9AzwvLDQOj4jIyMqYomxefNmUbiC0nOqtdZjdna2rvIiIZbk1so0c41YtyyZjllas2aN5jHRnA3X3t4um/8JGFsE2O12R7hG+uE77NHRUVWxFAgETLnhmOtUmg1bK2ZpxowZwueXX34Zv/nNb3DGGWcI2+QGVKlYCiVmKSUlBcXFxUhJSTEU4M1bboyKJXY9aad700034cwzzxRl+9aL1A2nZFniRR7vfghXzFJKSoqmG06PEIonsURuuNDgX3Yuuugi0b5QLUt68izpcWsx5Ga9RlIsvf/++7oEjJZYUhLxapYlOaH62WefYdu2beqV1sAKNyZ/jYQXS0oL5u7cuVMwzUYzF4PX65V1DwJAbm4uOjo6Ilwj/fADjZZlaXR0VNOylJ2dLRIyPExwGLUsHX744aLPP/jBD0TXkBss1SxLegZXPbPhlFxiTHTYbDahLKOWJatmlyQlJQVZluTE0scffyzqOD/66CPhs5GYJSD8liW54+IFsiyFBj9h49vf/rZoXzjdcPy9rHW8XB6mcMcs8c/s+++/r6u+/PMsVz+9Ak9NmH355Ze47bbbcMopp5gK/g6l/QKBAF5++WUh91M8JKtlmH59qqurw4oVK0Tbent70dXVhc7OTuzYsQOFhYW6cvFYjVZaA4fDoSv1gZTBwUHRoMREod/vD9vMA971IsfIyIhImcsFg2dlZeH111/HqlWrcNddd4n2TZw4EX6/P8gnbrPZUFxcjI6ODkyZMgV+v1/0FsO28fXkOyapqwkQiyObzSY6RilZW3Z2tihOiL+mXMJO/lge9lulp6cLD7sRy5Lf71eMzwgFXtR98MEHsiLvjTfeEAV7fvjhh5g3bx78fn+QWFJKX3Dw4EH4/X5FMcXuXVa+lhjSK4Lkfv94gM9ZJf2fEMM/zw6HA7/97W+FBaGHhobg9/sNp9UYGBjQbG/+WRkcHDQ0a1TPOWbh69ff3y9yGyp9P2kMoxTWnlLkXPpK5/Azwfft26crQeq2bdswa9Ys0W8tHe/48UepXW+//XY89NBD+P73v4/GxkbRfTEyMhKVZ0xvmabFUldXV9C2nJwcUe6dVatWRUUshYsVK1Zg+fLlQdvb2tp0LRMRCv/6179UOxyv1yusj5ScnIz29vagN4rh4WG88sor+Pjjj4PO9/v9WL9+fVDcTUdHBxYvXoxnn30WF110EdavXy96KHbt2oX169eLzuHf0D/66KOg/XxGWbvdjpdffln4WymA/aSTTsI//vEPoU58W8i94bz77ruy12G5Wex2u1Avubc8Pnv1hAkTMDAwgJGREaxfv141TYYRent7RZbN559/Xva4d999V+Qubm9vxzHHHAMgeBD4z3/+I3sNt9uNKVOm4NNPP5Xdz37/999/H8BYPIjSbwGMdcZ61nL68ssvg37/WIW/p+Tu21h22UeShx9+GG63W7hf+ftk586dOOuss3DjjTfiscceE/qHTZs2GSrD4/HgoosuwmWXXaaYcJIfe7744gvV+0zu5emdd94JWr3ASqQW6A8++ED47PP5ZOvL52+Tu9/+/ve/y36XQCAguh4/W/fAgQOifY8//rjw+YUXXsDUqVNVv8eWLVtw991346STTsJ9990nbO/r6xNdl6+X3Hc7ePAgHnroIQDA//7v/+LSSy8VvXhu3bo1Kn2FXk+BabGkJ2272+3GT37yE7NFGYZPYSBHKFYlYMyadvvttwt/9/X1YebMmZg3b57ibAWznHjiiaqLRR555JG46KKLcNttt2Hy5MmYP39+kIjIzs7G/Pnz0dPTgyeeeEK0r6CgAPPnz8eePXtE+84++2yUlJTguuuuEx3f2NiIvr4+LFq0SLTd7/fj2WefFf4+5ZRTMH/+fNExb775piCQsrKyRPt//vOfB323Z555Bnv37hXE0llnnRVk7pdy/PHHq+7Pzc0VypUTuJmZmcKDPGnSJOzZsweBQAAXXnihZa7bnJwcnHTSSYr7jz/+eHR3d+Orr74SWXIGBgYwf/58+P1+PPzww6JzpDFnjHPOOQf5+fn4wx/+ILt/dHQU3/nOd4T1tVjgPu/y48nOzlZ0wfMceeSRQb9/rMKL5mOOOUaot9/vh9vthsvlIvfc/3PuuecKYmnmzJnCxI7S0lKUlpYKLxRTp07F/PnzRSJAD7t378bu3bvx73//W/HlhLeQZGZmqt5nciEMJSUlOP/88w3VywjSGEY+Ma/dbpet75/+9Cfh87nnnhu0v6SkRNbwkJSUhPnz5+O6667Dp59+KurTbDabqCw2UYRdT60PAsaEDTC20C1/nbS0NNHfP/vZz4TPct/tlVdeET6fffbZmD9/vihv1PHHHx+VvkKvp8CwWFq6dCm8Xq+gyDs6OkRTmaV0dHTE7PT87u5uxeBvNdLS0mR94FpxHmZISkpSvbbNZkN+fj6amppw5JFHyh7L6nfNNdfgvffew6OPPirsmzRpEux2u+iBBg7NkpOi9pvyHVNqamrQ+bygnDBhgmpgbX5+Pq644grRIC89BxjLm/X+++9j1qxZQTlN5MjMzBSuIWeG5cXS0UcfLSzXII0zMkMgEFCcQVJbW4vFixejoqICX331lcjit3PnTuzatQtHHHFE0PlKdbPZbLDb7aoxAjabTbheWlqa6Hd87rnncMUVVwhWRbvdriubczifCavh7z255y2evku44QdjPgZx6tSpohxdo6OjsNvtIU8R3717t2Kb8y+Dfr9fV3JMKeH8PaUBy7zlcnBwULNsuckuSuPA0NAQXnrpJcWlY/hz+HoMDAxo1oNvO/7YkZER0d98feWuyZfL7gve2BIIBKLyfOkt07BYqq+vBzCm6quqqmCz2RSj2B0OB+rr64OsD5HE4XAo5nny+XwoLi6OcI1CQxrALbcfgGpbs8EgJSUFjzzyCN58803BlcIEjJ7ZcFrwQlJrNpx02rHSbDj+OnJ1evXVV+H3+3HppZcC0F56hK+DnLDi63XqqacKi32OjIxYFuA9OjqqGCTJnjPWlryJ+9NPP4XT6cTChQt1B3hr5VkCDsVkAcGz4Zg4Yu2aiAHeRhc3Hc/wzyD/rLDJHex3HxkZwf333x9kAbUCI9Pj2f3f3NyM+vp6eDyesAcXS+vE9xv8vmeffRZbt27F0NCQyFpmNMD78ssv16zH8PCwqM/Qs2SJ1qQQhtZsNr5PZnXirxGJJKFmCNkNV15ejsLCQixduhTNzc1W1slSFi5cqOo20bMQcDS54oor8Morr+Cqq67CU089pXicnjc36eDGD4ZZWVkA9M2G04K/rtZsOKl1QiqW2APIX0dukGYZzll92YOZlJQk2zZ8HeQ6A/7NmV/Ha3h42HCwqhKbN2/GLbfconoMax85U3Fzc3OQqV5tNhzv/pSDD9iUWlFSU1NFif2UklJKidcp+LE+jTna8L8rLzJZ/8H2f/HFF2ELwQhFLJ100kkiIRdOpHWSWpYYCxYskD3fquVEmEBKSkoKqpPSRBgepSBuNXHzxhtv4JxzzhHdG7xYZN9fah184403MG3aNCEmM5YwlTrA6XSisrLSqrqEhYqKCng8nqD4pPb2dgBAWVlZFGqln9bWVvT09GDy5Mmqg5Nc5y4dqKR/84MhEw+RtixpiaVf/vKXQddRG4CZUOMH9ZKSEtU6yIkl/m35tNNOEz6PjIxYJpb0wNpHqVMzkjpALgZCeoySWJK63ZSWO5EST5YlHhJL6vD3Bh/gzZImst/digDq119/XVYkKM3+koMJD94iGmmxxLeTntQIVq69xtpK2k5GLUvS2Ws8vDA677zzcMcdd4j2a1mW/vWvf+G8887Dscceiy+//DLmsnubzrOkpIp5+IDfSFNWVoby8vKg9AYNDQ1oaWnRDAKPNsxiAqiLBLkb6/XXXxcF4crFYDCstCzxYklusOStNlI3HP/Abd68WQhG1CuWpJallJQU/O1vf8Pvf/970XFaYol/CzrxxBOFz8PDw2FbxVsO1pZay0kwtJJSqiG1LEndcNI1/BItKSUPiSV1+H6BF/Ls+WXPqBXZo88//3xR4DMjFMsSf9/qcfuYuQ/U3HB+vx+BQED1+nLPbKgCz4xY4lMY8OJXq/0eeOAB0ffTEkudnZ3C58MPPzyq4TtymBZLemhsbIxEMYq0tLQgLy8PtbW1WLlyJaqrq1FdXR2zy5wooSdmiefMM88UrW6tZlliYikrK0skJMy64YxalviHi193zqhliXVMycnJyMzMxJlnnqlYBzmxdPbZZ2PmzJkoLy8XiQTeshQOIZCamipYPQHjliUlq5eegWF4eFi0kK6aZSkRY5Z4SCypw//2cn0P+92NrpeoREtLS9A2s2JJS3gMDAzgxBNPFGYCa61xumXLFhx77LFYt24dgOB+RfpsDg0NqVrE5J7ZUK0tu3btwtNPPx3Uj+hxw/GpXoyIJQCi5a+03HBSHnvsMc3rRxJLAgpWrVol3CBSfD4fvF6vFcWYoqamJtpVMI1RN5wUtZglJh5sNhuOOuooISdIpAO8+c5AKU5JrU5yliW5ejBxKC2TMXHiROzYsUOoA8u7xMcs5ebmihLNWcG6detEOcqMiiU+vwqPHtehnpglRiIud8JDYkkd/t5YsmQJ/vOf/+DKK68UtrHfXS1XV6jlMYxk8A5FLL300kvYtm0btm3bhksvvRTl5eVYtmyZKNcQz0033YRPPvkEV155JSorK1UtS6zOahMurHTDnX766bKpXvRYlvSKJblnxuPxCN4NLctSrCd9NS2Wli5diqamJhQXF8tOw9+3b5+mIif0ofYmr6dz1+OGA8byplgllowGePPw52oFeDOkMUtKYklurSge6cLETCzxs+EmTZpkuViSukG13HB614bTk1NMz2w4hl7LEgV4Jyb87zp58uSgBKVsfzjFktLacDt27MCECRNwxBFHCNtCiVniB+8f/OAHAIBf/epXimJJaiXSEkuDg4OC5c1utyMQCIjEA6uf0+nEzJkz8cYbb4TshmP9x1/+8hfRdi2xNDg4KHpR4/sRaV3krGT895ELcOf3h5r3MFKY7sm8Xq+mGFq4cKHZYghYb1mSC/AGxsSS3DF6MeKGk1qW+O/Bf18+8NtIzJLS+m9aYkna1ikpKRgaGsK7774rDABqq60DwKxZs/DLX/4SN9xwg+5psVKxpGVZYh1WVlaWasenpyNSsyyNt5glfgYkEQwfWyj3PFo940wuqaR0sdhAIACfzye8tPN9SSgxS0r5hRh33303PvroIzz11FNISkoK6lP0WJZYP8WWX5ITS7wV12zQ84wZM0Q527TccLxVCRBblqRjjpzLlReccqkT+DaWq8uTTz6J5cuX47nnntNMnhluTMcsyc00ktLQ0GC2GALGY5YYTGhIZ/6pWZbkjtELb1mSy5ukN2ZJ6ftGyrLEw+pSWVmJBx54IOgYubqed955+N73vmdIMEjXaWJtqRRUzjrXgoIC1esatSyN19lwHo8Hq1ev1jVxZbxzww034Pzzz5fNVWf1765lWQLGBma5Rab5QGojbjgtsfSLX/wC69atwxtvvAEgeCUAufpJ97PnOiMjI0gQ8gKP9TVmxZL0ZUzuBWvlypU466yzcODAAVWxJEVOLClZluTccHJcc801+OSTT3DTTTepHhcJImIj37FjB/Lz8yNRVEITqlj6+OOP8be//Q3XX3+9aLuSZYkPrA6l0+MfermHgS9LOuCaFUtWWZakYok/n9XxpJNOEta8Sk1NDYoLYselpKRoTm1mSFM3aGXJ5sXSli1bFI/TM4VbbTacNGYpUS1Ls2fP1rWMC3FoGQw5rHa/yl1PTnxIXXPp6eki0ROqWOLp6+vDtm3bhL/Zi4iWWJIyODgoEkvSPpw927xY4uukx5uQlpYm6nukblE5a05tbS0AYO3atSgqKhLtU+tHtNxw/AsfS8ir1+IeyXQtSpi2LFVVVWHVqlWqi4tGezZcohCqG66goAA333yz7pil4447TvgcSlZjfpCV6zD4TkWPRQwQf79oWpZ4iouL8dxzz6Gjo0P2zVMuoaYWUnGkJZZY52nUsiS3XI9Vs+H4eynexBJhDUq/+6WXXiqKJdKLlhsOGBus5eKYpGJJb8yS0v7TTz8dp59+uvA3s6houeHk6s+74aTPJC+WWHv29/ejuLgYP/7xj3W5OPlnEQgWR2qu+4GBgSBBatSypJSjCRj7/nrFUiz0I4bk/wUXXCAbnxQIBFBbWwuHw4Hc3FxR7qJYmQ2XCKjdMHziRL3wNyr/oJ9xxhn47//+75A6NUAsTOQ6DL5TkIqlUALV5fbt2rVLVBcrLUuM9PR0XHbZZQDGcjGxhX4ZvGVJi3POOQerVq0K2i4nanhYB3X00UerHicVSxkZGcJgwoLXectSSkqKSCgbccNlZmYKnXC8BngT5lD63efMmYPnn3/e8EuYHjecdCq+klgKJWaJh7cqAWPP1nXXXSdaqkROaEiRWpakYoN9v5SUFKE/am5uRmdnJzo7O4NyB8qRlZUlWtBXSyzxbZKamhrUxkru/EAgYMgNBxgTS7HQjxiqwcaNGxVnvSmtsUaz4axDSSzV19fjxz/+seHr8Q+CdJXq3/72t4avp1UGf32GtEPSI5aU1o8DDr2BsllqSm446RuXVhlybc+7pZ566incfvvtuPPOO4NyOul5K7rkkktk4/+0LEusAzrqqKNU12mUdnL8cRMmTMA333wTFLPEYyTAmxeisfBGSESeUFzoasidJyeWeDeT3IwrI244vQP5fffdh88//1y0rbe3V5dliRdLUuHCv8ywOqtl0JZD6hrUyrPU09MjfE5JSdFtWdKzfpzcbMCEtSw5nU60tbUZLoRmw1mD3A0zZcoUwcdsFP4tLFw3o1aHEYplSe2tVDrIh9OyxHdE+fn5eO6550T7jViWlKxlesVSenq6yKIjRSqWDjvsMGFbWloavvnmG3z44YfC8y2tT2pqKrnhCN2EIpYeeughDA0N4fbbbw/aJycM5GKW+MFfy7KkJyklQ61fkgolYOx5M+qGk7oa2flKMUt6hIb0miyFQHZ2Nvbv3x/UX/CGjYMHDwa1sTTmia05p5R8NJEsS4ZiltasWRNSITQbzhr0BDkawYqlCJRgbqHvfOc7qsfptSzpzXsj7RxYxygVP7xYYjMteMuOUcuSGmYWnNVyw7G3tdTUVFVrGRNGJ5xwAtrb20UzHlkZS5YsEbZpWZZSUlIU60ZiiVC6n9Vc6CeccAImT54su0+un5OLWeIHfzmxxLuPrRJLcuixLD344IPCNH652XByYkkutYAa0muy89kaflKxJHXZaaU/YL+LklhSSh0AGBNLb7/9Nk499VRhQk00MCSW5GaJPPvss4oJ8xhsJpzWcYQ6cgNPrIqlTZs24auvvhINynLotSzxM/TUkHbGrM1sNpuoA+fF0urVq9HZ2Yl77rkn6DylvwFtsRQJyxJrP7vdrmotY+bzc889F6WlpaLy5ESPtD42my3IsjR9+nTZssgNR4RiWVJz7W7cuBHr168X/h4ZGRHufXZf6rEs8clma2tr8fTTTyvWx8wMLDWxxCzSzz33nLCyhFyANx+zJBdnpSQ0pLNY5eDFEt/n8palAwcOBI0vSq5CtQW8GdL2NOKGA4D3338/rGOWFqZmw+3YsQNPP/20ZgB3b28vSkpK4HA4kJycjPvvv99MseMWuY7EyM0mJZw3XlpaGqZMmaJ5nPQYJbF03HHHobm5GW+99Zbq9aSdA/82w7fVxIkThc8pKSkoLCwUnavHDae1CLOR2XBKYknLssTQsiyxWAR2DP/99IglALLLnfDtyODroRZfRiQuesWS9N5TOu+f//wnLrroInR0dAAQvySygV9LLPFLFzGuuuqqoD7nkUceQWVlpSg+x2h+IzU3HKsvj5xliY9ZknPDKQmUm266Cenp6SgrK9MUS6OjoyIRw1uWnnzyyaBlzKRiiX3HcLvhGHr7w3BgqifLz89HZWUlTjvtNPT19YkWzeNh69H09PRgZGQE//znP7FhwwYzRY9LrHbD6c39Ew6efvppXHTRRfjZz34m2q5m7q6oqMBZZ52lel3pIM8HLPJIkz8CyhnDpfsYerM8h9OyxJ8vJ5ZYvZkbTq9Ykquz1LIEjLlOpPCWJbNJ9Ij4RMsNd8YZZwAAbrnlFmHfyMiI5ovFK6+8AkD8osfu6aGhIVU3nJxYAsQCIRAI4NZbb0VzczOam5uF7UpJYZXo7e1V7JvlLMB63XC8WFKyfM2cORN79+7Fyy+/rNiv8H0FLzB5y9LevXvx0ksvATj0cikNCGdtrOWGCwQCQhvylkA5sfSrX/1KMdeZkviLBKZf+7q7u7FkyRJMmjQJRUVFSE5OxrJly0THtLa2oqGhATk5OQDGpj/KrSJNqCPXkZgZjKJp0qysrMSLL74YlITRLNKHSU4sfetb35I914hYuvjiizWtJkYsS0qDixGxJNcJFxYWiv5mb5RaYknu3pDOhgOA448/Pug4/s3ZjOWTiF+0LEvr16/Hs88+K1pnTY9YYhm6+fuT3ffSdcyUxJL0Wfvkk0+Ez11dXcJnflKE0TXu1Nxwci81WgHecm44JbGUnp6OjIyMoPUdefj+gheYvHDkYWO3kmVJyw3H72cW+cHBQdm4q2uvvRannHKK7PXi1rIEjN1QXV1dWL16Ndra2rBu3Tps375dyBmzY8cO2Gy2oNQCcukHCHWsjv+IplhSwuwCptI3KbmHWDq1n8F3olIhxD/UZ599Nh577DHNulgRsyTtHJRmAiq54aSZ89kx/HXkOqCRkZGgsqR5lgDg6quvVi3TqrXBiPhCSyxNmjQJl19+uUiADw8PIykpSXW2q1QsJScniywVoViWtm/fLnxmbj6+DMB4/FIobji1mCXWH6kleWRIF8CWw263C/XgFwNXSvPDxJJUNOp1w/F1ZddSsixlZ2fLthEQ52Kpq6sLbW1tWLRoEUpLS1FeXo7m5uagNWWksQ2hZIYe71gtlh588EEAwE9/+lNLr2sGs2JJzUzLrCyLFy+W3a9mWeI7qfb2dl3xWFbELEktS2pvilKx9K1vfSvounJiiS/juOOOQ1VVFebPnx/0W8hZlubNm4cXX3wRl156qbCPfxEiN9z4JJTZcPzCsUps3LgRixYtEp7H1NRUYQANJWYJEFuWOjs7Fcs2woEDBxTFktzYp9cNx4sSPQJO6WWF7y++/e1vC4HmWmJJKcBbyw3HXHApKSmCRUtJLGVlZSmKpbh2wyllDmY3hJ4FPAl9XHTRRQDGzM7PPfccJk2aJPiUQ6GsrAz79+/HL3/5S6uqaBqzYkkthmvDhg3YsWMHjj32WNn9esWS3sWFwzEbjsV6SJFaltxuN954442ggUErZumuu+5CY2OjbH3kYpaAsfuSXxF81qxZwmeyLI1PjMyGKy0txeGHH46zzz5b8RietWvXCjOr+TUMrbAsKcXdGuXAgQP44IMPZPfJxT/JueG0klLqEUtK7jFetADAr3/9awCHVj6QIjeRA9B2wz355JO48MILhVxU/Kw/udlwGRkZihNHgOhalkxnenrvvffwt7/9Deeff76wbcOGDYJIUlKq/A1K6OPcc8/FP/7xDxQUFGDy5Mm49NJLTVvotDJZxxtqaxfl5OQIb0hy6BVLRtvcyjxLt9xyCy655BKsWrUKe/bsEbZLY5ZmzZqFCRMmBF2XdUJKYumoo44SPh922GGic6Wz4Xh4UURiiVC6n+W2u91u0aBpt9s1hQCzBqWmpgoiI9SYJX42t1LMjlHWrl2ruE9OLOm1LOkRS/wLp5KI4d1wjC+++ALvvvuu7PFK/Sa7PqtLfn4+LrnkEiQlJeGBBx4AALz66qt48803AYhFoZxlidUpNzdXtry4dsPV19djwYIFSE5OxjHHHIPk5GRUVFQAGMvBVFtbi9LSUtG6V6tWrYLL5TJb9Ljk9NNPFxK3JaIr06xlyYz4U4tZCmXWYTgsSxkZGbj99tuDLExSy5JS5nImlqTrvjF4sXTTTTfh+uuvx1NPPRVUF6kA5IUb39FRgPf4ROkFQe4+t9lsou1Kz0tqair+67/+C8Ch2CXeDadXLEnrxlujjAZyh4Jcln21mCXesqQnZolHyT0m57Zvbm7G6Ogopk6dGnS8llhi7Xbsscfit7/9bdAsWVbXGTNm6BJLcrOVgTh3wzmdTuzYsQMrVqxAaWkpVq9eDa/Xi0cffRQ5OTloaWlBW1sbtm/fjpKSEpSUlMDtdmPBggVW1J9IMMyKpUsuuQQ/+MEP8N3vfhcAcP311+s+l+9EpR2qGbGkx7Kk9L2lHSjrLKQiStr5aYklXgzyb9PTpk0TlfX444/jqquuAiAfs8Tgl3xQW/uPGB+YWRtOzcrKJg/wYondlwcPHgzJDcdbX4ymCADGkkuqzVrlhWBqaqpscma12XB8gDePkljiy1OzLElnz957770A5GM6ldxi7LuwdmPXVPoNv/3tb6u64dQsS0rtECksWXAlJydHCBDjKS0tFT6vXr0ar732WtB2guAxK5aSk5PxP//zPwDGJh/wLiE95zKssCwxK42eAULp+koB3vz2pKQkJCUliTo/VqbU1C4nlniho1ZXudlwDKU3TwrwHp8YccNJUZsZyiYPbNu2TTiW3ZdtbW2isA+9Yom3vhi1LC1YsACXXXYZUlNTFa04mZmZQliKnFhKSUnBSSedhH/84x+i7XJJKXnkxNIpp5yC733ve0HXkJKSkhLUN+zduxennnoq7rjjDvziF78Q7UtLS0NaWlrQ9dhELtZuesQSs0IbtSxF0wUHWGBZ0gMLmistLSWhRKhiVizxFBQUGJpBaJUb7pVXXsENN9yAn//85wD0WZaU3FXStyk2kPDChQkofmFfVqb0DU1OLO3evVuzfoC6Zen+++/HhRdeGLTQNlmWxidmLEtK5/KWpY8//hiA2LIknezCBnb2bCnFLPECwKhYYgO42qQP/rmx2+2CW/vBBx9Ed3c3du/eDafTqSvAm0cqlrKzs7FlyxbRS5OaZUkuZOGyyy6TnYlmt9tF/QtDKpbYMUrt8e1vfzvkmKVxIZZWrFgRiWKIBMBKsWQUqyxLF1xwAf73f/9XECZ6BgglUWGz2USdBOtopB2wdBufz4Y/jl2L/35qQfE8ajFLs2bNwssvvyzEIrLYkoULF+q6NpFYGIlZksLfm/zabbxliSWb5cUSgyWdNWpZGh4eFlxfRpYZ4v+XgxcZqampuOqqq9Db24vbbrsNkyZNEtKQ6Anw5pGKJbkYViVrV0pKiqxYysrKkr2OXDsDh8SSlhvu6KOPRllZGfLz83W54eQsS9GMVwIMuuFWrVqFjo4O0Q3Mr9Quh8/n01w7jiAYsSqWzAQq67EsKWWsBcZECusY5dxwcttYZ8W/oU2cOFHoCPkOccmSJXj00UdRVVWlWke12XBSNm7ciK+//hozZsxQPY5ITMy44Xh4K0daWhpmz56N5ORkQQClpqYGuaqPPfZYvPfee4Zjlnir0pw5c/C3v/1Ns35s4DcilgD5GCBpn2M0ZklO5ChNeFGyLLH23rVrF0455RTBrWm321XFktSyJP2dt27dGtRWapYlubLiyrL06KOPBi1T0tXVJSx1Ivdv7ty5mguOEgQjVsWSGbTigP71r3/hyCOPVDzm8MMPFz6rWZbkxBJvWeI7aP77/eY3v0FbW5uQpFQJvrPSEoBpaWkklMYxZtxw/KDP3+dpaWmYPn06KisrhW1ygzi756ViiZUtrdvQ0BACgYAw4CclJSkuiSRFj1ji+zQ1y5pU7MjNhuORiiW5Puvpp5/GSSedhD/+8Y+i7SkpKbJuNSagpk+fLlp1Q68bTsmyxLcPb1mSWtSVklHy50ULQzLf4/EE5U0qLi4OilOQorSYKUHEEvwDbmVaBjVhMWPGDJx88smq5xcXF+PDDz8EcKjT4eMSzIqlCRMm6ErlYbPZMGHCBAwMDFieTZ5ILMy44XhxwQ/QbLC89dZbhbif5ORk3WJJKWaJHcsP+HySVTX0xCzxbm45wcGQih2jAd5yfdZpp52GrVu3YnBwENdee62w3W63y7r+eWsTL3CMuuH49rDb7aK66bEsyRFtN5yh1+ecnJygtaYaGxs1z2toaDBWK2LcEiuWJSvFEt85v/DCC8JUfL3l8K5u1gnJdWpyMUW8G47viEL9fqwMo+4UYnyhdH8ZvW+kliVAvBD2hx9+GHTfs/tcrxuOHcsP+Mcdd5yu+umxLPGWYbnFrhnhsCwxpGLObrfLBn9L3Z788XJC76uvvsLChQvR3NwMQN4NJ22bUMVStC1Lpn0NHo9HZBaVQyqwCEKJWBFLVtaDv+6JJ54ovBUD+tx9/OAgJ5ZYx1RQUBB0XWmAt5Fy5WCDF1mWiFAw6oaTc9/w1/jqq69EgopP7jg4OIi+vj5dYmlgYEBkWSoqKsI555yjWVctsfSDH/xAtNi0mlgKR8wSf21p4k+5Po7vV6RiSc6ytGPHDlFojpwbTipy1AK81ZIKx71YamxsRGdnp2yiLYIwSjTFEv+AW1kP/rrSjlqPaDn11FOFz+zNi+90WUc9adIkbNu2DTt27Ag6HhDPuCPLEhEN9K6rKHc8P1jecMMNAIDZs2criqVXXnkFOTk5WLp0KQBtyxIfpJycnIw33ngDd955p2r9tGbD3X///aJ6h2pZMiuWpHV0OBy48cYbRVYvQN0Nx1uWlL6HnBtO2jbsmG+++WZ8WZZcLhe2b9+umOETgGipE4KIVfhO1MpkimqB43osNBMmTEBHRwfeeustoTPjOzW+YzruuONEiTj5DpTvmMiyREQDPSJbKSCaHywfeeQR/PrXv0Zra6vIDceLJbb0ydatWwEYi1liaMXJsOdI6XvZ7XbRPjXLiVTsGM2zpCWW+DQCeXl5mDJlCj7//HNRzKSaG44XpXl5ebJ9iB43HNMKvb29qmLp9ddfx/Tp0xWvE2lMi6WysjKsWrVK1bK0ceNGs8UQ44REdMOZtSwBQFFREc466yzhbyMdOoO3LIUqltjARGKJCAWjFkklsZSeno6f/OQncDqdipYlKUbccHLlqyEnVCZMmACbzSb6zqG44ayyLPHPP59vTUnM8f2KNGbJbrcLa5Ty6HHDsWz/fX19qmLp3HPPxZYtW4S/o23NNl16c3MzfD4f8vPz4XQ6kZubK0oV4PP50N7ebrYYYpwQzeUx1MTSU089he9///tCIGOo15V2eqG6w5QsS2rwHVOo5bKBKdodFxGf6LlXlRZ51rPYdKhiSY9lic/vxGD9hFIiR2m9jbjhWGb9adOmCbPOeP75z3+K/g71BYjvc5VilqSz4ZKSknDkkUfiq6++El1LLoO3mmWJj6kEgt1w/LnRXjjekpglt9uNoqIiTJo0CYFAAD09PcK/aFoKiPgjmvcL39lIRdtVV12FAwcO4JJLLjF8XbUFekPt4EIRS/ybYajlsgzKanmhCEIJM5YlJWumFZYl6Ww4Bj9Yyy3BwZAbyOVcdEbEEuOiiy7SZcnVKyakrkBeACpZ8qRuuOTkZGHtSx45y5JULOm1LEnPjbZYMv166HQ60dHRoXoMLXlAxAP8wygn2kK1pvDXDSVmSQ6+09W6xu9//3vce++9eOSRRxTroZdHH30Ud9xxh2rGcYLgsdvtwnJBZsSS0j1rVCzJ1UGPG06aKRxQtywx9IqlOXPmBG07/fTTMW3aNF3Pq95nmi2vwlCy5kvdcNKktHrFkvT34C1LvCsxIyMDRxxxhOhYK2bwWoXp0vXkUKqrqzNbDDFOiJVs71ZauPiONByWJa1162688Ubs2LFDFMgZ6ltaZmYmTj311Ki/5RHxAx/bYvR+12NZstoNJ7c8ifSzFLnngW3TK5Zmz56Nt956C8ccc4ywjeVY09Nuep9JvWJJ6objv39SUpJoIglDjxuOWZa+/vpruN1uAMB9992Hzs7OIMsS/72j3eeYFkulpaWax8yePdtsMcQ44emnn8app56Kv/zlL1GtR7jEklUxS3xHxgJBjRDttzRi/CAdnI2gJ6u+1LIklw8ICC3AW0ss6ekn9MYsAcBZZ50lyktoZPap3r5Emi5Aj1jSa1kyMhuOJz8/H8cff7xqveNSLO3cuRNr167Fhg0bgvYtWbJEmOp45ZVXRj3/ks/nQ21tLaqrq+FyuVBUVITW1tao1olQ5uSTT8bmzZtDig2ykli3LPHXDEUsqU1hJggruO666zBp0iTccsstIV9D7UWDwYujCRMmKGbfDiXAW81CAhwSB1ZYlqT1BA5Zzaxww917772YOnUqfv3rX4u2yy17AgS74aTB7nJiiX1XNTecXC4lo8lKo4HhIIy1a9eiurpa+NvlcuGVV14BMGYy7OzsFPY1NzfD4/Hg448/tqCqxmFCqaGhQXDveDweFBUVoby8PGhRYIJgREosWTEFPxSx9KMf/QgvvfQSxRMSYePxxx/H8PAwXnrpJUuupyQGpMudZGVlYcqUKfj6669Fx5nNsyQVS+eccw5uuukmAOEXS1ZYlu666y789Kc/xfDwsLDWJKAslqRuOD2WJYY0WzgP+40OHDggbNNjVIm2WDL0Wrtp0yZUVVVh9uzZuPPOO3HFFVegra0Ny5Ytw7PPPouenh643W6Mjo5idHQUzc3N2Lt3L+6///5w1V+VFStWiIQSABQWFqKhoQGtra2U0oBQJFJuOCvcYXJrPGnhcDiwceNGzQzFBBEqNpsNdrsd5557LgDoXm9NCT2WJfasnXHGGUHHsUGbFx5s28DAALZv3w5gbKo+Qyl1wbHHHos33ngjpixLesSE3DF6A7ylMUsOhwM1NTWydeO/s9x+aT30ZOeOK7HExEdHRwfq6+vR0tKC7du3o7m5GW63G52dnaIYpvLycrS3t+Ppp5+2vOJ6aG1tRVFRUdD2srIyACDLEqGIlfme1CxLVnQAWgHeBBFNJk2ahP379+Nf//qXqesoCQa5uKaHHnpI0YrLb2fu6N7eXiEBIr8Wo9Sy9IMf/AAA8Ktf/Up0bbXn2EjMkrR+VoslOb73ve8BGDMk8GjFLAFjE7z279+PM844Az/72c+E/VpiibcqLVu2DBUVFZr1jCux5PP5gt5EnU4n6uvrARyKcucpLCyMWu4cp9OJ7u7uoO3M0iS3jyCA+LIsheKGI4hIkpWVZXhdOClG4naOPPJIIakjQ04ssWDjf/7zn/D7/ZgyZYrIvcSLpRkzZuDBBx/EF198gQULFuiud6y44ZS455578PTTT6OtrU20XW02HF+fjIwMvPvuu/jFL34hbNOa8s/3r/fdd58uy1K0J6UYKl1pWrfL5RKS1cmhti+cuN1u9PT0BG33eDwADk3LJAgGm15/0UUXWXZNvhOTdmgklghCjNKgb9S6Ih2A5WKWWLDxm2++CWDMqsRf4/jjj0daWhpOPvlkrFixAjabDYcddpiuukXCDcfHHoXal6SlpaGyshJ5eXmi7Xy95dxwamhZlkIh2pYlQwHeSpXNyclR/SJqmU+jQWNjIxwOB6qqqkI6f3BwUBQnwoLT/H7/uHeJsO8fr+3w3nvv4cCBA3A4HJZ9B96lJ71mUlKS6XKGhobitr3jgXi/p+MF1r5Kz0sgEND8DfhjpGOSzWaD3+8XWTWYeNm3bx8AoLi4WFTGEUccgV27diE7O1v1WZWzRLPy+O+Tmpqq+R2ky71Ir8HD54QC9N2jeu9naZnSZZuMPA9qx+q9jp7fPxT0XtOQWFJzTajti7Yi5Glvb0d7eztaWlpCToC4YsUKLF++PGh7W1tb0M07XmHJxghgx44dwuf169eL9u3duzdom1EOHjxo+hqENnRPRwY+qzN/X+/cuVPzPv/000+FY6SD/eeff47169fD6/UK26STI9LT00N6lqQz74Cxl5j169cLgeMA8PbbbyvmgWJ88cUXwud//etfSE1NxdatW2WPfffdd4XP/f39huqudT8zDwwAbNiwAR999JHwd3d3t+6yvvjiC8Vjc3JydF9nz549Yenn2DI3WlhiWdLaF0uxQRUVFWhsbER5eXnI16irq8Ptt98u/N3X14eZM2di3rx5sgm3xhN+vx9utxsul8t0jEKi8M477wif58+fL9p32GGHBW3TywMPPIAf//jH+OMf/xjyNQht6J6ODKydeTExf/58ZGRkoL+/H7fccgvOPvts1WvMnTtX9CykpKQI64/l5+dj/vz5eP/994X9+fn5QmB3Tk4OfvjDH4a0rNGaNWuCtk2YMAHz588XCZrLL79c0y313HPPCZ/PPvtsnHfeeejt7ZU99vLLLxc+Z2dn6+oH9N7PfDt897vfFSWzPPzww3X3OTNmzAg69uWXX8ayZcvw6KOPBgWWK3HmmWeGpZ/TmwvS0F3R0tKC3t5eWSuS1+tVVKqhTtEvKCgwJLRyc3PhdrsVY6QqKipQV1cXsvuNkZaWJhuQZrfbqTP9f6gtDiH1/fM4nc6Q2+nWW2/F9OnTcdlll1FbRwC6pyMDLybsdjt27tyJrq4u2XQAjFdffRWvvvoqbrnllqDFYJlYYr8fL8b4SUlz587VtPoooRSAbbfbRakI9AQy8/VngfFKS61I44iM3J9a97N0XTzea5KcnKy7rJSUlKBjL7zwQlx44YW6zm9ra8MzzzyDn/70p2F5/nR/D6MXlkbM83R1dcluD9UNp3S9UKitrUVJSQlqamosuyZB6EHu/m9ra8Mf//hH3HfffaaurbZeFUEkAlOmTNFcMmXevHmYN29e0Pa0tDQh2aTcbDg+m7TL5Qq5jmpj3NFHH40//vGPsoHhcoQ6G87q2WJquZCMJNM1Wy+Xy2Xqt7EKQ9+isLBQSDhp5F+014ZrampCXl5ekFBqamqKUo2I8YRcR+pyufDEE0/EzMLBBBErWJHVniE33Z2/Pv9smhmQ2WLx11xzDU455RQAQGVlpbD/mmuu0X39SOdZUuKEE04Q/a2UOkCLaE/5twpDlqVQUwBEK3UAMOYC9Pl8shYln88X+QoR445YmuBAELHOE088gUsuuUTI32cGOWsIP3jzgdkFBQUhl3P66aejp6cHOTk52Lt3L1599VVDuZh4tCxLxx13HM4//3xcfPHFAMbyFP30pz/Fww8/HHL95ZgxYwY2btwovNDxbWlEAI1LsdTc3BxSIaGeZxav14vq6mqUlZWhtrYWwCGBxPYRRLghsUQQ+ikqKsKXX35pyXMjJ5b4mNv//u//xltvvYU77rjDdHlMVEyZMgXXXHNNyNfRsixlZ2fj0UcfFf5etmwZfvKTn4TFJV9cXCx8JstSAuNyueD1ehXdbQ0NDRGuETEeIbFEEMaw6pmRE0v8wrGnnHIKdu3aZUlZVsGLOTmxJBeQHInYxVDLILEUB1gZIE4QoUJiiSCig5xYkiaJjDX4JIlybrhozcjk25LNMNQDiSWCIHRBYokgogM/wLMUHgUFBTjssMPgcDgsDSa3Cl4ssfrzgiNaAo8vl8QSQRCWQ2KJIKIDP8Cztc/sdjv+85//IDk5OSafTV4sMaGh5YaLBHy5vCtTi2OOOSYc1Yk4JJYIIszEYodMEOMB3rI0depU4XMsut8YcmuVxYIbju/H9FiWNmzYgFdffRWLFy8OZ7UiBoklgggzJJYIIjrwYklvUshoIyeWYsENx6NHLJ1//vk4//zzI1CbyJAYzkSCiGGOOOKIaFeBIMYliSKWYsGyxGPEDZcokGWJIMLM9ddfj46ODpSWlka7KgQxruBnvsWzWIqFmCUeIwHeiQKJJYIIMykpKaIkcgRBRAZ+lQZ+0dxYZmhoKGgbiaXoQ244giAIIiHhxVK8xA5queHiJWYp0SCxRBAEQSQk8bj+Zzy44cZjzBKJJYIgCCIhIbEUHsiyRBAEQRAJwrHHHhvtKhhm2rRpQdvIDRd9SCwRBEEQCckf//hHfO9738OmTZuiXRXd/M///A8uvvhivPLKK8I2sixFH5oNRxAEQSQkTqcTf/rTn6JdDUNMnz4dL7zwgmhbrIklilkiCIIgCCKmiLWklOPRskRiiSAIgiBimHhc7iTRILFEEARBEDFMrFmWyA1HEARBEERMESsxS/PnzwcA3HLLLVGrQ7SgAG+CIAiCiGFiRSy1traio6MDc+bMiVodogWJJYIgCIKIYWIlz1J6ejrOPvvsqJUfTcgNRxAEQRAxTKxYlsYzJJYIgiAIIoYhsRR9SCwRBEEQRAwTa7PhxiMklgiCIAgihuEtSykpFGocDUgsEQRBEEQMw1uWeOFERA5qdYIgCIKIYXiBxAsnInKQWCIIgiCIGIbEUvQhsUQQBEEQMQwvkEgsRQcSSwRBEAQRw5BlKfqQWCIIgiCIGIbEUvQhsUQQBEEQMQy54aIPiSWCIAiCiGFsNpvwmcRSdCCxRBAEQRAxDC+Wpk2bFsWajF8oFShBEARBxDjvvPMODhw4gMMPPzzaVRmXkFgiCIIgiBhnzpw50a7CuIbccARBEARBECqMS7FUVFQU7SoQBEEQBBEnjDuxVF1dDY/HE+1qEARBEAQRJ4wrseTxeNDR0RHtahAEQRAEEUeMK7G0bt06VFZWRrsaBEEQBEHEEeNGLK1cuRJ1dXXRrgZBEARBEHHGuBBLHo8HTqcTDocj2lUhCIIgCCLOGBd5ltatW4eGhgbLrjc4OIjBwUHh776+PgCA3++H3++3rJx4hH3/8d4OkYDaOjJQO0cGaufIQO0sRm87JLxYCof7bcWKFVi+fHnQ9ra2NmRkZFhaVrzidrujXYVxA7V1ZKB2jgzUzpGB2nmM/v5+XccltFgKl/utrq4Ot99+u/B3X18fZs6ciXnz5mHixImWlhVv+P1+uN1uuFwu2O32aFcnoaG2jgzUzpGB2jkyUDuLYZ4hLWJaLBUUFKC7u1v38bm5uXC73XA6nQCsd78x0tLSkJaWFrTdbrfTzff/UFtEDmrryEDtHBmonSMDtfMYetsgpsVSV1dXyOe2trbC4/GgurpatJ3lWWLbGxoaKPCbIAiCIAhFYlosmaG8vBzl5eVB21kG78bGxijUiiAIgiCIeCNhxVIkCQQCAPT7PhMZv9+P/v5+9PX1kYk3zFBbRwZq58hA7RwZqJ3FsHGbjeNKjDuxZCQGSi/79+8HAMycOdPyaxMEQRAEEV7279+PnJwcxf22gJacShCamprgdrvR2toKACgsLERxcbEl7rjR0VHs3r0b2dnZsNlspq8Xz7CZgZ999tm4nxkYbqitIwO1c2Sgdo4M1M5iAoEA9u/fj2nTpiEpSTlP97gRS0Rk6OvrQ05ODnp7e+lBDDPU1pGB2jkyUDtHBmrn0BgXy50QBEEQBEGECoklgiAIgiAIFUgsEZaSlpaGu+++WzZpJ2Et1NaRgdo5MlA7RwZq59CgmCWCIAiCIAgVyLJEEARBEAShAoklgiAIgiAIFUgsEQRBEARBqEBiiSAIgiAIQgUSSwRBEARBECqQWCIIgiAIglCBxBJBEARBEIQKJJYIgiAIgiBUILFEEARBEAShAoklgiAIgiAIFVKiXQEiPvH5fFixYgV8Ph+8Xi+6u7tRV1eH8vJy2eNXrlyJffv2IS8vD11dXXC5XJYcO57w+XxYtGgRKisrVduD2toaqG1Ch+7V8EN9cIQJEIRBenp6AlVVVYGenh5hW2dnZwBAoLy8POj4qqqqQENDg2hbWVlZoLGx0dSx44Xy8vJAVVVVoKqqKgAg0NLSongstbU1UNuEBt2rkYH64MhDYokwTE1NjeghZTQ0NAQABNxut7CNPcBS5LYbOXY80tXVpToAUVtbA7WNeeheDS/UB0ceilkiDNPa2oqioqKg7WVlZQCAlpYWYVtjYyMKCwuDjmXbWltbQzqWCIba2hqobcIP3avmoD448pBYIgzjdDrR3d0dtN3hcACAaF97ezucTqfsdRwOB9xud0jHEsFQW1sDtU34oXvVHNQHRx4SS4Rh3G43enp6grZ7PB4AQElJibDN6/UiNzdX9jq5ubno6OgI6VgiGGpra6C2CT90r5qD+uDIQ2KJsIzGxkY4HA5UVVUBGJutoYbD4RCOMXIsEQy1tTVQ24QfulfDB/XB4YNSBxCW0N7ejvb2drS0tAimYIIgCCIyUB8cXsiyRFhCRUUFGhsbRbk4tB5Y/i3FyLFEMNTW1kBtE37oXg0P1AeHF7IsjUMKCgpkgwOVyM3NhdvtVgz8q6ioQF1dnWD61Ut3d7fiNc0cG0tY3dahMh7aOhJQ24QfuleNQ31w+CGxNA7p6uqy7Fq1tbUoKSlBTU2N7H6Hw6EoFnw+H4qLi0M6Nl6wsq21GO9tbRXUNuGH7lXroD44MpAbjgiZpqYm5OXlBT2kTU1NwueFCxfC6/UqXsPlcoV0LBEMtbU1UNuEH7pXrYH64MhBYokIifb2dvh8Ptm3Gd63XVFRAY/HE+Tvbm9vB3AoiZrRY4lgqK2tgdom/NC9ah7qgyNMtFOIE/FHV1dXwOl0BqqqqgI1NTWBmpoaYT2osrKyoCUOysvLAzU1NaJtcscZPXa8wZYcUFujidraGqhtzEH3anihPjjy2AKBQCDago2ILwoKClRNtZ2dnUEp82nF69Cpra2F1+uFx+OB1+uFw+FAWVkZcnNz0djYGHQ8tbU1UNsYh+7VyEB9cOQhsUQQBEEQBKECxSwRBEEQBEGoQGKJIAiCIAhCBRJLBEEQBEEQKpBYIgiCIAiCUIHEEkEQBEEQhAoklgiCIAiCIFQgsUQQBEEQBKECLaRLEARBGMbr9QqJJhsaGqJcG4IIL2RZIgiCIAzR3t4uZOkmiPEAiSWCIIgQGa9ioaysDOXl5cjNzZXd7/F4IlwjgggvJJYIgiBCwOPxyK53RoxRW1sb7SoQhGWQWCKIBKWpqQkVFRWw2Wyw2WwoKChAdXW17LGtra0oKCgQjtMz0Hm9XhQVFWHSpEkoKiqyuvpB+Hw+FBQUYOXKlWEvS09damtrhVid9vZ2oS0qKiqiXLvoU1hYiLy8PDQ1NUW7KgRhCSSWCCJBqaqqQktLC6qqqgCMBeEqWULKy8vhdrvhdDrR1dWlK2DX6XSis7MTxcXFltZbie7ubni9XmzcuFF2fyQtGRUVFaK2LCsrQ2dnJ5xOZ8TqEOvU1NSgsbFx3LoqicSCxBJBJDhMRGi5jLxeb0izmiIlEJxOJwKBAFpaWmT3R2pQbm9vF+ojRSmGZ7xSV1dH7jgiIaDUAQSR4DidThQWFqK9vR0+nw8Oh0P2uJaWlriNwWltbYXP54tIWY2NjQknAFpbW+F2u1WPKSoqEqyUeikvL8eiRYtU7zuCiAdILBHEOKC6uhrV1dVoampCTU1N0P54Hsy8Xi8WLVoUMXdge3u7onUrXikvL0d5eXlYrl1WVqZ43xFEvEBiiSDGAVVVVaiurkZjY6PsoNXU1BQU/M2CmAsKCrBv3z54vV7U1dWhsLBQd7l8MPa+fftQUFCgaJ1oampCZ2enINpcLhfKysoAjAmi6upqdHR0CLFSwJhFZN26dQCAjo4OIbja6XSioaEBTU1NaGxshMfjgcPhQFVVleBqbG1tFY4vKyvTtKywc4y4HT0eD0pLS+Hz+eB0OlFbW4vi4mIsWrQIXq8XZWVlWLNmDZqamuBwOOB2u5Gbm4vGxkb4fD4hQHrjxo3CdzIKSx5ZUFAgWN+6urpQXV1t6LcMFZfLhZaWFhJLRHwTIAhiXFBVVRUAEHC73UH7ysvLRX93dXUFHA5HoLOzU3Ubu25hYWHQNQsLCwMtLS2ibTU1NUFlsfKrqqqCjpXWtaysTLassrKyQFlZWdB2htPplN1fU1MTaGxsVDxP7viamhrF/WVlZUHfr7CwULYMVueGhgbRdofDEaipqZHdLt2mB6X2kv6ORujs7Aw0NDQEnE5nwOl0BhoaGhSv19nZGaChhoh36A4miHECG7Skg7nb7Q4SNXKDfiAwJmqkokNOLFVVVSmKF4fDISqvsbExACDQ09MjOs7pdAbVQUmYaYmlhoYG2QFbTfjIUV5eripYpO2mJiLKy8sDAAJdXV1B11DaLvfd1ejs7Aw4nc6g7S0tLabEkhF6enpkf1+CiCdoNhxBjBMKCwtRWFgYFAzd0tIiilfx+Xxob29HSUlJ0DVcLhc6Ojo0y2I5nuQoKysTBUjX1taivLw8KGaqsLAQlZWVmmXpgbn++Lw/Pp8PeXl5hq5jJLaruroaZWVliq6u3NxcOByOILee0+lU3G40iN3pdMLr9cLlcgmz+ICxGKVIuOAACO3V3d0dkfIIIhyQWCKIcQSLS2KiQW7wZ2Koq6sLTU1Non+A9qKpbKkLpWn0ubm5wjR/n88nxPNIkYo4MzgcDpSVlYlm+zU3Nxue3eX1ejXTA7BYr+bmZs3ZhWptZAUOhwMtLS3o6OiAy+WCzWZDUVGRSDhFikjNViSIcEAB3gQxjpAGessFdrOB2uVyhSRWjFgQ2LFGLTxayInA6upqVFRUwOv1Csk3jc4AdDgcmt/P6/WipaUFJSUlqKioQEVFhRCoLnc9I9tDgc10a29vh9vtRmtrK1wuF9xut2K9wgEl7CTiGbIsEcQ4o6qqCl6vFx6PR5hlxcPcM6EmeWRT+JXO7+7uFspk/3d1dYVUlhIdHR1B5TPhx7JKy7kZtdDjCisrK4PD4RBEitISM5Ggvb1dsCKVlZWhoaEBXV1dKC8vj1j6A9Ze8ZqagiAAEksEMe5gg3dFRQVcLpfsMWypCrXzlWBCgU3pl9La2iqKWWJWDzn0ri0mHYh9Pp+sK6uqqgpNTU1obW0NyWqWm5uLffv26T5+zZo16O7ujmoSSzlRVF1dHbEYIq/XS0KJiHtILBHEOKOwsFAI/FWK2WloaIDD4QhatJbPTcQjtbasWbNGlCeIwYKe+XLXrFkDIHhtN4/HIyt45Cw70sBzpQG6uroaPp/PkODh0RPvw4sQh8OBuro6rFy5Uojl4o+T+y5q20OhqakpyMrmdrstC57XoqOjI2IJQwkiXNgCgUAg2pUgCCKytLa2YuPGjZrB2kzAsJgifnaX1+tFbW2tsIwKcznxcTD8+V1dXapLZtTW1sLn86GgoECYDcYnpZSWJU2QWVtbC4/HA5fLhcLCQsV4nIKCArS0tIQ0G8zr9aKgoADSbtPj8WDFihVobW0VgskbGhqQm5uLoqIiQbyx7XLtxpJWam0vKytDdXW1LstYe3u7kJCTF2BOpzNsGbulVFdXh7RUCkHEEiSWCIIYV9TW1oaUCZtRVFSEhoaGiAZHxzMFBQWWx6QRRKQhNxxBEOMGZrkyQ11dXdwuOBxpWltbI5bPiSDCCYklgiASFo/HI4qFam5uxsKFC01ds7y8HD6fL+TZguOJxsZGU1Y8gogVSCwRBJGwrFu3Dq2trQAOzZCzYmZWY2NjVFMCxANNTU1wuVyUX4lICChmiSCIhMXn82HFihVCgHpNTY1l12ZJHslyEozH48G6deuobYiEgcQSQRBEiHg8HorJkYHahUg0SCwRBEEQBEGoQDFLBEEQBEEQKpBYIgiCIAiCUIHEEkEQBEEQhAoklgiCIAiCIFQgsUQQBEEQBKECiSWCIAiCIAgVSCwRBEEQBEGoQGKJIAiCIAhCBRJLBEEQBEEQKvwfkrr35Zu9VQUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from bayes_spec import SpecData\n",
    "from bayes_spec.models import GaussNoiseModel\n",
    "\n",
    "# Generate dummy data format for simulation\n",
    "velocity_axis = np.linspace(-250.0, 250.0, 501) # km/s\n",
    "noise = 1.0 # K\n",
    "brightness_data = noise * np.random.randn(len(velocity_axis)) # K\n",
    "\n",
    "observation = SpecData(\n",
    "    velocity_axis,\n",
    "    brightness_data,\n",
    "    noise,\n",
    "    xlabel=r\"Velocity (km s$^{-1}$)\",\n",
    "    ylabel=\"Brightness Temperature (K)\",\n",
    ")\n",
    "dummy_data = {\"observation\": observation}\n",
    "\n",
    "# Initialize and define the model\n",
    "model = GaussNoiseModel(dummy_data, n_clouds=3, baseline_degree=2, seed=1234, verbose=True)\n",
    "model.add_priors(\n",
    "    prior_line_area = 500.0, # mode of k=2 gamma distribution prior on line area (K km s-1)\n",
    "    prior_fwhm = 25.0, # mode of k=2 gamma distribution prior on FWHM line width (km s-1)\n",
    "    prior_velocity = [0.0, 50.0], # mean and width of normal distribution prior on centroid velocity (km s-1)\n",
    "    prior_baseline_coeffs = [1.0, 1.0, 1.0], # width of normal distribution prior on normalized baseline coefficients\n",
    "    prior_rms = 1.0, # width of half-normal distribution prior on spectral rms (K)\n",
    ")\n",
    "model.add_likelihood()\n",
    "\n",
    "# Simulate observation\n",
    "sim_brightness = model.model.observation.eval({\n",
    "    \"fwhm\": [25.0, 40.0, 35.0], # FWHM line width (km/s)\n",
    "    \"line_area\": [250.0, 125.0, 175.0], # line area (K km/s)\n",
    "    \"velocity\": [-35.0, 10.0, 55.0], # velocity (km/s)\n",
    "    \"baseline_observation_norm\": [-0.5, -2.0, 3.0], # normalized baseline coefficients\n",
    "    \"rms_observation\": noise, # spectral rms (K)\n",
    "})\n",
    "\n",
    "# Plot the simulated data\n",
    "plt.plot(dummy_data[\"observation\"].spectral, sim_brightness, 'k-')\n",
    "plt.xlabel(dummy_data[\"observation\"].xlabel)\n",
    "_ = plt.ylabel(dummy_data[\"observation\"].ylabel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "290a1ff7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:10.597143Z",
     "start_time": "2024-08-05T19:59:10.593885Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pack simulated data into SpecData\n",
    "observation = SpecData(\n",
    "    velocity_axis,\n",
    "    sim_brightness,\n",
    "    noise,\n",
    "    xlabel=r\"Velocity (km s$^{-1}$)\",\n",
    "    ylabel=\"Brightness Temperature (K)\",\n",
    ")\n",
    "data = {\"observation\": observation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20a2225f",
   "metadata": {},
   "source": [
    "## `Optimize`\n",
    "\n",
    "We use the `Optimize` class for optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ec488d1d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:10.725568Z",
     "start_time": "2024-08-05T19:59:10.599441Z"
    }
   },
   "outputs": [],
   "source": [
    "from bayes_spec import Optimize\n",
    "\n",
    "# Initialize optimizer\n",
    "opt = Optimize(\n",
    "    GaussNoiseModel,  # model definition\n",
    "    data,  # data dictionary\n",
    "    max_n_clouds=8,  # maximum number of clouds\n",
    "    baseline_degree=2,  # polynomial baseline degree\n",
    "    seed=1234,  # random seed\n",
    "    verbose=True,  # verbosity\n",
    ")\n",
    "\n",
    "# Define each model\n",
    "opt.add_priors(\n",
    "    prior_line_area = 200.0, # mode of k=2 gamma distribution prior on line area (K km s-1)\n",
    "    prior_fwhm = 30.0, # mode of k=2 gamma distribution prior on FWHM line width (km s-1)\n",
    "    prior_velocity = [0.0, 50.0], # mean and width of normal distribution prior on centroid velocity (km s-1)\n",
    "    prior_baseline_coeffs = [1.0, 1.0, 1.0], # width of normal distribution prior on normalized baseline coefficients\n",
    "    prior_rms = 2.0, # width of half-normal distribution prior on spectral rms (K)\n",
    ")\n",
    "opt.add_likelihood()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4b0d06",
   "metadata": {},
   "source": [
    "`Optimize` has created `max_n_clouds` models, where `opt.models[1]` has `n_clouds=1`, `opt.models[2]` has `n_clouds=2`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e66910c3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:10.730059Z",
     "start_time": "2024-08-05T19:59:10.727121Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bayes_spec.models.gauss_noise_model.GaussNoiseModel object at 0x7eade3784e90>\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "print(opt.models[4])\n",
    "print(opt.models[4].n_clouds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cc17df3",
   "metadata": {},
   "source": [
    "By default (`approx=True`), the optimization algorithm first loops over every model and approximates the posterior distribution using variational inference. We can supply arguments to `fit` and `sample` via dictionaries. Whichever model is the first to have a BIC within `bic_threshold` of the minimum BIC is the \"best\" model, and is then sampled with MCMC. The algorithm will stop early if two successive models fail to converge to a unique solution with an improving BIC. Note the `start_spread` parameter, which is like passing `start = {\"velocity_norm\": np.linspace(-3.0, 3.0, n_clouds)}` to `model.fit` or to `model.sample` within `init_kwargs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cf08c08c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:59.141293Z",
     "start_time": "2024-08-05T19:59:10.731559Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Null hypothesis BIC = 2.895e+03\n",
      "Approximating n_cloud = 1 posterior...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd929842bf6f4f26ac8e65b2862e6451",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 23000\n",
      "Interrupted at 22,999 [2%]: Average Loss = 1,139.8\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "137a453e78f341918215ad9cee18483a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cloud = 1 BIC = 1.965e+03\n",
      "\n",
      "Approximating n_cloud = 2 posterior...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "75f3161a91bb4059a0be51c20e3e47db",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 28600\n",
      "Interrupted at 28,599 [2%]: Average Loss = 1,075.5\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eec0d889dad74be792c71f79767bc6e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cloud = 2 BIC = 1.533e+03\n",
      "\n",
      "Approximating n_cloud = 3 posterior...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c754644267eb4c4283e08a4bbfa4e8a0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 37200\n",
      "Interrupted at 37,199 [3%]: Average Loss = 962.58\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ca9fe5207cb04180bb97e2c1d77031c0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cloud = 3 BIC = 1.507e+03\n",
      "\n",
      "Approximating n_cloud = 4 posterior...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ce822c3544d54ac7a986f59260c60793",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 29700\n",
      "Interrupted at 29,699 [2%]: Average Loss = 1,026\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "58a94e53ae934a9089a72848c0d1b6b9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cloud = 4 BIC = 1.522e+03\n",
      "\n",
      "Stopping criteria met.\n",
      "Approximating n_cloud = 5 posterior...\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9d95a1bcf1034bfca25e4ce08f4f6878",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 21200\n",
      "Interrupted at 21,199 [2%]: Average Loss = 1,241.4\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9bbe1e505d9c45b99d190823150068e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cloud = 5 BIC = 1.546e+03\n",
      "\n",
      "Stopping criteria met.\n",
      "Stopping early.\n",
      "Sampling best model (n_cloud = 3)...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "aca61a0e47934693b26e7228ab4d5812",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 37200\n",
      "Interrupted at 37,199 [3%]: Average Loss = 962.58\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e47a76f7cbd4a68b6273b2f0a93af77",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "945c97f2f38c42adbb14ecfbfa65d68e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GMM converged to unique solution\n"
     ]
    }
   ],
   "source": [
    "fit_kwargs = {\n",
    "    \"rel_tolerance\": 0.01,\n",
    "    \"abs_tolerance\": 0.01,\n",
    "    \"learning_rate\": 0.001,\n",
    "}\n",
    "sample_kwargs = {\n",
    "    \"chains\": 4,\n",
    "    \"cores\": 4,\n",
    "    \"init_kwargs\": fit_kwargs,\n",
    "    \"nuts_kwargs\": {\"target_accept\": 0.8},\n",
    "}\n",
    "opt.optimize(\n",
    "    bic_threshold=10.0,\n",
    "    kl_div_threshold=0.1,\n",
    "    sample_kwargs=sample_kwargs,\n",
    "    fit_kwargs=fit_kwargs,\n",
    "    start_spread = {\"velocity_norm\": [-3.0, 3.0]},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02ec3edd",
   "metadata": {},
   "source": [
    "The \"best\" model is saved in `opt.best_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3d5d558e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T19:59:59.535930Z",
     "start_time": "2024-08-05T19:59:59.143289Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best model has n_clouds = 3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[0]</th>\n",
       "      <td>-0.380</td>\n",
       "      <td>0.071</td>\n",
       "      <td>-0.514</td>\n",
       "      <td>-0.249</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2770.0</td>\n",
       "      <td>2446.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[1]</th>\n",
       "      <td>-1.870</td>\n",
       "      <td>0.156</td>\n",
       "      <td>-2.177</td>\n",
       "      <td>-1.594</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.003</td>\n",
       "      <td>4185.0</td>\n",
       "      <td>3016.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[2]</th>\n",
       "      <td>1.309</td>\n",
       "      <td>0.884</td>\n",
       "      <td>-0.289</td>\n",
       "      <td>3.045</td>\n",
       "      <td>0.017</td>\n",
       "      <td>0.013</td>\n",
       "      <td>2826.0</td>\n",
       "      <td>2878.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[0]</th>\n",
       "      <td>0.229</td>\n",
       "      <td>0.142</td>\n",
       "      <td>-0.025</td>\n",
       "      <td>0.509</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1076.0</td>\n",
       "      <td>1820.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[1]</th>\n",
       "      <td>-0.706</td>\n",
       "      <td>0.010</td>\n",
       "      <td>-0.724</td>\n",
       "      <td>-0.689</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2519.0</td>\n",
       "      <td>2386.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[2]</th>\n",
       "      <td>1.129</td>\n",
       "      <td>0.041</td>\n",
       "      <td>1.049</td>\n",
       "      <td>1.203</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1302.0</td>\n",
       "      <td>1911.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[0]</th>\n",
       "      <td>0.825</td>\n",
       "      <td>0.221</td>\n",
       "      <td>0.418</td>\n",
       "      <td>1.225</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.003</td>\n",
       "      <td>855.0</td>\n",
       "      <td>1528.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[1]</th>\n",
       "      <td>1.116</td>\n",
       "      <td>0.081</td>\n",
       "      <td>0.972</td>\n",
       "      <td>1.277</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1518.0</td>\n",
       "      <td>1885.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[2]</th>\n",
       "      <td>0.724</td>\n",
       "      <td>0.177</td>\n",
       "      <td>0.398</td>\n",
       "      <td>1.040</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.003</td>\n",
       "      <td>811.0</td>\n",
       "      <td>1767.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[0]</th>\n",
       "      <td>2.093</td>\n",
       "      <td>0.512</td>\n",
       "      <td>1.153</td>\n",
       "      <td>3.020</td>\n",
       "      <td>0.017</td>\n",
       "      <td>0.008</td>\n",
       "      <td>907.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[1]</th>\n",
       "      <td>0.793</td>\n",
       "      <td>0.040</td>\n",
       "      <td>0.719</td>\n",
       "      <td>0.868</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2154.0</td>\n",
       "      <td>2492.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[2]</th>\n",
       "      <td>1.188</td>\n",
       "      <td>0.172</td>\n",
       "      <td>0.851</td>\n",
       "      <td>1.505</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>2049.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rms_observation_norm</th>\n",
       "      <td>0.501</td>\n",
       "      <td>0.016</td>\n",
       "      <td>0.471</td>\n",
       "      <td>0.531</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>4012.0</td>\n",
       "      <td>2812.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[0]</th>\n",
       "      <td>164.955</td>\n",
       "      <td>44.245</td>\n",
       "      <td>83.563</td>\n",
       "      <td>245.096</td>\n",
       "      <td>1.509</td>\n",
       "      <td>0.690</td>\n",
       "      <td>855.0</td>\n",
       "      <td>1528.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[1]</th>\n",
       "      <td>223.113</td>\n",
       "      <td>16.157</td>\n",
       "      <td>194.479</td>\n",
       "      <td>255.482</td>\n",
       "      <td>0.414</td>\n",
       "      <td>0.278</td>\n",
       "      <td>1518.0</td>\n",
       "      <td>1885.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[2]</th>\n",
       "      <td>144.807</td>\n",
       "      <td>35.410</td>\n",
       "      <td>79.565</td>\n",
       "      <td>207.919</td>\n",
       "      <td>1.254</td>\n",
       "      <td>0.548</td>\n",
       "      <td>811.0</td>\n",
       "      <td>1767.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[0]</th>\n",
       "      <td>62.782</td>\n",
       "      <td>15.372</td>\n",
       "      <td>34.597</td>\n",
       "      <td>90.608</td>\n",
       "      <td>0.510</td>\n",
       "      <td>0.236</td>\n",
       "      <td>907.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[1]</th>\n",
       "      <td>23.777</td>\n",
       "      <td>1.197</td>\n",
       "      <td>21.566</td>\n",
       "      <td>26.031</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>2154.0</td>\n",
       "      <td>2492.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[2]</th>\n",
       "      <td>35.653</td>\n",
       "      <td>5.162</td>\n",
       "      <td>25.516</td>\n",
       "      <td>45.154</td>\n",
       "      <td>0.158</td>\n",
       "      <td>0.089</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>2049.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[0]</th>\n",
       "      <td>11.450</td>\n",
       "      <td>7.103</td>\n",
       "      <td>-1.238</td>\n",
       "      <td>25.439</td>\n",
       "      <td>0.220</td>\n",
       "      <td>0.122</td>\n",
       "      <td>1076.0</td>\n",
       "      <td>1820.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[1]</th>\n",
       "      <td>-35.296</td>\n",
       "      <td>0.478</td>\n",
       "      <td>-36.211</td>\n",
       "      <td>-34.451</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>2519.0</td>\n",
       "      <td>2386.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[2]</th>\n",
       "      <td>56.427</td>\n",
       "      <td>2.063</td>\n",
       "      <td>52.458</td>\n",
       "      <td>60.155</td>\n",
       "      <td>0.058</td>\n",
       "      <td>0.036</td>\n",
       "      <td>1302.0</td>\n",
       "      <td>1911.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[0]</th>\n",
       "      <td>2.467</td>\n",
       "      <td>0.240</td>\n",
       "      <td>2.036</td>\n",
       "      <td>2.923</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3168.0</td>\n",
       "      <td>3381.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[1]</th>\n",
       "      <td>8.812</td>\n",
       "      <td>0.404</td>\n",
       "      <td>8.084</td>\n",
       "      <td>9.619</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.006</td>\n",
       "      <td>2156.0</td>\n",
       "      <td>2960.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[2]</th>\n",
       "      <td>3.774</td>\n",
       "      <td>0.545</td>\n",
       "      <td>2.766</td>\n",
       "      <td>4.680</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.007</td>\n",
       "      <td>1009.0</td>\n",
       "      <td>2302.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rms_observation</th>\n",
       "      <td>1.002</td>\n",
       "      <td>0.032</td>\n",
       "      <td>0.942</td>\n",
       "      <td>1.063</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4012.0</td>\n",
       "      <td>2812.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                 mean      sd   hdi_3%  hdi_97%  mcse_mean  \\\n",
       "baseline_observation_norm[0]   -0.380   0.071   -0.514   -0.249      0.001   \n",
       "baseline_observation_norm[1]   -1.870   0.156   -2.177   -1.594      0.002   \n",
       "baseline_observation_norm[2]    1.309   0.884   -0.289    3.045      0.017   \n",
       "velocity_norm[0]                0.229   0.142   -0.025    0.509      0.004   \n",
       "velocity_norm[1]               -0.706   0.010   -0.724   -0.689      0.000   \n",
       "velocity_norm[2]                1.129   0.041    1.049    1.203      0.001   \n",
       "line_area_norm[0]               0.825   0.221    0.418    1.225      0.008   \n",
       "line_area_norm[1]               1.116   0.081    0.972    1.277      0.002   \n",
       "line_area_norm[2]               0.724   0.177    0.398    1.040      0.006   \n",
       "fwhm_norm[0]                    2.093   0.512    1.153    3.020      0.017   \n",
       "fwhm_norm[1]                    0.793   0.040    0.719    0.868      0.001   \n",
       "fwhm_norm[2]                    1.188   0.172    0.851    1.505      0.005   \n",
       "rms_observation_norm            0.501   0.016    0.471    0.531      0.000   \n",
       "line_area[0]                  164.955  44.245   83.563  245.096      1.509   \n",
       "line_area[1]                  223.113  16.157  194.479  255.482      0.414   \n",
       "line_area[2]                  144.807  35.410   79.565  207.919      1.254   \n",
       "fwhm[0]                        62.782  15.372   34.597   90.608      0.510   \n",
       "fwhm[1]                        23.777   1.197   21.566   26.031      0.026   \n",
       "fwhm[2]                        35.653   5.162   25.516   45.154      0.158   \n",
       "velocity[0]                    11.450   7.103   -1.238   25.439      0.220   \n",
       "velocity[1]                   -35.296   0.478  -36.211  -34.451      0.010   \n",
       "velocity[2]                    56.427   2.063   52.458   60.155      0.058   \n",
       "amplitude[0]                    2.467   0.240    2.036    2.923      0.004   \n",
       "amplitude[1]                    8.812   0.404    8.084    9.619      0.009   \n",
       "amplitude[2]                    3.774   0.545    2.766    4.680      0.018   \n",
       "rms_observation                 1.002   0.032    0.942    1.063      0.001   \n",
       "\n",
       "                              mcse_sd  ess_bulk  ess_tail  r_hat  \n",
       "baseline_observation_norm[0]    0.001    2770.0    2446.0   1.00  \n",
       "baseline_observation_norm[1]    0.003    4185.0    3016.0   1.00  \n",
       "baseline_observation_norm[2]    0.013    2826.0    2878.0   1.00  \n",
       "velocity_norm[0]                0.002    1076.0    1820.0   1.00  \n",
       "velocity_norm[1]                0.000    2519.0    2386.0   1.00  \n",
       "velocity_norm[2]                0.001    1302.0    1911.0   1.00  \n",
       "line_area_norm[0]               0.003     855.0    1528.0   1.01  \n",
       "line_area_norm[1]               0.001    1518.0    1885.0   1.00  \n",
       "line_area_norm[2]               0.003     811.0    1767.0   1.00  \n",
       "fwhm_norm[0]                    0.008     907.0    1559.0   1.01  \n",
       "fwhm_norm[1]                    0.001    2154.0    2492.0   1.00  \n",
       "fwhm_norm[2]                    0.003    1065.0    2049.0   1.00  \n",
       "rms_observation_norm            0.000    4012.0    2812.0   1.00  \n",
       "line_area[0]                    0.690     855.0    1528.0   1.01  \n",
       "line_area[1]                    0.278    1518.0    1885.0   1.00  \n",
       "line_area[2]                    0.548     811.0    1767.0   1.00  \n",
       "fwhm[0]                         0.236     907.0    1559.0   1.01  \n",
       "fwhm[1]                         0.018    2154.0    2492.0   1.00  \n",
       "fwhm[2]                         0.089    1065.0    2049.0   1.00  \n",
       "velocity[0]                     0.122    1076.0    1820.0   1.00  \n",
       "velocity[1]                     0.008    2519.0    2386.0   1.00  \n",
       "velocity[2]                     0.036    1302.0    1911.0   1.00  \n",
       "amplitude[0]                    0.003    3168.0    3381.0   1.00  \n",
       "amplitude[1]                    0.006    2156.0    2960.0   1.00  \n",
       "amplitude[2]                    0.007    1009.0    2302.0   1.00  \n",
       "rms_observation                 0.001    4012.0    2812.0   1.00  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(f\"Best model has n_clouds = {opt.best_model.n_clouds}\")\n",
    "az.summary(opt.best_model.trace.solution_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8480130",
   "metadata": {},
   "source": [
    "With `approx=False`, the optimization algorithm samples each model with MCMC (or SMC if `smc=True`) and determines which is the first model with a BIC within `bic_threshold` of the minimum across all models. This is more accurate, but slower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5bc57929",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T20:02:48.194751Z",
     "start_time": "2024-08-05T19:59:59.537776Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Null hypothesis BIC = 2.895e+03\n",
      "Sampling n_cloud = 1 posterior...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0d558086573e485ca9a59dbc7efb2c09",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 23000\n",
      "Interrupted at 22,999 [2%]: Average Loss = 1,139.8\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c0e9098348514bcabe22a4f1af32f3af",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6782c611b0244d95b3add18faff1070b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GMM converged to unique solution\n",
      "n_cloud = 1 solution = 0 BIC = 1.964e+03\n",
      "\n",
      "Sampling n_cloud = 2 posterior...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a59c2866f6264f589142c48c49a0fc57",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 28600\n",
      "Interrupted at 28,599 [2%]: Average Loss = 1,075.5\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8d9ded7e330b49c194d843bbb7120799",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ad944ca92a804299a91c5b9fcc0ebffa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GMM converged to unique solution\n",
      "n_cloud = 2 solution = 0 BIC = 1.531e+03\n",
      "\n",
      "Sampling n_cloud = 3 posterior...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a9117b6e9c644caacac22691e6233e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 37200\n",
      "Interrupted at 37,199 [3%]: Average Loss = 962.58\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c3b59ee68114763a78e61ebcc1821b9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 3 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0845e4b3f265413db46330a03b3bc6db",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GMM converged to unique solution\n",
      "n_cloud = 3 solution = 0 BIC = 1.504e+03\n",
      "\n",
      "Sampling n_cloud = 4 posterior...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "28f38b25414741fcb262604fb08fe55e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 29700\n",
      "Interrupted at 29,699 [2%]: Average Loss = 1,026\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "39a067d8f1c3476086fd6483e38fb044",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 9 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "92604b2616494bd3982593a4cbc85b7c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There were 668 divergences in converged chains.\n",
      "No solution found!\n",
      "0 of 4 chains appear converged.\n",
      "\n",
      "Stopping criteria met.\n",
      "Sampling n_cloud = 5 posterior...\n",
      "Initializing NUTS using custom advi+adapt_diag strategy\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "22e266450f584c1284f45fe0a57b4fb1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Convergence achieved at 21200\n",
      "Interrupted at 21,199 [2%]: Average Loss = 1,241.4\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [baseline_observation_norm, line_area_norm, fwhm_norm, velocity_norm, rms_observation_norm]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "24841d62527744a2bf0e8137d6ab7bd6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 12 seconds.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1e7f5b2f38b44946bc4d02c59d701736",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding log-likelihood to trace\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<pre style=\"white-space:pre;overflow-x:auto;line-height:normal;font-family:Menlo,'DejaVu Sans Mono',consolas,'Courier New',monospace\"></pre>\n"
      ],
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There were 2 divergences in converged chains.\n",
      "No solution found!\n",
      "0 of 4 chains appear converged.\n",
      "\n",
      "Stopping criteria met.\n",
      "Stopping early.\n"
     ]
    }
   ],
   "source": [
    "opt.optimize(\n",
    "    bic_threshold=10.0,\n",
    "    kl_div_threshold=0.1,\n",
    "    sample_kwargs=sample_kwargs,\n",
    "    fit_kwargs=fit_kwargs,\n",
    "    start_spread = {\"velocity_norm\": [-3.0, 3.0]},\n",
    "    approx=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "752176c2",
   "metadata": {},
   "source": [
    "Notice that models with more complexity than is present in the data (i.e., `n_clouds` > the true number of clouds) tend to have more divergences and difficulty converging. Eventually the stopping criteria terminate the algorithm early."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "eb0f2a29",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T20:02:48.370691Z",
     "start_time": "2024-08-05T20:02:48.198400Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best model has n_clouds = 3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[0]</th>\n",
       "      <td>-0.380</td>\n",
       "      <td>0.071</td>\n",
       "      <td>-0.514</td>\n",
       "      <td>-0.249</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2770.0</td>\n",
       "      <td>2446.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[1]</th>\n",
       "      <td>-1.870</td>\n",
       "      <td>0.156</td>\n",
       "      <td>-2.177</td>\n",
       "      <td>-1.594</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.003</td>\n",
       "      <td>4185.0</td>\n",
       "      <td>3016.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>baseline_observation_norm[2]</th>\n",
       "      <td>1.309</td>\n",
       "      <td>0.884</td>\n",
       "      <td>-0.289</td>\n",
       "      <td>3.045</td>\n",
       "      <td>0.017</td>\n",
       "      <td>0.013</td>\n",
       "      <td>2826.0</td>\n",
       "      <td>2878.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[0]</th>\n",
       "      <td>0.229</td>\n",
       "      <td>0.142</td>\n",
       "      <td>-0.025</td>\n",
       "      <td>0.509</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1076.0</td>\n",
       "      <td>1820.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[1]</th>\n",
       "      <td>-0.706</td>\n",
       "      <td>0.010</td>\n",
       "      <td>-0.724</td>\n",
       "      <td>-0.689</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2519.0</td>\n",
       "      <td>2386.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity_norm[2]</th>\n",
       "      <td>1.129</td>\n",
       "      <td>0.041</td>\n",
       "      <td>1.049</td>\n",
       "      <td>1.203</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1302.0</td>\n",
       "      <td>1911.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[0]</th>\n",
       "      <td>0.825</td>\n",
       "      <td>0.221</td>\n",
       "      <td>0.418</td>\n",
       "      <td>1.225</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.003</td>\n",
       "      <td>855.0</td>\n",
       "      <td>1528.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[1]</th>\n",
       "      <td>1.116</td>\n",
       "      <td>0.081</td>\n",
       "      <td>0.972</td>\n",
       "      <td>1.277</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1518.0</td>\n",
       "      <td>1885.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area_norm[2]</th>\n",
       "      <td>0.724</td>\n",
       "      <td>0.177</td>\n",
       "      <td>0.398</td>\n",
       "      <td>1.040</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.003</td>\n",
       "      <td>811.0</td>\n",
       "      <td>1767.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[0]</th>\n",
       "      <td>2.093</td>\n",
       "      <td>0.512</td>\n",
       "      <td>1.153</td>\n",
       "      <td>3.020</td>\n",
       "      <td>0.017</td>\n",
       "      <td>0.008</td>\n",
       "      <td>907.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[1]</th>\n",
       "      <td>0.793</td>\n",
       "      <td>0.040</td>\n",
       "      <td>0.719</td>\n",
       "      <td>0.868</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2154.0</td>\n",
       "      <td>2492.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm_norm[2]</th>\n",
       "      <td>1.188</td>\n",
       "      <td>0.172</td>\n",
       "      <td>0.851</td>\n",
       "      <td>1.505</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>2049.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rms_observation_norm</th>\n",
       "      <td>0.501</td>\n",
       "      <td>0.016</td>\n",
       "      <td>0.471</td>\n",
       "      <td>0.531</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>4012.0</td>\n",
       "      <td>2812.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[0]</th>\n",
       "      <td>164.955</td>\n",
       "      <td>44.245</td>\n",
       "      <td>83.563</td>\n",
       "      <td>245.096</td>\n",
       "      <td>1.509</td>\n",
       "      <td>0.690</td>\n",
       "      <td>855.0</td>\n",
       "      <td>1528.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[1]</th>\n",
       "      <td>223.113</td>\n",
       "      <td>16.157</td>\n",
       "      <td>194.479</td>\n",
       "      <td>255.482</td>\n",
       "      <td>0.414</td>\n",
       "      <td>0.278</td>\n",
       "      <td>1518.0</td>\n",
       "      <td>1885.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>line_area[2]</th>\n",
       "      <td>144.807</td>\n",
       "      <td>35.410</td>\n",
       "      <td>79.565</td>\n",
       "      <td>207.919</td>\n",
       "      <td>1.254</td>\n",
       "      <td>0.548</td>\n",
       "      <td>811.0</td>\n",
       "      <td>1767.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[0]</th>\n",
       "      <td>62.782</td>\n",
       "      <td>15.372</td>\n",
       "      <td>34.597</td>\n",
       "      <td>90.608</td>\n",
       "      <td>0.510</td>\n",
       "      <td>0.236</td>\n",
       "      <td>907.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[1]</th>\n",
       "      <td>23.777</td>\n",
       "      <td>1.197</td>\n",
       "      <td>21.566</td>\n",
       "      <td>26.031</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>2154.0</td>\n",
       "      <td>2492.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fwhm[2]</th>\n",
       "      <td>35.653</td>\n",
       "      <td>5.162</td>\n",
       "      <td>25.516</td>\n",
       "      <td>45.154</td>\n",
       "      <td>0.158</td>\n",
       "      <td>0.089</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>2049.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[0]</th>\n",
       "      <td>11.450</td>\n",
       "      <td>7.103</td>\n",
       "      <td>-1.238</td>\n",
       "      <td>25.439</td>\n",
       "      <td>0.220</td>\n",
       "      <td>0.122</td>\n",
       "      <td>1076.0</td>\n",
       "      <td>1820.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[1]</th>\n",
       "      <td>-35.296</td>\n",
       "      <td>0.478</td>\n",
       "      <td>-36.211</td>\n",
       "      <td>-34.451</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>2519.0</td>\n",
       "      <td>2386.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity[2]</th>\n",
       "      <td>56.427</td>\n",
       "      <td>2.063</td>\n",
       "      <td>52.458</td>\n",
       "      <td>60.155</td>\n",
       "      <td>0.058</td>\n",
       "      <td>0.036</td>\n",
       "      <td>1302.0</td>\n",
       "      <td>1911.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[0]</th>\n",
       "      <td>2.467</td>\n",
       "      <td>0.240</td>\n",
       "      <td>2.036</td>\n",
       "      <td>2.923</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>3168.0</td>\n",
       "      <td>3381.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[1]</th>\n",
       "      <td>8.812</td>\n",
       "      <td>0.404</td>\n",
       "      <td>8.084</td>\n",
       "      <td>9.619</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.006</td>\n",
       "      <td>2156.0</td>\n",
       "      <td>2960.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>amplitude[2]</th>\n",
       "      <td>3.774</td>\n",
       "      <td>0.545</td>\n",
       "      <td>2.766</td>\n",
       "      <td>4.680</td>\n",
       "      <td>0.018</td>\n",
       "      <td>0.007</td>\n",
       "      <td>1009.0</td>\n",
       "      <td>2302.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rms_observation</th>\n",
       "      <td>1.002</td>\n",
       "      <td>0.032</td>\n",
       "      <td>0.942</td>\n",
       "      <td>1.063</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4012.0</td>\n",
       "      <td>2812.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                 mean      sd   hdi_3%  hdi_97%  mcse_mean  \\\n",
       "baseline_observation_norm[0]   -0.380   0.071   -0.514   -0.249      0.001   \n",
       "baseline_observation_norm[1]   -1.870   0.156   -2.177   -1.594      0.002   \n",
       "baseline_observation_norm[2]    1.309   0.884   -0.289    3.045      0.017   \n",
       "velocity_norm[0]                0.229   0.142   -0.025    0.509      0.004   \n",
       "velocity_norm[1]               -0.706   0.010   -0.724   -0.689      0.000   \n",
       "velocity_norm[2]                1.129   0.041    1.049    1.203      0.001   \n",
       "line_area_norm[0]               0.825   0.221    0.418    1.225      0.008   \n",
       "line_area_norm[1]               1.116   0.081    0.972    1.277      0.002   \n",
       "line_area_norm[2]               0.724   0.177    0.398    1.040      0.006   \n",
       "fwhm_norm[0]                    2.093   0.512    1.153    3.020      0.017   \n",
       "fwhm_norm[1]                    0.793   0.040    0.719    0.868      0.001   \n",
       "fwhm_norm[2]                    1.188   0.172    0.851    1.505      0.005   \n",
       "rms_observation_norm            0.501   0.016    0.471    0.531      0.000   \n",
       "line_area[0]                  164.955  44.245   83.563  245.096      1.509   \n",
       "line_area[1]                  223.113  16.157  194.479  255.482      0.414   \n",
       "line_area[2]                  144.807  35.410   79.565  207.919      1.254   \n",
       "fwhm[0]                        62.782  15.372   34.597   90.608      0.510   \n",
       "fwhm[1]                        23.777   1.197   21.566   26.031      0.026   \n",
       "fwhm[2]                        35.653   5.162   25.516   45.154      0.158   \n",
       "velocity[0]                    11.450   7.103   -1.238   25.439      0.220   \n",
       "velocity[1]                   -35.296   0.478  -36.211  -34.451      0.010   \n",
       "velocity[2]                    56.427   2.063   52.458   60.155      0.058   \n",
       "amplitude[0]                    2.467   0.240    2.036    2.923      0.004   \n",
       "amplitude[1]                    8.812   0.404    8.084    9.619      0.009   \n",
       "amplitude[2]                    3.774   0.545    2.766    4.680      0.018   \n",
       "rms_observation                 1.002   0.032    0.942    1.063      0.001   \n",
       "\n",
       "                              mcse_sd  ess_bulk  ess_tail  r_hat  \n",
       "baseline_observation_norm[0]    0.001    2770.0    2446.0   1.00  \n",
       "baseline_observation_norm[1]    0.003    4185.0    3016.0   1.00  \n",
       "baseline_observation_norm[2]    0.013    2826.0    2878.0   1.00  \n",
       "velocity_norm[0]                0.002    1076.0    1820.0   1.00  \n",
       "velocity_norm[1]                0.000    2519.0    2386.0   1.00  \n",
       "velocity_norm[2]                0.001    1302.0    1911.0   1.00  \n",
       "line_area_norm[0]               0.003     855.0    1528.0   1.01  \n",
       "line_area_norm[1]               0.001    1518.0    1885.0   1.00  \n",
       "line_area_norm[2]               0.003     811.0    1767.0   1.00  \n",
       "fwhm_norm[0]                    0.008     907.0    1559.0   1.01  \n",
       "fwhm_norm[1]                    0.001    2154.0    2492.0   1.00  \n",
       "fwhm_norm[2]                    0.003    1065.0    2049.0   1.00  \n",
       "rms_observation_norm            0.000    4012.0    2812.0   1.00  \n",
       "line_area[0]                    0.690     855.0    1528.0   1.01  \n",
       "line_area[1]                    0.278    1518.0    1885.0   1.00  \n",
       "line_area[2]                    0.548     811.0    1767.0   1.00  \n",
       "fwhm[0]                         0.236     907.0    1559.0   1.01  \n",
       "fwhm[1]                         0.018    2154.0    2492.0   1.00  \n",
       "fwhm[2]                         0.089    1065.0    2049.0   1.00  \n",
       "velocity[0]                     0.122    1076.0    1820.0   1.00  \n",
       "velocity[1]                     0.008    2519.0    2386.0   1.00  \n",
       "velocity[2]                     0.036    1302.0    1911.0   1.00  \n",
       "amplitude[0]                    0.003    3168.0    3381.0   1.00  \n",
       "amplitude[1]                    0.006    2156.0    2960.0   1.00  \n",
       "amplitude[2]                    0.007    1009.0    2302.0   1.00  \n",
       "rms_observation                 0.001    4012.0    2812.0   1.00  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(f\"Best model has n_clouds = {opt.best_model.n_clouds}\")\n",
    "az.summary(opt.best_model.trace.solution_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a73fe34",
   "metadata": {},
   "source": [
    "With posteriors sampled with MCMC, we can also use leave-one-out cross-validation to perform model comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9157b48e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-08-05T20:02:53.439889Z",
     "start_time": "2024-08-05T20:02:51.123757Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rank</th>\n",
       "      <th>elpd_loo</th>\n",
       "      <th>p_loo</th>\n",
       "      <th>elpd_diff</th>\n",
       "      <th>weight</th>\n",
       "      <th>se</th>\n",
       "      <th>dse</th>\n",
       "      <th>warning</th>\n",
       "      <th>scale</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>-717.831101</td>\n",
       "      <td>11.490170</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>9.386881e-01</td>\n",
       "      <td>14.599520</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>False</td>\n",
       "      <td>log</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>-739.711946</td>\n",
       "      <td>9.757259</td>\n",
       "      <td>21.880846</td>\n",
       "      <td>6.131194e-02</td>\n",
       "      <td>14.725176</td>\n",
       "      <td>7.180182</td>\n",
       "      <td>False</td>\n",
       "      <td>log</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>-964.569764</td>\n",
       "      <td>8.267954</td>\n",
       "      <td>246.738664</td>\n",
       "      <td>2.671086e-12</td>\n",
       "      <td>22.670923</td>\n",
       "      <td>23.981153</td>\n",
       "      <td>False</td>\n",
       "      <td>log</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   rank    elpd_loo      p_loo   elpd_diff        weight         se  \\\n",
       "3     0 -717.831101  11.490170    0.000000  9.386881e-01  14.599520   \n",
       "2     1 -739.711946   9.757259   21.880846  6.131194e-02  14.725176   \n",
       "1     2 -964.569764   8.267954  246.738664  2.671086e-12  22.670923   \n",
       "\n",
       "         dse  warning scale  \n",
       "3   0.000000    False   log  \n",
       "2   7.180182    False   log  \n",
       "1  23.981153    False   log  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# leave-one-out cross validation\n",
    "loo = az.compare({\n",
    "    n_gauss: model.trace for n_gauss, model in opt.models.items()\n",
    "    if model.solutions # ignore models with no valid solutions\n",
    "})\n",
    "loo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a9508ba",
   "metadata": {},
   "source": [
    "The model with the largest `elpd_loo` is preferred."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}