{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# `qp` practical example\n",
    "\n",
    "_Alex Malz, Phil Marshall, Eric Charles_\n",
    "\n",
    "In this notebook we use the `qp` module to study some photo-Z PDFs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup, reading ensemble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports\n",
    "\n",
    "First let's import the packages we will need for this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import qp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data files\n",
    "\n",
    "Now lets download the data files we will need, if we haven't already"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "  0     0    0     0    0     0      0      0 --:--:--  0:00:01 --:--:--     0\n",
      "curl: (60) SSL certificate problem: unable to get local issuer certificate\n",
      "More details here: https://curl.se/docs/sslcerts.html\n",
      "\n",
      "curl failed to verify the legitimacy of the server and therefore could not\n",
      "establish a secure connection to it. To learn more about this situation and\n",
      "how to fix it, please visit the web page mentioned above.\n",
      "Warning: Got more output options than URLs\n"
     ]
    }
   ],
   "source": [
    "# base_url = 'https://slac.stanford.edu/~echarles/qp_example'\n",
    "# if not os.path.exists('qp_test_ensemble.hf5'):\n",
    "#     os.system('curl -o %s -OL %s/%s' % ('qp_test_ensemble.hf5', base_url, 'qp_test_ensemble.hf5'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading ensemble\n",
    "\n",
    "Now we read the ensemble, note that we only need the name of the data file, the name of the metadata file is assumed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ens = qp.read('qp_test_ensemble.hf5')\n",
    "QP_DIR = os.path.abspath(os.path.dirname(qp.__file__))\n",
    "data_file = os.path.join(QP_DIR, 'data', 'test.hdf5')\n",
    "ens = qp.read(data_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploration\n",
    "\n",
    "This will show use that the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ensemble =  <qp.ensemble.Ensemble object at 0x7fc2772aa510>\n",
      "Rep      =  mixmod\n",
      "NPDF     =  100\n",
      "Metadata =  {'pdf_name': array([b'mixmod'], dtype='|S6'), 'pdf_version': array([0])}\n"
     ]
    }
   ],
   "source": [
    "# Confirm that we have read the ensembles\n",
    "print(\"Ensemble = \", ens)\n",
    "# Print some simple information about the ensemble\n",
    "print(\"Rep      = \", ens.gen_class.name)\n",
    "print(\"NPDF     = \", ens.npdf)\n",
    "print(\"Metadata = \", ens.metadata())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting Example\n",
    "\n",
    "Now we are going to plot some PDFs from the ensemble\n",
    "Note that the first call to the `plot` specifies the x-axis limits, but does not specify the key (i.e., which PDF in the ensemble), so that defaults to 0 (i.e., the first PDF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 1300 is out of bounds for axis 0 with size 100",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m axes \u001b[38;5;241m=\u001b[39m ens\u001b[38;5;241m.\u001b[39mplot(xlim\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m0.\u001b[39m, \u001b[38;5;241m0.5\u001b[39m), label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPDF 0\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m      2\u001b[0m _ \u001b[38;5;241m=\u001b[39m ens\u001b[38;5;241m.\u001b[39mplot(key\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m, axes\u001b[38;5;241m=\u001b[39maxes, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPDF 1\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m----> 3\u001b[0m _ \u001b[38;5;241m=\u001b[39m \u001b[43mens\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1300\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxes\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43maxes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlabel\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mPDF 1300\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m      4\u001b[0m legend \u001b[38;5;241m=\u001b[39m axes\u001b[38;5;241m.\u001b[39mfigure\u001b[38;5;241m.\u001b[39mlegend()\n",
      "File \u001b[0;32m~/code/desc-rail/forked/qp/src/qp/ensemble.py:620\u001b[0m, in \u001b[0;36mEnsemble.plot\u001b[0;34m(self, key, **kwargs)\u001b[0m\n\u001b[1;32m    612\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mplot\u001b[39m(\u001b[38;5;28mself\u001b[39m, key\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[1;32m    613\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Plot the pdf as a curve\u001b[39;00m\n\u001b[1;32m    614\u001b[0m \n\u001b[1;32m    615\u001b[0m \u001b[38;5;124;03m    Parameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    618\u001b[0m \u001b[38;5;124;03m        Which PDF or PDFs from this ensemble to plot\u001b[39;00m\n\u001b[1;32m    619\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 620\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_gen_class\u001b[38;5;241m.\u001b[39mplot(\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n",
      "File \u001b[0;32m~/code/desc-rail/forked/qp/src/qp/ensemble.py:65\u001b[0m, in \u001b[0;36mEnsemble.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m     63\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m k, v \u001b[38;5;129;01min\u001b[39;00m md\u001b[38;5;241m.\u001b[39mitems():\n\u001b[1;32m     64\u001b[0m     red_data[k] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39msqueeze(v)\n\u001b[0;32m---> 65\u001b[0m dd \u001b[38;5;241m=\u001b[39m \u001b[43mslice_dict\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mobjdata\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     66\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m k, v \u001b[38;5;129;01min\u001b[39;00m dd\u001b[38;5;241m.\u001b[39mitems():\n\u001b[1;32m     67\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(np\u001b[38;5;241m.\u001b[39mshape(v)) \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m2\u001b[39m:\n",
      "File \u001b[0;32m~/code/desc-rail/forked/qp/src/qp/dict_utils.py:130\u001b[0m, in \u001b[0;36mslice_dict\u001b[0;34m(in_dict, subslice)\u001b[0m\n\u001b[1;32m    128\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m key, val \u001b[38;5;129;01min\u001b[39;00m in_dict\u001b[38;5;241m.\u001b[39mitems():\n\u001b[1;32m    129\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m--> 130\u001b[0m         out_dict[key] \u001b[38;5;241m=\u001b[39m \u001b[43mval\u001b[49m\u001b[43m[\u001b[49m\u001b[43msubslice\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m    131\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m (\u001b[38;5;167;01mKeyError\u001b[39;00m, \u001b[38;5;167;01mTypeError\u001b[39;00m):\n\u001b[1;32m    132\u001b[0m         out_dict[key] \u001b[38;5;241m=\u001b[39m val\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1300 is out of bounds for axis 0 with size 100"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAG5CAYAAABIqqroAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJFElEQVR4nO39eZxU5Z33/7+ruqurequi927oBpqtm31plE1Ek4iazWVyizFB85u4zZ1FwmOSGx86CWSSMGbibTQRo5mo0YxIjDo6d5JvIBMXECSKjYIgIFuz9L7vS9X5/XGql6IKpPv0crr79XzMmeo+ddWpqziBevu5rnMdh2EYhgAAANBnzqHuAAAAwHBHoAIAALCIQAUAAGARgQoAAMAiAhUAAIBFBCoAAACLCFQAAAAWRQ91B4arQCCgs2fPKjExUQ6HY6i7AwAALoJhGKqvr9fYsWPldPZfXYlA1Udnz55VTk7OUHcDAAD0walTp5Sdnd1vxyNQ9VFiYqIk84R4vd4h7g0AALgYdXV1ysnJ6foe7y8Eqj7qHObzer0EKgAAhpn+nq7DpHQAAACLCFQAAAAWEagAAAAsIlABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWMTNkTGsBAKGOgKG/AFDfsN8DAR/NgzJMAwFDMmQ+SiZ+87ldDjkcEgOOeR0mDfJdDqkKKdDDodDUU6Hop3dj/19E00AwMhCoEKvGYah5na/Glo71NTqV1ObX01tHV2Pze1+NbcF1NLuV3O7Xy1dm7mvtcN8bPMH1NoeUJs/oLYOc2v3m7+3+wNq9xtq9wfU4TeDU3sgoAjZaFA4HVK006noKIdcUU65ohyKdjrlijZ/j4lyyh3tlCvKKbfLKXd0lNzRzuAWJbfLKY8rKrg55YmOUmxMlGJd5mNc8Oe4mGjFu7sfY11RhDkAGAYIVKOMYRhqaO1QbXN711bX3KG6lnbVNberrqVD9S3tqm/pUENLhxpaO1Tf2qGGlnY1tHaosdWvxraOIQs2F+JwBCtPMh9l/l8YI/j/DJlVrYDRXc06n4AhM/j5Jcnf310/L4dDig+Gq3h3tBJ6bp5oeT0uJXqig5tLXo9L3tho+WI7f3bJ64lWdBSj+wAwkAhUw1iHP6DqpnZVNbapqrFN1U3mY01Tm6oa21XT1Kaa5nZVN7WptqldNcEA5f+k9HCRHA4pzhWlOHe04mLMqkpcsNricZnVFU+wMhPripLbZVZtPK7u6k1MdHc1Jya4RTuDVZ8eP0dHmcNvLqdTzh7DcU5H56MsV3IMo3soMRCQOaToN9QRCMgfHGrs8Btq8wfUEQh0/dze0V1la/cH1NrRY+tRketZpWvpCKg5WM1ravOruc2s5jW2+tXc1qFGM7nJMKSGVjPYSq19/myJ7mh5Y13yxbo0Js6lpLiYsMfk+BglxccoOS5GSfEuJbijqY4BwEUiUNlMS7tfFQ2tKq83t4qGNlU0tKqioVWVDW0qb2hVZUOrGZya2/tcKYqJdsoX/IL1xZpVjs7qRmJn1cNt/txZDUlwRyvebVZLEtzR8kRHyekcOV+4DodD0VEOW/ylCASMYMAyw1VjMFQ1tHSosa1D9S2dm1k57Py5Z7Wxtrm9K5jVByuNZ2qaL7oPMVFOJcebQSslIUYp8TFKSXArJSFGqQlupQV/Tkt0KyXerZhoqmAARi87fHeMCi3tfpXVtaq0vkUltS0qrWtRWX2ryoKPpXUtKq9vVV1LR6+O63BIY2JdPSoL3Y9m5cGlMXExGhMbfIwzA5THFTVAnxT9wel0BMOrtb+iHf6A6lo6VNPUptrmYJWyyaxa1jSZVczq4O/VTW2qbjQrns3BOW4ldS0qqWu5qPcaE+dSWoJbaYnmlp7oVobXE/zZowyvW5k+j+Ji+GcHwMgzLP5l27Rpk/793/9dxcXFmjlzpn7+859r+fLln/i6t956SytWrNCsWbO0d+/ekOdefPFF/cu//IuOHj2qyZMn68c//rFuuOGGPvWvpd2v4toWFdc062xti0pqOx9bVBz8vbqp/aKPFxPlVFqiW6mJbqUFqwGpPSoDZrXAbQ7RxLmYH4Pziu5RZeqN5ja/KhvNSmhlY5uqGtpU2WhWSTurppWNraqoN3/uCBjBgNauI2UNFzx2ojta6cFwleH1KMvnUabXo0xfrPmzz6OU+BiGGwEMK7YPVFu2bNGaNWu0adMmLVu2TI8//riuvfZaHThwQOPHjz/v62pra3Xrrbfq05/+tEpLS0Oe27Vrl1atWqV//dd/1Q033KCXX35ZN910k3bs2KFFixb1qn+XP/A31fhdF9U2JtppfnF4PUr3mv/13vlf8eld/1XvkTeWuSsYWrExUcqOiVN2Utwntg0EDNU2t6s8OFRdVm9WW82KbGgVtqnNbw4/lnfoaHnjeY8ZE+3sClpjx8Rq7JjOx1iNCz4mWKzeAUB/chiRFumxkUWLFmnBggV67LHHuvZNnz5d119/vTZu3Hje1918882aOnWqoqKi9F//9V8hFapVq1aprq5Of/7zn7v2XXPNNUpKStLmzZsjHq+1tVWtrd2Tguvq6pSTk6OcNb+X0x2nWFdU1z/65n9lm4+d/8Wd6fXIF+siKGFUa2jtUGldi0prW7qGE0s7K7l15mNFQ+tFzQ30xbo0bkysspNiNS7JDFo5yXHKTjIfvZ6L+w8dAKNLXV2dfD6famtr5fV6++24tv5PvLa2Nu3Zs0fr1q0L2b9y5Urt3LnzvK976qmndPToUf3ud7/Tj370o7Dnd+3ape985zsh+66++mr9/Oc/P+8xN27cqA0bNoTt/8PdSzQtJ52wBFyEBHe0EtISNDkt4bxt2joCKg2Gq+LaZp2paVZxTYvO1pg/n61pVl1L99IfB4rrIh7HF+tSdlKsxifHaXxynLKDj+OT4zRuTCyT6AH0K1sHqoqKCvn9fmVkZITsz8jIUElJScTXHDlyROvWrdP27dsVHR3545WUlPTqmJJ07733au3atV2/d1ao8rO88sb1bn4KgPOLiXYqJzlOOcnnH26sb2nXmZpmnak2Q9bp6madrm7SmepmnapuVlVjW1fg+vBseOByOqQsX6wmpMRpQkqcxifHa2JKnCakxGtCSpzliwEAjD7D4l+Ncys/hmFErAb5/X7dcsst2rBhg6ZNm9Yvx+zkdrvldrt70WsAAyXR41J+pkv5mZHL9Y2tHTpd3axTVU06Vd2koqomnaoyfy+qalJzu98MZDXN2nm0Muz1aYluTUyJ08SUeE1MjVduanzw5ziuUgQQka3/ZUhNTVVUVFRY5aisrCyswiRJ9fX1evfdd1VYWKhvfvObkqRAICDDMBQdHa2tW7fqU5/6lDIzMy/6mACGn3h3tPIyE5WXmRj2nGEYKm9oVVGlGa5OVjbpZGWjTgZ/rmps61oH7p0T1WGvz/J5lBsMWZPSEjQpNV6T0uKVnRSnqBG0LhuA3rF1oIqJiVFBQYG2bdsWsqTBtm3bdN1114W193q92rdvX8i+TZs26W9/+5v+8Ic/KDc3V5K0ZMkSbdu2LWQe1datW7V06dIB+iQA7MLhcCg90aP0RI8WTkwOe762uV1FlU06XtmoExXmdqyiUScqG1XT1B6c29USVtmKiXZqYkqcJgfniE1Oj9eUtERNSotnCBEYBWz/t3zt2rVavXq1Fi5cqCVLluiJJ55QUVGR7r77bknm3KYzZ87omWeekdPp1KxZs0Jen56eLo/HE7L/nnvu0eWXX64HHnhA1113nV555RX99a9/1Y4dOwb1swGwH1+sS7OzfZqd7Qt7rrqxTccrG3WsvFHHKxp0vML8+VhFo9o6Ajpc2qDDpeHrcI31eTQ5PUFT0hM0NT1RUzMSNCUtQUm9XB8MgH3ZPlCtWrVKlZWV+uEPf6ji4mLNmjVLf/rTnzRhwgRJUnFxsYqKinp1zKVLl+r555/X/fffr3/5l3/R5MmTtWXLll6vQQVgdEkK3u9wwfikkP3+gKGzNc36uLxBR8sadKyiUR+XNehYeYMqGtp0trZFZ2tbtP1IRcjrUhNiugLW1IxETUtP0LSMRIIWMAzZfh0quxqodSwAjCzVjW06Wt6go+UNOlLaoCNlDfq4rOGC91VMS3QrLyNR0zISlZdphqxpGYkMHQL9YKC+vwlUfUSgAmBFY2uHjpabQ4RHSut1uLReh0svHLRykmOVl+FVfnDC/fSsRE1Mief2U0AvEKhshkAFYCA0tnboSFmDDpfU61AwaB0qqVdZfWvE9jHRTk3LSFB+phm0ZmR5lZ/l7fX9G4HRgkBlMwQqAIOpqrFNh0rqdaikTodK6/VRiRm0mtr8EdtneN2anuXVjCyvZow1HyemxMvJ0g4Y5QhUNkOgAjDUAgFDp6qbdLC4XgeL6/RRSZ0OFterqKopYvu4mCizijXWq5ljfZo51qtpGYnyuKIGuefA0CFQ2QyBCoBd1be061CJGbIOFNfrQHGdPiquU2tHIKxttNOhKekJmjnWp1njvJo1zqfpWV4lMAEeIxSBymYIVACGkw5/QCcqG/Xh2TodOFunD8/W6cOztapuag9r63BIuSnxmjXOp9njfJo1zqeZ47zyelxD0HOgfxGobIZABWC4MwxDxbUtXeFq/xnzsbi2JWL73FQzZM0ZZy58OnOsV4mELAwzBCqbIVABGKkqGlq1/0ytPjxbp32na7XvTG3E5RwcDjNkzRnn05zsMZqb49PMsT7mZMHWCFQ2Q6ACMJpUNbZp35la7T9Tqw9O12jf6VqdjVDJinI6NC0jUXOzfZqbM0Zzs8doWkYCa2XBNghUNkOgAjDaVTS0at/pWn1w2gxZ75+uVUVD+HpZHpdTs8f5NDd7jObmjNG8nDHKToqVw8ESDhh8BCqbIVABQCjDMFRS16L3T3UGrBp9cKpW9a0dYW1TE2I0N9sMV/PHJ2lOjo9J7xgUBCqbIVABwCcLBAwdr2zU3iIzYO09VaODxXVq94d+9Tgc0pS0hK6AtWDCGE1NT1QUC5GinxGobIZABQB909Lu14HiOu0tMgNW4alqnaoKn/Se4I7W3Byf5ueYAWt+TpKSuKUOLCJQ2QyBCgD6T0VDq/YWmeGqsKhG75+qUWOE2+pMSo3X/PFJKphgblPTE7idDnqFQGUzBCoAGDj+gKHDpfUqLKrRe0XVeq+oWsfKG8PaJXqizYAVDFnzxo9hlXdcEIHKZghUADC4qhvbtPdUjfacNAPW3lM1YTeHdjqk/EyvFk40A9YlE5M1dkzsEPUYdkSgshkCFQAMrQ5/QB+V1Ou9omrtOVmtd09UR1yAdNyY2GC4SlLBhGTlZTLZfTQjUNkMgQoA7KektkXvnqzSuyfMKtaHZ+vkD4R+zSW6o7VgQpIuzU3WwglJmpszhtXdRxEClc0QqADA/hpbO/T+qRq9c6Ja756s0nsnq8Mmu8dEOTU726dLJiZrUW6yFkxIki+WNbFGKgKVzRCoAGD46RwmfOeEWcX6+4kqldeHru7uCM7DWpSbrEtzk3XJxGSlJbqHqMfobwQqmyFQAcDwZxiGiqqatPt4ld45XqV3TlTpRGVTWLtJafFdAevS3BSNY6L7sEWgshkCFQCMTGV1LXrnRLXeOVGl3cer9FFJnc79psxOitWi3BQtyk3WoknJGp8cx70JhwkClc0QqABgdKhtate7J81wtft4lfafqQ2b6J7l82hRbrIWT0rRokkpmphCwLIrApXNEKgAYHRqaO3QnpPV+vvxSu0+VqX3T9eE3Zsww+vW4kkpXRsByz4IVDZDoAIASFJzm1/vFVVr97FKvX28SnuLatTmD4S06QxYSyalaMnkFIYIhxCBymYIVACASFra/SosqtHbxyq161hlxIA11ufR4sndASs7KW6Iejv6EKhshkAFALgYLe1+vXeyujtgnQofIhyfHKclk1K0dIoZstK9niHq7chHoLIZAhUAoC+a2sw5WLuOmgHrg9Phk9ynpCdo6eQULZ1szsEaExczRL0deQhUNkOgAgD0h4bWDr1zvEq7jlVq59EKfXg2dJkGh0OakeXVsimpWjI5RZdOTFa8O3roOjzMEahshkAFABgItU3tevt4pXYdrdRbH1foSFlDyPPRTofmjx+jpZNTtWxKqubljFFMtHOIejv8EKhshkAFABgMZfUt2nW0Ujs/rtRbRyt0uro55Pm4mChdmpusy6aYASsvI1FOJ1cQng+BymYIVACAoVBU2aS3jlborY8rtOtopSob20KeT02I0dLJqbpsSqqWTuEKwnMRqGyGQAUAGGqBgKGPSuq182iFdnxcob8fr1JTmz+kzcSUOF021QxYSyanyhfrGqLe2gOBymYIVAAAu2nrCGjvqRrt+NisYO09VRNyBaHTIc3OHqPlU1J12dRULRifNOrmXxGobIZABQCwu/qWdr19rEpvfVyh7UfKdbS8MeT5uJgoLcpN1mVT03T51FRNSU8Y8Su4E6hshkAFABhuimubteNIRVcFq6IhdP5Vptejy6amanlwiDAlwT1EPR04BCqbIVABAIazQMDQwZK6roD19+NVau0IvUXOzLFeLQ9WrwomJskdHTVEve0/BCqbIVABAEaSlna/3jlRpe1HKrT9SIUOFteFPB/ritLiSclmwJqWpslp8cNyeJBAZTMEKgDASFZe36odH5dr++EKvXmkQhUNrSHPjxsTq+VTU3X5tDQtm5wqX9zwuHqQQGUzBCoAwGhhGObyDG8eLtf2I+bwYJu/e3jQ6ZDm5YzR5dPM6tXc7DGKsuniogQqmyFQAQBGq+Y2v94+XhmsXpXr43Nuj+OLdemyqalaERwezPR5hqin4QhUNkOgAgDAdKamWW8eLtebh8u14+MK1bd0hDyfl5GoFXlpWjEtTQuHeHI7gcpmCFQAAITr8Af0/ukavXG4Qm8cLtcHp2vUM2nEuqK0dHKKrshL04pp6RqfMri3xiFQ2QyBCgCAT1bV2KYdH1fo9UNlevNw+OT23NR4rZiWphV5aVqcm6LYmIGtXhGobIZABQBA73SuffXG4XK9cahce05Wq6PHrXFiop1alJusK/LSdUVemial9v/SDAQqmyFQAQBgTX1Lu976uDIYsMp0trYl5Pmc5FhdMc0MV0smpyguJtryexKobIZABQBA/zEMQx+XNej1Q+V6/XCZ3jleHbI0Q2f1asW0NF2Zn97n6hWBymYIVAAADJzG1g7tOlqp1w+X6bWPynWmpjnk+ZzkWF2Zl64r89K1eNLFz70iUNkMgQoAgMFhGIaOlgerV4fKwxYWdUc7tXhSiq7MM6tXE1Liz3ssApXNEKgAABgaja0d2nm0Uq8fKtPrh8KrV5NS43VFXrquzE/TpbnJIeteEahshkAFAMDQ65x79dohc2jwnRNVIVcOxrqitGxKqq7MT9OVeelKcLYTqOyEQAUAgP2YVw5W6LWPyvXaoTKV1YeuezVljFP/c+9n+/372/r1hwAAADaR6HHpmllZumZWlgzD0IHiOr1+qFx/+6hMhUXVOlza8MkH6QMqVH1EhQoAgOGlurFN/1/hMd2yfHq/f387++1IAAAANpYUH6PPzRk7IMcmUAEAAFhEoAIAALCIQAUAAGARgQoAAMCiYRGoNm3apNzcXHk8HhUUFGj79u3nbbtjxw4tW7ZMKSkpio2NVX5+vh566KGQNk8//bQcDkfY1tLScp6jAgAAnJ/t16HasmWL1qxZo02bNmnZsmV6/PHHde211+rAgQMaP358WPv4+Hh985vf1Jw5cxQfH68dO3borrvuUnx8vO68886udl6vV4cOHQp5rcfjGfDPAwAARh7br0O1aNEiLViwQI899ljXvunTp+v666/Xxo0bL+oYN954o+Lj4/Xss89KMitUa9asUU1NTZ/7xTpUAAAMPwP1/W3rIb+2tjbt2bNHK1euDNm/cuVK7dy586KOUVhYqJ07d2rFihUh+xsaGjRhwgRlZ2fr85//vAoLCy94nNbWVtXV1YVsAAAAks0DVUVFhfx+vzIyMkL2Z2RkqKSk5IKvzc7Oltvt1sKFC/WNb3xDt99+e9dz+fn5evrpp/Xqq69q8+bN8ng8WrZsmY4cOXLe423cuFE+n69ry8nJsfbhAADAiGH7OVSS5HA4Qn43DCNs37m2b9+uhoYGvf3221q3bp2mTJmiL3/5y5KkxYsXa/HixV1tly1bpgULFugXv/iFHnnkkYjHu/fee7V27dqu3+vq6ghVAABAks0DVWpqqqKiosKqUWVlZWFVq3Pl5uZKkmbPnq3S0lKtX7++K1Cdy+l06pJLLrlghcrtdsvtdvfyEwAAgNHA1kN+MTExKigo0LZt20L2b9u2TUuXLr3o4xiGodbW1gs+v3fvXmVlZfW5rwAAYPSydYVKktauXavVq1dr4cKFWrJkiZ544gkVFRXp7rvvlmQOxZ05c0bPPPOMJOnRRx/V+PHjlZ+fL8lcl+pnP/uZvvWtb3Udc8OGDVq8eLGmTp2quro6PfLII9q7d68effTRwf+AAABg2LN9oFq1apUqKyv1wx/+UMXFxZo1a5b+9Kc/acKECZKk4uJiFRUVdbUPBAK69957dfz4cUVHR2vy5Mn6t3/7N911111dbWpqanTnnXeqpKREPp9P8+fP15tvvqlLL7100D8fAAAY/my/DpVdsQ4VAADDz6hchwoAAGA4IFABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYBGBCgAAwCICFQAAgEUEKgAAAIsIVAAAABYRqAAAACwiUAEAAFhEoAIAALCIQAUAAGARgQoAAMAiAhUAAIBFBCoAAACLCFQAAAAWEagAAAAsIlABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYBGBCgAAwCICFQAAgEUEKgAAAIsIVAAAABYRqAAAACwiUAEAAFhEoAIAALCIQAUAAGARgQoAAMAiAhUAAIBFBCoAAACLCFQAAAAWEagAAAAsIlABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYNGwCFSbNm1Sbm6uPB6PCgoKtH379vO23bFjh5YtW6aUlBTFxsYqPz9fDz30UFi7F198UTNmzJDb7daMGTP08ssvD+RHAAAAI5jtA9WWLVu0Zs0a3XfffSosLNTy5ct17bXXqqioKGL7+Ph4ffOb39Sbb76pgwcP6v7779f999+vJ554oqvNrl27tGrVKq1evVrvv/++Vq9erZtuukm7d+8erI8FAABGEIdhGMZQd+JCFi1apAULFuixxx7r2jd9+nRdf/312rhx40Ud48Ybb1R8fLyeffZZSdKqVatUV1enP//5z11trrnmGiUlJWnz5s0Rj9Ha2qrW1tau3+vq6pSTk6Pa2lp5vd6+fDQAADDI6urq5PP5+v3729YVqra2Nu3Zs0crV64M2b9y5Urt3Lnzoo5RWFionTt3asWKFV37du3aFXbMq6+++oLH3Lhxo3w+X9eWk5PTi08CAABGMlsHqoqKCvn9fmVkZITsz8jIUElJyQVfm52dLbfbrYULF+ob3/iGbr/99q7nSkpKen3Me++9V7W1tV3bqVOn+vCJAADASBQ91B24GA6HI+R3wzDC9p1r+/btamho0Ntvv61169ZpypQp+vKXv9znY7rdbrnd7j70HgAAjHT9FqgMw1BFRYXKy8vV3Nys1NRUpaWlKS4urs/HTE1NVVRUVFjlqKysLKzCdK7c3FxJ0uzZs1VaWqr169d3BarMzMw+HRMAACASS0N+R44c0Y9+9COtXLlSXq9XmZmZmj17ti699FJNmjRJiYmJys/P1x133KEXXnhB7e3tvTp+TEyMCgoKtG3btpD927Zt09KlSy/6OIZhhEwoX7JkSdgxt27d2qtjAgAAdOpTheqFF17QL3/5S+3YsUOSGVgkyel0yufzKTY2VlVVVWppadHhw4d1+PBhPfnkk0pOTtatt96qtWvXaty4cRf1XmvXrtXq1au1cOFCLVmyRE888YSKiop09913SzLnNp05c0bPPPOMJOnRRx/V+PHjlZ+fL8lcl+pnP/uZvvWtb3Ud85577tHll1+uBx54QNddd51eeeUV/fWvf+36PAAAAL1i9MJf//pXY+HChYbT6TQcDocxb9484/777zdeffVVo7i42AgEAiHtW1pajD179hiPPfaY8ZWvfMXwer2Gw+Ew4uLijHXr1hk1NTUX9b6PPvqoMWHCBCMmJsZYsGCB8cYbb3Q9d9tttxkrVqzo+v2RRx4xZs6cacTFxRler9eYP3++sWnTJsPv94cc84UXXjDy8vIMl8tl5OfnGy+++GJv/iiM2tpaQ5JRW1vbq9cBAIChM1Df371ah6qzAvVP//RPuu2225SXl9er8Nba2qr//u//1i9+8Qtt375d69ev1/e///1eRkB7GKh1LAAAwMAZqO/vXg35bdiwQd/+9rfl8/n69GZut1tf+tKX9KUvfUnbt29XTU1Nn44DAABgJ7ZfKd2uqFABADD8jMqV0gEAAIYDy+tQ/fCHP5QkffnLX9bUqVMv2Pbpp59WUVHRsJ03BQAAEInlIT+n0ymHwyGfz6fnn38+7B55PS1fvlw7d+6U3++38pa2wJAfAADDj62H/Dwej2pqavS5z31ODz74YH8cEgAAYNjol0BVUFCg5557Ti6XS9/73ve0evXqkJXJAQAARrJ+m5R+8803a8eOHRo3bpyee+45XX755Tp79mx/HR4AAMC2+vUqvwULFmjPnj1atmyZ3nnnHV1yySV6++23+/MtAAAAbKffl01IS0vT3/72N915550qLi7WlVdeqSeffLK/3wYAAMA2LC+bEPGg0dH61a9+pblz52rNmjW64447tHfvXnV0dAzE2wEAAAypAQlUnf7pn/5Js2bN0pe+9CU9+uijA/lWAAAAQ2bAV0pfvny53nnnHc2dO1fc5QYAAIxElitUgUDgE9uMHz9eO3fu1O7du62+HQAAgO0M6JBfTx6PRytWrBistwMAABg03BwZAADAol4Fqs9+9rP9MmzX2Niof/u3f9OmTZssHwsAAGCo9SpQbd++XUuXLtVnPvMZPfvss6qvr+/VmxUWFuqf//mfNWHCBN13333cngYAAIwIvZpDdezYMa1fv17/8R//oddee01ut1uXXXaZLr30UhUUFCgrK0vJyclyu92qqalRVVWVDh48qHfffVc7duzQ0aNHZRiGpk+frieffFJf/OIXB+pzAQAADBqH0Ye1DI4fP65f/epX+u1vf6uysjLzQA7HedsbhiGHw6FPfepTuvPOO/UP//APcjqH9/Sturo6+Xw+1dbWyuv1DnV3AADARRio7+8+BapOHR0deu211/Tmm29q586dOnnypCoqKtTS0qLk5GSlp6dr3rx5uuyyy3TVVVdpwoQJ/dbxoUagAgBg+LFloBrNCFQAAAw/A/X9PbzH3QAAAGyAQAUAAGBRv62U3traqueff15/+ctfdPjwYdXX1ysxMVHTpk3TypUrdfPNN8vj8fTX2wEAANhGv8yh2rlzp7761a/q5MmTEW+A7HA4NH78eP3ud7/TsmXLrL6dLTCHCgCA4Wegvr8tV6g+/PBDXXXVVWpublZmZqZuv/12TZ8+XRkZGSorK9PBgwf1m9/8RidPntTKlSu1e/duzZo1qz/6DgAAYAuWK1Q33HCDXnnlFX31q1/Vb37zG7lcrrA27e3tuv322/Xss8/q+uuv10svvWTlLW2BChUAAMOPbZdNSElJkd/vV0lJyQXnSLW0tCgzM1NOp1NVVVVW3tIWCFQAAAw/tl02oa2tTXl5eZ844dzj8SgvL0/t7e1W3xIAAMBWLAeq6dOn6/Tp0xfV9tSpU5o5c6bVtwQAALAVy4FqzZo1Ki4u1sMPP3zBdo888ohKSkq0Zs0aq28JAABgK5av8rvlllt05swZ/Z//83/0xhtv6H//7/+t6dOnKz09XeXl5Tp48KA2bdqkP/7xj/rpT3+qm2++uT/6DQAAYBuWJ6VHRUVZ74TDoY6ODsvHGUxMSgcAYPix7TpU/XFvZe7PDAAAhjPLgSoQCPRHPwAAAIYtbo4MAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYBGBCgAAwCICFQAAgEUEKgAAAIsIVAAAABYRqAAAACwiUAEAAFhEoAIAALCIQAUAAGARgQoAAMAiAhUAAIBFBCoAAACLCFQAAAAWEagAAAAsIlABAABYNCwC1aZNm5SbmyuPx6OCggJt3779vG1feuklXXXVVUpLS5PX69WSJUv0l7/8JaTN008/LYfDEba1tLQM9EcBAAAjkO0D1ZYtW7RmzRrdd999Kiws1PLly3XttdeqqKgoYvs333xTV111lf70pz9pz549uvLKK/WFL3xBhYWFIe28Xq+Ki4tDNo/HMxgfCQAAjDAOwzCMoe7EhSxatEgLFizQY4891rVv+vTpuv7667Vx48aLOsbMmTO1atUqff/735dkVqjWrFmjmpqaPverrq5OPp9PtbW18nq9fT4OAAAYPAP1/W3rClVbW5v27NmjlStXhuxfuXKldu7ceVHHCAQCqq+vV3Jycsj+hoYGTZgwQdnZ2fr85z8fVsE6V2trq+rq6kI2AAAAyeaBqqKiQn6/XxkZGSH7MzIyVFJSclHHePDBB9XY2Kibbrqpa19+fr6efvppvfrqq9q8ebM8Ho+WLVumI0eOnPc4GzdulM/n69pycnL69qEAAMCIY+tA1cnhcIT8bhhG2L5INm/erPXr12vLli1KT0/v2r948WJ99atf1dy5c7V8+XL9/ve/17Rp0/SLX/zivMe69957VVtb27WdOnWq7x8IAACMKNFD3YELSU1NVVRUVFg1qqysLKxqda4tW7bo61//ul544QV95jOfuWBbp9OpSy655IIVKrfbLbfbffGdBwAAo4atK1QxMTEqKCjQtm3bQvZv27ZNS5cuPe/rNm/erK997Wt67rnn9LnPfe4T38cwDO3du1dZWVmW+wwAAEYfW1eoJGnt2rVavXq1Fi5cqCVLluiJJ55QUVGR7r77bknmUNyZM2f0zDPPSDLD1K233qqHH35Yixcv7qpuxcbGyufzSZI2bNigxYsXa+rUqaqrq9MjjzyivXv36tFHHx2aDwkAAIY12weqVatWqbKyUj/84Q9VXFysWbNm6U9/+pMmTJggSSouLg5Zk+rxxx9XR0eHvvGNb+gb3/hG1/7bbrtNTz/9tCSppqZGd955p0pKSuTz+TR//ny9+eabuvTSSwf1swEAgJHB9utQ2RXrUAEAMPyMynWoAAAAhgMCFQAAgEUEKgAAAIsIVAAAABYRqAAAACwiUAEAAFhEoAIAALCIQAUAAGARgQoAAMAiAhUAAIBFBCoAAACLCFQAAAAWEaisOvm21Now1L0AAABDKHqoOzDsPfclyRMlpeZJ4xZIY+ebjxmzpGj3UPcOAAAMAgKVVQmZUnupVH7Q3Pb+p7nf6ZIyZnYHrLELpLR8KYo/cgAARhqHYRjGUHdiOKqrq5PP51Ntba28apTOFkpn35POvGc+NleHvyg6Vsqc3V3JGjtfSpkiOaMG/wMAADAKhXx/e739dlwCVR9d8IQYhlRzsjtcnd0rFb8vtdaFHygmQcqa2x2wsuZJyZMkJ9PbAADobwQqm+n1CQkEpKqjZiXrzHtS8V4zZLU3hbd1e3uErHnmY1Ku5HD098cAAGBUIVDZTL+ckIBfqjgcGrJK9kkdLeFtPT6zejV2XvBxvpQ0kZAFAEAvEKhsZqBOiPztUvlHwTlZheZwYel+yd8W3tYzJljJmtcdtqhkAQBwXgQqmxmwQBVJR5t5BeHZvWbIKt4rlX54npDlM0NW1rzuYcOkXOZkAQAgApXtDGqgiqSjTSo7YIars3svHLLcXilzTrCSFQxbKZO5uhAAMOoQqGxmyANVJB1t5nBhz5BVsl/yt4a3dcVLWXNCq1mp01gnCwAwohGobMaWgSoSf7tUfqhHyHo/OPG9ObxttMdc4b1rXtZcKW26FB0zyJ0GAGBgEKhsZtgEqkj8HVLlETNcdYWsD6S2CPckdLqkjBnBStZcKXOuuQJ8TNygdxsAAKsIVDYzrANVJIGAVHUsWMkqNANW8ftSS214W4fTvHdh1tzuYcPM2eaEeAAAbIxAZTMjLlBF0rnie/H7we0DM3A1lkdunzQxGK7mBOdlzZES0gexwwAAXBiBymZGRaCKxDCk+pLuClZn0Kotitw+IbO7kpU5x3wcM4G1sgAAQ4JAZTOjNlCdT1NVj5AVfKz8WFKE/3l5fMFwFaxmZc7mCkMAwKAgUNkMgeoitDaYa2MVvy+VBINW2UEp0B7eNtojpc8IVrJmM/kdADAgBur7m5IABo47QRq/yNw6da76XvyBWdEq2WdubQ3S2ffMrZPDKaVMNQNW55Bh5hwpPmXwPwsAABdAhaqPqFD1o0BAqj7evUZWyQdm4Gosi9zeO657qLCzosW8LADARWDIz2YIVIOgvqRHJStYzao6Frmtxydl9AhYmbPNpR1YlBQA0AOBymYIVEOkpU4q3R9ayTrfvCynS0rPN+djdYaszFmslwUAoxhzqABJ8nilCUvNrVPnPQw7g1ZxsJrVWts9R6unMROC4apHNcuXzZAhAKDPqFD1ERUqm+tclLRkf3eoKvlAqj0Vub1nTHjISsuTolyD2m0AwMBiyM9mCFTDVFNVjyHD4Fb+kRToCG8bFWOGqsw55k2jO4cMY5MGv98AgH5BoLIZAtUI0tFqhqqeIatkn9RaF7m9L6fHnKzZZthKmsiQIQAMAwQqmyFQjXCGIdUUhYes891ix+0NVrFmdYes9OmSK3Zw+w0AuCAClc0QqEap5mpz9feSHlcaln8k+dvC2zqc5i11QoLWbCkxY/D7DQCQRKCyHQIVuvjbpYoj3QGrZJ85T6upMnL7+HQzYGXMCk6Cn2WuCM+9DAFgwBGobIZAhQsyDHNh0tL9wZC13/y54ogi3jA6ym0OEWbOMqtYmbPNexnGjhnsngPAiMY6VMBw4nBI3ixzm3pV9/62JnMh0pIPgmErGLTaGqTivebWk298j2pWcNhwzETJ6RzEDwMA+CRUqPqIChX6TSAg1ZwIDhnu7w5a55sAH5NgVq86Q1bGbCljhhQTP6jdBoDhiCE/myFQYcD1nABfGgxbZQclf2uExg4peVKPIcNgNcs7juUcAKAHApXNEKgwJPwdUuWR0JBVul9qKI3cvnMF+J4VrbTpksszqN0GALsgUNkMgQq20lAWOierZL9UcSjyCvCOKCl1auiQYeYsKSGDahaAEY9AZTMEKtheR6tUfqhH0ApWtJqrIrePS+0xAT64OCn3MwQwwhCobIZAhWHJMKT64uDCpD2uNKz8WBGXc3C6pLT80CsNM2ZL8SmD3nUA6A8smwDAOodD8o41t2kru/d3LudQuj906LC1zqxsle4LPU5iVo95WcFqVsoUFicFMGpRoeojKlQY8TrvZ3jukGH18cjtoz09qlnBeVkZM6XYpMHtNwBcAEN+NkOgwqjVWi+VHgi9yrD0gNTeGLm9L6fHcGGwopWUy+KkAIYEgcpmCFRAD4GAWbnqvI9hZ9CqPRW5vSveXIz03MVJ3YmD228Aow6BymYIVMBFaK4OVrP2d4etsoNSR0vk9km55jBh57yszFnSmAks5wCg3xCobIZABfSRv0OqOnpONetDqf5s5PZub/itdtKnSzFxg9tvACMCgcpmCFRAP2usDL3KsGSfVP6RFGgPb+twSsmTzUpWz0nwiVlUswBcEIHKZghUwCDoaJMqDocOGZbsl5oqIrePTT7nKsNZ5pWH0TGD228AtkWgshkCFTBEDMO8d+G59zOsOCIZ/vD2zmgpNa9HNSt4pWF86uD3HcCQI1DZDIEKsJn2Fqn8YOj9DEv2Sa21kdsnZnXPy8qcbVa1UiZLzqjB7TeAQTVQ39/DYiGYTZs2KTc3Vx6PRwUFBdq+fft527700ku66qqrlJaWJq/XqyVLlugvf/lLWLsXX3xRM2bMkNvt1owZM/Tyyy8P5EcAMNBcHmnsfGnBaunaB6T/3x+ldSelNfukmzdLV94nTf+ieSWhZN6C5+Nt0o6HpD/8o/ToJdJPxklPXCm9+i3p77+WTu6SWuqG9nMBGBZsf5+ILVu2aM2aNdq0aZOWLVumxx9/XNdee60OHDig8ePHh7V/8803ddVVV+knP/mJxowZo6eeekpf+MIXtHv3bs2fP1+StGvXLq1atUr/+q//qhtuuEEvv/yybrrpJu3YsUOLFi0a7I8IYKA4HNKY8eaW/9nu/ecuTlqyTyo7ILU3SWffM7eekiYGq1lzuocNx4xnAjyALrYf8lu0aJEWLFigxx57rGvf9OnTdf3112vjxo0XdYyZM2dq1apV+v73vy9JWrVqlerq6vTnP/+5q80111yjpKQkbd68+aKOyZAfMMIE/FLV8dCbRpful+rORG7v8YXeyzBztjkB3uUZ3H4D6JVReXPktrY27dmzR+vWrQvZv3LlSu3cufOijhEIBFRfX6/k5OSufbt27dJ3vvOdkHZXX321fv7zn5/3OK2trWptbe36va6OYQBgRHFGSalTzG3Wjd37m6p6XGEYrGiVfyS11Eon3zK3To4oKS0vNGQxAR4YFWwdqCoqKuT3+5WRkRGyPyMjQyUlJRd1jAcffFCNjY266aabuvaVlJT0+pgbN27Uhg0betF7ACNCXLI0aYW5depokyoOBQPWvu7A1VxtDh2WHZC0pbt9YtY5IWuOlDyJ+xkCI4itA1UnxznzFAzDCNsXyebNm7V+/Xq98sorSk9Pt3TMe++9V2vXru36va6uTjk5ORfTfQAjTXRMd/Wpk2FIdWeD4WpfdzWr6qg5Ab6+WDqytbt95/0MO4+TOUdKn8EK8MAwZetAlZqaqqioqLDKUVlZWViF6VxbtmzR17/+db3wwgv6zGc+E/JcZmZmr4/pdrvldrt7+QkAjBoOh+QbZ25513Tvb20wb61T2rOadUBqb5ROv2NuXcfouQJ8MGRlzpYSL/zvHYChZ+tAFRMTo4KCAm3btk033HBD1/5t27bpuuuuO+/rNm/erH/8x3/U5s2b9bnPfS7s+SVLlmjbtm0h86i2bt2qpUuX9u8HAAB3gjR+kbl1Cvilyo/DhwwbSqXKI+b24Uvd7ePTQgNW5hzWzAJsxtaBSpLWrl2r1atXa+HChVqyZImeeOIJFRUV6e6775ZkDsWdOXNGzzzzjCQzTN166616+OGHtXjx4q5KVGxsrHw+nyTpnnvu0eWXX64HHnhA1113nV555RX99a9/1Y4dO4bmQwIYXZzByetpedLsL3Xvry8NrWSV7DNXgG8sl47+zdw6ueLMIcKe1ayMGVJM/OB/HgD2XzZBMhf2/OlPf6ri4mLNmjVLDz30kC6//HJJ0te+9jWdOHFCr7/+uiTpiiuu0BtvvBF2jNtuu01PP/101+9/+MMfdP/99+vYsWOaPHmyfvzjH+vGG28Me935sGwCgEHR1iiVHTSXc+iqZn1orpl1LodTSplyzpDhHCkhbfD7DdgUt56xGQIVgCET8EtVx0JDVvEHUmNZ5PadVxn2HDZMyuUqQ4xKBCqbIVABsJ360mDA+qD7sfKopAj/zMckdt/HsDNkpU+Xorn4BiMbgcpmCFQAhoXOqwx7hqzSA5K/NbytM1pKm26Gq6w53UHLw79xGDkIVDZDoAIwbPk7pIrDPapZH5hDhi01kdsnTTTDVdYcKXOu+ZiQwb0MMSwRqGyGQAVgRDEMqfZ0dyWrOBi0ak9Fbh+f1iNkzZGy5jIvC8MCgcpmCFQARoXOexl2VrFKPjCrW0YgvG3nvKysud1hKy1finINfr+B8yBQ2QyBCsCo1dZk3q+w+P3uoFV2QOpoCW8bFWOul9VVyZonZczkFjsYMgP1/W37hT0BADYTEydlLzS3Tl3zsnpUsoo/kFprpeK95tbJ4ZRSppqVrKzgcGHmbCk2abA/CdBvqFD1ERUqAPgEhiFVnwgPWQ0lkduPGR8MWXODk9/nch9D9DuG/GyGQAUAfVRfGgxXe82AVfy+VHMyctuEzNBKVtZcyZfDFYboMwKVzRCoAKAfNVcHry58vztkVRxWxEVJY5O6w1XWXHNeFlcY4iIRqGyGQAUAA6yt0VyU9OxeqeR9M2SVHZQCHeFtYxJDq1hZ86TUqeaNqIEeCFQ2Q6ACgCHQ0WpeUXh2b3DY8H2pZH/kld9dcVLGLGnsvO6gxTIOox5X+QEAEO2Wxs43t07+dnN4sGfIKv5Aam+UTv/d3DpFuaWMGWYFK2uuGbbSZ3APQ1hGhaqPqFABgI0F/OaNoYv3BgNWcGutC2/rdJk3hh47zwxaY+dJ6TMll2dw+4xBwZCfzRCoAGCYCQSk6uM9ri7ca1a1It3DsPNG0WOD87HGzjcXJHXFDmqX0f8IVDZDoAKAEcAwpJqi7krW2b3mz02V4W0dUeYcrLHzzICVNc+81Q4ha1ghUNkMgQoARqjOG0V3VrA6H5sqwtsSsoYdApXNEKgAYBQxDKnuTGjAKt4rNZaHt+0KWfO7g1bGLOZk2QSBymYIVAAwyhmGVHf2nEpWYeSQ1TUna1530MqYxdWFQ4BAZTMEKgBAmJCQVWgGrbOFkYcLna7uJRw6l4JInyFFxwxyp0cXApXNEKgAABclZE5Wj5DVXBXeNiomuBhpMGCNWyCl5klRLBvZXwhUNkOgAgD0Wc+rC88Wdm8tteFto2PN2+p0hqyxC6SUKdy7sI8IVDZDoAIA9CvDMNfJ6qxgnS08/2KkMYndK72PW2AGraRcyeEY7F4POwQqmyFQAQAGXCAgVR0NrWIVvy+1N4W39YzpHibsrGR5xxKyzkGgshkCFQBgSAT8UvmhYMB6z3ws2Sf528LbJmR0h6vOoBWfOvh9thEClc0QqAAAttHRJpUdMAPWmWDIKjsoGf7wtr7x0rj50rgCM2iNnSe5Ewe9y0OFQGUzBCoAgK21NZmVq66Q9Z5U+XGEhg4pdZpZweoMWZkjd40sApXNEKgAAMNOc033lYVngkGr7nR4O6fLvBn0uILuoJU6TXJGDXaP+x2BymYIVACAEaGhrLuC1fkY6ebQMQnmlYU9K1ljxg+7Se8EKpshUAEARiTDkGpO9ghZwasL2xvD28aldgeszpAVnzL4fe4FApXNEKgAAKNGwC9VHJbO7AkOFe6RSvdLgY7wtkkTuwPWuAIpc44UEzfoXT4fApXNEKgAAKNae4sZqk6/2z1cWHkkvJ0jyrxnYc+QlZY/ZPOxCFQ2Q6ACAOAczTXBCe/vmkOFZ96VGkrD27niuxchHVcgZS+UvOMGZT4WgcpmCFQAAHwCw5DqzgaHCt/tXiOrrSG8bUKGNG6hlN1jPpan/79fCVQ2Q6ACAKAPes7HOv1ucD7WhxEWIQ2uj5W9sHuoMGOmFOWy9PYEKpshUAEA0E/amqSSD3qErHelmqLwdtGx5tINnSEre6Hky+nVUCGBymYIVAAADKCG8u6hws6J7y214e3i00MD1icMFRKobIZABQDAIAoEpKqj3RWs0++eZ+kGh5SWFwxZC83HtOlSVLQkApXtEKgAABhi7c1S8QfBgPWOdHqPVBthqLDzqsLshaobM1O+S1f1+/d3dL8dCQAAYDC5YqXxi8ytU31pcC7WO91XFrY1SCd3mFvrwNSRCFQAAGDkSMyQ8j9rbpJ5VWH5oe4q1pHdkt7p97dlyK+PGPIDAGD4Gajvb2e/HQkAAGCUIlABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYBGBCgAAwCICFQAAgEXRQ92B4cowDElSXV3dEPcEAABcrM7v7c7v8f5CoOqjyspKSVJOTs4Q9wQAAPRWZWWlfD5fvx2PQNVHycnJkqSioqJ+PSHovbq6OuXk5OjUqVPyer1D3Z1Rj/NhH5wL++Bc2Edtba3Gjx/f9T3eXwhUfeR0mtPPfD4ffzlswuv1ci5shPNhH5wL++Bc2Efn93i/Ha9fjwYAADAKEagAAAAsIlD1kdvt1g9+8AO53e6h7sqox7mwF86HfXAu7INzYR8DdS4cRn9fNwgAADDKUKECAACwiEAFAABgEYEKAADAIgIVAACARQSqC9i0aZNyc3Pl8XhUUFCg7du3X7D9G2+8oYKCAnk8Hk2aNEm/+tWvBqmnI19vzkVxcbFuueUW5eXlyel0as2aNYPX0VGgN+fipZde0lVXXaW0tDR5vV4tWbJEf/nLXwaxtyNfb87Hjh07tGzZMqWkpCg2Nlb5+fl66KGHBrG3I1tvvzM6vfXWW4qOjta8efMGtoOjSG/Oxeuvvy6HwxG2ffTRR717UwMRPf/884bL5TJ+/etfGwcOHDDuueceIz4+3jh58mTE9seOHTPi4uKMe+65xzhw4IDx61//2nC5XMYf/vCHQe75yNPbc3H8+HHj29/+tvHb3/7WmDdvnnHPPfcMbodHsN6ei3vuucd44IEHjL///e/G4cOHjXvvvddwuVzGe++9N8g9H5l6ez7ee+8947nnnjP2799vHD9+3Hj22WeNuLg44/HHHx/kno88vT0XnWpqaoxJkyYZK1euNObOnTs4nR3hensuXnvtNUOScejQIaO4uLhr6+jo6NX7EqjO49JLLzXuvvvukH35+fnGunXrIrb/3ve+Z+Tn54fsu+uuu4zFixcPWB9Hi96ei55WrFhBoOpHVs5FpxkzZhgbNmzo766NSv1xPm644Qbjq1/9an93bdTp67lYtWqVcf/99xs/+MEPCFT9pLfnojNQVVdXW3pfhvwiaGtr0549e7Ry5cqQ/StXrtTOnTsjvmbXrl1h7a+++mq9++67am9vH7C+jnR9ORcYGP1xLgKBgOrr6/v9pqSjUX+cj8LCQu3cuVMrVqwYiC6OGn09F0899ZSOHj2qH/zgBwPdxVHDyt+L+fPnKysrS5/+9Kf12muv9fq9uTlyBBUVFfL7/crIyAjZn5GRoZKSkoivKSkpidi+o6NDFRUVysrKGrD+jmR9ORcYGP1xLh588EE1NjbqpptuGogujipWzkd2drbKy8vV0dGh9evX6/bbbx/Iro54fTkXR44c0bp167R9+3ZFR/NV3F/6ci6ysrL0xBNPqKCgQK2trXr22Wf16U9/Wq+//rouv/zyi35vzuIFOByOkN8Nwwjb90ntI+1H7/X2XGDg9PVcbN68WevXr9crr7yi9PT0gereqNOX87F9+3Y1NDTo7bff1rp16zRlyhR9+ctfHshujgoXey78fr9uueUWbdiwQdOmTRus7o0qvfl7kZeXp7y8vK7flyxZolOnTulnP/sZgcqq1NRURUVFhaXZsrKysNTbKTMzM2L76OhopaSkDFhfR7q+nAsMDCvnYsuWLfr617+uF154QZ/5zGcGspujhpXzkZubK0maPXu2SktLtX79egKVBb09F/X19Xr33XdVWFiob37zm5LM4XDDMBQdHa2tW7fqU5/61KD0faTpr++MxYsX63e/+12v3ps5VBHExMSooKBA27ZtC9m/bds2LV26NOJrlixZEtZ+69atWrhwoVwu14D1daTry7nAwOjrudi8ebO+9rWv6bnnntPnPve5ge7mqNFffzcMw1Bra2t/d29U6e258Hq92rdvn/bu3du13X333crLy9PevXu1aNGiwer6iNNffy8KCwt7P1XH0pT2Eazzssvf/OY3xoEDB4w1a9YY8fHxxokTJwzDMIx169YZq1ev7mrfuWzCd77zHePAgQPGb37zG5ZN6Ce9PReGYRiFhYVGYWGhUVBQYNxyyy1GYWGh8eGHHw5F90eU3p6L5557zoiOjjYeffTRkMuRa2pqhuojjCi9PR+//OUvjVdffdU4fPiwcfjwYePJJ580vF6vcd999w3VRxgx+vLvVE9c5dd/ensuHnroIePll182Dh8+bOzfv99Yt26dIcl48cUXe/W+BKoLePTRR40JEyYYMTExxoIFC4w33nij67nbbrvNWLFiRUj7119/3Zg/f74RExNjTJw40XjssccGuccjV2/PhaSwbcKECYPb6RGqN+dixYoVEc/FbbfdNvgdH6F6cz4eeeQRY+bMmUZcXJzh9XqN+fPnG5s2bTL8fv8Q9Hzk6e2/Uz0RqPpXb87FAw88YEyePNnweDxGUlKScdlllxl//OMfe/2eDsMIzpwGAABAnzCHCgAAwCICFQAAgEUEKgAAAIsIVAAAABYRqAAAACwiUAEAAFhEoAIAALCIQAUAAGARgQrAqLJ+/Xo5HA6tX7/e0nEmTpwoh8OhEydO9Op1V1xxhRwOh15//fWw55qamvTP//zPys3NlcvlksPh0Ne+9jVL/QQwOKKHugMAANMdd9yh5557TnFxcZo3b57cbremTZsmSV0B0GoQBDAwCFQAMIjGjx+vvLw8xcXFheyvrq7W888/r7i4OH300UfKyckJeX7Dhg2SCFSAXRGoAGAQPfPMMxH3HzlyRIFAQLNmzQoLUwDsjzlUAGADzc3NkqTY2Ngh7gmAviBQARgSDodDDodDkvTiiy/q8ssv15gxY8ImeldVVem+++7TrFmzFB8fr8TERC1evFi//vWvFQgEIh67o6NDP/3pT5Wfny+Px6Nx48bpjjvuUGlp6Xn7YxiGnnnmma5+xMTEKDMzUwUFBfre976n06dPn/e1b7/9tq699lolJSUpPj5ey5cv19/+9reIbc+dlH7ixAk5HA5dccUVkqQ33nij68+mc1J6559Tzz+3zq23k+IBDAyG/AAMqQceeEDr1q1TRkaGpk2bFhIQPvzwQ1199dU6c+aMYmJiNGXKFLW2turvf/+7du/era1bt+r3v/99SODw+/268cYb9d///d+SpGnTpik2NlZPPfWUtm7dqi9+8YsR+/Hd735XDz74oCRzntO0adNUUVGh/fv367333tPSpUuVnZ0d9rr/9//+n9auXSuv16vJkyfr448/1o4dO3T11Vdr27ZtXUHpfDwej5YtW6ba2lrt379fXq9Xs2fP7np+2rRpWrZsmd566y1J0rJly8JeD8AGDAAYApIMSUZMTIzxxBNPGIFAwDAMw2hvbzfa29uNhoYGY/LkyYYk49vf/rZRW1vb9doPP/zQmDlzpiHJ+OUvfxly3IcfftiQZCQlJRnbt2/v2n/8+HFj1qxZhsvlMiQZP/jBD7qeKysrM5xOp+Hz+YwdO3aEHK+5udnYvHmz8f7774fsnzBhgiHJcLlcxsaNG42Ojg7DMAyjra3N+MpXvmJIMhYtWhT2uVesWGFIMl577bWQ/a+99pohyVixYsUF/7wA2BNDfgCG1F133aU77rijq8oUHR2t6OhoPfnkkzp69KhuuOEGPfzww/J6vV2vmTFjhp577jk5HA793//7f7v2G4bRVWX60Y9+pMsuu6zruYkTJ+q3v/2t2tvbw/pw9OhRBQIBfepTn4pYAbr55ps1Z86ciP2/5pprtG7dOkVFRUmSXC6Xfv7zn8vtdmv37t2qrq7u458MgOGEQAVgSN16660R97/00kuSpNtvvz3i83PmzNHEiRN17NixrvlNBw8eVFFRkTweT8QFMRcsWKDFixeH7e+8qm737t0qKirqVf8j9S81NVUTJ06UJB07dqxXxwMwPDGHCsCQmj59esT9+/btkyR9//vf109+8pOIbSoqKiRJZ86cUXZ2tg4fPixJmjBhQtg6Tz3f7+233w7ZN27cOP2v//W/9MILL2jKlCm68sordcUVV2j58uVavHixoqPP/0/l5MmTI+5PT0/XoUOH1NDQcN7XAhg5CFQAhlR8fHzE/bW1tZKkPXv2fOIxOpcc6AwvaWlp522bkZERcf8zzzyjGTNm6D/+4z+0detWbd26tetY3/ve97R27Vo5neFF/fP1v7OtYRif2H8Awx9DfgBsKSEhQZK54KVhGBfcOq+k63xNeXn5eY9bVlYWcb/H49H69et1+vRpHTx4UI8//ri+8IUvqLKyUt/97ndD5moBwLkIVABsacaMGZKk/fv3X/RrOu97V1RUpKampohtDh48+InHyc/P15133qlXX31VmzZtkiT9+te/vuh+ABh9CFQAbOnGG2+UJD3yyCMXPWyWn5+vnJwcNTc3R7zFy969e7Vr165e9aNzEvvZs2d79br+1rmCeufwJgB7IVABsKW77rpLkyZN0muvvaavfOUrKi4uDnm+oaFBv//977V27dqufU6ns+v3++67Tzt37ux67uTJk7rtttvkcrnC3ut//ud/9N3vflcHDhwIe49///d/l2ReITiUJk2aJMlcSR2A/RCoANhSQkKC/vjHPyo3N1ebN29Wdna2ZsyYocWLFysvL09jxozRqlWrQkKTJH3rW9/SZz/7WVVVVWnZsmWaPn265s+frylTpqimpkZ33XVX2HvV19frZz/7mWbOnKn09HRdcsklmjdvnjIyMvSf//mf8vl8euihhwbro0e0atUqSdLnP/95LViwQFdccYWuuOIKlZSUDGm/AJi4yg+AbeXn5+v999/Xpk2b9PLLL+vgwYM6duyYsrKytGLFCn32s5/VP/zDP4S8JioqSv/1X/+lBx98UE899ZSOHTumlJQU3Xbbbfrxj3+sxx57LOx9li9frkceeUTbtm3T/v37deDAAblcLk2ZMkXXXHONvvOd7ygzM3OwPnZE69atk9/v1/PPP68DBw6otbVVktTS0jKk/QJgchhc0wsAAGAJQ34AAAAWEagAAAAsIlABAABYRKACAACwiEAFAABgEYEKAADAIgIVAACARQQqAAAAiwhUAAAAFhGoAAAALCJQAQAAWESgAgAAsIhABQAAYNH/H/fW+2ahDHT6AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "axes = ens.plot(xlim=(0., 0.5), label=\"PDF 0\")\n",
    "_ = ens.plot(key=1, axes=axes, label=\"PDF 1\")\n",
    "_ = ens.plot(key=1300, axes=axes, label=\"PDF 1300\")\n",
    "legend = axes.figure.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timing benchmarks\n",
    "\n",
    "Now we are going to extract some information from the ensemble and time how long it takes to do so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are the grid points and quantiles at which we will extract values\n",
    "test_xvals = ens.gen_obj.xvals\n",
    "test_quantiles = np.linspace(0, 1, 51)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Time the pdf evaluation for 20000 PDFs\n",
    "pdfs = ens.pdf(test_xvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Time the cdf (Cumulative distribution function) evaluation for 20000 PDFs\n",
    "cdfs = ens.cdf(test_xvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ppfs = ens.ppf(test_quantiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Time the sf (survival fraction, 1-cdf) evaluation for 20000 PDFs\n",
    "sfs = ens.sf(test_xvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Time the isf (inverse survival fraction) evaluation for 20000 PDFs\n",
    "isfs = ens.isf(test_quantiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Time the generation of 100 samples for each of the 20000 PDFs\n",
    "samples = ens.rvs(size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changing the grid for the representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a grid using 51 grid points\n",
    "ens_g51 = qp.convert(ens, 'interp', xvals=np.linspace(0, 3, 51))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a grid using 21 grid points\n",
    "ens_g21 = qp.convert(ens, 'interp', xvals=np.linspace(0, 1, 21))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 0\n",
    "axes_g = ens.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_g51.plot(key, axes=axes_g, label=\"g51\")\n",
    "_ = ens_g21.plot(key, axes=axes_g, label=\"g21\")\n",
    "leg_g = axes_g.figure.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conversion to other represenations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  quantile representaion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert using 51 quantiles\n",
    "ens_q51 = qp.convert(ens, 'quant', quants=np.linspace(0.01, 0.99, 51))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert using 21 quantiles\n",
    "ens_q21 = qp.convert(ens, 'quant', quants=np.linspace(0.01, 0.99, 21))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 0\n",
    "axes_q = ens.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_q51.plot(key, axes=axes_q, label=\"q51\")\n",
    "_ = ens_q21.plot(key, axes=axes_q, label=\"q21\")\n",
    "leg_q = axes_q.figure.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Histogram representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a histogram using 51 bins\n",
    "ens_h51 = qp.convert(ens, 'hist', bins=np.linspace(0, 3.0, 51))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a histogram using 21 bins\n",
    "ens_h21 = qp.convert(ens, 'hist', bins=np.linspace(0, 3.0, 21))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 0\n",
    "axes_h = ens.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_h51.plot(key, axes=axes_h, label=\"h51\")\n",
    "_ = ens_h21.plot(key, axes=axes_h, label=\"h21\")\n",
    "leg_h = axes_h.figure.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other representations\n",
    "\n",
    "`qp` also includes spline-based and Gaussian mixture represenations.   The conversion to these forms much slower, so we will first reduce the base ensemble from 20000 PDFs to 100 PDFs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ens_red = ens[np.arange(100)]\n",
    "print(\"Reduced ensemble has %i PDFs\" % (ens_red.npdf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spline representation\n",
    "\n",
    "We can convert to the spline representation a few different ways.  This particular method specifies that we should evaluate each PDF at a grid of points and then use those to construct the spline represenation.  We do this for 2 different grids.  Note how much slower this conversion function is that the ones above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a histogram using 51 grid points\n",
    "ens_s51 = qp.convert(ens_red, 'spline', xvals=np.linspace(0, 3.0, 51), method=\"xy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ens_s21 = qp.convert(ens_red, 'spline', xvals=np.linspace(0, 3.0, 21), method=\"xy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 0\n",
    "axes_s = ens_red.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_s51.plot(key, axes=axes_s, label=\"s51\")\n",
    "_ = ens_s21.plot(key, axes=axes_s, label=\"s21\")\n",
    "leg_s = axes_s.figure.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gaussian mixture model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a gaussian mixture using 301 sample points and 3 components\n",
    "ens_m3 = qp.convert(ens_red, 'mixmod', xvals=np.linspace(0, 3.0, 301), ncomps=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Convert to a gaussian mixture using 301 sample points and 3 components\n",
    "ens_m5 = qp.convert(ens_red, 'mixmod', xvals=np.linspace(0, 3.0, 301), ncomps=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 0\n",
    "axes_m = ens_red.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_m3.plot(key, axes=axes_m, label=\"m3\")\n",
    "_ = ens_m5.plot(key, axes=axes_m, label=\"m5\")\n",
    "leg_m = axes_m.figure.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "key = 1\n",
    "axes_m1 = ens_red.plot(key, xlim=(0, 0.5), label=\"orig\")\n",
    "_ = ens_m3.plot(key, axes=axes_m1, label=\"m3\")\n",
    "_ = ens_m5.plot(key, axes=axes_m1, label=\"m5\")\n",
    "leg_m1 = axes_m1.figure.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qp",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}