Repository: stan-dev/cmdstanpy
Pull request #40 (2019-05-27 00:42:01)

Feature/30 extract adaptation info
#### Submission Checklist

- [x] Run unit tests
- [x] Declare copyright holder and open-source license: see below

#### Summary
Added logic to parse adaptation information from stan_csv files:

 + check metric consistency
 + parse stepsize and metric from stan_csv files
 + make available as attributes on `PosteriorSample` object

#### Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):  Columbia University



By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)




Patch:
diff --git a/cmdstanpy/cmds.py b/cmdstanpy/cmds.py
index 0e6eb6e..793d0f4 100644
--- a/cmdstanpy/cmds.py
+++ b/cmdstanpy/cmds.py
@@ -148,11 +148,13 @@ def sample(
             cores, cpu_count()))
         cores = cpu_count()
     runset = RunSet(args=args, chains=chains)
-    tp = ThreadPool(cores)
-    for i in range(chains):
-        tp.apply_async(do_sample, (runset, i))
-    tp.close()
-    tp.join()
+    try:
+        tp = ThreadPool(cores)
+        for i in range(chains):
+            tp.apply_async(do_sample, (runset, i))
+    finally:
+        tp.close()
+        tp.join()
     if not runset.check_retcodes():
         msg = 'Error during sampling'
         for i in range(chains):
diff --git a/cmdstanpy/lib.py b/cmdstanpy/lib.py
index a1c0678..3d8d5ca 100644
--- a/cmdstanpy/lib.py
+++ b/cmdstanpy/lib.py
@@ -394,11 +394,10 @@ class PosteriorSample(object):
     def __init__(self, run: Dict = None, csv_files: Tuple[str] = None) -> None:
         """Initialize object."""
         self._run = run
-        """sampler run info."""
         self._csv_files = csv_files
-        """sampler output csv files."""
         self._sample = None
-        """assembled draws across all chains, stored column major."""
+        self._metric = None
+        self._stepsize = None
         if run is None:
             raise ValueError('missing sampler run info')
         if csv_files is None:
@@ -409,6 +408,8 @@ class PosteriorSample(object):
         self._chains = run['chains']
         self._draws = run['draws']
         self._column_names = run['column_names']
