cmdstanpy_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CmdStanPy Tutorial\n",
    "\n",
    "\n",
    "## Workflow Outline\n",
    "\n",
    "Given a dataset and a model specification written as a Stan program, the CmdStanPy workflow is:\n",
    "\n",
    "### Assemble input data as either:\n",
    "  + A Python `dict` object consisting of key-value pairs where the key corresponds\n",
    " to Stan data variables and the value is of the correct type and shape.\n",
    "  + An existing data file on disk in either JSON or Rdump format.\n",
    "\n",
    "### Compile the model\n",
    "  + The `Model` object stores the filepath of the Stan program.\n",
    "  + The `compile` method translates the Stan program to C++ then calls the C++ compiler which creates the executable program.\n",
    "\n",
    "### Fit the model conditioned on the data to obtain an estimate of the posterior\n",
    "\n",
    "#### Run the NUTS-HMC sampler to obtain the full posterior density\n",
    "\n",
    "  + The `Model` class method `sample` invokes Stan's NUTS-HMC sampler which conditions the model on the input data and returns a `StanFit` object. The `StanFit` object contains the set of draws from the posterior plus metadata.\n",
    "  + Runs any number of chains - default is 4 chains.\n",
    "  + The output of each chain is stored on disk as a Stan csv file.\n",
    "\n",
    "#### Estimate the mode of the posterior density via optimization\n",
    "\n",
    "  + The `Model` class method `optimize` runs one of Stan's [optimization](https://mc-stan.org/docs/reference-manual/optimization-algorithms-chapter.html) algorithms to obtains an estimate of the mode of the posterior density.\n",
    "\n",
    "### Summarize and check the sample\n",
    "\n",
    "   + The `StanFit` class method `summary` invokes CmdStan's `stansummary` utility. Returns a Pandas DataFrame with estimates of posterior means, standard deviations, Monte-Carlo standard error, effective sample size, and convergence diagnostic statistic for all parameters in the model.\n",
    "   + The `StanFit` class method `diagnose` invokes CmdStan's `diagnose` utility which checks for the following problems:\n",
    "    + transitions that hit the maximum treedepth\n",
    "    + divergent transitions\n",
    "    + low E-BFMI values (sampler transitions HMC potential energy)\n",
    "    + low effective sample sizes\n",
    "    + high R-hat values\n",
    "   + See the Stan reference manual section on [posterior analysis](https://mc-stan.org/docs/reference-manual/analysis-chapter.html) for further details.\n",
    "\n",
    "###  Assemble the sample in-memory\n",
    "  + The resulting sample is accessed via the `StanFit` object:\n",
    "    + `sample`  - all draws from all chains, stored as a 3-D numpy.ndarray.\n",
    "    + `chains` - number of chains run by sampler\n",
    "    + `draws` - draws per chain\n",
    "    + `column_names` - names of the parameters, transformed parameters, and generated quantities variables returned in each draw\n",
    "    + `csv_files` - list of Stan csv output files which comprise the sample\n",
    "  + The method `get_drawset` flattens the 3-D sample array into a 2-D pandas.DataFrame for downstream analysis.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "* Install Python package from PyPI or directly from GitHub:\n",
    "\n",
    "  + `pip install --upgrade cmdstanpy`\n",
    "  + `pip install -e git+https://github.com/stan-dev/cmdstanpy`\n",
    "\n",
    "\n",
    "* CmdStanPy uses CmdStan directly to compile and run Stan programs, therefore CmdStan must be installed locally.\n",
    "\n",
    "  + if you have a working installation of CmdStan, set environment variable `CMDSTAN` to the full path to the top-level CmdStan directory.\n",
    "  \n",
    "  + if you don't already have CmdStan installed, run Python script `install_cmdstan` which downloads and compiles the latest release from https://github.com/stan-dev/cmdstan/releases.  By default this installs the latest version of CmdStan in the location `~/.cmdstanpy`.  Flags -d and -v are used to specify the directory and version, respectively.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hello World:  Compile and run CmdStan example model `bernoulli.stan`\n",
    "\n",
    "Every CmdStan release has an `examples/bernoulli` directory which contains the Stan model and test data files.  In this example we compile the model and run the sampler on the model and data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import relevant classes and methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import os.path\n",
    "from cmdstanpy import Model, StanFit, cmdstan_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CmdStanPy directory `examples/bernoulli` contains the model and data files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bernoulli_dir = os.path.join(cmdstan_path(), 'examples', 'bernoulli')\n",
    "bernoulli_path = os.path.join(bernoulli_dir, 'bernoulli.stan')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create a `Model` object for the specified Stan program file, then we call the `compile` method to create the corresponding C++ executable program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:cmdstanpy:compiling c++\n",
      "INFO:cmdstanpy:compiled model file: /Users/mitzi/.cmdstanpy/cmdstan-2.20.0/examples/bernoulli/bernoulli\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model(name=bernoulli,  stan_file=\"/Users/mitzi/.cmdstanpy/cmdstan-2.20.0/examples/bernoulli/bernoulli.stan\", exe_file=\"/Users/mitzi/.cmdstanpy/cmdstan-2.20.0/examples/bernoulli/bernoulli\")\n"
     ]
    }
   ],
   "source": [
    "bernoulli_model = Model(stan_file=bernoulli_path)\n",
    "bernoulli_model.compile()\n",
    "print(bernoulli_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the HMC-NUTS sampler on the model and data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Input data is either a Python `Dict` with entries corresponding to input data values, or it can be a file in JSON or Rdump format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bern_json = os.path.join(bernoulli_dir, 'bernoulli.data.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a `Dict` is specified, CmdStanPy writes it to a temp file in JSON format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bern_data = { \"N\" : 10, \"y\" : [0,1,0,0,0,0,0,0,0,1] }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Model` method `sample` runs the NUTS-HMC sampler and returns a `StanFit` object.  By default, `sample` runs 4 sampler chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:cmdstanpy:start chain 1\n",
      "INFO:cmdstanpy:finish chain 1\n",
      "INFO:cmdstanpy:start chain 2\n",
      "INFO:cmdstanpy:finish chain 2\n",
      "INFO:cmdstanpy:start chain 3\n",
      "INFO:cmdstanpy:finish chain 3\n",
      "INFO:cmdstanpy:start chain 4\n",
      "INFO:cmdstanpy:finish chain 4\n"
     ]
    }
   ],
   "source": [
    "bern_fit = bernoulli_model.sample(data=bern_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Access the sample: the `StanFit` object attributes and methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summarize the results\n",
    "\n",
    "The `StanFit` method `summary` returns output of CmdStan bin/stansummary as pandas.DataFrame.  The `summary` report provides estimates of both the parameter value and the goodness of fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "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>Mean</th>\n",
       "      <th>MCSE</th>\n",
       "      <th>StdDev</th>\n",
       "      <th>5%</th>\n",
       "      <th>50%</th>\n",
       "      <th>95%</th>\n",
       "      <th>N_Eff</th>\n",
       "      <th>N_Eff/s</th>\n",
       "      <th>R_hat</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>lp__</th>\n",
       "      <td>-7.265550</td>\n",
       "      <td>0.022104</td>\n",
       "      <td>0.776810</td>\n",
       "      <td>-8.746650</td>\n",
       "      <td>-6.980130</td>\n",
       "      <td>-6.750090</td>\n",
       "      <td>1235.07</td>\n",
       "      <td>9644.46</td>\n",
       "      <td>1.00074</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>theta</th>\n",
       "      <td>0.253735</td>\n",
       "      <td>0.003018</td>\n",
       "      <td>0.119649</td>\n",
       "      <td>0.084138</td>\n",
       "      <td>0.239654</td>\n",
       "      <td>0.471913</td>\n",
       "      <td>1571.63</td>\n",
       "      <td>12272.60</td>\n",
       "      <td>1.00056</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Mean      MCSE    StdDev        5%       50%       95%    N_Eff  \\\n",
       "name                                                                         \n",
       "lp__  -7.265550  0.022104  0.776810 -8.746650 -6.980130 -6.750090  1235.07   \n",
       "theta  0.253735  0.003018  0.119649  0.084138  0.239654  0.471913  1571.63   \n",
       "\n",
       "        N_Eff/s    R_hat  \n",
       "name                      \n",
       "lp__    9644.46  1.00074  \n",
       "theta  12272.60  1.00056  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run sampler diagnostics\n",
    "\n",
    "The `StanFit` method `diagnose` prints the output of CmdStan bin/stansummary to the console.  This is useful when the summary report shows `R_hat` values outside of the range of .99 to 1.01, or the number of effective samples (`N_eff` value) is below 5% of the total draws in the sample. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:cmdstanpy:Processing csv files: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-1-mgz2phsb.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-2-85nxjmmr.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-3-_504d8s7.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-4-ajyoemgm.csv\n",
      "\n",
      "Checking sampler transitions treedepth.\n",
      "Treedepth satisfactory for all transitions.\n",
      "\n",
      "Checking sampler transitions for divergences.\n",
      "No divergent transitions found.\n",
      "\n",
      "Checking E-BFMI - sampler transitions HMC potential energy.\n",
      "E-BFMI satisfactory for all transitions.\n",
      "\n",
      "Effective sample size satisfactory.\n",
      "\n",
      "Split R-hat values satisfactory all parameters.\n",
      "\n",
      "Processing complete, no problems detected.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'Processing csv files: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-1-mgz2phsb.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-2-85nxjmmr.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-3-_504d8s7.csv, /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp52kam1bg/stan-bernoulli-draws-4-ajyoemgm.csv\\n\\nChecking sampler transitions treedepth.\\nTreedepth satisfactory for all transitions.\\n\\nChecking sampler transitions for divergences.\\nNo divergent transitions found.\\n\\nChecking E-BFMI - sampler transitions HMC potential energy.\\nE-BFMI satisfactory for all transitions.\\n\\nEffective sample size satisfactory.\\n\\nSplit R-hat values satisfactory all parameters.\\n\\nProcessing complete, no problems detected.'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.diagnose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assemble all draws from all chains as an in-memory pandas DataFrame\n",
    "\n",
    "The `StanFit` method `get_drawset` returns a pandas.DataFrame, one draw per row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bern_drawset = bern_fit.get_drawset()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, `get_drawset` returns a DataFrame which contains all columns from the sampler's csv output file, i.e., it contains both the sampler state and the values for all parameter, transformed parameter, and generated quantities variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((4000, 8),\n",
       " Index(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__',\n",
       "        'divergent__', 'energy__', 'theta'],\n",
       "       dtype='object'))"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_drawset.shape, bern_drawset.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_drawset` function argument `params` takes a list of parameter or column names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4000, 1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thetas = bern_fit.get_drawset(params=['theta'])\n",
    "thetas.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "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>theta</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.418464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.418464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.267071</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      theta\n",
       "0  0.418464\n",
       "1  0.418464\n",
       "2  0.267071"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thetas[0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x112972518>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD8CAYAAACCRVh7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8XOV97/HPT7tlbZYl2bK8yDZeMGAbY8wOhiRsYcnNVkjThDYJISRpk3SjyX2RNG1vmvSVLoTcNpCN9GYhISk1BEJIgg0kYGyIF2xsY7xJ3iRZq7XYkuZ3/5gRlW2NNJZ15syMvu8X8/LozJkzX+QZ/+Y8z3Oex9wdERGRoWSFHUBERFKXioSIiMSlIiEiInGpSIiISFwqEiIiEpeKhIiIxKUiISIicalIiIhIXCoSIiISV07YAU5XRUWF19bWhh1DRCStvPzyy03uXnm6z0u7IlFbW8v69evDjiEiklbMbO9onqfmJhERiUtFQkRE4lKREBGRuFQkREQkLhUJERGJS0VCRETiUpEQEZG40u46CZGhdPT08otXD9HQcYwVs8u5sLY87EgiGUFFQtLezzcd5HOPbqa1q/fNbbevmME/vOM8srIsxGQi6U9FQtLaN9a8wZee3Mb5M8u496ZFzJtSzNd+8zrfWLOLquICPv22+WFHFElrKhKSliIR50tPvsaDz+3mpsXVfPW9S8jPyQbgnusXcrith68/s5Mbz6tmwdTikNOKpC91XEvaOdbXz589vIEHn9vNBy+ZxX23nf9mgQAwM+69+RwK87L5l6d3hJhUJP3pTEJSVvfxfrYcaAOgfGIexQW5bKhr5au/3M62Qx3cc8NCPnrlHMxO7Xcon5jHBy6p5eurd/JG41HmVhYlO75IRlCRkJS0Zkcjn354A82dx095bFppAd/64HLecvaUYY9xx2W1fOPZN/jh2n3875sWBRVVJKOpSEjK2dlwlDu/t545lUV8+V2LycvJornzGG1dvdRWTOSysyrIzR65pbSiKJ9rFlbx6IYD3HPDQnISeI6InEhFQlLO3/98K/k5WTz0JxdSVVxwRsd657LpPLXlMM/vbGLlgqoxSigyfuirlaSU7Yc6WL29kY9eNfeMCwTAVfMrKczL5umth8cgncj4oyIhKeXRDfvJzjJuXzFzTI5XkJvNFfMq+PVrDbj7mBxTZDxRkZCU4e48ufkgl86dTPnEvDE77lvPnsKh9h62HGgfs2OKjBcqEpIyth/uYM+RLm44t3pMj3vNwirM4FevqclJ5HSpSEjKeOGNIwCsXFA5psedXJTPudNK+V3s+CKSOBUJSRkv7W5m+qQJTCubMObHvnTuZDbsa6X7eP+YH1skkwVWJMyswMxeMrONZrbFzP52iH3yzexhM9tpZmvNrDaoPJLa3J11e5pZMTuYKb4vnjuZ4/0RXt7bEsjxRTJVkGcSx4Br3H0JsBS43swuPmmfDwEt7n4W8C/AlwPMIylsz5Eumo4eD2wdiAtry8nJMl7Y1RTI8UUyVWBFwqOOxn7Mjd1OHoN4K/BQ7P4jwFtsqIl4JONtjY08Oq+mNJDjF+XnsGRGmfolRE5ToH0SZpZtZhuABuBpd1970i41QB2Au/cBbcDkIDNJatp6sI2cLGPelOAm4rtkzmQ21bdx9FhfYK8hkmkCLRLu3u/uS4HpwAozO3c0xzGzO81svZmtb2xsHNuQkhK2HmjnrKqiE6b8HmsXz5lMf8R5Rf0SIglLyugmd28FngGuP+mh/cAMADPLAUqBU9oD3P0Bd1/u7ssrK8d2eKSkhtcOdnB2dUmgr7FsVhk5Wcba3WpyEklUkKObKs2sLHZ/AvA2YNtJu60CPhi7/27gN665E8ad5s7jHGrvYVHARaIwL4dza0p5aXdzoK8jkkmCPJOoBp4xs03AOqJ9Eo+b2RfN7JbYPt8CJpvZTuAzwD0B5pEU9drBaKd10GcSABfNLmdjXRs9vbpeQiQRgU0V7u6bgPOH2H7voPs9wHuCyiDpYcfhDgDmTw1+9bgVs8v5xrO72FDXysVzNEZCZCS64lpCt6epk+L8HCqL8gN/reWzyjFDTU4iCVKRkNDtauqktmLikGtVj7XSwlwWTi1RkRBJkIqEhG53UyezKyYm7fUuml3Oy3tb6O2PJO01RdKVioSEqqe3n/2t3UkvEt29/Wze35a01xRJVyoSEqq65i7cYU5l8orEhbFJBNXkJDIyFQkJ1a6mTgBqJyevSFQU5TO3cqKKhEgCVCQkVLsHikQSm5sAVsyezLo9zfRHdO2myHBUJCRUuxs7qSjKo3RCblJf96LZ5XT09L15IZ+IDE1FQkK1+0hyRzYNWKF+CZGEqEhIqPY0dTIrif0RA6aVTWBG+QQVCZERqEhIaHp6+2noOMaMSYWhvP6K2sm8tKcZzSkpEp+KhITmQGs3ANMnTQjl9S+aXU5z53HeaDw68s4i45SKhISmviXcIjHQL7FWTU4icalISGjeLBLl4TQ3zZpcSFVxPmt3qUiIxKMiIaGpb+kiJ8uYUhz87K9DMTNWzC7npd3qlxCJR0VCQlPf0k11WQE52eG9DS+aXc6h9h7qmrtDyyCSylQkJDT1LV1MLwunqWnARbGFh7TutcjQVCQkNPUt3aF1Wg84q7KISYW5ul5CJA4VCQnFsb7oNRLTQ7pGYkBWlnFhbTkv7VGREBmKioSE4kBrDxDe8NfBVswuZ++RLg62qV9C5GQqEhKK+pYuIDWKxOXzKgB4bkdTyElEUk9gRcLMZpjZM2a21cy2mNmfDbHPSjNrM7MNsdu9QeWR1BL2NRKDLZhSzJSSfNbsaAw7ikjKyQnw2H3An7v7K2ZWDLxsZk+7+9aT9nvO3W8KMIekoLCvkRjMzLhyXiW/3HqYvv5IqENyRVJNYJ8Gdz/o7q/E7ncArwE1Qb2epJdUuEZisCvnV9LW3cvGeq17LTJYUj6hZlYLnA+sHeLhS8xso5k9aWbnJCOPhK++pTv0ayQGu/ysCrIMnlWTk8gJAi8SZlYE/BT4lLufvAzYK8Asd18CfA14NM4x7jSz9Wa2vrFRH+JMUN/SlRKd1gMmTcxj8fQy9UuInCTQImFmuUQLxPfd/WcnP+7u7e5+NHb/CSDXzCqG2O8Bd1/u7ssrKyuDjCxJcKyvn8Ptx6hJoSIBcNX8SjbWt9LSeTzsKCIpI8jRTQZ8C3jN3f85zj5TY/thZitieTQ/QoY7+OY1EqnT3ARw1YJK3OH5nRoKKzIgyNFNlwF/BGw2sw2xbZ8FZgK4+38A7wY+ZmZ9QDdwm2s6zow3sNjQtLKCkJOcaMn0Mkon5LJmRyM3L5kWdhyRlBBYkXD35wEbYZ/7gfuDyiCpqX5gRboU6rgGyM4yLp9XwbM7GnF3Yie5IuNaaow/lHFlf0s3ZjC1NLXOJACumldJQ8cxth3qCDuKSEpQkZCkO9DaTVVxPnk5qff2u3J+dGCERjmJRKXep1Qy3v7WbmrKUmtk04CppQUsnFrM6u0NYUcRSQkqEpJ0+1u7mZaiRQJg5YIq1u9poaOnN+woIqFTkZCkikScg609KXeNxGArF1TSF3F+q6GwIioSklxNR49xvD/C9BQ+k7hg1iSK83NYvV39EiIqEpJU+9+8RiJ1i0RudhaXz6tg9fboUFiR8UxFQpJqoEikcnMTRJucDrX3aCisjHsqEpJU+1tS/0wC4Kr5VQBqcpJxT0VCkupAazfFBTmUFOSGHWVYU0sLOLu6RENhZdxTkZCkSuVrJE62ckElL+9toV1DYWUcU5GQpNrf2pM+RWJ+bCjs6xoKK+OXioQk1f6WrpTvtB6wTENhRVQkJHk6enpp7+lL+U7rAbnZWVwxv4LVOxo0FFbGLRUJSZoDscWG0qW5CWDl/CoOtx/jtYMaCivjk4qEJM3+1i4g9Ye/DnbVguissKt3aJSTjE8qEpI0+99ctjR9isSUkgIWVZeoX0LGLRUJSZr9Ld3kZhuVRflhRzktGgor45mKhCTN/tZuqksnkJWVXsuCrlxQRX/EeV5DYWUcUpGQpDmQRhfSDbZsZhnFBTm6+lrGJRUJSZr9Lam92FA8OdlZXDmvkjU7NCusjD+BFQkzm2Fmz5jZVjPbYmZ/NsQ+Zmb3mdlOM9tkZsuCyiPh6u2PcLgjtRcbGs5VCyo1FFbGpSDPJPqAP3f3RcDFwMfNbNFJ+9wAzIvd7gT+PcA8EqJDbT24k9KLDQ1n5fzoUNhn1OQk40xgRcLdD7r7K7H7HcBrQM1Ju90KfM+jXgTKzKw6qEwSnvo0mSI8nqqSAs6ZVsIaDYWVcSYpfRJmVgucD6w96aEaoG7Qz/WcWkgkA6TLYkPDWbmgkpf3tdDWraGwMn4EXiTMrAj4KfApd28f5THuNLP1Zra+sVHf5NJRXXMXZjCtrCDsKKOmobAyHgVaJMwsl2iB+L67/2yIXfYDMwb9PD227QTu/oC7L3f35ZWVlcGElUDVtXQxtaSA/JzssKOM2vkzyijRUFgZZ4Ic3WTAt4DX3P2f4+y2CvhAbJTTxUCbux8MKpOEp765mxmTCsOOcUZysrO4Yr6Gwsr4EuSZxGXAHwHXmNmG2O1GM7vLzO6K7fMEsAvYCTwI3B1gHglRXUsX08vTtz9iwMr5lTR0HGPrwVG1nIqknZxEdjKznxE9K3jS3SOJPMfdnweGnX/Bo1/HPp7I8SR9Hevr51B7T9qfScCgWWG3N3LOtNKQ04gEL9Ezif8LvA943cz+0cwWBJhJMsz+lm7cYUZ5+heJquICzq0pUb+EjBsJFQl3/5W7/yGwDNgD/MrMfmdmfxzrnBaJqy52jcTMDCgSEF2I6JV9rbR1aSisZL6E+yTMbDJwB/Bh4PfAvxEtGk8HkkwyRl1zdLGhGRnQJwFw9cJK+iPOczs1HFsyX0JFwsz+C3gOKARudvdb3P1hd/8kUBRkQEl/dS1d5GVnMaU4fa+RGGzpjEmUTsjVQkQyLiTUcQ086O5PDN5gZvnufszdlweQSzJIfXM3NZPSbx2JeLKzjCvmVbBmRyORiGfM/5fIUBJtbvr7Iba9MJZBJHPta+5KqyVLE3H1gioaNRRWxoFhzyTMbCrRuZQmmNn5/M+Q1hKiTU8iI6pr6eK86Zk1b+OV8weGwjZwbo2GwkrmGqm56TqindXTgcFXTXcAnw0ok2SQjp5eWrt6M+IaicEqi/M5r6aU1dsb+cQ188KOIxKYYYuEuz8EPGRm73L3nyYpk2SQuubMGv462MoFlXz9mZ20dfVSWqiR4JKZhu2TMLP3x+7WmtlnTr4lIZ+kubqWzBr+OtjKBVVEHJ59XaOcJHON1HE9MfZnEVA8xE1kWG9eI5FhzU0AS2eUUVaoobCS2UZqbvpG7M+/TU4cyTR1zV0U5edQloHNMdGhsJUaCisZLdGL6b5iZiVmlmtmvzazxkFNUSJx1bV0M33SBKIzx2eeqxdU0nRUQ2ElcyV6ncS1sVXlbiI6d9NZwF8GFUoyx77mroyY2C+egaGwz2zThH+SmRItEgPNUm8HfuLubQHlkQzSH3H2HelidsXEkXdOUxVF+SyeXsrqHeqXkMyUaJF43My2ARcAvzazSqAnuFiSCQ6193C8P8KsyZl7JgHRUU6/39dCa9fxsKOIjLlEpwq/B7gUWO7uvUAncGuQwST97W3qBGD25Mw9k4Do9RIRh+debwo7isiYS3SCP4CFRK+XGPyc741xHskge45Eh7/OyuDmJoAl08uYVJjLb7Y1cPOSaWHHERlTiS5f+p/AXGAD0B/b7KhIyDD2HOkkLyeL6pLMmCI8nuws4y1nT+GXWw7R2x8hNzvIpeNFkivRM4nlwKLYmtQiCdnT1MnM8sJxcf3AdedM5ZGX63lx1xGumFcZdhyRMZPoV55XgalBBpHMs/dIF7UZ3h8x4Ip5FRTmZfPUlkNhRxEZU4kWiQpgq5k9ZWarBm7DPcHMvm1mDWb2apzHV5pZm5ltiN3uPd3wkroiEWdvcye1GT6yaUBBbjZXza/kl1sOE4nohFsyR6LNTV8YxbG/C9zP8P0Wz7n7TaM4tqS4wx099PRGMr7TerDrzpnKk68eYkN9K8tmTgo7jsiYSHQI7BqiV1rnxu6vA14Z4TnPAs1nGlDS056m6Mim8XImAXD1wipyskxNTpJREp276SPAI8A3YptqgEfH4PUvMbONZvakmZ0zBseTFLH3SPQaifHSJwFQOiGXS+ZO5qlXD6ExHpIpEu2T+DhwGdAO4O6vA1Vn+NqvALPcfQnwNYYpOmZ2p5mtN7P1jY2a/iAd7DnSRW62Ma0s89aRGM5150xlz5Eudhw+GnYUkTGRaJE45u5vzjkQu6DujL4quXu7ux+N3X8CyDWzijj7PuDuy919eWWlhhemg71HOplRXkj2OBj+Oti1i6ZghpqcJGMkWiTWmNlngQlm9jbgJ8BjZ/LCZjbVYvNHm9mKWJYjZ3JMSR27mzrHVVPTgKqSAi6YOYknNh8MO4rImEi0SNwDNAKbgY8CTwD/e7gnmNkPgReABWZWb2YfMrO7zOyu2C7vBl41s43AfcBtulgvM7j7uLpG4mRvX1zNtkMd7GzoCDuKyBlLaAisu0fM7FHgUXdPqFPA3W8f4fH7iQ6RlQzT2HGM7t5+aivGz8imwd5+XjVffHwrj208yKffplV+Jb0NeyZhUV8wsyZgO7A9tiqdLnyTuN6c2G+cnklUlRRw0exyHtt0QKOcJO2N1Nz0aaKjmi5093J3LwcuAi4zs08Hnk7S0p6mgeGv4/NMAuCmxdPY1djJawfV5CTpbaQi8UfA7e6+e2CDu+8C3g98IMhgkr72HOkkJ8uoGWfDXwe74dypZGcZj286EHYUkTMyUpHIdfdTVlKJ9UvkBhNJ0t3eI9F1rXPG8ZTZk4vyuXTuZB7fdFBNTpLWRvoUD7ceo9ZqlCHtOdKZ8UuWJuLmxdPY19zFpnotCS/pa6QiscTM2oe4dQDnJSOgpBd3Z3dTJ7PH0cR+8Vx3zlRys9XkJOlt2CLh7tnuXjLErdjd1dwkpzjcfoyu4/3MUZGgtDCXK+dV8vNNBzV9uKSt8dtoLIHY1RSds2hOZVHISVLDTUuqOdDWwyv7WsKOIjIqKhIypnY1Roe/qrkp6q1nTyE/J4vHN2maDklPKhIypnY3dTIhN5upJQVhR0kJxQW5XL2gip9vPki/mpwkDalIyJja1XiU2oqJZI2z2V+Hc8vSaTR2HOPFXZq/UtKPioSMqV1NncypVFPTYNcsrKIoP4dVGzTKSdKPioSMmeN9EeqauzSy6SQFudlce84Unnj1IMf6+sOOI3JaVCRkzOxr7iTi6ExiCLcsmUZHTx9rtmtlRUkvKhIyZv5nZJOGv57ssrMqKJ+Yx39vVJOTpBcVCRkzu5o0/DWe3Ows3n5eNb9+7TBHj/WFHUckYSoSMmZ2N3ZSUZRH6QRdjD+UW5ZOo6c3wtNbtf61pA8VCRkzu5qOMkdNTXFdMHMSNWUTNMpJ0oqKhIyZ3Rr+OqysLOOmJdU893oTzZ2aRFnSg4qEjIm27l6ajh5Xf8QIbl1SQ1/EeWKzpumQ9KAiIWNid6zTWhP7De/s6mLOqipSk5OkjcCKhJl928wazOzVOI+bmd1nZjvNbJOZLQsqiwRvV2N09ledSQzPzLh1yTRe2tPMgdbusOOIjCjIM4nvAtcP8/gNwLzY7U7g3wPMIgHb3dRJdpYxs1wr0o3k5iXTAHhM10xIGgisSLj7s0DzMLvcCnzPo14EysysOqg8EqxdjZ3MmDSBvBy1YI6ktmIiS2aUsUpFQtJAmJ/oGqBu0M/1sW2ShqIT+6k/IlG3LJnGlgPt7Gw4GnYUkWGlxdc+M7vTzNab2frGRs19k2oiEWd301H1R5yGmxdXY4bOJiTlhVkk9gMzBv08PbbtFO7+gLsvd/fllZWVSQkniTvU3kNPb0TXSJyGqpICLpkzmVUb9uOuxYgkdYVZJFYBH4iNcroYaHN3DR5PQ1qydHTesbSGPUe6tP61pLQgh8D+EHgBWGBm9Wb2ITO7y8zuiu3yBLAL2Ak8CNwdVBYJ1u6maLv6XPVJnJYbF1dTmJfNj9fVhx1FJK6coA7s7reP8LgDHw/q9SV53mjsZGJeNlXF+WFHSStF+TnceF41j286wOdvWURhXmAfR5FRS4uOa0ltu5s6mV05ETOta3263rt8Bp3H+3lis2aGldSkIiFn7I1Gzf46WhfWTqJ2ciE/WV838s4iIVCRkDPSdbyP/a3dnFWlIjEaZsZ7ls9g7e5m9h7pDDuOyClUJOSM7GrsxB0ViTPwzmU1ZBk88rI6sCX1qEjIGRm4YlhFYvSqSydw5fxKHnm5nv6IrpmQ1KIiIWdkZ8NRsrOM2sm6RuJMvOeCGRxs6+H5nU1hRxE5gYqEnJGdDUeZNblQE/udobcuqqKsMFcd2JJy9MmWM/J6Qwdn6SK6M5afk807ltbwyy2HtbSppBQVCRm13v4Ie490qT9ijNy2YgbH+yP87BV1YEvqUJGQUdt7pJO+iKtIjJGFU0tYOqOMH62r06R/kjJUJGTUNLJp7L1vxUx2Nhzl5b2a9E9Sg4qEjNpAkdDEfmPnpiXVFOXn8IOX9oUdRQRQkZAz8HrDUWrKJjAxXxPTjZXCvBxuWTqNn286SFtXb9hxRFQkZPR2Nhxlrpqaxtz7VszkWF+ERzcMuQaXSFKpSMioRCLOG41HNfw1AOfWlHJuTQk/fGmfOrAldCoSMir7W7vp6Y2o0zogt104k22HOthY3xZ2FBnnVCRkVF5v6ABg3hQViSDcunQaE3Kz+eFadWBLuFQkZFReOxgtEgumFoecJDMVF+Ry85JqHtt0gI4edWBLeFQkZFS2HeqgpmwCJQW5YUfJWLevmEnX8X5WbTwQdhQZx1QkZFS2HWzn7GqdRQRp6YwyFk4t5kcvadI/CY+KhJy2nt5+djV1snBqSdhRMpqZcfuKmWze38aGutaw48g4FWiRMLPrzWy7me00s3uGePwOM2s0sw2x24eDzCNjY2fDUfojzkKdSQTunctqKM7P4cHndoUdRcapwIqEmWUDXwduABYBt5vZoiF2fdjdl8Zu3wwqj4ydbYeindY6kwhecUEut180kyc3H6SuuSvsODIOBXkmsQLY6e673P048CPg1gBfT5Jk28F28nOyqJ1cGHaUceGOS2vJMuM7v90TdhQZh4IsEjXA4B63+ti2k73LzDaZ2SNmNiPAPDJGth3qYP6UYnKy1aWVDNPKJnDT4moeXrePtm4Nh5XkCvtT/hhQ6+6LgaeBh4bayczuNLP1Zra+sbExqQHlVNsOtbNQ10ck1UeunEPn8X6+q7MJSbIgi8R+YPCZwfTYtje5+xF3Pxb78ZvABUMdyN0fcPfl7r68srIykLCSmIaOHpqOHmdhtfojkumcaaVcu2gK33xuF61dWt5UkifIIrEOmGdms80sD7gNWDV4BzOrHvTjLcBrAeaRMbCpLjqX0OLppSEnGX8+c+18jh7v00gnSarAioS79wGfAJ4i+o//j919i5l90cxuie32p2a2xcw2An8K3BFUHhkbG+tbyc4yzpmmM4lkWzi1hLefV813fruHho6esOPIOBFon4S7P+Hu8919rrv/Q2zbve6+Knb/b9z9HHdf4u5Xu/u2IPPImdtY38a8qiIK87TQUBj+4toF9PZH+MovtocdRcaJsDuuJY24O5vqW1kyvSzsKONWbcVE/uTy2Tzycr2uwpakUJGQhNU1d9Pa1cviGeqPCNMnr5lHZXE+X1i1hUhEixJJsFQkJGEb6qPfXHUmEa6i/Bz++vqFbKhr5YfrtN6EBEtFQhK2sa6VvJwsrSGRAt61rIZL507mS09s40Brd9hxJIOpSEjC1u9pZun0MnJ1pXXozIx/fOdi+iPOZ/9rs9bClsDo0y4JOXqsj1cPtHPRnPKwo0jMzMmF/NX1C1i9vZEfrdOaExIMFQlJyPo9zfRHnItmTw47igzygUtqufysCr6wagtbD7SHHUcykIqEJGTt7mZysoxls9RpnUqys4x/vW0pZYW53P39l2nr0gSAMrZUJCQhL+46wnnTS3URXQqqKMrn/vct40BrD3d89yU6j/WFHUkyiIqEjKil8zgb61q5Yp4mV0xVF9aWc9/t57OxrpUPfPslmo4eG/lJIglQkZARrdnRSMThmoVVYUeRYVx/7lTuf98ythxo4+33PceTmw9q1JOcMRUJGdFvtjVQUZTH4hpdaZ3qbjyvmkfuupRJhXl87PuvcON9z/Nfv6+ntz8SdjRJUyoSMqze/ghrdjSyckEVWVkWdhxJwLk1pTz2ycv5yrsX09cf4dMPb+SKLz/Dv69+QyvbyWlTkZBhPbujkbbuXq47Z2rYUeQ05GZn8d7lM3jqU1fynT++kLOqivjyL7ax8p+e4eF1+9QMJQlTkZBh/ez3+ymfmMdV89VpnY6ysoyrF1Tx/z58EY9/8nLOqirir3+6mY98b71WuJOEqEhIXO09vTy99TA3L64mL0dvlXR3bk0pP/7oJXz+5kWs2dHIjf/2HK/sawk7lqQ4ffIlrp+sr+d4X4R3XTA97CgyRsyMP75sNj/92KVkZxvv/Y8X+M5vd6v5SeJSkZAh9fZH+Pbzu1lRW85iTQ2ecRZPL+PxT1zBygVV/O1jW/nED35PR486teVUKhIypB+tq2N/azd3rZwTdhQJSGlhLg9+4ALuuWEhv9hyiJu/9jwvvHEk7FiSYlQk5BQtncf516d3cNHscq5eoAvoMpmZcddVc/nBhy+i353bH3yRz/x4g9aokDepSMgJ3J17fraJ9p5ePn/zOZjp2ojx4KI5k/nlp67i7pVzeWzjAVb+02ru/e9X2a9iMe4FWiTM7Hoz225mO83sniEezzezh2OPrzWz2iDzyPDcnf/zxGs8teUwf3XdQhZNKwk7kiTRhLxs/ur6haz+y6t51wXT+cHafVzx5d/w4YfWs3p7A326antcsqBGNZhZNrADeBtQD6wDbnf3rYP2uRtY7O53mdltwP9y9z8Y7rjLly/39evXB5J5PGto7+Hzq7bw5KuH+OAls/jCLTqLGO/2t3bzg7V7eXhL5M53AAAIHklEQVRdHU1HjzOpMJdrFk7hqgWVLJleyszyQr1H0oiZvezuy0/7eQEWiUuAL7j7dbGf/wbA3b80aJ+nYvu8YGY5wCGg0ocJpSIxNrqP93OovYctB9pYs72RxzcdpD/ifOba+Xz0yjn68MubjvdF+PVrh3lqyyF+s62B9p7oVOSFedlMLSmgsjifgtxscrONnKwssrLAMGL/YWYYkJNllEzIjd4KciidkEvphFyKC3Ipjv1cXJBDUX4OOVoid8yNtkgEuThADTB4TcV64KJ4+7h7n5m1AZOBprEOs2ZHI3/3ePQkZqAGnVCJnFO2Da5V/ua2wU/xU7cNUd5GdZwhjzfoOMPkPXHbqRn6I07n8f43f56Yl80tS6Zx18q5zK6YeOr/gIxreTlZ3HBeNTecV01vf4TthzrYVN/GzoajHG7voaGjh5au4/T1O32RCBGPvt8cwKPvR3ent9/p6Oml41jfkJ+TwQrzsinIzSY7y8jJMrLMyMk2ss0wI+2/xIw2/R9cOIMPX5HcEYdpsYKMmd0J3Akwc+bMUR2jKD+HBVOKBx30hD8GXmeIbac85YQ3qJ1yJ/YtKs5zB2+3QVvf3HbCu+cMjnPCfie+JbPMqCjOY2pJAXMri1g0rYRcfXOTBORmZ3FuTSnnnsGMwP0R52hPH+09vbR199Le00t7d1+0gMS2d/T0cayvn/4I9Eci9EWcSMTpi/iIBSbVOaP/H6goyh/DJIkJskjsB2YM+nl6bNtQ+9THmptKgVMGarv7A8ADEG1uGk2YC2ZN4oJZk0bzVBEZQ9lZRmlhLqWFuSf8AyGpKcivj+uAeWY228zygNuAVSftswr4YOz+u4HfDNcfISIiyRXYmUSsj+ETwFNANvBtd99iZl8E1rv7KuBbwH+a2U6gmWghERGRFBFon4S7PwE8cdK2ewfd7wHeE2QGEREZPfVWiohIXCoSIiISl4qEiIjEpSIhIiJxqUiIiEhcgc3dFBQzawT2JrBrBQFM7zGGUj0fpH7GVM8HqZ8x1fNB6mdM9XwQzTjR3StP94lpVyQSZWbrRzOZVbKkej5I/Yypng9SP2Oq54PUz5jq+eDMMqq5SURE4lKREBGRuDK5SDwQdoARpHo+SP2MqZ4PUj9jqueD1M+Y6vngDDJmbJ+EiIicuUw+kxARkTOUMUXCzMrN7Gkzez32Z9zFI8ysxMzqzez+VMpnZkvN7AUz22Jmm8xs2PW+xzDb9Wa23cx2mtk9Qzyeb2YPxx5fa2a1ych1Gvk+Y2ZbY7+zX5vZrGTmSyTjoP3eZWZuZkkdDZNIPjN7b+z3uMXMfpDMfIlkNLOZZvaMmf0+9nd9Y5LzfdvMGszs1TiPm5ndF8u/ycyWpVi+P4zl2mxmvzOzJQkd2N0z4gZ8Bbgndv8e4MvD7PtvwA+A+1MpHzAfmBe7Pw04CJQFnCsbeAOYA+QBG4FFJ+1zN/Afsfu3AQ8n8feWSL6rgcLY/Y8lM1+iGWP7FQPPAi8Cy1MpHzAP+D0wKfZzVar9Dom2q38sdn8RsCfJGa8ElgGvxnn8RuBJogtDXgysTbF8lw76+70h0XwZcyYB3Ao8FLv/EPCOoXYyswuAKcAvk5RrwIj53H2Hu78eu38AaABO++KX07QC2Onuu9z9OPCjWNbBBmd/BHiLJW+R4RHzufsz7t4V+/FFoqsgJlMiv0OAvwO+DPQkMxyJ5fsI8HV3bwFw94YUzOhASex+KXAgiflw92eJrnsTz63A9zzqRaDMzKqTk27kfO7+u4G/X07jc5JJRWKKux+M3T9EtBCcwMyygK8Cf5HMYDEj5hvMzFYQ/Ub1RsC5aoC6QT/Xx7YNuY+79wFtwOSAc53y2jFD5RvsQ0S/zSXTiBljTQ8z3P3nyQwWk8jvcD4w38x+a2Yvmtn1SUsXlUjGLwDvN7N6ouvUfDI50RJ2uu/VMCX8OQl00aGxZma/AqYO8dDnBv/g7m5mQw3buht4wt3rg/giPAb5Bo5TDfwn8EF3j4xtysxlZu8HlgNXhZ1lsNiXk38G7gg5ynByiDY5rST6DfNZMzvP3VtDTXWi24HvuvtXzewSoqtanqvPyOkxs6uJFonLE9k/rYqEu7813mNmdtjMqt39YOwf2aFOly8BrjCzu4EiIM/Mjrp73I7GJOfDzEqAnwOfi52yBm0/nLAe/fTYtqH2qTezHKKn+keSkG3waw8YKh9m9laixfgqdz+WpGwDRspYDJwLrI59OZkKrDKzW9x9fQrkg+i33rXu3gvsNrMdRIvGuiTkg8Qyfgi4HsDdXzCzAqJzEiW7aSyehN6rYTKzxcA3gRvcPbHPcDI7VgLutPknTuwY/soI+99BcjuuR8xHtHnp18CnkpgrB9gFzOZ/OgzPOWmfj3Nix/WPUyzf+USb5eaF9N4bMeNJ+68muR3XifwOrwceit2vINpsMjnFMj4J3BG7fzbRPglL8t91LfE7ht/OiR3XL4XwXhwu30xgJ3DpaR0z2f8TAf5yJsf+gX0d+BVQHtu+HPjmEPsnu0iMmA94P9ALbBh0W5qEbDcCO2L/0H4utu2LwC2x+wXAT2JvsJeAOUn+ux0p36+Aw4N+Z6tCeP8Nm/GkfZNaJBL8HRrRJrGtwGbgtlT7HRId0fTbWAHZAFyb5Hw/JDrisJfomdeHgLuAuwb9Dr8ey785hL/jkfJ9E2gZ9DlZn8hxdcW1iIjElUmjm0REZIypSIiISFwqEiIiEpeKhIiIxKUiISIicalIiIhIXCoSIiISl4qEiIjE9f8BM4NA+RUJFJMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bern_drawset.theta.plot.density()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `StanFit` property `sample` is a 3-D numpy ndarray which contains all draws across all chains.  This array is created only as needed; therefore the first time that this property is accessed CmdStanPy will read in the contents of the sampler's csv output files.  Because the csv output files also contain stepsize and metric information, the `stepsize` and `metric` arrays will also be created.\n",
    "\n",
    "The numpy ndarray is stored column major format so that values for each parameter are stored contiguously in memory, likewise all draws from a chain are contiguous.  Thus the dimensions of the ndarray are arranged as follows:  (draws, chains, columns):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 4, 8)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.sample.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python's index slicing operations can be used to access the information by chain.\n",
    "For example, to select all draws and all output columns from the first chain,\n",
    "we specify the chain index (2nd index dimension).  As arrays indexing starts at 0,\n",
    "the index '0' corresponds to the first chain in the ``RunSet``.\n",
    "\n",
    "The following expression selects the first 3 draws from chain 1 for the parameter `theta`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('theta', array([0.418464, 0.418464, 0.267071]))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.column_names[7], bern_fit.sample[0:3,0,7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get HMC sampler tuning parameters\n",
    "\n",
    "##### stepsize\n",
    "\n",
    "The `StanFit` property `stepsize` property is a 1-D numpy ndarray which contains the stepsize used by the sampler for each chain.  This array is created at the same time as the `sample` and `metric` arrays are created.\n",
    "\n",
    "At the end of adaptation, the stepsize for the 4 chains in this example is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.09157 , 0.945093, 1.37647 , 0.882146])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.stepsize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### metric_type, metric\n",
    "\n",
    "The `metric` property is an numpy ndarray which contains the metric used by the sampler for each chain.  This array is created at the same time as the `sample` and `stepsize` arrays are created.\n",
    "\n",
    "At the end of adaptation, the metric for the 4 chains in this example is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('diag_e', array([[0.598581],\n",
       "        [0.380448],\n",
       "        [0.458512],\n",
       "        [0.466716]]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bern_fit.metric_type,  bern_fit.metric"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}