+        self._num_params = run['num_params']
+        self._metric_type = run['metric']
 
     def __repr__(self) -> str:
         return 'PosteriorSample(chains={},  draws={}, columns={})'.format(
@@ -465,7 +466,7 @@ class PosteriorSample(object):
                 if not (p in self._column_names or p in pnames_base):
                     raise ValueError('unknown parameter: {}'.format(p))
         if self._sample is None:
-            self._sample = self.get_sample()
+            self.parse_sample()
         data = self._sample.reshape(
             (self._draws * self._chains), len(self._column_names), order='A'
         )
@@ -491,7 +492,7 @@ class PosteriorSample(object):
 
     @property
     def chains(self) -> int:
-        """Number of chains"""
+        """Number of chains run by sampler"""
         return self._chains
 
     @property
@@ -518,6 +519,25 @@ class PosteriorSample(object):
         """
         return self._csv_files
 
+    @property
+    def metric_type(self) -> str:
+        """Metric type, either 'diag_e' or 'dense_e'"""
+        return self._metric_type
+
+    @property
+    def metric(self) -> np.ndarray:
+        """Array of per-chain metric used by sampler."""
+        if self._metric is None:
+            self.parse_sample()
+        return self._metric
+
+    @property
+    def stepsize(self) -> np.ndarray:
+        """Array or per-chain stepsize."""
+        if self._stepsize is None:
+            self.parse_sample()
+        return self._stepsize
+
     @property
     def sample(self) -> np.ndarray:
         """
@@ -527,27 +547,53 @@ class PosteriorSample(object):
         all draws from a chain are contiguous.
         """
         if self._sample is None:
-            self._sample = self.get_sample()
+            self.parse_sample()
         return self._sample
 
-    def get_sample(self) -> np.ndarray:
+    def parse_sample(self) -> None:
         """
-        Returns posterior sample.
-        The first time this function is called it assembles the sample
-        from the stan_csv files; subsequent calls to this function
-        return the assembled sample.
+        Allocates and populates the stepsize, metric, and drawset arrays
+        by parsing the validated stan_csv files if they don't already exist.
         """
-        sample = np.empty(
+        if not (self._stepsize is None and self._metric is None and
+                    self._sample is None):
+            return
+        self._stepsize = np.empty(self._chains, dtype=float)
+        if self._metric_type == 'diag_e':
+            self._metric = np.empty(
+                (self._chains, self._num_params), dtype=float)
+        else:
+            self._metric = np.empty(
+                (self._chains, self._num_params, self._num_params),
+                dtype=float)
+        self._sample = np.empty(
             (self._draws, self._chains, len(self._column_names)), dtype=float,
             order='F'
-        )
+            )
         for chain in range(self._chains):
-            draw = 0
-            with open(self._csv_files[chain], 'r') as fd:
-                for line in fd:
-                    if line.startswith('#') or line.startswith('lp__,'):
-                        continue
-                    vs = [float(x) for x in line.split(',')]
-                    sample[draw, chain, :] = vs
-                    draw += 1
-        return sample
+            with open(self._csv_files[chain], 'r') as fp:
+                # read past initial comments, column header
+                line = fp.readline().strip()
+                while len(line) > 0 and line.startswith('#'):
+                    line = fp.readline().strip()
+                line = fp.readline().strip()  # adaptation header
+                # stepsize
+                line = fp.readline().strip()
+                label, stepsize = line.split('=')
+                self._stepsize[chain] = float(stepsize.strip())
+                line = fp.readline().strip()  # metric header
+                # metric
+                if self._metric_type == 'diag_e':
+                    line = fp.readline().lstrip(' #\t')
+                    xs = line.split(',')
+                    self._metric[chain, :] = [float(x) for x in xs]
+                else:
+                    for i in range(self._num_params):
+                        line = fp.readline().lstrip(' #\t')
+                        xs = line.split(',')
+                        self._metric[chain, i, :] = [float(x) for x in xs]
+                # draws
+                for i in range(self._draws):
+                    line = fp.readline().lstrip(' #\t')
+                    xs = line.split(',')
+                    self._sample[i, chain, :] = [float(x) for x in xs]
diff --git a/cmdstanpy/utils.py b/cmdstanpy/utils.py
index 8a1efd6..ad9e3a6 100644
--- a/cmdstanpy/utils.py
+++ b/cmdstanpy/utils.py
@@ -6,7 +6,7 @@ import os.path
 import json
 import subprocess
 import numpy as np
-from typing import Dict
+from typing import Dict, TextIO
 
 
 def validate_cmdstan_path(path: str) -> None:
@@ -128,46 +128,121 @@ def check_csv(filename: str) -> Dict:
 def scan_stan_csv(filename: str) -> Dict:
     """Process stan_csv file line by line."""
     dict = {}
-    draws_found = 0
     lineno = 0
-    with open(filename) as fp:
-        line = fp.readline().strip()
+    with open(filename, 'r') as fp:
+        lineno = scan_config(fp, dict, lineno)
+        lineno = scan_column_names(fp, dict, lineno)
+        lineno = scan_metric(fp, dict, lineno)
+        lineno = scan_draws(fp, dict, lineno)
+    return dict
+
+
+def scan_config(fp: TextIO, config_dict: Dict, lineno: int) -> int:
+    """
+    Scan initial stan_csv file comments lines and
+    save non-default configuration information to config_dict.
+    """
+    cur_pos = fp.tell()
+    line = fp.readline().strip()
+    while len(line) > 0 and line.startswith('#'):
         lineno += 1
-        while len(line) > 0 and line.startswith('#'):
-            parse_header(line, dict)
-            line = fp.readline().strip()
-            lineno += 1
-        dict['column_names'] = tuple(line.split(','))
-        num_cols = len(dict['column_names'])
+        if not line.endswith('(Default)'):
+            line = line.lstrip(' #\t')
+            key_val = line.split('=')
+            if len(key_val) == 2:
+                if key_val[0].strip() == 'file' and not key_val[1].endswith(
+                        'csv'):
+                    config_dict['data_file'] = key_val[1].strip()
+                elif key_val[0].strip() != 'file':
+                    config_dict[key_val[0].strip()] = key_val[1].strip()
+        cur_pos = fp.tell()
         line = fp.readline().strip()
-        lineno += 1
-        while len(line) > 0 and line.startswith('#'):
-            line = fp.readline().strip()
+    fp.seek(cur_pos)
+    return lineno
+
+
+def scan_column_names(fp: TextIO, config_dict: Dict, lineno: int) -> int:
+    """
+    Process columns header, add to config_dict as 'column_names'
+    """
+    line = fp.readline().strip()
+    lineno += 1
+    names = line.split(',')
+    config_dict['column_names'] = tuple(names)
+    return lineno
+
+
+def scan_metric(fp: TextIO, config_dict: Dict, lineno: int) -> int:
+    """
+    Scan stepsize, metric from  stan_csv file comment lines,
+    set config_dict entries 'metric' and 'num_params'
+    """
+    if 'metric' not in config_dict:
+        config_dict['metric'] = 'diag_e'
+    metric = config_dict['metric']
+    line = fp.readline().strip()
+    lineno += 1
+    if not line == '# Adaptation terminated':
+        raise ValueError(
+            'line {}: expecting metric, '
+            'found:\n\t "{}"'.format(lineno, line))
+    line = fp.readline().strip()
+    lineno += 1
+    label, stepsize = line.split('=')
+    if not label.startswith('# Step size'):
+        raise ValueError(
+            'line {}: expecting stepsize, '
+            'found:\n\t "{}"'.format(lineno, line))
+    try:
+        float(stepsize.strip())
+    except ValueError:
+        raise ValueError(
+            'line {}: invalid stepsize: {}'.format(
+                lineno, stepsize))
+    line = fp.readline().strip()
+    lineno += 1
+    if not ((metric == 'diag_e' and
+                 line == '# Diagonal elements of inverse mass matrix:') or
+                 (metric == 'dense_e' and
+                      line == '# Elements of inverse mass matrix:')):
+        raise ValueError(
+            'line {}: invalid or missing mass matrix '
+            'specification'.format(lineno))
+    line = fp.readline().lstrip(' #\t')
+    lineno += 1
+    num_params = len(line.split(','))
+    config_dict['num_params'] = num_params
+    if metric == 'diag_e':
+        return lineno
+    else:
+        for i in range(1,num_params):
+            line = fp.readline().lstrip(' #\t')
             lineno += 1
-        while len(line) > 0 and not line.startswith('#'):
-            draws_found += 1
-            if len(line.split(',')) != num_cols:
+            if len(line.split(',')) != num_params:
                 raise ValueError(
-                    'file {}, at line {}: bad draw, expecting {} items, '
-                    'found {}'.format(
-                        filename, lineno, num_cols, len(line.split(',')))
-                )
-            line = fp.readline().strip()
-            lineno += 1
-    dict['draws'] = draws_found
-    return dict
+                    'line {}: invalid or missing mass matrix '
+                    'specification'.format(lineno))
+        return lineno
 
 
-def parse_header(line: str, dict: Dict) -> None:
+def scan_draws(fp: TextIO, config_dict: Dict, lineno: int) -> int:
     """
-    Parse initial stan_csv file comments lines and
-    save non-default configuration information to dict.
+    Parse draws, check elements per draw, save num draws to config_dict.
     """
-    if not line.endswith('(Default)'):
-        line = line.lstrip(' #\t')
-        key_val = line.split('=')
-        if len(key_val) == 2:
-            if key_val[0].strip() == 'file' and not key_val[1].endswith('csv'):
-                dict['data_file'] = key_val[1].strip()
-            elif key_val[0].strip() != 'file':
-                dict[key_val[0].strip()] = key_val[1].strip()
+    draws_found = 0
+    num_cols = len(config_dict['column_names'])
+    cur_pos = fp.tell()
+    line = fp.readline().strip()
+    while len(line) > 0 and not line.startswith('#'):
+        lineno += 1
+        draws_found += 1
+        if len(line.split(',')) != num_cols:
+            raise ValueError(
+                'line {}: bad draw, expecting {} items, found {}'.format(
+                    lineno, num_cols, len(line.split(',')))
+                )
+        cur_pos = fp.tell()
+        line = fp.readline().strip()
+    config_dict['draws'] = draws_found
+    fp.seek(cur_pos)
+    return lineno
diff --git a/cmdstanpy_tutorial.ipynb b/cmdstanpy_tutorial.ipynb
index 1055209..c19c347 100644
--- a/cmdstanpy_tutorial.ipynb
+++ b/cmdstanpy_tutorial.ipynb
@@ -37,8 +37,10 @@
     "\n",
     "#### Analyze resulting `PosteriorSample`\n",
     "  + `PosteriorSample` contains names of csv output files, names of csv column headers\n",
-    "  + attribute `sample` assembles in-memory sample from csv files\n",
-    "  + methods `summary` and `diagnose` invoke CmdStan tools `bin/stansummary` and `bin/diagnose` respectively\n"
+    "  + methods `summary` and `diagnose` invoke CmdStan tools `bin/stansummary` and `bin/diagnose` respectively\n",
+    "  + `sample` property holds the in-memory sample\n",
+    "  + `stepsize` property holds the stepsize used by sampler for each chain\n",
+    "  + `metric` property holds the metric used by sampler for each chain"
    ]
   },
   {
@@ -47,11 +49,13 @@
    "source": [
     "### Installation\n",
     "\n",
-    "* Download this repo from GitHub\n",
+    "* Clone or download CmdStanPy from GitHub repository https://github.com/stan-dev/cmdstanpy\n",
     "\n",
-    "* run script `make_cmdstan.sh` - (need equivalent for Windows) - this installs and builds cmdstan in subdirectory `releases`\n",
+    "* CmdStanPy uses CmdStan directly to compile and run Stan programs, therefore CmdStan must be installed locally.\n",
     "\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 script `make_cmdstan.sh` 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"
    ]
   },
   {
@@ -150,12 +154,12 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "input data tempfile: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmpnn3si19v/tmpvfm5ytyn.json\n",
+      "input data tempfile: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmp546biv20/tmpktez1i7q.json\n",
       "start chain 1.  start chain 2.  \n",
       "\n",
-      "finish chain 2.  \n",
-      "start chain 3.  \n",
       "finish chain 1.  \n",
+      "start chain 3.  \n",
+      "finish chain 2.  \n",
       "start chain 4.  \n",
       "finish chain 3.  \n",
       "finish chain 4.  \n"
@@ -246,42 +250,42 @@
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>lp__</th>\n",
-       "      <td>-7.288520</td>\n",
-       "      <td>0.017688</td>\n",
-       "      <td>0.738774</td>\n",
-       "      <td>-8.849610</td>\n",
-       "      <td>-7.008580</td>\n",
-       "      <td>-6.750170</td>\n",
-       "      <td>1744.39</td>\n",
-       "      <td>8684.48</td>\n",
-       "      <td>1.00096</td>\n",
+       "      <td>-7.251000</td>\n",
+       "      <td>0.01974</td>\n",
+       "      <td>0.705775</td>\n",
+       "      <td>-8.749310</td>\n",
+       "      <td>-6.97809</td>\n",
+       "      <td>-6.749950</td>\n",
+       "      <td>1278.29</td>\n",
+       "      <td>5581.30</td>\n",
+       "      <td>1.00199</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>theta</th>\n",
-       "      <td>0.249456</td>\n",
-       "      <td>0.003236</td>\n",
-       "      <td>0.120517</td>\n",
-       "      <td>0.073465</td>\n",
-       "      <td>0.235161</td>\n",
-       "      <td>0.464975</td>\n",
-       "      <td>1386.90</td>\n",
-       "      <td>6904.73</td>\n",
-       "      <td>1.00046</td>\n",
+       "      <td>0.246911</td>\n",
+       "      <td>0.00306</td>\n",
+       "      <td>0.115661</td>\n",
+       "      <td>0.077695</td>\n",
+       "      <td>0.23484</td>\n",
+       "      <td>0.451042</td>\n",
+       "      <td>1428.82</td>\n",
+       "      <td>6238.56</td>\n",
+       "      <td>1.00078</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "text/plain": [
-       "           Mean      MCSE    StdDev        5%       50%       95%    N_Eff  \\\n",
-       "name                                                                         \n",
-       "lp__  -7.288520  0.017688  0.738774 -8.849610 -7.008580 -6.750170  1744.39   \n",
-       "theta  0.249456  0.003236  0.120517  0.073465  0.235161  0.464975  1386.90   \n",
+       "           Mean     MCSE    StdDev        5%      50%       95%    N_Eff  \\\n",
+       "name                                                                       \n",
+       "lp__  -7.251000  0.01974  0.705775 -8.749310 -6.97809 -6.749950  1278.29   \n",
+       "theta  0.246911  0.00306  0.115661  0.077695  0.23484  0.451042  1428.82   \n",
        "\n",
        "       N_Eff/s    R_hat  \n",
        "name                     \n",
-       "lp__   8684.48  1.00096  \n",
-       "theta  6904.73  1.00046  "
+       "lp__   5581.30  1.00199  \n",
+       "theta  6238.56  1.00078  "
       ]
      },
      "execution_count": 7,
@@ -323,9 +327,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### PosteriorSample.sample\n",
+    "#### PosteriorSample.extract()\n",
     "\n",
-    "The `sample` property is a 3-D numpy ndarray which contains all draws across all chains.  The 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):"
+    "The `extract` function returns a pandas.DataFrame, one draw per row.  It 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."
    ]
   },
   {
@@ -336,7 +340,10 @@
     {
      "data": {
       "text/plain": [
-       "(1000, 4, 8)"
+       "((4000, 8),\n",
+       " Index(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__',\n",
+       "        'divergent__', 'energy__', 'theta'],\n",
+       "       dtype='object'))"
       ]
      },
      "execution_count": 9,
@@ -345,17 +352,8 @@
     }
    ],
    "source": [
-    "a_sample.sample\n",
-    "a_sample.sample.shape"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### PosteriorSample.extract()\n",
-    "\n",
-    "The `extract` function flattens the `sample` 3-dimensional ndarray (draws X chains X columns) to a  pandas.DataFrame, (i.e., a 2-dimensional ndarray), one draw per row.  "
+    "samp = a_sample.extract()\n",
+    "samp.shape, samp.columns"
    ]
   },
   {
@@ -366,11 +364,7 @@
     {
      "data": {
       "text/plain": [
-       "((4000, 8),\n",
-       " Index(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__',\n",
-       "        'divergent__', 'energy__', 'theta'],\n",
-       "       dtype='object'),\n",
-       " <matplotlib.axes._subplots.AxesSubplot at 0x11f2ed438>)"
+       "<matplotlib.axes._subplots.AxesSubplot at 0x11c0e72b0>"
       ]
      },
      "execution_count": 11,
@@ -379,7 +373,7 @@
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl43Fd97/H3V7ut1dpted/jfVHsLGQnwQmQhYTGoZTADZgQoBTofZ7c0sullD4tbaElkNKGJCWsCSQhMc1CQhKykNiJvO+WvMu2Ni+SRrK1nvvHjIQsaxnJ+s1vRvN5Pc88nuXMbz6WJX/1O+f8zjHnHCIiIgAJfgcQEZHooaIgIiLdVBRERKSbioKIiHRTURARkW4qCiIi0k1FQUREuqkoiIhINxUFERHpluR3gKHKz893U6dO9TuGiEhM2bBhQ51zrmCwdjFXFKZOnUpZWZnfMUREYoqZHQqnnbqPRESkm4qCiIh0U1EQEZFuKgoiItJNRUFERLqpKIiISDcVBRER6RZz1ylIfDrb1sEru2o4eKKJiePGcN1FRWSk6ttXZKTpp0qi3rbKej7/i40cPtnc/Vx+RgrfunUBqxaM9zGZyOij7iOJauXVjfz5w+vo6HT896cuZtc3V/HEmkuYkDOGe3+2kYff3O93RJFRRUVBolZLewdf/OUmUpISeHzNJVwzp5AxKYmsnJ7Hrz57KTcuKOZbz+3iyQ2VfkcVGTVUFCRqPfb2QXZXNfJPH1nEpNyx57yWlpzI91Yv5fKZedz/1FY2HzntU0qR0UVFQaJSw9k2HnxtH1fNLuD984r6bJOSlMB/fGw5RVlpfOnxTTSebYtwSpHRR0VBotKvyyqpP9PGV2+YPWC77LHJ/PvqJRw52czXn90RoXQio5eKgkSdzk7Hz9YdYvmUcSyamDNo+4un5vKFa2fxm01HeXF7VQQSioxeKgoSdd6qqONAXROfuHRK2O/54rUzmTc+i799Zjunmlo9TCcyuqkoSNR5dvMxstKSWLWgOOz3JCcm8C8fXcTp5lb+31p1I4kMl4qCRJWW9g5e2lnFDfOLSU1KHNJ750/I5gvXzmTtlmPqRhIZJhUFiSpvldfReLadDy4a3pXKn7+mqxtpGyfVjSQyZCoKElVe3F5FVloSl8/IH9b7kxMT+NePLuZ0c5u6kUSGwbOiYGZpZvaumW0xsx1m9nd9tEk1syfMrMLM1pvZVK/ySPRzzvFGeS1XzC4gJWn435rzJmTxxWtn8dstx3hx+/ERTCgy+nl5ptACXOucWwwsAVaZ2SW92twDnHLOzQT+Dfi2h3kkypXXBKhuaOHKWcM7S+jpvmtmMH9CcDaSupFEwudZUXBBgdDD5NDN9Wp2C/BY6P6TwHVmZl5lkuj2xt5aAN43q+CCj9XVjVR/Rt1IIkPh6ZiCmSWa2WagBnjZObe+V5MS4AiAc64dqAfyvMwk0euN8jpmFKRTkjNmRI530Xh1I4kMladFwTnX4ZxbAkwEVpjZguEcx8zWmFmZmZXV1taObEiJCmfbOli//wRXjMBZQk+fu3oGC0qy+NtndtDU0j6ixxYZjSIy+8g5dxp4DVjV66WjwCQAM0sCsoETfbz/IedcqXOutKBgZP/TkOiw5chpWto7uXzmhY8n9JScmMA3b1lAXaCFH2nvBZFBeTn7qMDMckL3xwDXA7t7NVsL3B26fwfwqnOu97iDxIGyQ6cAKJ0ybsSPvWzyOG5cUMxDb+znRKBlxI8vMpp4eaYwHnjNzLYC7xEcU/gfM/ummd0cavMIkGdmFcBXgPs9zCNRbMOhU8woSGdceoonx//K9bNpbu3gp+sOeXJ8kdHCsz2anXNbgaV9PP/1HvfPAh/1KoPEhs5Ox4ZDp1g1P/y1joZqVlEm180t5CfvHOKzV85gTMrQltAQiRe6oll8t682QP2ZNpZPHfmuo57WXDmdk02trN1y1NPPEYllKgriOy/HE3paMS2XmYUZ/KpMezqL9EdFQXy34dAp8tJTmJaf7unnmBkfXT6RDYdOsa82MPgbROKQioL4bmvlaRZPyiESF7PftqyExATjyQ06WxDpi4qC+Kq5tZ2KmgALSrIj8nmFmWlcNiOPF7YdR7OfRc6noiC+2nW8gU4HiyJUFABWLSjm4Ilmdlc1RuwzRWKFioL4amtlPQALJ0auKNwwrxgzeEG7s4mcR0VBfLXtaD0FmakUZaVF7DMLMlO5eGouL+1QURDpTUVBfLWtsj6iXUddrp1byO6qRqobzkb8s0WimYqC+Ka5tZ19tZEbZO7pqtnBhRVf36tVd0V6UlEQ3+w8FhxkXuhDUZhbnElhZqqKgkgvKgrim21HIz/I3MXMuGp2AW+V19He0RnxzxeJVioK4ptdxxvIz0iJ6CBzT1fNKaD+TBtbKk/78vki0UhFQXyzu6qRucVZvn3+ZTOCG/q8s++8fZ1E4paKgviio9Oxp6qROcWZvmXITU9hbnEm6/af9C2DSLRRURBfHDrRREt7J3N9LAoAl0zPY8OhU7S2a1xBBFQUxCddS0xcNN6/7iOAldNyOdPWwbajGlcQARUF8cnu4w0kGMwszPA1x4ppuQDqQhIJUVEQX+yuamRafjppyf5ui5mXkcqcokzW7ddgswioKIhP/J551NPK6blsOHSKNl2vIKKiIJEXaGnn8Mlm3weZu1wyPY/m1o7ui+lE4plnRcHMJpnZa2a208x2mNmX+mhztZnVm9nm0O3rXuWR6LG3OjjIPNfnQeYufxpXUBeSiJdnCu3AV51z84BLgM+b2bw+2r3pnFsSun3TwzwSJXYfDxWFKDlTyM9IZWZhBu8e0GCziGdFwTl33Dm3MXS/EdgFlHj1eRI79lQ1kJGaREnOGL+jdFs5LZeyg6e0DpLEvYiMKZjZVGApsL6Ply81sy1m9oKZzY9EHvHXrqpGZhdlkJBgfkfptnJ6HoGWdnYeb/A7ioivPC8KZpYBPAX8lXOu90/cRmCKc24x8H3gmX6OscbMysysrLZWSx3HMue6lreIjvGELitD4wrrdb2CxDlPi4KZJRMsCD93zj3d+3XnXINzLhC6/zyQbGb5fbR7yDlX6pwrLSgo8DKyeKw20EL9mTZm+XzRWm9FWWlMzRvLeo0rSJzzcvaRAY8Au5xz3+2nTXGoHWa2IpRHU0BGsYqaAOD/lcx9WTktj/cOnqSz0/kdRcQ3Xp4pXA78BXBtjymnN5nZvWZ2b6jNHcB2M9sCPACsds7pJ3IU2xfFRWHFtFzqz7R1r8skEo+SvDqwc+4tYMCRROfcD4AfeJVBos++2ibSUxIZn+3PxjoDWTk9OK7w7oETzJsQXWMeIpGiK5oloipqAswozCDUaxhVJo4bS0nOGI0rSFxTUZCIqqgJMLMg+rqOuqycnsu7B06iXkyJVyoKEjGNZ9uoajjLjCgcT+iyclouJ5pa2Vcb8DuKiC9UFCRi9tU2AdE5yNxl5bQ8QPsrSPxSUZCIiebpqF2m5I2lKCtV4woSt1QUJGIqagIkJxpTcsf6HaVfZsaKaXm8e+CExhUkLqkoSMRU1ASYmpdOUmJ0f9utnJZLdUMLh040+x1FJOKi+6dTRpV9tYGo7jrqcknoeoX1B3RxvcQfFQWJiJb2Dg6daIqJojCjIIO89BQtjidxSUVBIuJgXTOdLroHmbuYGZfNzOeN8jqtgyRxR0VBIqJr5tGMKL5wradr5hRQF2jR/goSd1QUJCIqagKYxU5RuHJ2cIn213bX+JxEJLJUFCQiKmoDlOSMYUxKot9RwpKfkcriidm8tkdFQeKLioJEREVNbMw86unqOYVsOnKaU02tfkcRiRgVBfFcR6djf210L4TXl6vnFOAcvFGuLWAlfqgoiOeOnjpDS3tnzJ0pLJqYQ35GCi/tqPY7ikjEqCiI5ypqgzuZxVpRSEwwPjC/mFd313CmtcPvOCIRoaIgnouFhfD688FF4znT1qEBZ4kbKgriuYqaAPkZKeSMTfE7ypCtnJZHfkYKz2097ncUkYhQURDPVdQEYub6hN4SE4wbF4znld3VBFra/Y4j4jkVBfGUcy4mp6P2dOvSEs62dfLc1mN+RxHxnIqCeKou0ErD2faYPVMAWDY5h1mFGTz+3hG/o4h4zrOiYGaTzOw1M9tpZjvM7Et9tDEze8DMKsxsq5kt8yqP+COWB5m7mBl3XjyJTYdPs6eq0e84Ip7y8kyhHfiqc24ecAnweTOb16vNjcCs0G0N8EMP84gPKmpjvygAfGTZRJITjcffO+x3FBFPeVYUnHPHnXMbQ/cbgV1ASa9mtwA/cUHrgBwzG+9VJom8fTUB0lMSGZ+d5neUC5KbnsIN84t5euNRmls14CyjV0TGFMxsKrAUWN/rpRKgZ0dtJecXDsxsjZmVmVlZba2WHIglFTUBZhRmYGZ+R7lgn7xsKvVn2nhq41G/o4h4xvOiYGYZwFPAXznnhrU4vXPuIedcqXOutKCgYGQDiqcqamJvzaP+lE4Zx+KJ2Tz61gFtviOjlqdFwcySCRaEnzvnnu6jyVFgUo/HE0PPySjQeLaNqoazzIjx8YQuZsb/et80DtQ16QpnGbW8nH1kwCPALufcd/tpthb4RGgW0iVAvXNOl46OEvtqm4DYH2Tu6aaF4xmfncYjbx3wO4qIJ7w8U7gc+AvgWjPbHLrdZGb3mtm9oTbPA/uBCuBHwH0e5pEIGw3TUXtLTkzg7sum8va+E+w4Vu93HJERl+TVgZ1zbwEDji465xzwea8yiL8qagIkJxpTcsf6HWVE3bViMg+8Us6jbx3kO3+22O84IiMqrDMFM3vazD5oZroCWsJWURNgal46SYmj69sme0wyH10+kbVbjlLTcNbvOCIjKtyf1v8APgaUm9k/mdkcDzPJKLGvNrbXPBrIpy6fRnun46frDvkdRWREhVUUnHO/d879ObAMOAj83szeNrNPhWYYiZyjpb2DQyeaRm1RmJqfzvsvKuLn6w9ztk0b8MjoEfZ5vZnlAZ8EPg1sAr5HsEi87EkyiWkH65rpdKNrkLm3e943jZNNrfxmk2ZRy+gR7pjCb4A3gbHAh51zNzvnnnDOfREYvT/1MmxdM49ieXXUwayclsv8CVk88tYBgnMmRGJfuGcKP3LOzXPO/WPXdQRmlgrgnCv1LJ3ErPKaRsxGd1EwM+553zQqagK8vlfLr8joEG5R+FYfz70zkkFkdCmvCTA5dyxjUhL9juKpDy2aQGFmKj9++6DfUURGxIDXKZhZMcEF6saY2VL+dN1BFsGuJJE+VVQHmDWKxxO6pCQlsHrFZL7/ajlHTjYzaZRdkyHxZ7AzhQ8A/0pwTaLvAt8J3b4C/I230SRWtXV0sr8uwMzCTL+jRMSdFweX7/pVmXZmk9g34JmCc+4x4DEzu90591SEMkmMO3SimbYOFxdnCgAlOWO4enYBT7x3hC9dN2vUXawn8WXA714z+3jo7lQz+0rvWwTySQyqqAluWTm7KD7OFAA+tnIKNY0tvLpbq6dKbBvsV5r00J8ZQGYfN5HzlFeHpqMWpg/ScvS4Zk4BRVmp/OJdbdcpsW2w7qP/Cv35d5GJI6PB3poAE8eNYWyKZ+stRp2kxATuLJ3E91+r4OjpM5TkjPE7ksiwhHvx2j+bWZaZJZvZK2ZW26NrSeQc5dWNcTOe0NMdyyfhHDyjK5wlhoU7InZDaCvNDxFc+2gm8L+9CiWxq72jk/11TXE1ntBlct5YLp46jqc3VuoKZ4lZ4RaFrn6ADwK/ds5pdxHp0+GTzbS2d47qNY8G8pFlE9lX28S2o/oRkdgUblH4HzPbDSwHXjGzAkALyct5ykNrHs2KwzMFCG7XmZKUwNMb1YUksSncpbPvBy4DSp1zbUATcIuXwSQ2jcYtOIcie0wy188rYu2WY7S2d/odR2TIhnKVzVzgTjP7BHAHcIM3kSSWlVc3UpIzhozU+Jl51Nvty0o42dSqRfIkJoX1k2tmPwVmAJuBrh1FHPATj3JJjNpbPXp3WwvXFbMKyM9I4emNlVw/r8jvOCJDEu6vc6XAPKcpFTKAjk7HvtoAl83I8zuKr5ITE/jw4gn8fN1h6s+0kT1GmxNK7Ai3+2g7UDyUA5vZo2ZWY2bb+3n9ajOrN7PNodvXh3J8iT6Vp5ppae+My+movd22tITWjk5e2Hbc7ygiQxLumUI+sNPM3gVaup50zt08wHt+DPyAgbuY3nTOfSjMDBLlupa3mFkU391HAAtLsplekM5vNh1l9YrJfscRCVu4ReEbQz2wc+4NM5s61PdJ7NobWggv3scUILgr221LSvjOy3u17IXElHCnpL5O8Erm5ND994CNI/D5l5rZFjN7wczmj8DxxEcV1QGKs9LISlMfOsAtS0oAeHazrlmQ2BHu2kefAZ4E/iv0VAnwzAV+9kZginNuMfD9gY5nZmvMrMzMymprNc0vWpXXBJilrqNuk/PGUjplHL/ZeFTLXkjMCHeg+fPA5UADgHOuHCi8kA92zjU45wKh+88DyWaW30/bh5xzpc650oKCggv5WPFIZ6ejoibArDjZbS1cty4tobwmwM7jDX5HEQlLuEWhxTnX2vXAzJIIXqcwbGZWbGYWur8ilOXEhRxT/HP09BnOtHXoTKGXDy4cT3KiaeVUiRnhFoXXzexvgDFmdj3wa+C3A73BzH4JvAPMMbNKM7vHzO41s3tDTe4AtpvZFuABYLWug4hd5aFB5nhcMnsg49JTuHpOIWu3HKOjU9/eEv3CnX10P3APsA34LPA88PBAb3DO3TXI6z8gOGVVRoGu6ajqPjrfbUtLeHlnNev2n+DymX32kIpEjbCKgnOu08yeAZ5xzmmkV86ztzpAYWYq2WM186i3a+cWkpmaxG82HVVRkKg3YPeRBX3DzOqAPcCe0K5ruvpYzlFR06jxhH6kJSdy48JiXtxexZnWjsHfIOKjwcYUvkxw1tHFzrlc51wusBK43My+7Hk6iQnOueB0VHUd9evWpSUEWtr5/a5qv6OIDGiwovAXwF3OuQNdTzjn9gMfBz7hZTCJHcfqz9LcqplHA7lkWh7js9M0C0mi3mBFIdk5V9f7ydC4gjqPBYC91V0zj3Sm0J+EBOPmJRN4fW8tJwItg79BxCeDFYXWYb4mcaSie+aRzhQGctvSEto7Hc9p5VSJYoMVhcVm1tDHrRFYGImAEv3KaxrJz0hhXHqK31Gi2tziLOYWZ6oLSaLagEXBOZfonMvq45bpnFP3kQBokHkIbltawsbDpzl0osnvKCJ9GsoezSLncc5RUa2F8MJ185IJmMEzm475HUWkTyoKckGqGs7S2NKu8YQwjc8ew6XT83hms1ZOleikoiAXpHt5C23BGbZbl5ZwoK6JLZX1fkcROY+KglyQ8hrNPBqqVQuKSU1K0ICzRCUVBbkg5dWN5KankJeR6neUmJGVlsz75xXx2y3HaOvo9DuOyDlUFOSClNcEtCfzMNy2pIQTTa28VX7etaEivlJRkGFzzlFe3chszTwasitnF5AzNplntH+zRBkVBRm22sYWGs626xqFYUhJSuCmheN5aUc1za3tfscR6aaiIMO2V8tbXJBbl5Rwpq2Dl3dq5VSJHioKMmxdW3DOVPfRsJROGceE7DSe3awL2SR6qCjIsJXXBMgZm0yBZh4NS0KC8eElE3hjby0nm7S+pEQHFQUZtvLqRmYVZmBmfkeJWbcu0cqpEl1UFGRYunZb03TUCzO3OJPZRRk8qwvZJEqoKMiw1AVaOd3cpplHF8jMuGVJCWWHTlF5qtnvOCLeFQUze9TMasxsez+vm5k9YGYVZrbVzJZ5lUVGXnlot7XZWvPogn1w4XgAfrdDs5DEf16eKfwYWDXA6zcCs0K3NcAPPcwiI6x7zSPNPLpgU/PTmVucyYvbNa4g/vOsKDjn3gBODtDkFuAnLmgdkGNm473KIyOrvKaRrLQkCjM182gkrFpQTNmhU9Q0nvU7isQ5P8cUSoAjPR5Xhp47j5mtMbMyMyurra2NSDgZ2N7qALOKMjXzaITcuGA8zsFL6kISn8XEQLNz7iHnXKlzrrSgoMDvOAJU1AR0JfMIml2UwbT8dF7cXuV3FIlzfhaFo8CkHo8nhp6TKFcXaOFkU6s21hlBZsaqBcW8s/8Ep5t1IZv4x8+isBb4RGgW0iVAvXNOI20xoFxrHnli1fxiOjqd1kISX3k5JfWXwDvAHDOrNLN7zOxeM7s31OR5YD9QAfwIuM+rLDKyKmo0HdULiyZmU5Izht/tUBeS+CfJqwM75+4a5HUHfN6rzxfv7K0OkJmaRFGWZh6NJDPj+nlF/OLdwzS1tJOe6tmPp0i/YmKgWaJLeU0jM4u05pEXPjC/mNb2Tl7fq1l24g8VBRmyipoAs7W8hScunjqOcWOT1YUkvlFRkCE52dRKXaBVVzJ7JCkxgevnFfHq7hpa2zv9jiNxSEVBhqRrzSNNR/XOB+YX03i2nXf2n/A7isQhFQUZkr01mo7qtctn5pOekqgL2cQXKgoyJBXVjWSkJjE+O83vKKNWWnIiV88p5OWd1XR0Or/jSJxRUZAh6dpYRzOPvHXD/CLqAi1sOnzK7ygSZ1QUZEj2VmvNo0i4Zm4hyYmmWUgScSoKErZTTa3UBVp0JXMEZKUlc/nMfH63o5rgdZ4ikaGiIGGrqA0OMs/UdNSI+MD8Yg6fbGZ3VaPfUSSOqChI2PZ2TUdV91FEvP+iIsxQF5JElIqChK28OkB6SiIlOWP8jhIXCjJTKZ0yjhe2VakLSSJGRUHCVl7TqJlHEXbz4gnsqW5kx7EGv6NInFBRkLCVh7bglMj58OIJpCQm8NTGSr+jSJxQUZCw1De3UdPYovGECMsZm8J1FxWydvMx2jq0FpJ4T0VBwlJe07XmkYpCpN2+bCInmlp5ZVeN31EkDqgoSFjKu9c8UvdRpF09p4AJ2Wn8+O0DfkeROKCiIGHZU9XIWM088kVSYgJ3XzaVdftPsuNYvd9xZJRTUZCw7KlqZE5xJgkJmnnkh9UrJjM2JZFH3tLZgnhLRUEG5Zxjd1UDc4vVdeSX7DHJ3HnxJJ7dfIz9oSvLRbygoiCDqmls4VRzG3OLs/yOEtfuu3omqUkJfOelvX5HkVHM06JgZqvMbI+ZVZjZ/X28/kkzqzWzzaHbp73MI8Oz63jwwimdKfirIDOVT18xnee2HWfzkdN+x5FRyrOiYGaJwIPAjcA84C4zm9dH0yecc0tCt4e9yiPD17Ugm84U/LfmyukUZqbyf57epusWxBNenimsACqcc/udc63A48AtHn6eeGRPVSMTstPIHpvsd5S4l5GaxLduXcCu4w385x/2+R1HRiEvi0IJcKTH48rQc73dbmZbzexJM5vkYR4Zpl3HG5ijrqOoccP8Yj60aDwPvFqundlkxPk90PxbYKpzbhHwMvBYX43MbI2ZlZlZWW1tbUQDxrvW9k721QaYO15dR9HkH25dSFFWGl/4xSZON7f6HUdGES+LwlGg52/+E0PPdXPOnXDOtYQePgws7+tAzrmHnHOlzrnSgoICT8JK3/bXBWjrcBpkjjLZY5N58GPLqGk8y1//eouW1pYR42VReA+YZWbTzCwFWA2s7dnAzMb3eHgzsMvDPDIMXTOPLtKZQtRZPCmHr910Eb/fVcOP3tzvdxwZJZK8OrBzrt3MvgD8DkgEHnXO7TCzbwJlzrm1wF+a2c1AO3AS+KRXeWR4tlbWMyY5ken56X5HkT7cfdlU3j14km+/uIflU8axfEqu35EkxlmsnXaWlpa6srIyv2PEjTt++DYAT37uMp+TSH8azrbx4e+/RWt7J8/95RXkpqf4HUmikJltcM6VDtbO74FmiWLtHZ3sONbAwonZfkeRAWSlBccXTgRa+cqvNtPZGVu/6El0UVGQfu2rbeJMWweLVBSi3oKSbP7vh+fxhz21PKTxBbkAKgrSry2VwaUUFk3M8TmJhOPjKydz44JivvPSHrYf1RLbMjwqCtKvbZX1ZKQmMS1Pg8yxwMz4x48sJC89lb98fBPNre1+R5IYpKIg/dp6tJ4FJVnaQyGG5IxN4bt/tpgDdU186znN8JahU1GQPp1p7WDH0XqWTBrndxQZostm5rPmyun8Yv1hfrejyu84EmNUFKRPm46cor3TsXKa5r3Hoq9eP4cFJVnc/9RWqhvO+h1HYoiKgvTp3QMnMYNlU3SmEItSkhL43uqlnGnr4K9/vUXTVCVsKgrSp/cOnuSi4iyyx2i57Fg1oyCDr39oPm+W1/HoH7W3s4RHRUHO09bRycZDp1mhrqOYd9eKSVw/r4h/fnEPO45pmqoMTkVBzrPh0CnOtHVwyfQ8v6PIBTIzvn37InLGJvOlxzfTcLbN70gS5VQU5Dyv7a4hOdF436x8v6PICMhNT+Hf71zCwbomPvXf79HUousXpH8qCnKe1/bUsGJaLhmpni2iKxF22cx8HrhrKZsOn+Ljj6ynLtAy+JskLqkoyDkqTzWztzrANXMK/Y4iI+ymheN58GPL2HmsgVsf/KO28pQ+qSjIOV7cHrzY6dq5Kgqj0Y0Lx/PEZy/FObjjP9/h317eS2t7p9+xJIqoKMg5ntl8lIUl2UwvyPA7inhkyaQcXvirK7hl8QS+90o5Nz3wJm9X1PkdS6KEioJ021vdyPajDdy6tMTvKOKxrLRkvnvnEh79ZCkt7R187OH1fPGXm6iq19XP8U4jidLtsbcPkpKUwK1LJvgdRSLk2rlFXDYjnx/+YR8/fH0fr+yq5tPvm8ZnrpxOZpouXIxHOlMQAE4EWnhqYyW3LSkhLyPV7zgSQWnJiXz5+tm8/OUruWZuIQ+8WsFV//IHHn5zPwFNX407KgoCwAOvlNPW4fjMldP9jiI+mZKXzoMfW8baL1zOvPFZfOu5XVz6j6/wD8/t5NCJJr/jSYSo+0jYcuQ0P1t/mNUXT2JmoQaY492iiTn87NMr2XzkNI+8dYBH/3iQH715gGWTc7hlSQk3LRxPQabOJkcrc8671RPNbBXwPSAReNg590+9Xk8FfgIsB04AdzrnDg50zNLSUldWVuZN4DhU29jC7T98m/aOTl740pVkj1U/spzreP0Zntld4JpgAAAIIUlEQVR0jGc3H2V3VSNmwcJx7ZxCrruokPkTsjDTRkzRzsw2OOdKB23nVVEws0RgL3A9UAm8B9zlnNvZo819wCLn3L1mthq4zTl350DHVVEYOXurG/nczzZw7PRZfv6ZlSybrGWyZWB7qhp5aUcVr+6pYfOR0zgH+RmprJyey8ppuVw8NZc5RZnarS8KhVsUvOw+WgFUOOf2hwI9DtwC7OzR5hbgG6H7TwI/MDNzXp6+xLH2jk6O159lx7F6XtxexW+3HidnTDI//tTFKggSljnFmcwpzuSL182iLtDC63tqeaO8lvcOnOS5rccBSEtOYEZBBrMKM5iUO5a89BTyMlLJTEsiNSmR1OQEUhITSEtOIDkxgaTEBJITjKTEBBITjOREIykhgeRE0xmID7wsCiXAkR6PK4GV/bVxzrWbWT2QB4z4lTSv763l7/9nJz3rjet1p2cl6mrnuh/3eC30bNdzfZWwIb3/nPedG6rvNgMcu9drPd/Y3NZBR2izlczUJD5x6RTuu3qm+odlWPIzUrl9+URuXz4RCC6R8u6Bk+w81kB5TYD3Dp5i7ZZjXMj+PokJFiwUoaKRnBh83FU0EhPiq3CsvngSn77C28kgMTHQbGZrgDUAkydPHtYxMlKTmFOUGTpgj2P/6TN6v0TX91rvNue0625j57zn3Ped2+ac53ofqM/PHeD9fWQ7v03wTnpqIiU5Y5hRmMGSSTkkJ2rymYyciePGMnHcWD6y7E/PdXQ6Tje3crKplYazbbS0d9IaurW0d9LW0Ul7h6OtM/RnRycdnY72Tnfeax39PBdP8iMwXdzLonAUmNTj8cTQc321qTSzJCCb4IDzOZxzDwEPQXBMYThhlk8Zx3JtLSkSUYkJRl5Gqq59iSFe/qr4HjDLzKaZWQqwGljbq81a4O7Q/TuAVzWeICLiH8/OFEJjBF8AfkdwSuqjzrkdZvZNoMw5txZ4BPipmVUAJwkWDhER8YmnYwrOueeB53s99/Ue988CH/Uyg4iIhE8jjSIi0k1FQUREuqkoiIhINxUFERHppqIgIiLdPF0l1QtmVgscAvLxYDkMjyirN5TVG8rqDb+zTnHOFQzWKOaKQhczKwtnxb9ooKzeUFZvKKs3YiWruo9ERKSbioKIiHSL5aLwkN8BhkBZvaGs3lBWb8RE1pgdUxARkZEXy2cKIiIywmKmKJhZrpm9bGbloT/73RzBzLLMrNLMfhDJjD0+f9CsZrbEzN4xsx1mttXMBtyb2oOMq8xsj5lVmNn9fbyeamZPhF5fb2ZTI5mvV5bBsn7FzHaGvo6vmNkUP3KGsgyYtUe7283MmZlvs1HCyWpmfxb62u4ws19EOmOPHIN9D0w2s9fMbFPo++Amn3I+amY1Zra9n9fNzB4I/T22mtmyvtr5yjkXEzfgn4H7Q/fvB749QNvvAb8AfhCtWYHZwKzQ/QnAcSAnQvkSgX3AdCAF2ALM69XmPuA/Q/dXA0/49LUMJ+s1wNjQ/c9Fc9ZQu0zgDWAdUBqtWYFZwCZgXOhxYRRnfQj4XOj+POCgT1mvBJYB2/t5/SbgBYKbJl4CrPcj50C3mDlTAG4BHgvdfwy4ta9GZrYcKAJeilCuvgya1Tm31zlXHrp/DKgBBr2wZISsACqcc/udc63A4wQz99Tz7/AkcJ35sxnuoFmdc68555pDD9cR3OXPD+F8XQH+Hvg2cDaS4XoJJ+tngAedc6cAnHM1Ec7YJZysDsgK3c8GjkUw359COPcGwb1h+nML8BMXtA7IMbPxkUkXnlgqCkXOueOh+1UE/+M/h5klAN8B/jqSwfowaNaezGwFwd+A9nkdLKQEONLjcWXouT7bOOfagXogLyLp+skR0lfWnu4h+JuYHwbNGuoumOScey6SwfoQztd1NjDbzP5oZuvMbFXE0p0rnKzfAD5uZpUE93D5YmSiDdlQv58jztNNdobKzH4PFPfx0td6PnDOOTPra9rUfcDzzrlKr3+pHYGsXccZD/wUuNs51zmyKeOLmX0cKAWu8jtLX0K/tHwX+KTPUcKVRLAL6WqCZ19vmNlC59xpX1P17S7gx86575jZpQR3dFygn6mhi6qi4Jx7f3+vmVm1mY13zh0P/Ufa16nspcAVZnYfkAGkmFnAOdfvgJ+PWTGzLOA54GuhU8lIOQpM6vF4Yui5vtpUmlkSwVPyE5GJ12eOLn1lxczeT7AgX+Wca4lQtt4Gy5oJLAD+EPqlpRhYa2Y3O+fKIpYyKJyvayXBPu824ICZ7SVYJN6LTMRu4WS9B1gF4Jx7x8zSCK415FeXV3/C+n72Uyx1H60F7g7dvxt4tncD59yfO+cmO+emEuxC+okXBSEMg2Y1sxTgNwQzPhnBbBD8oZ5lZtNCOVYTzNxTz7/DHcCrLjRSFmGDZjWzpcB/ATf72O8Ng2R1ztU75/Kdc1ND36PrCGaOdEEYNGvIMwTPEjCzfILdSfsjGTIknKyHgesAzOwiIA2ojWjK8KwFPhGahXQJUN+jqzk6+D3SHe6NYH/2K0A58HsgN/R8KfBwH+0/iX+zjwbNCnwcaAM297gtiWDGm4C9BMcxvhZ67psE/5OC4A/Vr4EK4F1guo//9oNl/T1Q3ePruDZas/Zq+wd8mn0U5tfVCHZ37QS2AaujOOs84I8EZyZtBm7wKecvCc4kbCN4pnUPcC9wb4+v6YOhv8c2P//9+7vpimYREekWS91HIiLiMRUFERHppqIgIiLdVBRERKSbioKIiHRTURARkW4qCiIi0k1FQUREuv1/t9Syt238cs8AAAAASUVORK5CYII=\n",
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl83XWd7/HXJ3vbpE2apAtt0nRJW1qQLqHsCLIIqDCOqDAjiIMy7oyO3uE69zozOo+rzsx1xasgKqKjoihaFVD2VbBpSwvd9z1blzR7k5zP/eOchDRNm5M0v/M7J+f9fDzOo2f55XfezSPJ53yX3/dr7o6IiAhARtgBREQkeagoiIhILxUFERHppaIgIiK9VBRERKSXioKIiPRSURARkV4qCiIi0ktFQUREemWFHWCoSkpKvKKiIuwYIiIpZeXKlQ3uXjrYcSlXFCoqKqiurg47hohISjGzXfEcp+4jERHppaIgIiK9VBRERKSXioKIiPRSURARkV4qCiIi0ktFQUREeqXcdQoija2dPLO5jt0HWykcm81l8yZRNnFs2LFERgUVBUkZkYhzz3PbufupLbQc6+593mwdNy8r53+/bQFjcjJDTCiS+lQUJCUc64rwsZ+u4vH1tVy9YDIfvXwOC6aOp6axnftf2skPX9rBppomHvi7ZYzL1Y+1yHBpTEGSnrvzmV+u4fH1tXz+7Qu455alLCorJCcrg/LisXz+HQu4++YlrN59mH/8xRrcPezIIilLRUGS3k9e2c3yNfv57Fvn8XcXz8TMTjjmbW+ayueuO5PH1tXwwxd3Jj6kyCihoiBJbffBVv7PHzZw6dxSPnrZ7FMee/vFM3nL/En81582sfdwa4ISiowuKgqS1L782AbM4CvvOnvAFkJfZsYX/+osAL74+/WJiCcy6qgoSNJavfswj7xWwx2XzmLqhDFxfc20wjF8+M2z+eO6WlbuOhxwQpHRR0VBktY3ntxC8bgcPnTJrCF93e0Xz6QkP5evPLZRg84iQ6SiIElpS20Tz2yq5/0XVgx5ium43CzuvGIOf9lxiGc31weUUGR0UlGQpHTf8zvIy87gfefPGNbXv/fccqYXjeGrj29Wa0FkCFQUJOkcajnGw6/u411LpjNxXM6wzpGTlcEn31LJ2r2NPLmhboQTioxeKgqSdH69ai/HuiLcekHFaZ3nnUumMaN4LF99fDORiFoLIvFQUZCk4u78onoP55QVMm9KwWmdKzszgzuvqGT9gaP8cV3NCCUUGd1UFCSprNnbyObaZt5TNX1EznfDomnMKh3H157YTLdaCyKDUlGQpPKL6j3kZWfwjnPOGJHzZWYY/3DlXDbXNvOH1w6MyDlFRrPAioKZ5ZnZX8xsjZmtM7N/G+CYXDN70My2mtkrZlYRVB5Jfh1d3fxuzX6uO2sq4/OyR+y8bz97KvMmF/D1JzbT1R0ZsfOKjEZBthQ6gLe4+znAIuAaMzu/3zG3A4fdfQ7wNeArAeaRJPfspnqa2ru4ftHItBJ6ZGQYn7qqku31Lfz21f0jem6R0SawouBRzbGH2bFb/07dG4Afxe4/BFxhgy1wI6PW79ceoGhsNhfNKRnxc1+9YAoLpo7nG09uoVOtBZGTCnRMwcwyzexVoA543N1f6XfINGAPgLt3AY1AcZCZJDm1HevmiQ21XHPWFLIzR/7HMiPD+PRVc9l9qJVfrdw74ucXGS0CLQru3u3ui4DpwDIzO2s45zGzO8ys2syq6+u1bMFo9PSmOlqPdfP2N41s11FfV5w5iXPKCvn6E1to67Odp4i8ISGzj9z9CPA0cE2/l/YBZQBmlgVMAA4O8PX3unuVu1eVlpYGHVdC8Pu1+ynJz+G8mRMDew8z43PXzqfmaDvff2F7YO8jksqCnH1UamaFsftjgKuAjf0OWw68P3b/RuAp10I1aaf1WBdPbazj2rOmkhVA11Ff580q5qoFk/nOM9uob+oI9L1EUlGQv4FTgafNbC2wguiYwu/N7Atmdn3smO8DxWa2Ffg0cFeAeSRJvbClgfbOCNecNSUh73fXtfNp74rwjSc3J+T9RFLJ0NYkHgJ3XwssHuD5z/e53w68O6gMkhqe3FBHQW4W51YE13XU1+zSfP5mWTk//ctubrtwJnMm5SfkfUVSga5ollBFIs5Tm+q4dF4pOVmJ+3G888pKxmRn8uVH+/doiqQ3FQUJ1Wv7Gqlv6uCK+ZMS+r4l+bn8/aWzeGJDLa/va0zoe4skMxUFCdWTG+vIMLhsXmKLAsCtF1YwLieT+57XTCSRHioKEqqnNtaypLxo2JvpnI4JY7K5aVk5v1t7gP1H2hL+/iLJSEVBQnOo5Riv7zvKm+eGd+3JbRdW0B1xHtJVziKAioKE6OXt0esULwxgraN4lU0cy0Vzinlo5V7tziaCioKE6MWtDeTnZnHO9Amh5nj30jJ2H2rllR2HQs0hkgxUFCQ0f952kGUzJwZ+FfNgrjlrCuNyMvndWi2rLaKiIKE40NjG9oYWLpwd/qK4edmZXDZ/En9aV6stOyXtqShIKP68LTaeMDu88YS+rlk4hYbmDlbtPhx2FJFQqShIKKp3HaYgL4v5UwrCjgLA5fMnkZOZwWOv14QdRSRUKgoSilW7DrOorJCMjOTYaC8/N4sLZhfzzKa6sKOIhEpFQRKuqb2TTbVNLJ1RFHaU41xSWcK2+hZdyCZpTUVBEm7NnkbcYUl5shWF6EV0L2xpCDmJSHhUFCThVu0+jBksKi8MO8px5k7OZ1JBLs9t0Zavkr5UFCThVu0+TOWkfMbnZYcd5ThmxsWVJby4tUFXN0vaUlGQhHJ3Vu8+knRdRz0uqSzhcGsn6w8cDTuKSChUFCSh9h5uo7Gtk7Omhbu0xcmcNzN6Md2KnVryQtKTioIk1Lr90U/gC88YH3KSgZ1ROIZphWNUFCRtqShIQq3f30iGwfwpyVkUAM6tKGLFzsO4a1xB0o+KgiTUuv1HmV2az5iczLCjnNS5MydS39TBroOtYUcRSbjAioKZlZnZ02a23szWmdmdAxxzmZk1mtmrsdvng8ojyWHd/qNJ23XUY1nFREDjCpKesgI8dxfwj+6+yswKgJVm9ri7r+933PPu/vYAc0iSONjcQc3RdhYkeVGYXZpP4dhsVuw8xLurysKOI5JQgbUU3P2Au6+K3W8CNgDTgno/SX5vDDIn58yjHhkZRtWMiazYqRVTJf0kZEzBzCqAxcArA7x8gZmtMbNHzWxhIvJIOHrm/i+YmtwtBYgONu9oaKGhuSPsKCIJFXhRMLN84FfAP7h7/yuCVgEz3P0c4FvAb05yjjvMrNrMquvrtQRBqtpS28ykglyKxuWEHWVQPYv1rdql1oKkl0CLgpllEy0I/+3uv+7/ursfdffm2P1HgGwzO2HXFXe/192r3L2qtLQ0yMgSoK11TVROzg87RlzOmjaB7ExjpTbdkTQT5OwjA74PbHD3r57kmCmx4zCzZbE8B4PKJOGJRJwtdc1UTkqOTXUGk5edyVnTJqilIGknyNlHFwG3AK+Z2aux5z4HlAO4+3eBG4GPmFkX0Abc5LpiaFTa39hG67HulGkpACwtL+KBl3dxrCtCTpYu6ZH0EFhRcPcXgFNuq+XudwN3B5VBkseWumaAlGkpQHRc4b4XdrBufyOLk3QBP5GRpo8/khBba3uKQuq0FJbEBptXqgtJ0oiKgiTE5tomSvJTY+ZRj8nj85heNIZVGmyWNKKiIAkRHWROnVZCj6Uzili5S4vjSfpQUZDAuTtb65qZm0KDzD2Wziii9mgH+460hR1FJCFUFCRwNUfbae7oYs7k1Blk7tGzQ5zGFSRdqChI4Lak4CBzj/lTChibk6nrFSRtqChI4HY0tAAwq3RcyEmGLiszg0VlhbqyWdKGioIEbkdDC+NyMinNzw07yrAsnVHEhgNNtHR0hR1FJHAqChK4nQdbqCgZR2xFk5SzZEYR3RFnzd4jYUcRCZyKggRuR0O0KKSqJWVaMVXSh4qCBKqzO8Lew23MSuGiMGFsNpWT8jUDSdKCioIEas+hVrojTkVx6hYFiI4rrNp9hEhEF7HJ6KaiIIHaeTA68yiVu48gOq7Q2NbJ9obmsKOIBEpFQQK1vT5aFGameFFYqsXxJE2oKEigdh5sYcKYbIrGZocd5bTMKhlH4dhsFQUZ9VQUJFA7G1pTejpqDzNjaXmRioKMeioKEqgdDS3MLB4bdowRsWRGEdvqWzjccizsKCKBUVGQwLR3drO/sS3lB5l79IwrrN6j1oKMXioKEpjdh1pxT/1B5h7nTC8kM8PUhSSjmoqCBKZnIbzRUhTG5GSy8IzxVO9UUZDRS0VBAtNTFEZL9xFE91dYs/cInd2RsKOIBCKwomBmZWb2tJmtN7N1ZnbnAMeYmX3TzLaa2VozWxJUHkm8nQ0tFI/LYXxeak9H7WvZzIm0d0ZYs0eL48noFGRLoQv4R3dfAJwPfMzMFvQ75lqgMna7A/hOgHkkwXY0tIyarqMeF80uIcPguS0NYUcRCURgRcHdD7j7qtj9JmADMK3fYTcAD3jUy0ChmU0NKpMkVqqvjjqQCWOzOaeskOc214cdRSQQCRlTMLMKYDHwSr+XpgF7+jzey4mFAzO7w8yqzay6vl6/jKmgpaOLuqaOUddSALikspS1e49wpFXXK8joE3hRMLN84FfAP7j70eGcw93vdfcqd68qLS0d2YASiF0HWwFSfnXUgbx5bgkRhxe2qgtJRp9Ai4KZZRMtCP/t7r8e4JB9QFmfx9Njz0mK61kddcYouZq5r3OmF1KQl8Xzm1UUZPQJcvaRAd8HNrj7V09y2HLg1tgspPOBRnc/EFQmSZzRsmT2QLIyM7h4TgnPbq7HXfsryOgSV1Ews1+b2dvMbChF5CLgFuAtZvZq7HadmX3YzD4cO+YRYDuwFfge8NGhhJfktbOhhZL8XPJzs8KOEogrzpxMzdF21uxtDDuKyIiK9zf2/wEfAL5pZr8Efujum071Be7+AnDKpTE9+jHrY3FmkBSy82ArM0tGX9dRj6vOnExWhvHo6wdYVFYYdhyRERPXJ393f8Ld/xZYAuwEnjCzl8zsA7FxA5Hj7GxoYcYoHGTuMWFsNhfMLuax12vUhSSjStzdQWZWDNwGfBBYDXyDaJF4PJBkkrJaj43e6ah9XXf2VHYdbGXDgaawo4iMmHjHFB4GngfGAu9w9+vd/UF3/wSQH2RAST07G6LTUUfjzKO+rl4wmcwMY/ma/WFHERkx8bYUvufuC9z9Sz2zg8wsF8DdqwJLJylpV8/Mo1HcfQRQnJ/L5fNK+dWqvXRpgTwZJeItCv8+wHN/HskgMnrsGMXTUft7T1UZ9U0dPL1JV9rL6HDK2UdmNoXoshNjzGwxb8wmGk+0K0nkBKN9Ompfl8+fRGlBLg+u2MNVCyaHHUfktA32W/tWooPL04G+F6A1AZ8LKJOkuNE+HbWv7MwMblw6nXuf286eQ62UTUyP/7eMXqfsPnL3H7n75cBt7n55n9v1J1m2QmTUT0ft7/0XVJBh8L3nt4cdReS0DdZ99D53/wlQYWaf7v/6KZavkDSVLtNR+5oyIY93Lp7GL6r3cOcVlRTn54YdSWTYBhto7vnNzgcKBriJHCddpqP2d8els+noivCDF3eEHUXktJyypeDu98T+/bfExJFUly7TUfubMymf686ayo9e2sUHL55F0bicsCOJDEu8F6/9h5mNN7NsM3vSzOrN7H1Bh5PUk07TUfv75BWVtBzr4vsvqLUgqSve6xSujm2Q83aiax/NAT4bVChJXbsaWtNmOmp/86YUcN3ZU7n/pZ0cbtGubJKa4i0KPb/hbwN+6e5aL1gGtONgCxVpNp7Q1yffotaCpLZ4i8LvzWwjsBR40sxKgfbgYkmq2nWwJS27jnr0tBZ++OIOtRYkJcW7dPZdwIVAlbt3Ai3ADUEGk9TTeqyL2qMdad1SALjzikpaO7u57wVdtyCpZygdv/OJXq/Q92seGOE8ksK210cHmWeVpvfCuXMnF/C2s6dy/4s7NRNJUk68s49+DPwXcDFwbuym1VHlONvqmwGYneZFAaIzkVo7u3WVs6SceFsKVcAC1xZTcgrb6lvIsPS7cG0gcycXcM3CKfzsL7u588pKcrMyw44kEpd4B5pfB6YEGURS37b6ZsomjiUvW38AAW5aVs7h1k6e3FAXdhSRuMXbUigB1pvZX4COnifd/fpAUklK2lbXrK6jPi6eU8LUCXn8onoP1509New4InGJtyj861BPbGY/IHqxW527nzXA65cBvwV6JnT/2t2/MNT3keQQiTg7Glq4eE5J2FGSRmaG8c7F0/jus9s42NyhhfIkJcQ7JfVZolcyZ8furwBWDfJl9wPXDHLM8+6+KHZTQUhh+4600dEVYfYktRT6uu7sqUQcnthQG3YUkbjEO/voQ8BDwD2xp6YBvznV17j7c8Ch00onKUMzjwa28IzxTCscwx/XqShIaoh3oPljwEXAUQB33wJMGoH3v8DM1pjZo2a28GQHmdkdZlZtZtX19doLNxlti12jMLs0fa9mHoiZ8daFU3hhSwNN7Z1hxxEZVLxFocPde6/Zj13AdrrTU1cBM9z9HOBbnKLl4e73unuVu1eVlpae5ttKELbVN1M4NpuJulDrBFctmMyx7ggvbTsYdhSRQcVbFJ41s88BY8zsKuCXwO9O543d/ai7N8fuPwJkm5lGKVNUz8wjMws7StJZOqOIMdmZvLi1IewoIoOKtyjcBdQDrwF/DzwC/K/TeWMzm2KxvyBmtiyWRR+lUtT2hhZ1HZ1ETlYG582ayAtbVBQk+cU1JdXdI2b2G+A37h5Xp76Z/Qy4DCgxs73AvwDZsfN9F7gR+IiZdQFtwE26Yjo1NbZ1Ut/UkfZrHp3KxXNK+PdNG9h3pI1phWPCjiNyUqcsCrFP8v8CfJxYq8LMuoFvDTaF1N1vHuT1u4G7h5RWktJ2zTwa1MWV0Z7RF7c28J6qspDTiJzcYN1HnyI66+hcd5/o7hOB84CLzOxTgaeTlKCZR4ObN7mAkvxcdSFJ0husKNwC3OzuvdtIuft24H3ArUEGk9Sxta6Z7EyjbKIWwjsZM+P8WRNZsVOX7khyG6woZLv7CR9tYuMK2cFEklSzpbaJ2aX5ZGfGO28hPVXNKOJAYzv7jrSFHUXkpAb7LT7VfoLaa1AA2FTbROXkgrBjJL2qiokAVKu1IElssKJwjpkdHeDWBJydiICS3Jo7uth7uI15kzXIPJj5UwoYm5NJ9c7DYUcROalTzj5ydy2ML6e0pbYJiG4qI6eWlZnB4vJCqnepKEjyUiewnJbNsaIwb4qKQjyWzpjIppqjWgdJkpaKgpyWzbXN5GVnUFakmUfxOLeiiIjD6t1Hwo4iMiAVBTktm2ubqJxUQEaG1jyKx+LyIjIMdSFJ0lJRkNOyqaZJ4wlDkJ+bxbwp41mloiBJSkVBhu1I6zHqmjqYN0Uzj4ZiSXkhr+45QndES31J8lFRkGHbXBtd80gthaFZUl5Ec0cXW+qawo4icgIVBRm2TZp5NCxLZhQBsGqXBpsl+agoyLBtrmmiIDeLKePzwo6SUiqKxzJxXA6rdmtcQZKPioIM26baJuZOKdBua0NkZiwpL9RgsyQlFQUZFndnc61mHg3X4vIitje0cLhFS4hJclFRkGGpb+7gSGun1jwapiXl0XGF1XvUWpDkoqIgw7K5JjbzSIPMw3JO2QQyM0yDzZJ0VBRkWDZpIbzTMjYnizOnFmiwWZKOioIMy6aaoxSPy6EkPzfsKClrSXkRa/Ycoas7EnYUkV4qCjIsm2qamD9VrYTTsaS8iJZj3b2tLpFkEFhRMLMfmFmdmb1+ktfNzL5pZlvNbK2ZLQkqi4ys7oizqbaJ+VPGhx0lpfUMNq/SiqmSRIJsKdwPXHOK168FKmO3O4DvBJhFRtCugy20d0aYr0Hm01I2cQwl+bms1vUKkkQCKwru/hxwqs1obwAe8KiXgUIzmxpUHhk5G2ui3R1nTlVL4XT0XsSmwWZJImGOKUwD9vR5vDf23AnM7A4zqzaz6vr6+oSEk5PbeOAoGQZzJukahdO1ZEYROw+2crC5I+woIkCKDDS7+73uXuXuVaWlpWHHSXsba5qYWTKOvGxt4X26ei9i07iCJIkwi8I+oKzP4+mx5yTJbaxpYr66jkbEm6ZPICvDWKkuJEkSYRaF5cCtsVlI5wON7n4gxDwSh+aOLnYfauVMDTKPiLzsTBaeoZ3YJHlkBXViM/sZcBlQYmZ7gX8BsgHc/bvAI8B1wFagFfhAUFlk5Gyq6dlDQS2FkbK4vIgHV+yhqztCVmZK9OjKKBZYUXD3mwd53YGPBfX+EoyeoqDpqCOnqqKI+1/ayWv7GlkcG2MQCYs+lsiQbKw5Sn5uFtOLxoQdZdS4YFYxAC9tOxhyEhEVBRmijQeamK+NdUZUcX4uZ04dz4tbG8KOIqKiIPFzdzbUHNWezAG4aHYx1bsO097ZHXYUSXMqChK3/Y3tNLV3aTpqAC6aU8KxrgjVOzULScKloiBx21RzFEDTUQOwbOZEsjKMF9SFJCFTUZC4bTgQ21hHRWHEjcvNoqqiiKc31oUdRdKcioLEbWNNE9OLxjA+LzvsKKPSlWdOZlNtE7sPtoYdRdKYioLEbeOBo7o+IUBXL5gCwJ/W14ScRNKZioLEpaOrm+0NLdpYJ0DlxWOZN7mAx9fXhh1F0piKgsRla10z3RHXFpwBu2rBZFbsPMShlmNhR5E0paIgcdl4oGd5C7UUgnTt2VOIOPxh7f6wo0iaUlGQuGysOUpOVgYVxWPDjjKqLTxjAmdOHc9DK/eGHUXSlIqCxGVjTRNzJ+drFc8EuHHpdNbsbWRLbVPYUSQN6Tdc4rLhQJO6jhLkhkVnkJVhai1IKFQUZFANzR00NHdoOmqClOTncuWZk/n5ij20HusKO46kGRUFGdQbeyiopZAoH7p0Jo1tnfyyWq0FSSwVBRnUhgPRNY80HTVxls6YyJLyQr73/HaOdUXCjiNpREVBBrVu/1EmFeRSkp8bdpS08okrKtl7uI2fvrIr7CiSRlQUZFDr9jey8Ax1HSXaZXNLuXB2Md98aitH2zvDjiNpQkVBTqm9s5tt9S0sPGNC2FHSjpnxP689k0Mtx/j641vCjiNpItCiYGbXmNkmM9tqZncN8PptZlZvZq/Gbh8MMo8M3aaaJrojrpZCSM6ePoFbzp/BD1/awcpd2oBHghdYUTCzTODbwLXAAuBmM1swwKEPuvui2O2+oPLI8KzbHx1kVkshPP907XzOmDCG//HQGm3XKYELsqWwDNjq7tvd/Rjwc+CGAN9PArBufyMFeVmUTRwTdpS0lZ+bxZf++my21bfw5Uc3hh1HRrkgi8I0YE+fx3tjz/X3LjNba2YPmVlZgHlkGNbtP8qCqeMxs7CjpLVL55bygYsquP+lnTyhpbUlQGEPNP8OqHD3NwGPAz8a6CAzu8PMqs2sur6+PqEB01l3xNlYc1RdR0nirmvns2DqeD770BpqGtvDjiOjVJBFYR/Q95P/9Nhzvdz9oLt3xB7eBywd6ETufq+7V7l7VWlpaSBh5UTb65tp74xokDlJ5GZl8q2/WUxHV4Q7f76a7oiHHUlGoSCLwgqg0sxmmlkOcBOwvO8BZja1z8PrgQ0B5pEh6h1knqaikCxml+bzb9cv5JUdh/j201vDjiOjUGBFwd27gI8DfyT6x/4X7r7OzL5gZtfHDvukma0zszXAJ4HbgsojQ7dufyM5WRnMLs0PO4r0cePS6dyw6Ay+/sRmVuw8FHYcGWXMPbWaoFVVVV5dXR12jLTw7u++RFfEefijF4UdRfppau/k7d96gc6uCI/ceQmFY3PCjiRJzsxWunvVYMeFPdAsSaqzO8Jr+xpZVFYYdhQZQEFeNt+6eTH1zR185pdrSbUPd5K8VBRkQJtqmmjvjLC4vCjsKHISb5peyD9dM58nNtRy/0s7w44jo4SKggzo1T1HAFislkJSu/3imVwxfxJfemQjr+9rDDuOjAIqCjKgV/ccoXhcDtOLdCVzMjMz/uvd51Ccn8PHf7qK5g7t1CanR0VBBrR692EWlRXqSuYUUDQuh2/ctJjdh1r554df0/iCnBYVBTlBY1sn2+pbNMicQpbNnMinrpzLb1/dry085bSoKMgJVuyIzn0/d+bEkJPIUHz08jlcNKeYzy9/nbV7j4QdR1KUioKc4OXtB8nJylBLIcVkZhjfuGkxxeNy+dAD1VofSYZFRUFO8PKOgywpLyQvOzPsKDJEJfm5fP+2Kprbu/jQA9W0HdP+CzI0KgpynMa2TtbtP8r5s4rDjiLDNH/KeL5582Je39/Ih3+yko4uFQaJn4qCHOfl7QdxR0UhxV1x5mS+9M6zeXZzPR/771Uc64qEHUlShIqCHOfpjXUU5GaxRFcyp7yblpXzhRsW8sSGOm79wSs0tnaGHUlSgIqC9IpEnCc31nHpvFJysvSjMRrcekEFX3vvOazadYR3fudFNhw4GnYkSXL6zZder+1rpL6pgyvPnBR2FBlB71w8nR/fvoym9i5uuPtF7nl2G53d6k6SgakoSK/H19eSYXDZXBWF0ea8WcU8duclXDavlC89upG3ffN5XtzaEHYsSUIqCgJEu45+8+o+Lq4spWic1uYfjYrzc7nnlqXce8tS2jq7+dv7XuGOB6rZWtcUdjRJIioKAsCKnYfYe7iNv148LewoEiAz4+qFU3j8U2/mM1fP5aVtB7n6a8/xTw+t5UBjW9jxJAmoKAgAD/x5FwV5WVy9cHLYUSQB8rIz+fhbKnn2s5dx24UzeXj1Pi77z2f40iMbONJ6LOx4EiIVBWHPoVYeff0Af3NeOWNzssKOIwlUnJ/L59+xgKc+82be9qap3Pv8di79j6f5zz9uVMshTakoCF99fDNZmRncdmFF2FEkJNOLxvLV9yzi0Tsv4YLZxXznmW1c/JWnueOBapav2a99GtKIPhamub/sOMTDq/fxkctmM3WCNtRJd/OnjOeeW6rYc6iVn7y8i1+v3sef1teSk5XBhbOLubSylEvnljC7NF91Sl2mAAAIHUlEQVR7bYxSFuSGHGZ2DfANIBO4z92/3O/1XOABYClwEHivu+881Tmrqqq8uro6mMBppq6pnb+6+0WyszL4wycvIT9XnxHkeN0RZ9Xuwzzy2gGe3VTP9oYWAKZOyOOSyhIurizlotnFFOfnhpxUBmNmK929arDjAvsrYGaZwLeBq4C9wAozW+7u6/scdjtw2N3nmNlNwFeA9waVSd6wta6JO368kiNtnTx4xwUqCDKgzAzj3IqJnFsxEd4RHX96YWsDz2+p57HXa/hFbEOfiuKxLCor5E3TC5lZMo7y4rFMLxpDbpZW2k01gbUUzOwC4F/d/a2xx/8TwN2/1OeYP8aO+bOZZQE1QKmfIpRaCsMTiTgHW47x2r4jPPZ6DQ+v3kdBXjb33LI0+gsvMkTdEee1fY28tK2BV3cfYfWeI9Q3dRx3zNicTCaMyWbCmGzG52VTkJcVu2X3+zer97jCsTmx47PIytSw50gJvaUATAP29Hm8FzjvZMe4e5eZNQLFwIhfavns5nq++PtoI6Wn5vRWnj4lqOdu37r0xnN9j/MTn+tXyoZ1jgHPdfxxfV864f/S97g+J2nvjHAstrRBblYGN51bzifeModJ4/MQGY7MDGNRWWHvZkzuTn1zB3sOtbLrYCv7DrfR2NZJY1snR9o6aWrvpOZoO1vqumhq76SpvYuuyKk/lObnZpGblUFmhpGVYWRmGlkZGZjBSI5ojOT4SJAjLe89t4wPXjIrwHdIkYFmM7sDuAOgvLx8WOfIz81i3uSCPic97p/jfijeeO6Ewwc8juOOs+O+tu8PyBvP2YnPHfeTZCc81z/TgOc47r2O/9HMzc7gjAljqJyUz5IZRdpAR0acmTGpII9JBXksnTF469Pdae+M0NTeydH2ThrbumhsOxYtJK3RQtLY1smxrggRd7q6ne6I0xVxukeyh2NETxXcGC1EN1EKWpBFYR9Q1ufx9NhzAx2zN9Z9NIHogPNx3P1e4F6Idh8NJ8zSGUUsnaHloEWShZkxJieTMTmZarEmkSA77FYAlWY208xygJuA5f2OWQ68P3b/RuCpU40niIhIsAJrKcTGCD4O/JHolNQfuPs6M/sCUO3uy4HvAz82s63AIaKFQ0REQhLomIK7PwI80u+5z/e53w68O8gMIiISP833EhGRXioKIiLSS0VBRER6qSiIiEgvFQUREekV6CqpQTCzemBX2DniUEIAy3UESHmDl2qZlTdYic47w91LBzso5YpCqjCz6ngWn0oWyhu8VMusvMFK1rzqPhIRkV4qCiIi0ktFITj3hh1giJQ3eKmWWXmDlZR5NaYgIiK91FIQEZFeKgojxMwmmtnjZrYl9u8JmzeY2SIz+7OZrTOztWaW8P2ozewaM9tkZlvN7K4BXs81swdjr79iZhWJztgvz2B5P21m62PfzyfNbEYYOfvkOWXePse9y8zczEKdfRJPXjN7T+x7vM7MfprojAPkGexnotzMnjaz1bGfi+vCyBnL8gMzqzOz10/yupnZN2P/l7VmtiTRGU/g7rqNwA34D+Cu2P27gK8McMxcoDJ2/wzgAFCYwIyZwDZgFpADrAEW9Dvmo8B3Y/dvAh4M8XsaT97LgbGx+x9J9ryx4wqA54CXgapkzgtUAquBotjjSWHlHULme4GPxO4vAHaGmPdSYAnw+klevw54lOjGiecDr4T5/XV3tRRG0A3Aj2L3fwT8Vf8D3H2zu2+J3d8P1AGDXkwygpYBW919u7sfA35ONHdfff8fDwFX2EhuYDs0g+Z196fdvTX28GWiO/yFJZ7vL8AXga8A7YkMN4B48n4I+La7HwZw97oEZ+wvnswOjI/dnwDsT2C+44O4P0d0r5iTuQF4wKNeBgrNbGpi0g1MRWHkTHb3A7H7NcDkUx1sZsuIftLZFnSwPqYBe/o83ht7bsBj3L0LaASKE5LuRPHk7et2op+6wjJo3lj3QJm7/yGRwU4inu/vXGCumb1oZi+b2TUJSzeweDL/K/A+M9tLdD+XTyQm2rAM9Wc8cIFusjPamNkTwJQBXvrnvg/c3c3spNO6Yp8Efgy8390jI5syPZnZ+4Aq4M1hZzkZM8sAvgrcFnKUocgi2oV0GdFW2HNmdra7Hwk11andDNzv7v/XzC4gurvjWfpdi4+KwhC4+5Une83Mas1sqrsfiP3RH7CZbWbjgT8A/xxrLibSPqCsz+PpsecGOmavmWURbX4fTEy8E8STFzO7kmhhfrO7dyQo20AGy1sAnAU8E+uRmwIsN7Pr3b06YSnfEM/3dy/Rfu5OYIeZbSZaJFYkJuIJ4sl8O3ANgLv/2czyiK4zFHbX10Di+hlPJHUfjZzlwPtj998P/Lb/AWaWAzxMtA/xoQRm67ECqDSzmbEsNxHN3Vff/8eNwFMeGxELwaB5zWwxcA9wfRL0d58yr7s3unuJu1e4ewXRMZCwCgLE9/PwG6KtBMyshGh30vZEhuwnnsy7gSsAzOxMIA+oT2jK+C0Hbo3NQjofaOzTDR2OsEe6R8uNaL/7k8AW4AlgYuz5KuC+2P33AZ3Aq31uixKc8zpgM9GxjH+OPfcFon+cIPoL9EtgK/AXYFbI39fB8j4B1Pb5fi5P5rz9jn2GEGcfxfn9NaJdXuuB14CbwswbZ+YFwItEZya9ClwdYtafEZ1l2Em01XU78GHgw32+v9+O/V9eC/vnwd11RbOIiLxB3UciItJLRUFERHqpKIiISC8VBRER6aWiICIivVQURESkl4qCiIj0UlEQEZFe/x82FaDPzi84FgAAAABJRU5ErkJggg==\n",
       "text/plain": [
        "<Figure size 432x288 with 1 Axes>"
       ]
@@ -391,8 +385,14 @@
     }
    ],
    "source": [
-    "samp = a_sample.extract()\n",
-    "samp.shape, samp.columns, samp.theta.plot.density()"
+    "samp.theta.plot.density()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `extract` function argument `params` takes a list of parameter or column names:"
    ]
   },
   {
@@ -418,10 +418,182 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 13,
    "metadata": {},
-   "outputs": [],
-   "source": []
+   "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.343484</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>0.232548</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>0.342342</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "      theta\n",
+       "0  0.343484\n",
+       "1  0.232548\n",
+       "2  0.342342"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "thetas[0:3]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### PosteriorSample.sample\n",
+    "\n",
+    "The `sample` property 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 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": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(1000, 4, 8)"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a_sample.sample\n",
+    "a_sample.sample.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "('theta', array([0.343484, 0.232548, 0.342342]))"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a_sample.column_names[7], a_sample.sample[0:3,0,7]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### PosteriorSample.stepsize\n",
+    "\n",
+    "The `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": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([0.850374, 0.789576, 0.826539, 1.10651 ])"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a_sample.stepsize"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### PosteriorSample.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": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.507758],\n",
+       "       [0.42248 ],\n",
+       "       [0.528419],\n",
+       "       [0.422747]])"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a_sample.metric"
+   ]
   }
  ],
  "metadata": {


Test patch:
diff --git a/test/data/output_bad_metric_1.csv b/test/data/output_bad_metric_1.csv
new file mode 100644
index 0000000..833101b
--- /dev/null
+++ b/test/data/output_bad_metric_1.csv
@@ -0,0 +1,53 @@
+# stan_version_major = 2
+# stan_version_minor = 19
+# stan_version_patch = 0
+# model = bernoulli_model
+# method = sample (Default)
+#   sample
+#     num_samples = 10
+#     num_warmup = 100
+#     save_warmup = 0 (Default)
+#     thin = 1 (Default)
+#     adapt
+#       engaged = 1 (Default)
+#       gamma = 0.050000000000000003 (Default)
+#       delta = 0.80000000000000004 (Default)
+#       kappa = 0.75 (Default)
+#       t0 = 10 (Default)
+#       init_buffer = 75 (Default)
+#       term_buffer = 50 (Default)
+#       window = 25 (Default)
+#     algorithm = hmc (Default)
+#       hmc
+#         engine = nuts (Default)
+#           nuts
+#             max_depth = 10 (Default)
+#         metric = diag_e
+#         metric_file =  (Default)
+#         stepsize = 1 (Default)
+#         stepsize_jitter = 0 (Default)
+# id = 1
+# data
+#   file = examples/bernoulli/bernoulli.data.json
+# init = 2 (Default)
+# random
+#   seed = 123456
+# output
+#   file = output_1.csv
+#   diagnostic_file =  (Default)
+#   refresh = 100 (Default)
+lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
+-6.76206,1,0.787025,1,1,0,6.81411,0.229458
+-6.81411,0.983499,0.787025,1,1,0,6.8147,0.20649
+-6.85511,0.994945,0.787025,2,3,0,6.85536,0.310589
+-6.85511,0.812189,0.787025,1,1,0,7.16517,0.310589
+-6.85511,0.310983,0.787025,1,1,0,8.35691,0.310589
+-10.0407,0.245981,0.787025,1,1,0,10.0474,0.614551
+-6.78719,1,0.787025,1,3,0,8.92907,0.21615
+-7.58503,0.845748,0.787025,1,3,0,7.63531,0.115185
+-8.08941,0.916692,0.787025,1,1,0,8.16887,0.0892886
+# 
+#  Elapsed Time: 0.001332 seconds (Warm-up)
+#                0.000249 seconds (Sampling)
+#                0.001581 seconds (Total)
+# 
diff --git a/test/data/output_bad_metric_2.csv b/test/data/output_bad_metric_2.csv
new file mode 100644
index 0000000..be05513
--- /dev/null
+++ b/test/data/output_bad_metric_2.csv
@@ -0,0 +1,57 @@
+# stan_version_major = 2
+# stan_version_minor = 19
+# stan_version_patch = 0
+# model = bernoulli_model
+# method = sample (Default)
+#   sample
+#     num_samples = 10
+#     num_warmup = 100
+#     save_warmup = 0 (Default)
+#     thin = 1 (Default)
+#     adapt
+#       engaged = 1 (Default)
+#       gamma = 0.050000000000000003 (Default)
+#       delta = 0.80000000000000004 (Default)
+#       kappa = 0.75 (Default)
+#       t0 = 10 (Default)
+#       init_buffer = 75 (Default)
+#       term_buffer = 50 (Default)
+#       window = 25 (Default)
+#     algorithm = hmc (Default)
+#       hmc
+#         engine = nuts (Default)
+#           nuts
+#             max_depth = 10 (Default)
+#         metric = diag_e
+#         metric_file =  (Default)
+#         stepsize = 1 (Default)
+#         stepsize_jitter = 0 (Default)
+# id = 1
+# data
+#   file = examples/bernoulli/bernoulli.data.json
+# init = 2 (Default)
+# random
+#   seed = 123456
+# output
+#   file = output_1.csv
+#   diagnostic_file =  (Default)
+#   refresh = 100 (Default)
+lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
+# Adaptation terminated
+# Step size = not a number
+# Diagonal elements of inverse mass matrix:
+# 1
+-6.76206,1,0.787025,1,1,0,6.81411,0.229458
+-6.81411,0.983499,0.787025,1,1,0,6.8147,0.20649
+-6.85511,0.994945,0.787025,2,3,0,6.85536,0.310589
+-6.85511,0.812189,0.787025,1,1,0,7.16517,0.310589
+-6.85511,0.310983,0.787025,1,1,0,8.35691,0.310589
+-10.0407,0.245981,0.787025,1,1,0,10.0474,0.614551
+-6.78719,1,0.787025,1,3,0,8.92907,0.21615
+-7.58503,0.845748,0.787025,1,3,0,7.63531,0.115185
+-8.08941,0.916692,0.787025,1,1,0,8.16887,0.0892886
+# 
+#  Elapsed Time: 0.001332 seconds (Warm-up)
+#                0.000249 seconds (Sampling)
+#                0.001581 seconds (Total)
+# 
diff --git a/test/data/output_bad_metric_3.csv b/test/data/output_bad_metric_3.csv
new file mode 100644
index 0000000..de2477d
--- /dev/null
+++ b/test/data/output_bad_metric_3.csv
@@ -0,0 +1,59 @@
+# stan_version_major = 2
+# stan_version_minor = 19
+# stan_version_patch = 0
+# model = bernoulli_model
+# method = sample (Default)
+#   sample
+#     num_samples = 10
+#     num_warmup = 100
+#     save_warmup = 0 (Default)
+#     thin = 1 (Default)
+#     adapt
+#       engaged = 1 (Default)
+#       gamma = 0.050000000000000003 (Default)
+#       delta = 0.80000000000000004 (Default)
+#       kappa = 0.75 (Default)
+#       t0 = 10 (Default)
+#       init_buffer = 75 (Default)
+#       term_buffer = 50 (Default)
+#       window = 25 (Default)
+#     algorithm = hmc (Default)
+#       hmc
+#         engine = nuts (Default)
+#           nuts
+#             max_depth = 10 (Default)
+#         metric = diag_e
+#         metric_file =  (Default)
+#         stepsize = 1 (Default)
+#         stepsize_jitter = 0 (Default)
+# id = 1
+# data
+#   file = examples/bernoulli/bernoulli.data.json
+# init = 2 (Default)
+# random
+#   seed = 123456
+# output
+#   file = output_1.csv
+#   diagnostic_file =  (Default)
+#   refresh = 100 (Default)
+lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
+# Adaptation terminated
+# Step size = 0.787025
+# Elements of inverse mass matrix:
+# 1, 0, 0
+# 0, 2, 0
+# 1, 0, 1
+-6.76206,1,0.787025,1,1,0,6.81411,0.229458
+-6.81411,0.983499,0.787025,1,1,0,6.8147,0.20649
+-6.85511,0.994945,0.787025,2,3,0,6.85536,0.310589
+-6.85511,0.812189,0.787025,1,1,0,7.16517,0.310589
+-6.85511,0.310983,0.787025,1,1,0,8.35691,0.310589
+-10.0407,0.245981,0.787025,1,1,0,10.0474,0.614551
+-6.78719,1,0.787025,1,3,0,8.92907,0.21615
+-7.58503,0.845748,0.787025,1,3,0,7.63531,0.115185
+-8.08941,0.916692,0.787025,1,1,0,8.16887,0.0892886
+# 
+#  Elapsed Time: 0.001332 seconds (Warm-up)
+#                0.000249 seconds (Sampling)
+#                0.001581 seconds (Total)
+# 
diff --git a/test/data/output_bad_metric_4.csv b/test/data/output_bad_metric_4.csv
new file mode 100644
index 0000000..9f07e4c
--- /dev/null
+++ b/test/data/output_bad_metric_4.csv
@@ -0,0 +1,59 @@
+# stan_version_major = 2
+# stan_version_minor = 19
+# stan_version_patch = 0
+# model = bernoulli_model
+# method = sample (Default)
+#   sample
+#     num_samples = 10
+#     num_warmup = 100
+#     save_warmup = 0 (Default)
+#     thin = 1 (Default)
+#     adapt
+#       engaged = 1 (Default)
+#       gamma = 0.050000000000000003 (Default)
+#       delta = 0.80000000000000004 (Default)
+#       kappa = 0.75 (Default)
+#       t0 = 10 (Default)
+#       init_buffer = 75 (Default)
+#       term_buffer = 50 (Default)
+#       window = 25 (Default)
+#     algorithm = hmc (Default)
+#       hmc
+#         engine = nuts (Default)
+#           nuts
+#             max_depth = 10 (Default)
+#         metric = dense_e
+#         metric_file =  (Default)
+#         stepsize = 1 (Default)
+#         stepsize_jitter = 0 (Default)
+# id = 1
+# data
+#   file = examples/bernoulli/bernoulli.data.json
+# init = 2 (Default)
+# random
+#   seed = 123456
+# output
+#   file = output_1.csv
+#   diagnostic_file =  (Default)
+#   refresh = 100 (Default)
+lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,theta
+# Adaptation terminated
+# Step size = 0.787025
+# Elements of inverse mass matrix:
+# 1, 0, 0
+# 0, 2, 0, 0
+# 1, 0, 1
+-6.76206,1,0.787025,1,1,0,6.81411,0.229458
+-6.81411,0.983499,0.787025,1,1,0,6.8147,0.20649
+-6.85511,0.994945,0.787025,2,3,0,6.85536,0.310589
+-6.85511,0.812189,0.787025,1,1,0,7.16517,0.310589
+-6.85511,0.310983,0.787025,1,1,0,8.35691,0.310589
+-10.0407,0.245981,0.787025,1,1,0,10.0474,0.614551
+-6.78719,1,0.787025,1,3,0,8.92907,0.21615
+-7.58503,0.845748,0.787025,1,3,0,7.63531,0.115185
+-8.08941,0.916692,0.787025,1,1,0,8.16887,0.0892886
+# 
+#  Elapsed Time: 0.001332 seconds (Warm-up)
+#                0.000249 seconds (Sampling)
+#                0.001581 seconds (Total)
+# 
diff --git a/test/data/runset-bad/bad-cols-bern-4.csv b/test/data/runset-bad/bad-cols-bern-4.csv
index 3812d91..fa5a256 100644
--- a/test/data/runset-bad/bad-cols-bern-4.csv
+++ b/test/data/runset-bad/bad-cols-bern-4.csv
@@ -40,7 +40,7 @@ lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,thet
 # Adaptation terminated
 # Step size = 0.730383
 # Diagonal elements of inverse mass matrix:
-# 0.508248
+# 0.508248, 0.508248
 -7.00518,1,0.730383,2,3,0,7.26237,0.345882
 -6.76715,0.968283,0.730383,1,3,0,7.334,0.226096
 -7.44199,0.932693,0.730383,2,3,0,7.44227,0.124863
diff --git a/test/test_posterior_sample.py b/test/test_posterior_sample.py
index 998a0d0..ec71db1 100644
--- a/test/test_posterior_sample.py
+++ b/test/test_posterior_sample.py
@@ -12,32 +12,8 @@ goodfiles_path = os.path.join(datafiles_path, 'runset-good')
 badfiles_path = os.path.join(datafiles_path, 'runset-bad')
 
 
-class PostSampleTest(unittest.TestCase):
-    def test_postsample_prefab_csv(self):
-        column_names = ['lp__','accept_stat__','stepsize__','treedepth__',
-                            'n_leapfrog__','divergent__','energy__', 'theta']
-        output = os.path.join(datafiles_path, 'runset-good')
-        csv_files = []
-        for i in range(4):
-            csv_files.append(os.path.join(output, 'bern-{}.csv'.format(i+1)))
-        run = {'draws': 100, 'chains': 4, 'column_names': tuple(column_names)}
-        post_sample = PosteriorSample(run, tuple(csv_files))
-        self.assertEqual(post_sample.chains,4)
-        self.assertEqual(post_sample.draws,100)
-        self.assertEqual(post_sample.column_names, tuple(column_names))
-        post_sample.sample
-        self.assertEqual(post_sample.sample.shape, (100, 4, 8))
-        df = post_sample.extract()
-        self.assertEqual(df.shape, (400, 8))
-
-    def test_postsample_bad_csv(self):
-        csv_files = []
-        csv_files = [os.path.join(datafiles_path, 'no-such-sample.csv')]
-        run = {'draws': 1, 'chains': 1, 'column_names': []}
-        with self.assertRaises(Exception):
-            PosteriorSample(run, tuple(csv_files))
-
-    def test_postsample_good(self):
+class PosteriorSampleTest(unittest.TestCase):
+    def test_sample_good(self):
         column_names = ['lp__','accept_stat__','stepsize__','treedepth__',
                             'n_leapfrog__','divergent__','energy__', 'theta']
         stan = os.path.join(datafiles_path, 'bernoulli.stan')
@@ -49,8 +25,12 @@ class PostSampleTest(unittest.TestCase):
         post_sample = sample(model, data_file=jdata)
         self.assertEqual(post_sample.chains,4)
         self.assertEqual(post_sample.draws,1000)
+        self.assertEqual(post_sample.columns, len(column_names))
         self.assertEqual(post_sample.column_names, tuple(column_names))
         post_sample.sample
+        self.assertEqual(post_sample.metric_type, 'diag_e')
+        self.assertEqual(post_sample.stepsize.shape,(4,))
+        self.assertEqual(post_sample.metric.shape,(4, 1))
         self.assertEqual(post_sample.sample.shape,(1000, 4, 8))
         df = post_sample.summary()
         self.assertTrue(df.shape == (2, 9))
@@ -60,9 +40,38 @@ class PostSampleTest(unittest.TestCase):
         sys.stdout = sys.__stdout__
         self.assertEqual(capturedOutput.getvalue(), 'No problems detected.\n')
 
-    def test_postsample_diagnose_divergences(self):
-        output = os.path.join(datafiles_path, 'diagnose-good',
-                                  'corr_gauss_depth8-1.csv')
+    def test_sample_bad_csv(self):
+        csv_files = []
+        csv_files = [os.path.join(datafiles_path, 'no-such-sample.csv')]
+        run = {'draws': 1, 'chains': 1, 'column_names': []}
+        with self.assertRaises(Exception):
+            PosteriorSample(run, tuple(csv_files))
+
+
+class SampleTest(unittest.TestCase):
+    def test_sample_prefab_runset(self):
+        column_names = ['lp__','accept_stat__','stepsize__','treedepth__',
+                            'n_leapfrog__','divergent__','energy__', 'theta']
+        num_params = 1
+        metric = 'diag_e'
+        run = {'draws': 100, 'chains': 4, 'column_names': tuple(column_names),
+                   'num_params' : num_params, 'metric' : metric}
+
+        output = os.path.join(datafiles_path, 'runset-good')
+        csv_files = []
+        for i in range(4):
+            csv_files.append(os.path.join(output, 'bern-{}.csv'.format(i+1)))
+
+        post_sample = PosteriorSample(run, tuple(csv_files))
+        self.assertEqual(post_sample.chains,4)
+        self.assertEqual(post_sample.draws,100)
+        self.assertEqual(post_sample.column_names, tuple(column_names))
+        post_sample.sample
+        self.assertEqual(post_sample.sample.shape, (100, 4, 8))
+        df = post_sample.extract()
+        self.assertEqual(df.shape, (400, 8))
+
+    def test_sample_diagnose_divergences(self):
         # TODO - use cmdstan test files instead
         expected = ''.join([
             '424 of 1000 (42%) transitions hit the maximum ',
@@ -71,8 +80,12 @@ class PostSampleTest(unittest.TestCase):
             'due to this limit will result in slow ',
             'exploration and you should increase the ',
             'limit to ensure optimal performance.\n'])
-        run = { 'draws': 10, 'chains': 1, 'column_names': ['a', 'b', 'c']}
+        run = { 'draws': 10, 'chains': 1, 'column_names': ['a', 'b', 'c'],
+                    'num_params' : 3, 'metric' : 'diag_e'}
+        output = os.path.join(datafiles_path, 'diagnose-good',
+                                  'corr_gauss_depth8-1.csv')
         csv_files = tuple([output])
+
         post_sample = PosteriorSample(run, csv_files)
         capturedOutput = io.StringIO()
         sys.stdout = capturedOutput
@@ -80,24 +93,36 @@ class PostSampleTest(unittest.TestCase):
         sys.stdout = sys.__stdout__
         self.assertEqual(capturedOutput.getvalue(), expected)
 
-    def test_postsample_big(self):
+    def test_sample_big(self):
         sampler_state = ['lp__','accept_stat__','stepsize__','treedepth__',
                              'n_leapfrog__','divergent__','energy__']
         phis = ['phi.{}'.format(str(x+1)) for x in range(2095)]
         column_names = sampler_state + phis
+        metric = 'diag_e'
+        run = { 'draws': 1000, 'chains': 2, 'column_names': column_names,
+                    'num_params' : len(phis), 'metric' : metric}
+
         output = os.path.join(datafiles_path, 'runset-big')
         csv_files = []
         for i in range(2):
             csv_files.append(os.path.join(
                 output, 'output_icar_nyc_{}.csv'.format(i+1)))
-        run = { 'draws': 1000, 'chains': 2, 'column_names': column_names}
+
         post_sample = PosteriorSample(run, tuple(csv_files))
+        self.assertEqual(post_sample.chains, 2)
+        self.assertEqual(post_sample.columns, len(column_names))
+        self.assertEqual(post_sample.column_names, column_names)
         post_sample.sample
+        self.assertEqual(post_sample.metric_type, 'diag_e')
+        self.assertEqual(post_sample.stepsize.shape,(2,))
+        self.assertEqual(post_sample.metric.shape,(2, 2095))
         self.assertEqual((1000,2,2102), post_sample.sample.shape)
         phis = post_sample.extract(params=['phi'])
         self.assertEqual((2000,2095), phis.shape)
         phi1 = post_sample.extract(params=['phi.1'])
         self.assertEqual((2000,1), phi1.shape)
+        mo_phis = post_sample.extract(params=['phi.1', 'phi.10', 'phi.100'])
+        self.assertEqual((2000,3), mo_phis.shape)
         phi2095 = post_sample.extract(params=['phi.2095'])
         self.assertEqual((2000,1), phi2095.shape)
         with self.assertRaises(Exception):
diff --git a/test/test_runset.py b/test/test_runset.py
index c7239b5..6aca584 100644
--- a/test/test_runset.py
+++ b/test/test_runset.py
@@ -123,7 +123,7 @@ class RunSetTest(unittest.TestCase):
         with self.assertRaisesRegex(ValueError, 'draws'):
             runset.validate_csv_files()
 
-    def test_validate_bad_cols(self):
+    def test_validate_extra_param(self):
         # construct runset using existing sampler output
         stan = os.path.join(datafiles_path, 'bernoulli.stan')
         exe = os.path.join(datafiles_path, 'bernoulli')
diff --git a/test/test_utils.py b/test/test_utils.py
index 21c3a09..c142dcc 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -87,6 +87,29 @@ class ReadStanCsvTest(unittest.TestCase):
         with self.assertRaisesRegex(Exception, 'found 9'):
             dict = check_csv(csv_bad)
 
+    def test_check_csv_metric_1(self):
+        csv_bad = os.path.join(datafiles_path, 'output_bad_metric_1.csv')
+        with self.assertRaisesRegex(Exception,
+                                   'expecting metric'):
+            dict = check_csv(csv_bad)
+
+    def test_check_csv_metric_2(self):
+        csv_bad = os.path.join(datafiles_path, 'output_bad_metric_2.csv')
+        with self.assertRaisesRegex(Exception, 'invalid stepsize'):
+            dict = check_csv(csv_bad)
+
+    def test_check_csv_metric_3(self):
+        csv_bad = os.path.join(datafiles_path, 'output_bad_metric_3.csv')
+        with self.assertRaisesRegex(Exception,
+                                    'invalid or missing mass matrix specification'):
+            dict = check_csv(csv_bad)
+
+    def test_check_csv_metric_3(self):
+        csv_bad = os.path.join(datafiles_path, 'output_bad_metric_4.csv')
+        with self.assertRaisesRegex(Exception,
+                                    'invalid or missing mass matrix specification'):
+            dict = check_csv(csv_bad)
+
 
 if __name__ == '__main__':
     unittest.main()