Part_1_Validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alignment of SMP and Snow Pits / Evaluation of Proksch (P15) coefficients\n",
    "*Josh King and Mike Brady, Environment and Climate Change Canada, 2020*\n",
    "\n",
    "This workbook introduces a snow on sea ice calibration procedure for SMP-derived estimates of density first introduced in [Proksch, et al., 2015](https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2014JF003266). Where indicated, the work modifies portions of the [SMP python package from SLF](https://github.com/slf-dot-ch/snowmicropyn) and uses a number of open source community packages to facilitate processing.\n",
    "\n",
    "I'm still not great at GIS in python so the maps in the publication were done in ESRI ArcMap.\n",
    "\n",
    "### ***Alignment takes a long time due to the large number of scaling candidates. If you want to skip that part and just load the result set, set `skip_alignment` below to `True`.***\n",
    "\n",
    "### Notes on settings and constants\n",
    "**CUTTER_SIZE** defines the half height in mm of the density cutter used as reference. Can be changed to accommodate different sampler sizes. No need to change this for ECCC data. \n",
    "\n",
    "**WINDOW_SIZE** defines the size of the rolling window used in SLF shot noise calculations. A 5 mm window was used in [Proksch, et al., 2015](https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2014JF003266) when there was separation between the SMP and density cutter. Increasing the window reduces sensitivity to sharp transitions and reduces resolution of the analysis. However, moving to something like 2.5 mm makes comparison difficult as some of the very fine structure resolved is very different over the ~10 cm separation between the SMP and density profiles. \n",
    "\n",
    "**NUM_TESTS** defines how many random scaling configurations to test against when attempting to align the SMP and snow pit data. We brute-force the alignment in our paper so `NUM_TESTS` must be large to ensure the test space searched is sufficient. A lower number of tests risks poor alignment and therefore poor calibration. In the paper we use 10k permutations.\n",
    "\n",
    "**MAX_STRETCH_LAYER** and **MAX_STRETCH_OVERALL** define how much an individual layer can be eroded or dilated, and the maximum change in total length of the SMP profile, respectively. We allow a rather large 70% change to individual layers to accommodate pinching out but restrict the total change to 10% to avoid overfitting.\n",
    "\n",
    "**H_RESAMPLE** and **L_RESAMPLE** define the resampled resolution of the SMP and the layer size used for matching profiles. These terms are interactive with the layer stretching and should be evaluated carefully if changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Community packages\n",
    "import os \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "from scipy import stats\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "\n",
    "# Local packages\n",
    "import smpfunc # SMP helper functions\n",
    "\n",
    "# Import SLF SMP Package\n",
    "# To ensure compatablility build from the environment.yml\n",
    "from snowmicropyn import Profile, proksch2015, loewe2012\n",
    "\n",
    "# Import field data\n",
    "pit_summary = pd.read_csv(\"./data/Pit/pit_summary.csv\")\n",
    "pit_density = pd.read_csv(\"./data/Pit/pit_density.csv\")\n",
    "input_data = os.path.abspath(\"./data/SMP/Calibration\")\n",
    "\n",
    "# Set constants\n",
    "CUTTER_SIZE = 15  # Half the height of the density cutter in mm\n",
    "WINDOW_SIZE = 5  # SMP analysis window in mm\n",
    "H_RESAMPLE = 1  # Delta height in mm for standardized SMP profiles\n",
    "L_RESAMPLE = 50  # Layer unit height in mm for SMP matching\n",
    "MAX_STRETCH_LAYER = 0.75  # Max layer change in % of height\n",
    "MAX_STRETCH_OVERALL = 0.15  # Max profile change in % of total height\n",
    "NUM_TESTS = 10000  # Number of scaling candidates to generate for alignment testing \n",
    "\n",
    "# Set conditions\n",
    "skip_alignment = True # Set as true to just load the results from a pickle instead of reprocessing\n",
    "paper_conditions = True # Set as true to reproduce the paper results with seeding\n",
    "\n",
    "# Small differences in comparison to the paper will occur if a seed is not set.\n",
    "# This is mainly because we use a brute-force approach for matching the SMP and \n",
    "# snow pit profiles with a modest search size (specified by NUM_TESTS).\n",
    "if paper_conditions:\n",
    "    np.random.seed(2019) \n",
    "\n",
    "os.makedirs('./output/figures', exist_ok=True)    \n",
    "\n",
    "# Helper functions\n",
    "def rmse(data):\n",
    "    return np.sqrt(np.mean(data**2))\n",
    "\n",
    "def load_smp(smp_file):\n",
    "    p = Profile.load(smp_file)\n",
    "    p = smpfunc.preprocess(p, smoothing = 0, noise_threshold = 0.01)\n",
    "    ground  = p.detect_ground()\n",
    "    surface  = p.detect_surface()\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the SMP calibration profiles at each snow pit (n = 25).  \n",
    "Preprocess each SMP file to remove noise and detect ground and air interfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SMP calibration profiles loaded: 25\n"
     ]
    }
   ],
   "source": [
    "file_list = []\n",
    "for file in os.listdir(input_data):\n",
    "    if file.endswith(\".pnt\"):\n",
    "        file_list.append(os.path.join(input_data, file))\n",
    "        \n",
    "smp_data = [load_smp(file) for file in file_list]\n",
    "\n",
    "print(\"SMP calibration profiles loaded: {}\".format(len(smp_data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary information on the snow pit properties\n",
    "\n",
    "**pit_summary** holds summary information on all now pits. Here is a breakdown of the variables within:\n",
    "  - **ID**: Site identification code  \n",
    "  - **IOP**: Campaign identity code where a = Alert and e = Eureka  \n",
    "  - **ICE**: Ice surface type where m = MYI and f = FYI  \n",
    "  - **DOY**: Day of year  \n",
    "  - **Date**: Date snow pit was completed  \n",
    "  - **SITEID**: Campaign-specific site identification code  \n",
    "  - **LAT** and **LON** decimal degree WGS-84 coordinates  \n",
    "  - **PITFILE**: Internal ECCC snow pit file names (Available on request)  \n",
    "  - **SMPF**: Last 4 digits of the corresponding SMP file name  \n",
    "  - **MPD** Average of nearby Magnaprobe measured snow thickness  \n",
    "  - **SD**: Total snow pit thickness  \n",
    "  - **PD**: Operator-recorded SMP penetration, better to check directly from SMPF  \n",
    "  - **OFF**: Offset between retracted SMP tip and the snow surface. When mounted on the support structure this generally does not change unless there is settling  \n",
    "  - **NOTESMP**: SMP operator notes  \n",
    "  - **NOTESPIT**: Snow pit operator notes  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>ID</th>\n",
       "      <th>IOP</th>\n",
       "      <th>ICE</th>\n",
       "      <th>DOY</th>\n",
       "      <th>Date</th>\n",
       "      <th>SITEID</th>\n",
       "      <th>LAT</th>\n",
       "      <th>LON</th>\n",
       "      <th>PITFILE</th>\n",
       "      <th>SMPF</th>\n",
       "      <th>MPD</th>\n",
       "      <th>SD</th>\n",
       "      <th>PD</th>\n",
       "      <th>OFF</th>\n",
       "      <th>NOTESMP</th>\n",
       "      <th>NOTESPIT</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>a</td>\n",
       "      <td>m</td>\n",
       "      <td>101</td>\n",
       "      <td>11/04/2017</td>\n",
       "      <td>3</td>\n",
       "      <td>83.98343</td>\n",
       "      <td>-66.35089</td>\n",
       "      <td>170411_Site3_Pit1</td>\n",
       "      <td>1848</td>\n",
       "      <td>0.52</td>\n",
       "      <td>0.45</td>\n",
       "      <td>0.6</td>\n",
       "      <td>0.16</td>\n",
       "      <td>Probe went down fine. Tilted forward a few deg...</td>\n",
       "      <td>Completed ~100 m away from the aircraft landin...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   ID IOP ICE  DOY        Date  SITEID       LAT       LON            PITFILE  \\\n",
       "0   1   a   m  101  11/04/2017       3  83.98343 -66.35089  170411_Site3_Pit1   \n",
       "\n",
       "   SMPF   MPD    SD   PD   OFF  \\\n",
       "0  1848  0.52  0.45  0.6  0.16   \n",
       "\n",
       "                                             NOTESMP  \\\n",
       "0  Probe went down fine. Tilted forward a few deg...   \n",
       "\n",
       "                                            NOTESPIT  \n",
       "0  Completed ~100 m away from the aircraft landin...  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pit_summary.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SMP profile total: 25\n",
      "Snow pit total: 25\n",
      "\n",
      "Density total: 219\n",
      "\n",
      "Snow pits by campaign:\n",
      "e    20\n",
      "a     5\n",
      "Name: IOP, dtype: int64\n",
      "\n",
      "Snow pits by ice type:\n",
      "f    17\n",
      "m     8\n",
      "Name: ICE, dtype: int64\n",
      "\n",
      "Density measurements by layer type:\n",
      "F    88\n",
      "H    86\n",
      "R    36\n",
      "I     6\n",
      "N     3\n",
      "Name: TYPE, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "print('SMP profile total: %i' % len(smp_data))\n",
    "print('Snow pit total: %i' % len(pit_summary))\n",
    "print('\\nDensity total: %i' % len(pit_density))\n",
    "print('\\nSnow pits by campaign:')\n",
    "print(pit_summary.IOP.value_counts())\n",
    "print('\\nSnow pits by ice type:')\n",
    "print(pit_summary.ICE.value_counts())\n",
    "\n",
    "# Note: The small number of recent snow measurements were rolled into rounded classification\n",
    "#       for the paper. In the future many more measurements are needed to chracterize recent snow\n",
    "\n",
    "print('\\nDensity measurements by layer type:')\n",
    "print(pit_density.TYPE.value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extraction of corresponding snow pit density and SMP properties\n",
    "This section relies on functions imported from `smpfunc.py`. We iterate through all the SMP profiles that are associated with snow pits and do the following:\n",
    "\n",
    "1. Restructure the snow pit density estimates into a dataframe (`density_df`)  \n",
    "2. Make first guess estimates of density using the P15 coefficients (`p2015`)  \n",
    "3. Resample the SMP and snow pit arrays to a common thickness axis  \n",
    "4. Generate scaling candidates (`random_tests`) according to `L_RESAMPLE`, `NUM_TESTS`, `MAX_STRETCH_OVERALL`, and `MAX_STRETCH_LAYER` \n",
    "5. Scale the profiles and extract the SMP estimates corresponding with each density cutter measurement (`compare_profiles`)  \n",
    "6. Evaluate retrieval skill according to correlation and rmse (`retrieved_skill`)  \n",
    "7. Save properties of the best scaling candidate for further analysis  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished with S34M0295\n",
      "Finished with S34M0309\n",
      "Finished with S34M0315\n",
      "Finished with S34M0335\n",
      "Finished with S34M0361\n",
      "Finished with S34M0439\n",
      "Finished with S34M0462\n",
      "Finished with S34M0509\n",
      "Finished with S34M0526\n",
      "Finished with S34M0541\n",
      "Finished with S34M0633\n",
      "Finished with S34M0649\n",
      "Finished with S34M0666\n",
      "Finished with S34M0735\n",
      "Finished with S34M0752\n",
      "Finished with S34M0776\n",
      "Finished with S34M0844\n",
      "Finished with S34M0869\n",
      "Finished with S34M0875\n",
      "Finished with S34M0956\n",
      "Finished with S34M1848\n",
      "Finished with S34M1856\n",
      "Finished with S34M1874\n",
      "Finished with S34M1894\n",
      "Finished with S34M1908\n",
      "Average scaling: 0.078\n"
     ]
    }
   ],
   "source": [
    "comparison_df = pd.DataFrame()\n",
    "min_scaling_coeff = []\n",
    "\n",
    "if not skip_alignment:\n",
    "    for smp in smp_data:\n",
    "        smp_file_num = int(smp.name[-4:])\n",
    "        pit_df  = pit_summary[pit_summary['SMPF'] == smp_file_num] # Select the matching pit\n",
    "        density_df = pit_density[pit_density['ID'] == pit_df['ID'].values[0]]\n",
    "        density_df = density_df.assign(relative_height=np.abs(((density_df['TOP']*10) - CUTTER_SIZE) - density_df['TOP'].max()*10).values)\n",
    "\n",
    "        # Make first guess at microstructure based on original profile\n",
    "        l2012 = loewe2012.calc(smp.samples_within_snowpack(), window=WINDOW_SIZE)\n",
    "        p2015 = proksch2015.calc(smp.samples_within_snowpack(), window=WINDOW_SIZE)\n",
    "\n",
    "        # Estimate offset of the snow depth and SMP profile\n",
    "        smp_profile_height = p2015.distance.max()\n",
    "        smp_height_diff = pit_df.MPD.values*1000 - smp_profile_height\n",
    "\n",
    "        # Create new SMP resampled arrays and determine the number of layers\n",
    "        depth_array = np.arange(0, p2015.distance.max() + smp_height_diff, H_RESAMPLE)\n",
    "        density_array = np.interp(depth_array,p2015.distance,p2015.P2015_density)\n",
    "        force_array = np.interp(depth_array,p2015.distance,l2012.force_median)\n",
    "        l_array = np.interp(depth_array,p2015.distance,l2012.L2012_L)\n",
    "\n",
    "        smp_df = pd.DataFrame({'distance': depth_array, \n",
    "                               'density': density_array,\n",
    "                               'force_median': force_array,\n",
    "                               'l': l_array})\n",
    "\n",
    "        # Generate a selection of random transformations to evaluate brute-force alignment\n",
    "        # We use a brute force approach because there was no gradient that could be used to optimize the relationship\n",
    "        num_sections = np.ceil(len(smp_df.index)/L_RESAMPLE).astype(int)\n",
    "        random_tests = [smpfunc.random_stretch(x, MAX_STRETCH_OVERALL, MAX_STRETCH_LAYER) for x in np.repeat(num_sections, NUM_TESTS)] \n",
    "\n",
    "        scaled_profiles = [smpfunc.scale_profile(test, smp_df.distance.values, smp_df.density.values, L_RESAMPLE, H_RESAMPLE) for test in random_tests]\n",
    "        compare_profiles = [smpfunc.extract_samples(dist, rho, density_df.relative_height.values, CUTTER_SIZE) for dist, rho in scaled_profiles]\n",
    "        compare_profiles = [pd.concat([profile, density_df.reset_index()], axis=1, sort=False) for profile in compare_profiles]\n",
    "        retrieved_skill = [smpfunc.calc_skill(profile, CUTTER_SIZE) for profile in compare_profiles]\n",
    "        retrieved_skill = pd.DataFrame(retrieved_skill,columns = ['r','rmse','rmse_corr','mae'])\n",
    "\n",
    "        min_scaling_idx = retrieved_skill.sort_values(['r', 'rmse_corr'], ascending=[False, True]).head(1).index.values\n",
    "        min_scaling_coeff.append(random_tests[int(min_scaling_idx)])\n",
    "        \n",
    "        dist, scaled_l =  smpfunc.scale_profile(min_scaling_coeff[-1], smp_df.distance.values, smp_df.l.values, L_RESAMPLE, H_RESAMPLE)\n",
    "        dist, scaled_force_median = smpfunc.scale_profile(min_scaling_coeff[-1], smp_df.distance.values, smp_df.force_median.values, L_RESAMPLE, H_RESAMPLE)\n",
    "\n",
    "        result = compare_profiles[int(min_scaling_idx)].assign(l=smpfunc.extract_samples(dist, scaled_l, density_df.relative_height.values, CUTTER_SIZE).mean_samp,\n",
    "                                                  force_median=smpfunc.extract_samples(dist, scaled_force_median, density_df.relative_height.values, CUTTER_SIZE).mean_samp)\n",
    "        comparison_df = comparison_df.append(result, ignore_index=True)\n",
    "        print(\"Finished with {}\".format(smp.name))\n",
    "    \n",
    "    # Summary of scaling\n",
    "    smp_thickness = [p.samples_within_snowpack().distance.max() for p in smp_data]\n",
    "    scaling_total = [s.sum() for s in min_scaling_coeff]\n",
    "    scaling_mean_abs =  np.round(np.abs(np.array(scaling_total)).mean(), 3)\n",
    "    print('Average scaling: %0.3f' % scaling_mean_abs) # in %\n",
    "\n",
    "    # save the results to a local file since tbe brute-force method takes a while to compute\n",
    "    comparison_df.to_pickle(\"./output/smp_pit_comparison.pkl\")\n",
    "else:\n",
    "    # Load previously-generated results if desired\n",
    "    comparison_df = pd.read_pickle('./output/smp_pit_comparison.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filtering of results\n",
    "Apply a few rules to clean up the data:\n",
    "1. Drop comparisons where any field is not a number (`NaN`), this can occur when a density sample overlaps with no SMP data (ie. depth hoar omissions)\n",
    "2. Drop comparisons where the density cutter is not fully within the SMP profile (identified by the SMP count)\n",
    "3. Drop comparisons where there is new snow or ice in the sample (Not enough samples to support robust analysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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>count_samp</th>\n",
       "      <th>mean_samp</th>\n",
       "      <th>median_samp</th>\n",
       "      <th>stdev_samp</th>\n",
       "      <th>index</th>\n",
       "      <th>ID</th>\n",
       "      <th>TOP</th>\n",
       "      <th>BOTTOM</th>\n",
       "      <th>RHO</th>\n",
       "      <th>NOTES</th>\n",
       "      <th>TYPE</th>\n",
       "      <th>relative_height</th>\n",
       "      <th>l</th>\n",
       "      <th>force_median</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>31</td>\n",
       "      <td>584.201003</td>\n",
       "      <td>617.731699</td>\n",
       "      <td>6.604018e+01</td>\n",
       "      <td>209</td>\n",
       "      <td>24</td>\n",
       "      <td>18.0</td>\n",
       "      <td>11.5</td>\n",
       "      <td>442.5</td>\n",
       "      <td>-9999</td>\n",
       "      <td>R</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.172577</td>\n",
       "      <td>10.598897</td>\n",
       "      <td>141.701003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>31</td>\n",
       "      <td>464.755238</td>\n",
       "      <td>536.770771</td>\n",
       "      <td>1.224404e+02</td>\n",
       "      <td>210</td>\n",
       "      <td>24</td>\n",
       "      <td>11.5</td>\n",
       "      <td>9.0</td>\n",
       "      <td>385.5</td>\n",
       "      <td>-9999</td>\n",
       "      <td>F</td>\n",
       "      <td>80.0</td>\n",
       "      <td>0.370917</td>\n",
       "      <td>6.120971</td>\n",
       "      <td>79.255238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>31</td>\n",
       "      <td>276.171140</td>\n",
       "      <td>281.557287</td>\n",
       "      <td>1.966974e+01</td>\n",
       "      <td>211</td>\n",
       "      <td>24</td>\n",
       "      <td>9.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>199.0</td>\n",
       "      <td>-9999</td>\n",
       "      <td>H</td>\n",
       "      <td>105.0</td>\n",
       "      <td>0.813759</td>\n",
       "      <td>0.446445</td>\n",
       "      <td>77.171140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>31</td>\n",
       "      <td>281.557287</td>\n",
       "      <td>281.557287</td>\n",
       "      <td>5.684342e-14</td>\n",
       "      <td>212</td>\n",
       "      <td>24</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>358.0</td>\n",
       "      <td>-9999</td>\n",
       "      <td>F</td>\n",
       "      <td>175.0</td>\n",
       "      <td>0.685753</td>\n",
       "      <td>0.315619</td>\n",
       "      <td>-76.442713</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>31</td>\n",
       "      <td>681.181976</td>\n",
       "      <td>690.728469</td>\n",
       "      <td>6.660636e+01</td>\n",
       "      <td>213</td>\n",
       "      <td>25</td>\n",
       "      <td>22.5</td>\n",
       "      <td>11.5</td>\n",
       "      <td>480.0</td>\n",
       "      <td>-9999</td>\n",
       "      <td>R</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.178738</td>\n",
       "      <td>36.294175</td>\n",
       "      <td>201.181976</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   count_samp   mean_samp  median_samp    stdev_samp  index  ID   TOP  BOTTOM  \\\n",
       "0          31  584.201003   617.731699  6.604018e+01    209  24  18.0    11.5   \n",
       "1          31  464.755238   536.770771  1.224404e+02    210  24  11.5     9.0   \n",
       "2          31  276.171140   281.557287  1.966974e+01    211  24   9.0     2.0   \n",
       "3          31  281.557287   281.557287  5.684342e-14    212  24   2.0     0.0   \n",
       "4          31  681.181976   690.728469  6.660636e+01    213  25  22.5    11.5   \n",
       "\n",
       "     RHO  NOTES TYPE  relative_height         l  force_median       error  \n",
       "0  442.5  -9999    R             15.0  0.172577     10.598897  141.701003  \n",
       "1  385.5  -9999    F             80.0  0.370917      6.120971   79.255238  \n",
       "2  199.0  -9999    H            105.0  0.813759      0.446445   77.171140  \n",
       "3  358.0  -9999    F            175.0  0.685753      0.315619  -76.442713  \n",
       "4  480.0  -9999    R             15.0  0.178738     36.294175  201.181976  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = comparison_df.dropna() \n",
    "result = result[result['count_samp']>=CUTTER_SIZE*2] # Remove comparisons outside the profile\n",
    "result = result[~result['TYPE'].isin(['N', 'I'])] # Remove new snow and ice because we don't have enough samples\n",
    "result['error'] = result['mean_samp']-result['RHO']\n",
    "result.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare extracted SMP density against *in situ*\n",
    "SMP density estimates using the Proksch et. al. 2015 coefficients are compared here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proksch et al. 2015 Eval.\n",
      "N: 196\n",
      "RMSE: 125.0\n",
      "bias: 104.0\n",
      "r^2: 0.72\n"
     ]
    }
   ],
   "source": [
    "# P2015 evaluation stats\n",
    "p2015_rmse = np.sqrt(np.mean(result['error']**2))\n",
    "p2015_bias = (result['error']).mean()\n",
    "p2015_r2 = np.ma.corrcoef(result['mean_samp'],result['RHO'])[0, 1]**2\n",
    "p2015_n = len(result['mean_samp'])\n",
    "p2015_p = stats.pearsonr(result['mean_samp'],result['RHO'])[1]\n",
    "\n",
    "print('Proksch et al. 2015 Eval.')\n",
    "print('N: %i' % p2015_n)\n",
    "print('RMSE: %0.1f' % np.round(p2015_rmse))\n",
    "print('bias: %0.1f' % np.round(p2015_bias))\n",
    "print('r^2: %0.2f' % p2015_r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAFUCAYAAAB2uskDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXzU5bn//9edyTIJkBA2USKy1CgqiJTiGhfWhNWWQku1rTzo8ZzWY885tdoj2Pb8VKA9be059Wi/eor8vr8uViNtRSUJiwtxFxDBNUgoGgRZEjMhyWSZ3L8/JolDyCSTmc8sSd7Px4MHyZ177rmS1nDNPdd93cZai4iIiIiIJIakeAcgIiIiIiKfU4IuIiIiIpJAlKCLiIiIiCQQJegiIiIiIglECbqIiIiISAJRgi4iIiIikkCUoIuISEIzxlxmjHnSGLOwB49JMsb8xhhz0hhTZIxJjWaMIiJOUoIuIiIJyxiTBZwLLKRn/2ZNAv4P8AMgHzjP+ehERKIjOd4BiIiIBGOtrTbGlLZ9boxJAx4B3gcOW2t/F+Rxu1vnjwGKgLejHqyIiEOUoIuISG+yAhgN/ANwyBhTB9zcYc4Ra+3XjTHjgF8DOcBEYE9MIxURCZMSdBER6U0mAtn4S162AU9ba//U2URrbbkxJg/YC8xBCbqI9BKqQRcRkd7kIDAS2AR8D6jvbJIx5mJjzGBr7VHgBeDl2IUoIhKZhEnQjTHjjDG2w59fGmNuNMbcb4z5izFmTsD8TsdFRKTvMMZkAPNaP80D/oh/J3wn8FVrbVOQh34beMUY84/AX621L0U9WBERhxhrbbxjAPwJN/BDYB/gA2YA84EX8Z/Gvwb4Of7awyFAWcdxa2117CMXEREREXFOItWgv2WtnQztOyZvAJfjj7EG8AKZ+BP30UHG/xL7sEVEREREnJMwCbq1dm/Ap/Pxt8U6t/Xz5oCvjQfGBBk/jTHmZlpP+A8YMOCL559/vhPhiojExIABA7jooot48MEHj1trh8c7nlgaNmyYHTNmTLzDEBHpsZ07d0b0OzthEvQOvg78DLihk68ZTk3MA8dPY619GHgYYOrUqXbHjh1OxSgiElWHDh1iy5YtpKSk8OCDDx6MdzyxNmbMGPQ7W0R6I2NMRL+zEy5BN8YMAi6w1r5ujLm0dTjwMOt+Pj+133FcRKRP2L9/P8899xxZWVkUFBTEOxwREYmhhOniEuB64KnWj/+Gf7c8G0gBTgBbuhgXEen19uzZw7Zt2zjjjDNYuHAhAwcOjHdIIiISQwm3gw58DbgbwFr7cWt3l9uAgcDXrLUewBNkXESk17LW8uqrr7J3717Gjh3L9OnTcblc8Q5LRERiLOESdGvt/A6fPwY81sm8TsdFRHojn8/H888/z/79+7nwwgu54oorMKbTozUiItLHJVyCLiLS3zQ2NrJ582Y++eQTLr30UiZNmtSnk3NjzDBgKf5bQfdYaz+Oc0giEqLGxkZuu+02Xn/9dXw+X7zDiRuXy8W0adP41a9+RWpqquPrK0EX6SfKysooLS0lLy+P3NzceIcjrWpraykuLqaqqorrrruOc889t/sH9WLGmAvwd9ZaYK2tinc8ItIzv/vd70hNTWX79u2kpaXFO5y4aWho4M477+R3v/sd3/ve9xxfPxEPiYqIw8rKypgyZQrf//73mTJlCmVlZfEOSYCqqiqefPJJPB4P+fn5/SE5TwKewN8W96/GmLXGGG0UifQif/vb3/jBD37Qr5NzgLS0NH7wgx/w5JNPRmV9/WIU6QdKS0ux1lJXV0dGRgalpaXaRY+zTz/9lOLiYpKSkliwYAHDhg2Ld0ixMAuYAFwHnAH8GfgYeLBtQuDlcqNHj45DiCLSlaqqKkaOHBnvMBLCmWeeSVVVdN4IVIIu0g/k5eVhjCEjIwNjDHl5efEOqV/7+9//zrZt2xg4cCBz585l0KBB8Q4pVia2/u0CtrV+fC0BCXrHy+ViGZyIdM9ae0p3qcrKSjye6DXSy8zMZMiQIVFbPxIul4uWlpaorK0EXaQfyM3NZdeuXapBTwDvvvsuL730EsOHDyc/Px+32x3vkGLJ2/r3cKBt26kpTrGISIQqKysZN24c1dXVUXuOrKwsysvLEzZJjxYl6CL9RG5urhLzOLLWsmPHDt58801Gjx7NzJkzSU7ud7+CNwG/As4DMlrHngo+XWKltrYWr9fb7Ty3282AAQNiEJH0Bh6Ph+rqapYvX05mZmZU1l+/fj0ej6fbBH337t1861vfYu/evXz729/m6NGj1NbWsnr1ajZt2sTWrVt5/fXXAXjxxRc5cuQIX/3qVx2P2Sn97l8HEZFYa2lpobS0lA8++IDzzz+fq666iqSk/ndG31pbHnDJXBbwc3SfRdzV1tZyy4/u4khNXbdzRw7K4IGf36skXU6RmZlJdnZ2XGOYPHkyU6ZMYe/evdx0001MnDiRYcOGsXTpUr73ve+1v4B47733WLZsGb/4xS/iGm93lKCLiERRU1MTW7du5eOPP+aLX/wiU6ZM6dM9zrtjrS0ECuMdh3zO6/VypKaOwTfcTlpm8CSrwVPFkT/+Aq/XqwRdEt7hw4cBGDhwIOvXr6e8vJzGxkZuv/12Kioq+O1vf4vH4+Hmm2+Oc6SdU4IuIhIl9fX1FBcXc/z4cfLy8pgwYUK8QxIJKi0zG3f20HiHIRKxP/zhDzQ2NnLXXXfxr//6r9xwww2Ul5eTmprKV7/6VZ555hm++93v8vWvfz3eoQalBF1EJAo8Hg9FRUXU1tYye/ZszjnnnHiHJCLSL9x4441ce+218Q4jIkrQRUQcduzYMYqLi7HWMm/ePM4444x4hyQiIr1I/zulJCISRR9//DFPP/00ycnJLFy4UMm5iEiMWRv8CoXU1FQAfD4fx48fj1VIPaYddBERh5SVlbF9+3aGDBlCfn4+GRkZ3T9IRKSXi9ZFRT1Zd/fu3ezatQuAv/3tb+1dXA4fPsy+ffsAePXVV5kxYwYTJkzg/vvv58ILL0zYW5yVoIuIRMhay+7du3njjTcYNWoUs2bNat+lERHpqzIzM8nKymL9+vVRe46srKyQeqxPnjyZPXv2nDZ+5plnsn///lPG3n33XcfiixYl6CIiEbDW8tJLL/Huu+9y7rnncvXVV59yDbaISF81ZMgQysvLo7aDDv4XAf3tFlFQgi4iErbm5maeffZZ/v73v3PxxRczbdq0ft3jXET6nyFDhvTLBDralKCLiIShoaGB4uJijh49yhVXXMFFF10U75BERKSPUIIuItJDJ0+epKioCI/Hw4wZMxg3bly8QxIRkT5ECbqISA+cOHGC4uJimpqamDt3LmeeeWa8QxIRkT5GfdBFYqSsrIx169ZRVlYW71AkTJ988glPPfUUAAsXLlRyLiIiUaEddJEYKCsrY8qUKVhrMcawa9cucnNz4x2W9MD+/ft5/vnnyczMpKCggIEDB8Y7JBER6aOUoIvEQGlpKdZa6urqyMjIoLS0VAl6L7J3715eeeUVRo4cyZw5c0hLS4t3SCIiCaG2thav1xu19d1uNwMGDIja+olKCbpIDOTl5WGMISMjA2MMeXl58Q5JQmCt5bXXXmPPnj2MHTuW6dOnq8e5iEir2tpabvnRXRypqYvac4wclMEDP7+33yXpStBFYiA3N5ddu3ZRWlpKXl6eds97AZ/Px/PPP8/+/fu58MILueKKK9TjXEQkgNfr5UhNHYNvuJ20zGzH12/wVHHkj7/A6/WGnKDff//97N+/n2HDhnH8+HFGjBjBxIkT+cY3vsHJkyf58Y9/zN13390+/6c//Sl33303Z511Fn/6058488wz+d73vse2bdu4/vrrSUlJ4cCBA9x1110sWrTI8e8xGCXoImEqKyvrUcKdm5urxLyXaGxsZMuWLRw6dIhp06Zx8cUXKzkXEQkiLTMbd/bQeIfBzp07+fGPf0xVVRXGGMrLy/nxj3/MypUrufHGG3nkkUd45JFH+OlPf4rL5cLn81FdXQ3AueeeyzXXXAPAVVddxbZt21i0aBHf/OY3GTZsGEuWLKGysjJm54/UxUUkDG2HPr///e8zZcoUdWbpQ2pra3nqqac4fPgw1113HZMnT1ZyLiLSCxw9epTq6mpuuOEGDh8+zLhx47jpppsAOOOMM1i6dCmHDh3imWeeAWDLli0sWLCgyzWrqqrwer0MHz48puePlKCLhCHw0Ke1ltLS0niHJA747LPP2LhxIx6Ph/z8fM4999x4hyQiIiG67rrruPLKK3n00UcZO3Yst912G5dddln712+99VYAHn74YQAOHTrE+PHjg6739NNP86Mf/YhbbrmFl19+mZSUlOh+AwGUoIuEwalDn+qNnjg+/fRTnnzySZqbm1mwYAE5OTnxDklERHqgqamJrVu3snr1ajIyMrjvvvu4+uqraW5uBmDatGlMmzaN4uJi3nzzTc4+++wu15s/fz7r1q3jl7/8Jeecc04svoV2StBFwtB26PM3v/lNlz3Nu0rAA8tkJk+ezOrVq5Wox8nf//53nn76adxuN4sWLWLYsGHxDklERHroySefpLa2lpUrV3LgwAG+8pWvsHv3bvbs2dM+59Zbb8Xn83Hbbbcxc+bMOEbbNR0SFQlTd4c+u7ucKLBMBuDuu+9m7dq1cbnEqKcHXvuS9957jxdffJHhw4czZ84c0tPT4x2SiESgsrISj8fT7bzMzEyGDBkSg4j6vgZPVUKs29zczE9+8hMeeOABsrKyWLVqFVu2bGHcuHHtc5YuXcoPf/hD5syZQ1JS1/vU1tqw4naCEnSRKOnucqK2MpnU1FQaGxtpbGwkOTk55pcY9ddbTq217Ny5k127djF69GhmzJgR0/pCEXFeZWUl48aNa+/M0ZWsrCzKy8uVpEfA7XYzclAGR/74i6g9x8hBGbjd7pDmpqSkkJyczB133MHIkSN59913eeqpp9i0aRN/+MMfmD17NldccQW33norK1aswOPx8PjjjwOwb98+tm/fzsiRI3nxxRcBKC4uZtasWXEpeVSCLhIl3dWpt5XJFBYWsnr1aowxcbnEqD/ectrS0kJpaSkffPAB5513Hnl5ed3upIhI4vN4PFRXV7N8+XIyMzO7nLd+/Xo8Ho8S9AgMGDCAB35+b8LcJHrDDTdwww03dPq1b3zjG+0fr1q1qv3jO+64gzvuuOOUuVu3bg0jUmcpQReJklAuJ8rNzWXVqlUsWbIkbiUm/e2W07ZDRB9//DFf/OIXmTJlitooivQxmZmZZGc7f3GOnG7AgAH97pbPWFCCLhJFoV5O5MQlRuHWkfenW07r6+spLi7m+PHj5OXlMWHChHiHJCIichol6CJ9QKR15P3hllOPx0NRURG1tbXMnj075i2zREREQpVwCboxZhiwFDgI7AGuAS4FRgEPWWtLWufd2Nm4SH/UH+vIe+LYsWMUFxdjrWXevHmcccYZ8Q5JRKRXMsbg8/lwuVzxDiXufD5f1M4vJVSCboy5AHgYWGCtrTLGjAXWA5PwJ+qPG2NGA0M6G7fWdn9sW6QP6m915D3x8ccfs3XrVtxuNwUFBQwePDjeIYmI9FrZ2dkcOXKEUaNGxTuUuDt8+HDUzjokTNsCY0wS8ARggL8aY9YCi/G/iKgBvEAmMANYFGRcJKYS5SbQUC9O6m/KysooKSkhMzOTRYsWKTkXEYnQ9ddfz3333UdDQ0O8Q4mrhoYG7rvvPhYtWhSV9RNpB30WMAG4DjgD+HPA15oDPh4PjAkyLhIzodZ9x+oSoP5QRx4qay27d+/mjTfeYNSoUcyaNYvU1NR4hyUi0ut95zvf4bbbbuPqq6/G5/PFO5y4cblcTJs2je985ztRWT+REvSJrX+7gG2tH28FOt7Dajg1MQ8cP33QmJuBmwFGjx4deZTS77Ul3EeOHOm27ru/XgIUT9ZaXn75Zd555x2+8IUvcM0116hWUkTEIampqdx///3xDqPPS6QEva3L/XCg7W7XKa1/B5bi7Afqg4yfxlr7MP66dqZOnRq/O1ulTwhMuNuuAO6q7luHN2OrubmZ5557jgMHDnDxxRczbdo09TgXEZFeJ5ES9E3Ar4DzgIzWsX8FHgGygRTgBLAFyALu62RcJKo6JtwrV67scr4Ob8ZOQ0MDJSUlfPrpp1x++eVMnDix+weJiIgkoIRJ0K215a2tE2/Dn4D/HPgD0Ng6NhD4mrXWA3gC5gaOi0RVx4R76tSpLF68GGsta9euPa2EpT9dAhRPJ0+epKioCI/Hw/Tp0xk/XkdSRESk90qYBB3AWlsIFHYYfqz1T8e5nY6L9ERPD3B2TLhDKWGJ9PBmsBhjdfg00VVWVlJUVERTUxMFBQWcddZZ8Q5JgjDGLMf/rmibpa2/90VEJEBCJegisdCW2Obk5LTvfnc8wNkx+e34edtY20HRaJWwlJSUsGjRIpKSkkhKSmqPUYdP/T755BM2b95MSkoKCxcuZMiQIfEOSbp2NbCu9eMW4Kk4xiIikrCUoEu/EpjY+nw+jDF4vd5Tdr87Jr8bNmxg8eLF+Hw+Wlpa2LhxI2PHjm2fA7By5UqWLFniaJJcVlbGwoULaWxsBCA9Pb09Rh0+hf379/P888+TmZlJQUEBAwcOjHdI0gVjzDjgBuBjoBT4qbXW2/WjRET6p4S5qEgkFgIT26SkJFpaWk7Z/S4rK2PNmjX4fD7q6uqw1vLoo4/i8/nwer00NjayaNEiCgsL29cxxjBy5EjHE+TS0tJTrhBuamoiJycH0OHTvXv38uyzzzJ8+HAWLlyo5Lx3uBx/t62zgG8DrxljRnScZIy52Rizwxiz49ixY7GOUUQkIShBlz4p2A2fgYltUlISGzdubL99E2DKlCk8/vjjeL1e3G43xhiWLVtGS0tL+xptSXO0E+S8vDxcLlf7BTsul4vFixdTVlbWb28Otdby6quv8sorrzBmzBjmzZtHWlpavMOSEFhr/2itnQCMBH4CjACu72Tew9baqdbaqcOHD491mCIiCUElLtLndFWf3VVXlXXr1mGtpb6+nvT0dJYuXcrKlSvJzc1l48aNp9SCL1myhKlTp/Loo4+ybNmyqCTIbbGuWbOGxx9/nPr6elwuV3s5S3+7OdTn8/HCCy/w4YcfcuGFF3LFFVeox3kvZK2tBu4xxkzA3yZXREQ6UIIufU539dnBEtuOZSNtyTnAnDlz2LNnT3tiD7QfMH3iiSeitoudm5vLypUreeKJJ/ptOQtAY2MjW7Zs4dChQ0ybNo2LL75YyXkvY4y5FpgB/Ke1tgY4BPwlrkGJiCQoJejS54Rbn91dz/LAxL5tt72nhzTDaY3Y33up19XVUVRURFVVFddee22/+/77kHTgW8ANxpjfAo9Zaw/HOSYRkYSkBF36nEgS2lDLRsJ5ERBJa8T+Vs7S5rPPPqOoqAiv18ucOXM4++yz4x2ShMlaWwScE+84JHYqKyvxeLq/QzAzM1MtUkU6UIIufVK0E9pwXgSoNWLPfPrppxQXF5OUlMT8+fPRgUGR3qOyspJx48ZRXV3d7dysrCzKy8uVpIsECCtBN8aM7uFDTlprK8N5LpF4Kysro7CwkOPHjzNs2DCWLFkCQGFhzy5A7O+tEXvi4MGDbN26lYEDB1JQUEBmZma8QxKRHvB4PFRXV7N8+fIu//v1eDysX78ej8ejBF0kQLg76H8HbA/m/xF/7aFIr1JSUnLKZUEA99xzDwANDQ0ArF69mt27d3e7G97fa8lD9d577/Hiiy8ybNgw8vPzSU9Pj3dIIhKmzMxMsrOz4x2GSK8TboK+AXg6xLkpQE933EXahXOw0qnn7Zicgz8xD7xAyOfznVauEizmnpbexOt7jwdrLTt37mTXrl2cffbZzJw5k5QUdeETEZH+J9wE/afW2ndDnWyMuSDM55E+LJTkM9SDldFIZDve5Bmo7RZS8F8gFFiuEslh0EBOrdMbtLS0UFpaygcffMB5551HXl5e0J+9iIhIXxdWgt6T5Dyc+dL3hZp8hnKwMlqJbNtNnm63m+bmZsB/e2hycjJ//etf2bFjBwBLliw55fmcOgzaXw6VNjU1sW3bNj766COmTJnCF7/4RfU4FxGRfq1HCboxZiRwJbDDWnswOiFJfxBq8hnKwcpga0W6q96xZrztudrWmzNnTqePc+owaH84VFpfX09JSQnHjh0jLy+PCRMmxDskERGRuOs2QTfGrAH+EfgtsALIAt4yxtxrrX0myvFJHxVq8hnKwcrO1nJqV71jzXgoazh1GLSvHyr1eDwUFRVx8uRJZs2axZgxY+IdkoiISEIIZQf9VuDXwBH8nVsmWWs/NMYsj2pk0qd1lXx23Pnu7mBlZ2v19KZPp2vYnerD3lcvKDp27BjFxcW0tLQwf/58zjjjjHiHJCIikjBCSdA3AeuttQeMMS+1JufpwC3A+uiGJ31ZZ8lnuDvfHdfqSXlIfzqMmQgqKirYsmULbrebgoICBg8eHO+QREREEkooCfo3gSnAAaDtKr8LgLeiFZT0X04djOzJDn0iH8bsa20Wy8rK2L59O9nZ2eTn5zNgwIB4hyTSp9XW1uL1erucU1VVFaNoRCRU3Sbo1tpG4NXWT+8Atlprd+KvRxdxlJMHI0PdoQ9Wwx7vxLgv7exba3nrrbd4/fXXGTVqFLNmzSI1NTXeYYn0WpWVlXg8ni7n1NXVcc9//Q9VDb5u16tLycCV5nYqPBGJUE/bLKr3mURVtA9GdrZbvmLFitO6tSRCYpzIO/s9Ya3l5Zdf5p133mH8+PFce+21uFyueIcl0mtVVlYybtw4qquru52bPGocC/7PRgYNG9nlPFeam5QMvaMlkih6mqDbqEQhEiCaByOD7dAHPmdPD5hGS19os+jz+Xj22Wc5cOAAkyZN4tJLL1WPc5EIeTweqqurWb58OZmZmUHnHT16lEe3vYRxD8CdPTSGEYpIpMK9SVSkVwq3bWM89PY2iw0NDZSUlHDkyBEuv/xyJk6cGO+QRPqUzMxMsrOzg369vr4+htGIiJN6mqB/EpUoRGKoJ20bc3JyKC0tbR+Ptd7aZvHkyZMUFRXh8XiYMWMG48ePj3dIIiIivUaPEnRr7U2hzjXGjLDWHu1xRCIxFOwwaNvHwWrRE+EQaaKqrKykqKiIpqYmCgoKOOuss+IdkoiISK8ScYmLMWYk8M/AGUDbyS8DTAYuiXR9kWjprktKsEOafam7itMOHz5MSUkJKSkpLFiwgKFDVfcqIiLSU07UoL8AnNvJuA6USkIrLS2lpaWF+vp60tPTTzsMGqwWva90V3FaeXk5zz33HJmZmRQUFDBw4MB4hyQiItIrOZGge4GvAIE3HbiALzuwtkjU5OTktB+i6uwwVbBDmolyiDSRvP3227zyyiuMGDGC/Px80tLS4h2SiIhIr+VEgv4Y8Ja19kDgoDGm1oG1RaJmx44duFwufD7/JR633HJLeyIeWGO+YsWpd3L19u4qTrLW8vrrr/PWW28xZswYpk+fTnKymkOJiIhEwol/ST8CdhtjTgSMGfw16RkOrN/v6ABi9JWVlbF69er25BwgKSmpvWNLdzXmvbW7ipN8Ph8vvPACH374IRdccAFXXnmlepyLiIg4wIkE/b+BQa1/AqkGPQw6gOisYC92CgsLT0nOU1JSSEpKIi8vTzXmIWhsbGTLli0cOnSIL33pS0yePFnJuUiCqqurpaampss5qakppKW5YxSRiHTHiQR9L3C9tfazwEFjzCIH1u53lBw6J9iLnbbd88bGRgDcbjd33XUXS5Ysaf9ZO11j3pfeFamrq6OoqIjKykquvfbaXv/9iPRVaWlpJCW7KC4uhvSuD22npqaybNkyJekiCcKJBH0LcI0x5s2AMQOMcWDtfkcHEJ0T7MVOaWlp+25vamoqd911F6tWrWp/nNM15n3pXZHPPvuMoqIivF4v+fn5nH322fEOSUSCSE9PZ/To0bgnXEDywKyg8xqbGin74AMaG5uUoIskCCcS9HsIXs7y3w6s36/oAGJ4OtuhDvZip+P4kiVLTlvPyRrzvvKuyKeffkpJSQnGGObPn8/w4cPjHZKIdMPlcpGWmkqyOiuJ9CpOJOjHgfeAloCxJKD3ZSAJQgcQg+ssEQ+2Q93xxQ7AunXryMvLi+mLoL7wrsjBgwfZtm0bGRkZzJ07l8zMzHiHJCIi0mc5kaBfaa3d13HQGHO+A2uLtAuWiHe1Q92WqHf22I7tE6Olt78r8t577/Hiiy8ybNgw8vPzSU9Pj3dIIiIifVrECXpnyXnr+PuRri0SKFgiHsoOdbzLTHrjuyLWWnbt2sXOnTs5++yzmTlzJikpKfEOS0REpM9LqBtFjDHLgUcChpYCacClwCjgIWttSevcGzsbl74rWCLetkNdWFh42mPaSmJycnJ6fZlJLLW0tPDiiy/y/vvvc95555GXl0dSUlK8wxIREekXEipBB64G1rV+3AK8DewBJgHXAI8bY0YDQ4D1HcettdWxD1lipbtSkbVr12KtZe3atezatQs49cKhDRs2UFFR0SvLTGKpqamJbdu28dFHH3HJJZcwdepU9TgXERGJoYRJ0I0x44AbgI+BUuCnwJfxx1gDeIFMYAYwOsj4X2IeuMRUsFKRzkpYgFPGKioqYlZ33lvV19dTUlLCsWPHuOqqq7jgggviHZKIiEi/E/F71saY+1v/HtJh/GJjTKEx5tfGmIwQlroc2A+cBXwbeA04r/VrzQHzxgPnBhnvLL6bjTE7jDE7jh07FkIY0ht1Vv7SF7qnxFJNTQ0bN27kxIkTzJo1S8m5iIhInDixg77EGPMVYKQxpgT4srW2AX8t+UXAfwCrWv8EZa39I/BHY0wW8H3g/wFu6mSq4dTEPHC8s3UfBh4GmDp1arB+7dLLdVaH3lWbRZW4nOr48eMUFRXR0tLCvHnzGDlyZLxDEumXKisr8Xg8Xc6pqKiIynPX1p7s8usNJ2uor6/no48+oqampsu50Ygxnj8bkVhzIkEfAZQBB/An5P8C/CdwPrDXWrvWGPPjUBdrrSO/xxgzAf8u+n9x6k7/fqC+9eOO49LPdVaHDnDgwAEWL14c9DbPzgtkpeUAACAASURBVPqr9xcVFRVs2bKFtLQ05s+fT3Z2drxDkj7OGLMMuNlae128Y0kklZWVjBs3jurq7o9TDRgwgDSHLh9yuVwkJ6ewcePGrifWn4QXXmBK4e9DWtfJGOP1sxGJFycS9F9ba29r+yQgGU/n80S62xNmxphr8deR/6e1tgY4hL+m/JdANpACnAC2AFnAfZ2MSy/lRILcsQ69sLCwPWH3+XxYa2lsbMTtdp/SZjFYf/V4fA+xtm/fPl544QWys7PJz89nwIAB8Q5J+rjWzZeHgZ3xjiXReDweqqurWb58ebeXgaWlpZGREUr1aPeSk5O58MIL8fl8Xc5rPllN3ZhzyL8+tPsQnIwxXj8bkXhxIkGfZoz5IVCJ//DmdGNMW52Bu/XvUP7VTwe+BdxgjPkt8Ji19uPWdoq3AQOBr1lrPYAnyLj0Qk4lyB1rzuHzQ6JpaWk0NjYC4PV6ycnJaX9csAOmPUm2nfoeYsVay549e3jttdc466yzmD17NqmpqfEOS/o4Y8wA4Bv4N1UkiMzMzJi/k5WcnExyctcpgasxlaaUFAYPHhy3BDgePxuReHAiQX8C+HXrxyeBxcD21s8bjTFPAY3dLWKtLQLO6WT8MeCxUMel93HqEqHOas7Xrl1LRkYGzc3NpKWl0dDQQHp6+il1ih0T+5ycnB4n2/G+CKknrLW8/PLLvPPOO4wfP55rr70Wl8sV77Ckf7gd/7ui34x3ICIiicyJBP3PwHv4e5I/ba193xhzPnALsAZYATzjwPNIHxVqt5VQSkg6tmFsS9hzcnJYvHgxLpfrtOfomNiHk2xH0jEmlqUxPp+PZ599lgMHDjBp0iQuvfRS9TiXmGi9iG4b4MNf9ugyxqRaaxs7zLsZuBlg9OjRMY9TRCQROJGg/4e19rvAZgBjzGD8h0SXW2vvobWDikgw3V1ABOGXkAQm7F09R8fEvqfJdijfQ2diWRrT0NDA5s2bOXz4MJdddhmTJk2KyvOIBPFt/BfLtRkNrMTf6audOm+JiDiToH/dGFNhrV3dWhf+S/ydXfSLVUIW7AKiNqWlpfh8Prxe72mHPJ16jsB5HVs2Orl+oFiVxpw8eZKioiI8Hg8zZsxg/PhOrw0QiaZb8R/wBygEjuBvxysiIh04kaBPBzDGlAJX4O9Rvhb4mwNriwCQk5OD1+sFTj/kGS0dWzZGI3GOxWVKlZWVFBUV0dTUREFBAWeddZbjzyHSHWvt3raPjTENQLW19qM4hiQikrAiTtCttW9Ce33hRuDr1to9xpj5ka4t0qaiooL09HTq6+tPO+QZDbHa2Q63NCZUhw8fpqSkhOTkZBYsWMDQoUMdXV8kHNbaMfGOQUQkkYWVoBtjumqW+mbroTMb7voiHeXl5ZGUlBSzQ5ix2NluE05pTCjKy8t57rnnyMzMJD8/n0GDBjn+HCLS/9TX13c/CX/rRrVvFQlPuAl0M3CY4HXmScCwMNcWOU2sD2FGe2c72t5++21eeeUVRowYQX5+vm7VE5GImZQ0mtKz2LTrnZDmDzQtzLv6ipCS9O7eFY32u6YiiSbcBP0Ka22Xt8AZYy4Nc22RTsX6EGa0drajyVrL66+/zltvvcWYMWOYPn16t5ePiIiEwuVOZ8iyH2CbGrqd21xbw8m/3E9zc3OXCXpaWhoDBgwI6V3KAQMGaLNB+o1w/+Xe190Ea+1rbR8bYwZaa0+G+VwiYXOyVCWW/crD4fP52L59O/v27eOCCy7gyiuvVI9zEXGUy50O7nTH1svIyOCWW26hoaH7pD8tLS1uN5iKxFq4CfqNwIOhTDT+DGEu8HiYzyUSdnLsVKlKLPuVh6OpqYnNmzdz6NAhvvSlLzF58mQl5yLSK2RkZCjxFukg3AT9f4wx9/dg/v9FCbqEKdLkONJSlbKyMtasWdPehz2aXV3CUVdXR3FxMSdOnOCaa67hvPPOi3dIIiIiEoFwE/S76dlFRLvDfB6RmLU87Ezbi4OWlpb2S5Ki3dWlJ6qrq9m0aRNer5f8/HzOPvvseIckIiIiEQorQbfW/ofDcYgEFcuWhx21vTho67++dOlSVq5cmRC7559++iklJSUYY5g/fz7Dhw+Pd0giIiLiALV3kF7hzjvv5Pjx4wwbNowDBw7E7LBmxxcHiZKcHzx4kG3btpGRkUFBQQFZWVndP0hERER6BSXoktDaSkza6r/T0tJoaGggPT2dpKQkNmzYQEVFRdSS9UTsh/7+++9TWlrKsGHDyM/PJz3duY4KIiIiEn9K0CVuQunM0lZi4vV6AdpbcdXX15OSksKCBQtISUmJameVROmHbq1l165d7Ny5k7PPPpuZM2eSkpIS77BERETEYUrQJS5C7czSVmLidrtP2UEHf2vBtr/DOTya6H3NA7W0tPDiiy/y/vvvk5uby9VXX01SUlK8wxIREZEoiDhBN8b8BLjHWtuTri7Sz4XamSWwxCQnJ4eKigrefvttHnzwQRobGwFITU3t8eHRRO9rHqi5uZlt27Zx8OBBLrnkEqZOnaoe5yIi0qvU1ta2vxveFbfbzYABA2IQUWJzYgf9CuAeY8xnwDPW2vccWFP6uJ50ZulYYlJWVsb//u//kpycjLWWVatWsWTJkh4l2PFs3dgTXq+X4uJijh07xlVXXcUFF1wQ75BERER6pLa2llt+dBdHauq6nTtyUAYP/Pzefp+kO5Ggf8NaW2mMSQEWG2N+DTwD/Nlae8yB9SVOolkCEsnhSycObsazdWOoampq2LRpEydPnmTmzJmMHTs23iGJiIj0mNfr5UhNHYNvuJ20zOyg8xo8VRz54y/wer1K0B1YY6QxZgjwD8ByYBgwG7jJGPMo8I61tsiB55EYikUJSCSHLyM9uJmI3VkCHT9+nOLiYnw+H/PmzWPkyJHxDklERCQiaZnZuLOHxjuMXsGJBP3N1nUM0Aj8GfittbYUwBhzgTHm19baf3PguSRGeksJSCQSpTtLRxUVFWzZsoW0tDTmzZtHdnbw3QYRERHpe5xI0FOACuC3wDpr7dG2LxhjMoHNwEBACXov0htKQPqiffv28cILL5CdnU1+fn6/f4tPRESkP3IiQT8A5FprfW0DxpgJwBGgGngIeN+B55EYSvQSkL7GWsuePXt47bXXOOuss5g9ezapqanxDktERETiwIkE/YPA5LzVQaDQWjsPuMeB55A4SNQSkN7UvzwU1lpeeeUV3n77bcaPH8+1116Ly+WKd1giIiISJ2El6MaYycCfgPNaP++YoAN4IohLeoloJ8sd14/k8GoiJvY+n49nn32WAwcOMHHiRC677DL1OBcREennwkrQrbW7jTGXAP8fsAD4NPDLQB3w88jDk0QW7U4vna0f7uHVRLyYqKGhgc2bN3P48GEuu+wyJk2aFNd4REREJDGEfVe4tbYBuAH4Z2vt2IA/46y1F1lrf+9cmJKIApNlay2lpaVRXz/cw6uBazU1NVFYWOhorD1VW1vLxo0b+fTTT5k+fbqScxEREWkXVoJujEk2xgyw1jZbax8JMmdxZKFJoot2p5fO1m87vPqb3/ymR7vgeXl5WGsBaGxsZPXq1ZSVlTkab6iqqqr429/+xsmTJykoKOALX/hCXOIQERGRxBTuIdGdwBhjzIX4D4R2xkawvvQC0e70Emz9cA6v5ubmsmrVKu6++24aGxsxxsSlt/vhw4fZvHkzLpeLhQsXMnSoLmyQxGaMeRaoBYqAv1prD8c5JBGRPi/cBPoIkAU0AM3AYfwJeZskYERkoUlv4HSnl44HOZ1cf8mSJaxdu5bk5OS49HYvLy/nueeeY9CgQRQUFDBo0KCYPr9ImKYCl1hr98c7EBGR/iLcQ6Jz2j42xsyw1r7YcY4x5ppIApP+J9oHOePZ2/2dd97h5ZdfZsSIEcyZMwe32x2z5xaJ0Ma25NwY8xXgZmAWUAI8bK39WzyDExHpiyIuQeksOW8dfyHStaV/CbdDS0/Eure7tZY33niD3bt3c8455zBjxgySk1X5Jb3K8bYPrLV/McZsBN4EFlprmwGMMVdYa1+OV4AiIn1NxJmCMWYl0AisB/KAHwOVwD/pLVHpiXAOnSZib/M2Pp+P7du3s2/fPiZMmMBVV12lHufSG11qjPlJh7FqYGXr/58N/h31q2IdmIhIX+XEVt4/A9cAE4FC/PXnxcD/Rb+wpYfuvPNOwF8v3l3CnYi9zds0NTWxZcsWKioqmDp1KpdccomSc+mtLgWmdTJ+eawDERHpL5xI0J/GfzHRo/iT8/+y1t5mjFnnwNrST3RMtpcsWdLtY0pLS2lpaaG+vp709PS4dGXpTF1dHcXFxZw4cYJrrrmG8847L94hiUTiM2Av0BLk60nARbELR0Sk73MiQT8T+ADIALYDdxhj5gBfA1aEs6AxZhlws7X2OmPMjfh3cEYBD1lrS1rndDouvVM49ec5OTnU19cDUF9fT05OTixC7VJ1dTWbNm2ivr6eOXPmMHr06HiHJBKpf7TWdnmzlzGm+1fUIiISMicS9KXAjfh3V36PP2E+A7g7nMWMMROAh4Gdxpix+GvbJ+Evo3ncGDMaGNLZuLW2OsLvReIknPrziooK3G43Xq8Xt9tNRUVFDCIN7ujRoxQXFwMwf/58RoxQp1HpE/7i0BwREQmREwn6CGAsMAy4snXMAJcA/9mThYwxA4BvACdahxa1xlgDeIFMYAYwOsi4/pHopcJpgZiXl4fL5YraTaY98dFHH7F161YyMjIoKCggKysrbrGIOMla6zP+AxT/BHwLOA9IBd7HvynzP9ZaXxxDFBHpc5xI0Ivx/8LuyHYy1p3bgV8C32z9/NzWv5sD5owHxgQZP40x5mb8fXtVbuCwYB1Uwu2s0tMWiPHsax7o/fffp7S0lGHDhpGfn096enpc4hCJBmNMMv5bRKfj33xpMwX/RswcY8w8a204v/Olj2srQ3RCcnIyqampjq0nkautrcXr9XY7r6qqKgbR9C1OJOjN+He6PQFjScCXe7KIMWY5sA3w4f9HwNW6zmlTOTUxDxw/jbX2YfwlM0ydOlX/gDgkWAeVWHdWiXVf80DWWnbt2sXOnTvJyclh1qxZpKSkxCUWkSi6Hf87lJ8Az+G/SboB/7um04B84F+BX3e3kDHm6/jfWc0A7rDWPhKlmCXOTEoaTelZbNr1jmNrDjQtzLv6CiXpCaK2tpZbfnQXR2rqQppfl5KBK02X9IXKiQT9IeCQtXZX4KAxZmAP1/k2/nryNqOBC1o/DkzU9wP1QcYlRoId6nTisqFE7m3epqWlhZdeeon33nuP3Nxcrr76apKSOns9KdLrLQW+a619qLMvGmPy8d9/0WWCboyZCJwD/Aj4n9b5StD7KJc7nSHLfoBtanBkvebaGk7+5X6am5uVoCcIr9fLkZo6Bt9wO2mZ2d3Od6W5SckYEIPI+gYnEvTfALaTHs+2h+vfCrQV7hbi36W5HvgQyAZS8Nemb2mdd18n4xIjwQ51hnPYM1BPduDjlcg3Nzezbds2Dh48yCWXXMLUqVPV41z6sveDJecA1tpiY8ziENZ521q7F9o7dTmTuUnCcrnTwa2Sv74uLTMbd/bQeIfR5ziRoDfhT6YDe+Qm4T88GrK2X9wAxpgGoNpae7C1neJtwEDga9ZaD+AJMi4xEqz+uyd14Z0l2KHuwLcl8j6fj5aWFjZu3MicOXOi880G8Hq9lJSUcPToUa688kouvPDCqD+nSJylGWPyCFJG2Dp+TneLWGutMcYNfAeYA9zpXIgiIn2LEwl6nrX29Y6DxpjLwl3QWjsm4OPHgMc6mdPpuMROsPrvUOrCg+2Uh7oDX1pais/naz+csmjRIvbs2RPVnfSamhqKioqoqalh5syZjB07NmrPJZJArsd/zsgJjfjvzfgE+JUx5lVr7cuBE3SwX0Sk80OYPdJZct46/mqka0vfFbhTbq2ltLQU+HwH/je/+U2X5S15eXm0tHz+pk1SUlL7GtFw/PhxnnzySerr65k3b56Sc+lvTDd/QmKtbbHWbgGWtQ4N6WTOw9baqdbaqcOHD484cBGR3ijiBN0Y4zbGrDLG/Efr52cZY66PODLp07raKc/NzWXFihVd7obn5uayceNG0tLSSE9PJykpKWp90A8dOsRTTz1FUlISCxcuZOTIkVF5HpEEdYe1Ngn/2R+3tTap4x/85YZdMsZcb4yZ2/rpB0AF/tunRUSkAydKXB4BZgMvAVhrPzHGnDDG/Mxa++8OrC8OSLTOKE70MJ8zZw579uyJ6ve1b98+XnjhBQYPHkxBQQEDBugEuvQ7640xRfh/zzcbYx601v5bhzm/D2Gd84F/N8Y8if/w/yydHRIR6ZwTCfoRa+0wY8yPAsbG4L91Tgl6Aoh1b/JQOdHDPFp90K217Nmzh9dee40zzzyT2bNnk5aW5vjziPQCv8J/qBP8XbO+b4zZZa1tT8qttSc6fWQAa+3PgJ9FJ0QRkb7FicbNZxhj/gGYZIy5zRizAVgHVDuwtjggWL13oisrK2PdunWUlZXF9Hmttbzyyiu89tprjBs3jrlz55KWlha3eETi7Drgn4GLgKnAg8CCwAnGGF0CICLiICd20NcAxcAoPj/40wT80IG1xQGR9iaPh3jt+vt8Pp577jnKy8uZOHEil112GcaYhH0XQiQGtltrHwz4fJcx5ncd5lwP/CWGMYmI9GkRJ+jW2neMMbnAfGAs/kuDngF0n2uCcKLeO9Y664feNt72PThdV9/Q0MDmzZs5fPgwl112GZMmTeoynt7wcxRxwLXGmGc7jI0PGEsCLkQJuoiIY8JK0I0x3V3PPAeoAVaEs750LZzENFq12tHScdc/JyfnlB3sDRs2sHjxYsd2tGtraykqKuKzzz5j+vTpfOELX+gynt7wLoSIQ0a1/uno7ICPbYxiERHpF8LdQb8phDnvh7m2dKG/lFp03PXvuIP96KOPOrajXVVVxaZNm2hsbKSgoIBRo0ad9iKoN74LIeKQEuBxTr0tOlASsDR24YiI9H3hJugPA3/Gv2vyLeBtYFfA12cD9ZGFJp3pT6UWHXf9jTGkp6fj8/nIy8vjiSeeiHhH+/Dhw2zevBmXy8XChQsZOnRo0BdBve1dCBGH/Lu19q2uJhhjdnX1dRER6ZlwE/QfWmtPAhhj7rTWnlLKYox5CX/Cfm+E8UkH8Sq1iHcf9dzcXDZs2MDChQtJSkriX/7lX9iwYQMVFRVhx3TgwAGeffZZBg0aREFBAYMGDQL614sgke50l5yHOkdEREIXVoLelpy3Gtd6icVuoAEYDMwAhkUennQUj1KLRCmrqaioIDk5uT1xrqioYMWK8I45vPPOO7z88suMGDGCOXPm4HZ/fqZZ9eYiIiIST060WbwVeIzPL7IAf5vFmxxYWzoR7VKLjrvlhYWFNDQ00NzcjNvtjmhHOZSd+GBznEicrbW88cYb7N69m3POOYcZM2aQnHzqfwaqNxcREZF4cqLNYokxZiwwD8gBjgHPWGuPRLq2xF7H3fINGzZw77330tzcDIDX6yUnJ6fTx4WSeHe3E9/VnEgT55aWFl544QX27dvHhAkTuPLKK0lK6vx+FdWbi4iISLw4sYOOtbYK+IMTa0l8ddYtpaXl8+YNKSkpVFRUnPKYUEtgQqnt7m5OuIlzU1MTW7ZsoaKigqlTp3LJJZdgjOnxOiIiIiLRpuuZ5RQdy0iWLVuGy+Vq/3pycvJppSWBSbW1tv1Soe7W7qxEJRr13/X19Tz99NMcOnSIq6++milTpig5FxERkYTlyA669B2dlZHs3r2bwsJCAJYsWQLA6tWr2z8PNakOXDsnJ6c9ke+4Q+5k/Xd1dTWbNm2ivr6eOXPmMHr06IjWExEREYk2Jej9UHf14h3LSHJzc1m1alX7YydPnkx9vb/N/erVq9m9e3enSXVnz9P2d1clMU7Vfx89epTi4mIA5s+fz4gRIyJeU0QkUGVlJR6Pp9t5mZmZDBkyJAYRSV9SW1uL1+vtdp7b7WbAgAExfd6qqirHnk9OF1aCboy5HnjKWuvrMHYu8La1tsih+MRhkbZMLC0txedr/58dn89HaWkpK1asOGWdrp4nFn3GP/roI7Zu3Up6ejpz584lKyvL0fVFRCorKxk3bhzV1dXdzs3KyqK8vFxJuoSstraWW350F0dq6rqdO3JQBg/8/F5HkvSePG9dSgauNHe386Tnwt1B/wvwFnCJMcYFbAJmtn5ttzFmBXBTh37pkgAiTY7z8vJOqUl3uVydlrR09TzR7jP+/vvvU1paytChQykoKCA9PT3steJ9QZOIJC6Px0N1dTXLly8nMzOzy3nr16/H4/EoQZeQeb1ejtTUMfiG20nLzA46r8FTxZE//gKv1+tIgh7q8wK40tykZDi3cy+fi6TEpW1L8lZgFv5Ligqstc8bYwYC/wT8MsL4xGGRJsed1aR3lrh29TzR6jNureXNN99kx44d5OTkMGvWLFJSUsJeL1EuaBKRxJaZmUl2dteJjEi40jKzcWcP7TfPK35O1KB/A7DA/2utfR78N40aYwpQgp5wcnNz2bBhA48++ijLli0LK+EMrEnvak5XSXhbnXlZWRnr1q2LOFFvaWnhpZde4r333uPcc8/lmmuuCdrjPFSxKMURERER6SjcBL0ZcBl/r7rzW8f+1PZFY8xNwOWRhSadibTkoqysjMWLF2Ot5YknnojqrnB3SbhTO9TNzc1s27aNgwcPMnnyZL70pS850kYx2qU4IiIiIp0JN0EfDiwE/jdgjQMBX78deDuCuKQTTiS0Tu4KR3p7qBOxeL1eSkpKOHr0KFdeeSUXXnhhWN9LZ6JViiMiIiLSlbASdGttNfB74PfGmH8F5gNHA6bcARyMPDwJ5ERC69SusBO3h0YaS01NDUVFRdTU1DBz5kzGjh0b1vfSFadaPoqIiIiEKuIa9NZOLX/uMPZMpOvK6ZxKru+8804g+AHPUHT3YqFtdz0nJycqh0VPnDhBUVERPp+PuXPncuaZZ4b1fYiIiIgkmqhdVGSMmW+tfTpa6/dHkZZcdNz1brsVNBxdvVjo+DwbNmygoqKiy8OiPXHo0CE2b95MamoqCxcuVPcEERER6VPCvajofqCrxpfJwLWAEnSHRVJy4WT9eduLhcLCQo4fP05hYWF7wr9mzRp8Ph9erxe3282jjz7KypUrHSkV+fDDD3n++ecZPHgwBQUFjt6cJiIiIpIIIuniclM3c2yYa0uURKMryerVq6mvrwfgnnvuaW9t6PV6SU1Nxev18vjjj5/SMSacTjTWWvbs2cNrr73GmWeeyezZs0lLS4s4fhER6T3a/r3pSnJyMqmpqTGIJrFUVVV1O8ftdmtjq5cIN0H/JbAIWAs0dvL1FGBZuEFJdDjdlaS0tBSfz9f+eXNzMy6Xi8bGRtLT07nkkkt48803qa+vb9+xB3rcicZay6uvvsrevXsZN24c11133Sm3mYqISN9mUtJoSs9i0653up070LQw7+or+k2S7kpzU5eSwT/f84tu544clMEDP79XSXovEG4Xl0PGmHuttY8Em2OM+TD8sCRanOxKkpeXd0qinJycTFJSUvsO/U9+8hMWL158yo59T8tsfD4fzz33HOXl5Vx00UVcfvnljvQ4FxGR3sPlTmfIsh9gmxq6nNdcW8PJv9xPc3Nzv0nQUzIG8IVb7sXX4O1yXoOniiN//AVer1cJei8Q9iFRa+0jxpiBQKO19rRd9LZbRaV366ocJTc3l927d1NYWAjA1KlT2bFjB/B5h5i2HfucnJxuu7p01NDQwObNmzl8+DCXXnopkyZNUnIuItJPudzp4E6PdxgJKSVjACkZSrr7knAPiQ4F/gDMBnzGmAestf/maGQSd6H0Os/NzWXVqlVBO8S0ze/Y1aUtkQ+mtraWoqIiPvvsM6677jrOPffc6HyTIiIiIgkmKczH/QqYAxj8Sf73jTHfdCwqiYqysjLWrVtHWVlZSPMDy1Gste015D2dW1hYSFNTU/vXduzYwdq1a1mzZg1Tpkw5LZ6qqiqefPJJampqKCgoUHIuIiIi/Uq4JS7XAf8MPA+4geXAAvy3iwJgjEmy1rZEGqA4I9SbPwN11+s8sPQl2NyysjJWr15NY6O/CqqpqYnjx48HrUM/cuQIJSUluFwuFixYwLBhw6L0ExERERFJTOEm6NuttQ8GfL7LGPO7DnOuB/4S5voSoY4JdDg90HNzc9mwYQOPPvooy5Yta58fLNnvrENMaWnpaXXjDz30EAAZGRlYazly5AhlZWWkpKTw7LPPMnDgQObOncugQYOi8JMRERERSWzhJujXGmOe7TA2PmAsCbiQHiboxpivA/8JZAB3tB5EvRG4FBgFPGStLWmd2+l4X9TTvuGdJdDh9EAvKytj8eLFWGtP6WMeLNnvrENM2/OmpqbS2NhIU1MTKSkprFy5EvD3UV+zZg2rV69m5cqVTJw4kTlz5uB2u8P7YYmIiIj0cuEm6KNa/3R0dsDHPbqoyBgzETgH+BHwP8CvjTHPAeuBScA1wOPGmNHAkM7GrbXVPf1GEl04pSmdJdArVqzocQ/0YIl4qMl+2wuL//7v/+bpp5+muLiYpKSk9kOkbbvrdXV1pKamcuLECebNm0d5ebljvdpFREREeptwE/QS4HEgWI15ErC0h2u+ba3dC2CMWQY04L8MKRmoAbxAJjADGB1kvM+V1IRTmhIsge64w93dznxn67Q9ZsOGDVRUVAR9bNsLC5/Ph9frJS0tjZaWFv7pn/6J7373u+Tm5tLS0kJLSwupqakkJSXxj//4j5SXl/f4BYmIiIhIXxJugv7v1tq3uppgjNnVkwWttdYY4wa+g79DzJ1AW/uO5oCp44ExQcb7nHBKU0K5MTTUFoq7du1q73N+4MABvvzlL+Pz+XC5XOzevTto8tz2wsLr9V+c0NDgkrHUGAAAIABJREFUv1zit7/9Ld/97ndpamqivLyclStX0tjYyI033sh5553HunXrevyCRER6B2PMD4F/x78B8w/W2k1xDknEMbW1te3/5nXH7XbrsiDpUrg3iXaZnIc6pxONwAfAJ/hbOW7tZI7h1MQ8cPz0QWNuBm4GGD16dBghxVcoyXawxwWbW1ZWxpo1a2hpaaG+vr7LRPjAgQPcfffdJCUl4fP5aGpqav9aYWEhq1at6vQ52l5YuN3uU35hJSUlsW3bNkaOHMnx48dZsmQJ559//mmP68kLEhFJfMaYPGAK/nc6v4O/NHGYtTa0jEYkgdXW1nLLj+7iSE1dSPNHDsrggZ/fqyRdggr7JtFoaG3LuKW1xOUVYCMwk1P7te8H6ls/7jje2ZoPAw8DTJ06tUd18YkisCNK4Oeh6FjG0rH0JD09PWgiXFZWxqJFi9pbJLpcrlO+fvz48aDPE/jCAuCWW25prz/3er1UVVUxe/ZszjnnnNO+13BekIhIwvvQWvsNAGNMDjAd8MU3JBFneL1ejtTUMfiG20nLzO5yboOniiN//AVer1cJugSVMAm6MeZ6oLH1Lc8PgAr8Cfp9QDaQApwAtgBZQcb7pHAOigZ7XGDpidvtZunSpaxcubLT9UpLS0lK6vouq3Xr1pGTk9Pe7SUwvsBd/Ly8PJ555hmampoYOnQo+fn5jBgxotM1u9r9F5HeyVp7OODTFGC9tbap47ze/q6n9G9pmdm4s4fGOwzpAyJK0I0xA/En1Y0OxHI+8O/GmCeBD4FZ1tqDre0UbwMGAl+z1noAT5DxPimcg6IdH5eamkphYSFLliw5pYQkWHIO/qQ6KSmJ9PR0WlpaeOCBB7j11lvx+XwYY3jooYcwxtDc3IzL5eqyXMbtdpOVlUV6ejpz584lKyvLkZ+NiPQuxphzgMHA7Z19vS+86ykiEqmwEnRjzFDgD8BswGeMecBa+2+RBGKt/Rnws07GHwMeC3W8Lwq3LjsvLw9r/f++NTY2snr1apYsWdK+k56Tk9Nl2Uxn5SZ5eXmUlpZy5MgR1qxZQ11dHW63m5aWlvb4cnJyWLduXftjPvjgA7Zv396+c56RkeHcD0dEeg1jTBL+xHy+tfakMeYsa+0n8Y5LRCTRhLuD/iv8nVba1vi+MWaXtfb3zoTVO/T0AqFwRXJQdNWqVdx99900NjZijGnviQ6EVDbTsdyk7fOSkhKam5txu924XK72tosdy11+//vfc+zYMUaNGsXs2bNJSUlx5ociIr3Rz/B35/qVMWYI8Fv8TQFE/v/27jw+6ure//jrQ/YQEkCQTUQ2KaKCgCheWdSqRIQqFhA3ikXQul7p9def1KVXRCtKvfV6r9D+tKitRRARqIkWEEWqlAoUFQsKIqCgLIEQsofz++M7iZMwk3XILHk/H495TObM+X7ncyZw8pkz53uOiPipb4J+EXAHsApIBiYBo4CKBN3Mmvku+oxJ9Z0XXpfzV73osqbzB/rAMHbsWB577DHi4+MrrWUeaBUXIOCHgEAXml5zzTXExcVx7NgxXnvtNS6/3Pu85r9MYnJyMq+//jqTJk1i6NChx11kGi6N9cFKRL5nZhM5flrLlHDEIiIS6eqboL/nnPsfv8frzez3VepcRQxuHFSe3O3du/eErdddXfIfLLkMdkzV0Xcg4Coup5xySsDjq7vQtDy53717d0Uc5dNxyjcmGj16NMOHD8cs4CqYje5Ef7ASkcCcc/OAeeGOQ0QkGtQ3QR9uZiurlHX3K2sG9CHGEnT/5K58bveJWK872EWhW7dupV+/fgE3ClqwYAElJSUUFxcf94HBf/S9fIS76iouwV4zUHl1c+JPPfVUnnjiCdauXcu4ceMYOXJkyN6XUKjvBbciIiIijaW+CXon362qzn4/x9zV91WTu/vvv5/27duHfKpEsAR4wYIFFBQUVNQr3yho69atPProoxXrlTvngn5g8D+3c46ePXtWKk9JSaGsrIxTTjklaCzB5sQfOXKErKwskpKSeOihh+jWrVvI3pNQ0UZIIiIiEunqm6C/BbwKBJtj3gwYV89zR6yqyd3YsWNPyOirfwLsv9JKMKtXr66YQhIXF8fUqVMrRtyrJtHl516wYAGPPvooM2fO5LHHHmP9+vW89tprjB49mmbNmnHNNddUTP8IlIxXnRN/4MABsrKyKCsr44orrqBDhw4hf19CQRshiUgsKC4uprQ00Kba3/Mf0BEpl5OT06DnpXHUN0H/hXPun9VVMLP19Tx3xKpNcheqCxDLj/WfL/3aa69VjHDHxcUxduxYoPJyimVlZcyZM4cRI0YE3Dyo/Nzt27fHzCpN9QCIj48/bvpHTReofv3117z99tskJiYyevRoWrWqfhe1cNNGSCISzYqLi/nLe38jz1W/kRxASUoGlpDUCFFFJ/9rqIJJT0+ndevWjRDNiRWXlEx+Qip3PDKrxrr5CanEJSU3QlQSTH0T9I9DVCfqVJfchfoCxKpTanbv3s3GjRsDjmZXXU7xlVdeqXaudbCpHnWd/vHFF1+watUqMjIyyMzMJC0trd7tFRGRmpWWlpLnmpE85k7im7eotq4lJBGXnNJIkUWPpKQkmjdvXqu/cxkZGWzfvj3qk/SE1Ob0uH0GZUWFNdaNS0omIbV5I0QlwdQ3QR8DLAxBnYjTkBHwUFyA6P/6weZ/Bzpn1eUUJ0yYwMKFCwMm2+WvUb52uX9b6zL9Y9OmTXz44Yd06NCByy67jKQkjdKIiDSW+OYtiG/RMtxhRKXU1FRuv/12ioqKqq2Xm5vLCy+8QG5ubtQn6OAl6Uq8o0N9E/S5ZnYH1c9BP4soS9AbOgLe0AsQA71+dQlz1Q8TVeuWzzWvSxv9PwD4nx++XyO9Z8+efPjhh3z88cd07dqViy++OGLWOBcREamN1NRU7WwtEau+CXpLYAhQ3eLWUbeKy+rVqyvWB09OTq7zCHhDL0AMNAL/05/+NOhc92Drnvt77LHHcM5VXAha21H+QEtKmhlmxrPPPktRURF9+vThggsuiJg1zkVERERiQX0T9O+Af+HtJBooEW8GZNbz3GFzyimnUFjozc0qLCysWGqwLhpyAWJdRuBrk2hXt4Z51eUUqzt/YmIi4F2YlJSUxIoVK5g2bRpnn322knMRERGREKtvgn4KMAoYBuwEXnbOfedfwcyWNzC2Rrd7925SUlIoKCggJSWlVld3h1JdRuBrk8wHm8MebDnFYMeWj6CXzzG/4YYb6Nu3bwhbLiIiIiLl6pWgO+dKgdeB182sIzDJzNoAK4Fs56l+8e4INGTIEJo1axbWTWxqOwJfm2Q+WJ3du3cHXE4x0LELFixg9+7d7Nixg5SUFO655x6GDh0amsaKiIiIyHHqO4JewTn3jZk9CVwN/A4oM7PngRecczsbev7GFG2b2NQmmQ9Up7ZTab788kt+9atfUVJSAkBycjKPP/54aIIXERERkYAalKCb2anALcDNQHu+v2h0HLAfeLZB0YVBU9jEprYbLo0ePboiOQc4duxYvZaOFBEREZHaq1eCbmajgKnA5XgXhBpQhLes4hzn3PshizBChWrH0HAJ9EHEv03z58+vmHteLi4uLizTfkRERESakvqOoL+Bt3qLAVuAucA859zB8gpm1s05t73hIUaeUO4YGimJftVlFW+55Rbi4uKIi4ujrKyM22+/ndtuuy0qP4yIiESCmhYeaOyFCWJZQUFBjXXi4+MrVikLhZKSEvLyjlASX/05i/KOUFJSQk5OTrX1anpeYltDprh8AMwB3sNL1tPMrHyf93jgV8CNDQsvMoVix1AIbaLfUFWXVUxOTmbDhg2sWbMm7B8eRESiWV22lW/evLl2ZW4AS0iiJCWDN9d/WmPdNDvGyKEXhCRJz8nJ4b333qW05WuQklZ95ZIi4jZ+ytQHZ5KQkFBt1fyEVOKSkhscn0Sf+ibo2cACvOktwzh+w6JmQM8GxBXRGrpjaLlQJfqhcP7553Ps2DESExOJi4vj5ptvplevXvzgBz8ISzwiIrGittvKg5fMa3fL+otLTqH1hHtxJdW/16VHj5C36BlKS0tDkqAfOXKE0pJSuvfoQUqrttXWLS4pZivQcvRo0tJaVFs3LimZhNTmDY5Pok99E/TpzrkNVQvNrAcQ55zbYmYfNyy0yBWq1V5Cleg3VEFBAVu3bmX69OkAjBs3TiPmIiIhpG3lG09ccgokp4TltRMSEmv3DUhCEkktTyK5RfUJujRd9U3Qh5rZbcAOvPnnR/BG1EcCmNmHwFWhCDBShWK1l0hY1vHw4cNkZWWRn5/PjTfeSJcuXRo9BhERERH5Xn0T9DhgF/CYc86Z2X8BV+LNRX8WyAF+jbf8olQjnMs67tu3j+zsbJxzjBw5knbt2oUlDhERERH5XrN6HtfVOfeILzk/HbgNLzmf5Zy7yzn3EHAgZFFKyO3atYulS5cSHx/P6NGjlZxLo1u5ciUdOnTAzJgzZw4A+fn5zJ49m9GjR/P111/XeI7du3dz3333ce2111aU7d+/n6uvvppf/epXDB06lH/961+Vjvnkk0947rnnePfdd7VKgoiIRKT6Juj+V188jjcSvw+Y4Vfesr5ByYm1ZcsWsrOzadmyJT/60Y9o2VK/Kml8F198MTfccANxcXHcddddrFmzhtTUVO6++26uvvpqOnXqVO3xe/fuZfHixcyaNYtDhw5VlE+bNo1//vOfPPTQQ5x88smMHTu24rn58+czc+ZMpk6dyrBhw2jVqtUJa5+IiEh91TdB72Zm15rZbLy55g643zmXB94a6MDY6k4gjc85x4YNG3j33Xfp2LEjV155pS5akrDq06cPTz31FMXFxVxzzTXs3r2buLg4zLyFofbs2UOHDh148MEHjzu2ffv23HHHHceV//3vf69YLaNHjx588sknbN68me3btzNp0iQKCgoYOnQoL7744oltnIiISD3VN0GfDjwK3IOXnM9wzj1vZm3MbCqwFtClyRHEOceaNWtYt24dPXv2ZMSIESHdoEGkvu6++25uvfVWvv32W8aMGUNhYWHFc8XFxeTm5lYaIa9J79692bNnD19++WXFB9D8/Hx+97vfUVBQwOuvv86ZZ57JpEmT+OSTT0LeHhERkYaq10WizrnPzOwHQB9gv3OufPuzjsC/gB+HKD4JgdLSUlauXMmOHTvo27cvgwYNqhihFIkEzzzzDF988QXLly9n6tSpXHTRRQB06dKFgwcP1mnjlv/6r/+iWbNmTJ06la+++orExER69uzJxx97K7+WlZVxySWX8Nxzz/Hee+9x5plnnpA2iYiI1Fe9dxJ1zpUAG6uUbWpwRBJSRUVFZGdn891333HBBRcoGZGIFB8fz4IFCxg8eDAvvvhiRYJeH507d2bhwoUUFhbSrl07rr32WjIyMkhO9nbj27dvX8Xc85p28RMRAW+n0BbVrFl++PDhRoxGmoJ6J+gS+Y4cOUJWVhZHjhzhkksuoVu3buEOSaSSkpKSip9btmzJsmXLOO+88yrKdu7cyRlnnMHkyZN5+umnjzveOVfp3t9DDz1Ely5dKo676qqreO2119iyZQv5+fkkJiZy2WWXhbpJIhJDSkpKIK0l9zz2G1JSgm9+VFBQABltaJagqaMSGkrQY9SBAwfIysqitLSUK664gg4dOoQ7JJFK3n//fZ555hkyMjIYN24cAN27d2fRokUVSyzGx8eTlpZGWlracccfPHiQRYsWAfDpp5/yzjvvcNFFF7Fu3Tpef/11MjIyWLt2bcUf1euvv57PPvuMxx9/nHbt2jFv3jxtzCUi1SotLYXUdNIn3EvrDp2D1svLOwJLltAsKTw7mErsUYIeg77++mv++te/kpCQwOjRo2ndunW4QxI5zoUXXsimTcfPihs6dGjFzx07dmTv3r0Bj2/dujWTJ09m8uTJlcrPPfdczj333OPqmxmPPvpoA6MWkaYosUVLkludFPT5kvhESKj9tTIiNVGCHmO2bdvGO++8Q0ZGBpmZmQFHHkVEREQkcilBjyGbNm3iww8/pH379lx++eV1WvlCRERERCKDEvQY4Jzjww8/5OOPP6Zr165cfPHFxMXFhTssEREREakHJehRrqysjFWrVrFt2zb69OnDBRdcoDXORURERKJYfXcSlQhQXFxMVlYW27ZtY9CgQUrOJaps3LiRs88+GzPjtttu4/7776d3797ceuutFBYWUlxczEsvvURcXBxmxvPPP1/p+IkTJ2Jm9O3bl48++gjwNjy65557mDFjBvfccw8zZ84EYPHixaSmpmJmzJgxgwcffJDBgwdz1VVX1SrW3bt3c99993HttddWlBUXF3PjjTeSlpZG//792b59+3HHTZkyhYcffrie75CIiDRVEZWgm9nPzWy/mX1tZlf4ym4ws2fMbJGZXe5XN2B5U3H06FGWLl3K3r17ueiii+jXr5+Sc4kq/fr1o3///gCMHz+emTNncvPNNzNnzhz+8Ic/kJiYyI033sioUaMAmDt3bsWxhw4d4tixYwCcc845DBgwgI8++ogHHniA3/zmN/zyl7/krrvu4tNPPwW8NdBPPvlkAH75y18yZcoUPvjgAy688MIa49y7dy+LFy9m1qxZHDp0qKL8mWeeoU2bNlxwwQVs2LCBW265pdJxL7/8Mr/73e8a8A6JiEhTFTEJupkNAfoDi4AOwKtm1ht4Afgf4G1fWYaZdQ1UHp7IG19OTg5vvPEGubm5jBgxgp49e4Y7JJGQ2LFjBwCdOnWqKOvXrx8XXngha9eurViWcenSpYwdO7bSsd999x2HDx/m+uuvZ8+ePXTr1o2f/OQnAV9nzpw5APz7v/87AHv27KFDhw48+OCDx9Vt3749d9xxx3Hl1113Hb/5zW/4y1/+QkZG5e7n008/ZcuWLbVrdBNjZs3N7E4z2xzuWEREIlXEJOjAF86565xzU4BsvPnxV/rujwCFQDpwCfCjIOUx79tvv2XJkiWUlZUxatQoTjnllHCHJNJg8+bNY8CAAbz44ossXry4YtS83J133gl8P4peXFxMenp6pToXXXQR//Zv/8Yrr7xC165dmTZtGueff/5xr3XzzTfz3//93wAVF1MXFxeTm5tbaYS8Jv6bfxUXFzNlyhQA8vLy+POf/8xPf/rTWp+rqTCzRGAkMB3Q1sYSdgUFBeTn5we9FRUVhTvEJiknJ4cDBw5Uezt69Gi4wzyhIuYiUefcHr+HCXgj5Kf5Hpf6Pde9mvLjmNkUYArAqaeeGoJIw2fHjh2sWLGCtLQ0MjMzj0tQRKLVTTfdxPTp08nLy/N27qtizJgxdOrUiZdffpkbbriBQYMGceDAgUp1SkpKWL58ObNnz+bJJ59k9uzZrFy5knXr1hEf/31X9/zzz1NWVlbp2C5dunDw4MF6LU366quvct111zF+/HgAHn30Ue68886KZL+4uJiysjKtrAQ454rxvvH8GXD8pyeRRmIJSZSkZPDm+k+rrVdSUgIZbYhLSm6kyJq2uKRk8hNSueORWTXWbd8ilWd/PYPmzZs3QmSNL2IS9HJm1gVoCfwHEGjbP6NyYu5ffhzn3FxgLsDAgQNdiMJsdJs3b2bNmjW0bduWyy+/vGL7cpFYYGb89re/5bzzzuPee+/liiuuqPRvPD4+nltvvZUHHniA5557jj/84Q+sWrWq0jneeOMNMjMzuf/++7n99tu5+eabWbRoEZs2baqY617uxz/+MeDNZd+3b1+9p4nt3buXtWvXMmfOHAoKCigqKuLxxx/n8ccfr6jz2GOPcdlllzF8+PB6vUZTE0uDKhK54pJTaD3hXlxJ9SPkRcXF8PnnxKfEZhIYaRJSm9Pj9hmUFRVWW68oN4e9f5xFYWFhzCbokTTFBTNrhpeYX+mcywMO+p7yj3Mb8EWQ8pjjnGPdunW8//77dO7cmZEjRyo5l5jjnGPgwIFMmjSJnTt3Vqy+4m/KlCkkJiaSmZkZ8BylpaUVc8gzMjKYPn06LVq0oFu3yjMpnHMVU2geffRRUlNT2blzJyeddBL33HNP0Pj878EbGZ8wYQL79u1j4sSJDB06lGbNmrF69WpWr17NggULAJg0aRLnnHNOHd+Rpss5N9c5N9A5N7Bt27bhDkdiWFxyCvEtWlZ/S8uABG3615gSUpuT3Oqkam9J6a3CHeYJF2kj6I8DPYGnzKw18AreaHkrvGkvB4C/AhnA7ADlMeXYsWOsXr2aLVu20KtXL4YMGUKzZhH1mUqk3jZu3Mj69esBbxnEs846i5kzZ7Jw4UJmzZrFueeey44dO1i6dCljxozh7LPP5u6772bMmDHs27ePJUuWALBhwwY++ugjEhISiI+P57777qN9+/Zs3ryZpUuX0rJlS5YsWcJ3330HwN13302rVq1Ys2YNa9as4YknnmDPnj2kpaWRlpZ2XJwHDx5k0aJFgHfx5zvvvMNFF13E1KlTK43id+/enfT09IqVYcoveD311FOPu4hURESkOhGToJvZRLzRc39TgBuAaUAaMN45lwvkmlmg8phRPp92165dDBgwgP79+2sZRYkp/fr1q1iVxV/VCzXvuuuuip+feOIJANq2bcvs2bOZPXt2xXMDBgzg+uuvD/hao0ePJj8/P2gsHTt2ZO/evQGfa926NZMnT2by5MmVyl944QVeeOGFoOc87bTTKo24SyVGkGmJIiISQQm6c24eMC/AU/N9t6r1A5bHgoKCArKzs9m/fz9Dhgyhd+/e4Q5JRKTBzCwOuBQ4A0g0s5uAV51z1U84FRFpYiImQRdPbm4uWVlZHD16lMsuu4wuXbqEOyQRkZBwzpXhLaOryeUiItVQgh5B9u3bR3Z2Ns45Ro4cSbt27cIdkoiIiIg0MiXoEWLXrl0sX76c5ORkMjMzadmyZbhDEhEREZEwUIIeAbZu3cp7771Hq1atyMzMJDU1NdwhiYiIiEiYaM2+MHLOsWHDBlatWkWHDh0YNWqUknMRPx988AGDBw9m2LBhTJgwwdvVT0REJMZpBD1MnHOsWbOGzZs306NHD4YNG6atwEWq6NKlCytXriQlJYXp06fzxhtvVOwCKlLu4MGD5OaGbqXd9PR0WrduHbLzSfQ7ejSvQc+L1JUS9DAoLS1l5cqV7Nixg759+zJo0CCtcS4SQMeOHSt+jo+P10ZdcpyDBw/SrVs3Dh8+HLJzZmRksH37diXpQlxcHPHxCRUbo1UnPj5BA20SMkrQG1lRURHZ2dl89913DB48mLPOOivcIYlEvC+//JKsrCymT58e7lAkwuTm5nL48GEmTZpEenp6SM73wgsvkJubqwRdiI+Pp0+fPpSVldVY10vmlVZJaGg4qhHl5eWxZMkS9u/fz8UXX6zkXJq8jRs3cvbZZ2Nm3Hbbbdx///307t2bW2+9lcJCb++a3NxcJk6cyEsvvURiYmKtzvvyyy9z5513MmbMGN56663jnt++fTtmVun285//HIAnn3ySNm3a0KlTJ958883QNVZOqPT0dFq1atXgWyiSfIkt8fHxJCUl1XhTci6hpH9NjeTAgQNkZ2dTUlJCZmZmpa/uRZqqfv360b9/fz7++GPGjx/P8OHDadWqFffddx/9+vVj8uTJTJgwgYcffphevXrV6pxffvklkyZNYtOmTbz77ruMGzeOnTt3kpGRUVHnb3/7G3379qVnz57ExcWxYsUKxo4dy+rVq1m/fj1jxozh97//PePGjWP//v0kJyefqLdARETkOBpBbwTffPMNS5cuBWD06NFKzkWqsWPHDgA6derEK6+8wtq1a/nP//xPhg8fzvz58yvqDR06lGHDhh13/BtvvEFpaSktWrQgOTmZ3NxcVqxYUalO37592bhxIwsWLOD555/n5JNP5rzzzqNHjx786U9/Yu7cuYwYMYLS0lLNKRURkUanEfQTbNu2baxatYr09HQyMzNJS0sLd0giEWnevHlMmzaNrVu3snjxYkaNGgXAjTfeGLD+wYMHA140+vnnnwNU+rp527Ztler4Ty9btmwZmZmZAHTo0KGivKSkhEmTJpGQkFDPFomIiNSPRtBPoI8//pgVK1bQtm1bRo8ereRcpBo33XQTSUlJ5OXlUVpaWmP9jRs3sn79+uPKA80Ddc4FPc+f//xnxo0bV6nsq6++4tChQ8yaNasWkYuIiISWEvQTwDnHhx9+yAcffEDXrl0ZOXIkSUlJ4Q5LJKKZGb/97W9p1qwZ9957LwUFBTXWD7Q8aY8ePQA4duxYRVn37t0BbxUlf0eOHGHz5s0MGjSoouzYsWPMmjWLZcuWkZaWxjfffFPvNomIiNSHEvQQKysrY+XKlWzatIk+ffrwwx/+UHNYRWrBOcfAgQOZNGkSO3fuZObMmdXWP+eccxgwYMBx5VdddRXx8fHk5ORQUlLCSSedxKWXXsprr71GamoqixYtqqjrP5Wm3C9+8Qs+//xzpk2bxhVXXMFHH30UmgaKiIjUkuagh1BxcTFvv/0233zzDYMGDaJv377agEikGv7TVBYvXsxZZ53FzJkzWbhwIbNmzWLw4MFcccUVAY9NT08POAe9c+fOvPzyyzz11FPk5eUxf/580tPTSU5Opnnz5qSkpFTUnT9/Pg8++GDF43nz5h03rWXu3LmhaKqIiEitKUEPkaNHj5KdnU1OTg7Dhw/n9NNPD3dIIhGvX79+bNq06bjyQ4cO1Xjs+++/H/S58ePHM378+EplI0eOPG47+GXLllV6PHHiRCZOnFjja4uINIaiokKKi0tqrJeYmEBSUtNbDjYnJ6fGOuWDM9FGCXoIHDp0iKysLAoLC7n88svp3LlzuEMSERGRKFZUVMgrr7xCcXFxjXUTExOZMGFCk0nS45KSyU9I5Y5Har6Qv32LVJ799YyoS9KVoDfQt99+S3Z2Ns2aNePKK6+kbdu24Q5JREREolxxcQnFxcVFMNWYAAASXUlEQVSc3qsXiQnBd1EuLilm65YtFBeXNJkEPSG1OT1un0FZUWG19Ypyc9j7x1kUFhYqQW9KduzYwYoVK0hLSyMzM1NbRIuIiEhIJSYkaiW4ABJSm5OQGl1Jd10oQa+nzz77jPfff582bdowYsSISheeiYiIiIjUlxL0OnLO8Y9//IMNGzbQuXNnfvjDH2qnQREREREJGSXodXDs2DFWr17Nli1b6NWrF0OGDAm4zFuobd26ldWrVzNkyBCtDiMiIiIS45Sg11JJSQnLly9n165d9O/fnwEDBjTKGudbt26lf//+OOcwM9avX68kXURERCSGaSfRWigoKGDZsmXs3r2bIUOGMHDgwEbbgGj16tU458jPz8c5x+rVqxvldUVEREQkPDSCXoPc3FyysrLIy8vj0ksv5bTTTmvU1x8yZAhmRmpqKmbGkCFDGvX1RURERKRxKUGvxr59+8jOzubYsWNceeWVtGvXrtFjOP3001m/fr3moIuIiIg0EUrQg9i1axfLly8nOTmZzMxMWrZsGbZYTj/9dCXmIiLVKCgoqHGt6Pj4eBITg2/44m/37t0hqSMSzNGjeQ16XmovJyenxjrJyckRtZmREvQAtm7dynvvvUerVq0YMWJERP3CRETke/n5+ZDWkuyNn9W45G2aHWPk0AuqTdKTkpJo3rx5racTNm/eXJvISJ3ExcURH5/AkiVLaqwbH59AXFxcI0QVm+KSkslPSOWOR2bVWLd9i1Se/fWMiMn5lKD7cc6xceNG1q1bR6dOnbj00ktrPdoiIiKNr6ioCFLTSbr6DlJbnRS0XunRI+QteobS0tJq+/XU1FRuv/1277y1kJSURGpqap3jlqYrPj6ePn36UFZWVmNdL5lXqlZfCanN6XH7DMqKCqutV5Sbw94/zqKwsFAJeqRxzrFmzRo2b95M9+7dGT58uD61iohEibjmLYhvEZqpiKmpqUq65YSKj49X4t1IElKbk5AaGUl3XehfB1BaWso777zDl19+ydlnn815553XaMsoioiIiIj4a/IJelFREW+99RZ79+5l8ODBnHXWWeEOSURERESasCadoOfl5ZGVlUVubi6XXHIJ3bt3D3dIIiIiItLENdkE/eDBg2RlZVFSUkJmZiYdO3YMd0giIiIiIk0zQf/mm294++23SUhIYNSoUZx0UvAr/0VEREREGlPEJehm1hy4GbjNOXeGr+wG4DygEzDHOfdWdeXV2bZtG6tWrSI9PZ3MzEzS0tJOVFNERMRPffpsEZGmKKISdDNLBEYC04GWvrKuwAvA2cAw4FUzOxVoHajcOXc42Pk7duzIypUrOfnkkxkxYoQ2lxARaSTB+vLq+mwRkaYqohJ051wxXqf9M+B8X/GP8OI8AhQC6cAlwKlByhcFO3/37t3p0qULF198sdYfFRFpXMH68qB9dn5+Pi+99FK1J922bRsAhTn7qq1XlneYkpISvvrqq6gfnCkqKqKkpIT8g/uIKy4OdzgSZiXFRVCQx66tn5HaXOv310dxbg5HDx5k/vz5tGjRotq6Z555Juecc84JjykastSevvtSv7LuwGlByisxsynAFN/DoqeffvqTUAcYJm2A/eEOIgRipR2gtkSiWGkHQK9wB9BAwfrySqr22TfddFOt+uwdD02sVRA7/16raidC6P8t/n11SE9XS7H0fyqm2vL+/8ZOWwjT72XdsoWhPF2D+uxoSNBLA5RZNeWVOOfmAnMBzOwfzrmBoQ0vPGKlLbHSDlBbIlGstAO8toQ7hgZSnx0DbYmVdoDaEqlipS0N7bOjIUH/wnffzK9sG1AQpFxERCJPsL5cRESqaFZzlbAwvh9ZWYw38tIKSAAOAH+tplxERCKP+mwRkVqKqBF0M4sDLgXOABLN7CbgVeAGYBqQBox3zuUCub4lu6qWV2fuCQu+8cVKW2KlHaC2RKJYaQdEeVucc7vUZ8eEWGkHqC2RKlba0qB2mHMuVIGIiIiIiEgDReoUFxERERGRJkkJuoiIiIhIBIm5BN3MmpvZnWa22a/sBjN7xswWmdnlNZVHAjP7uZntN7OvzewKX1k0tuNaM9vpa8vNvrKoa4c/M5tgZu/4fo7KtpjZJDNzfrex0doWADNrY2Y/M7ORZtY5GttiZt2q/E6cmT0ZjW2pC/XZkdUOiL1+W3125FGfXYu2OOdi5gYkAuOAvUChr6wrUAL0Bm4FDgMZwcrD3QZfzEOAP+FdYHAMyPPFGW3tOAv4P8AEvBUbDkfj76NKm3rj7YS4Kprbgrfl+u99t7nR+O/Lry1nAO8DrXyPo/L3gncx/EZgAfBnYB9wXjS2pQ5tVp8dQe3wtSWm+m312RHZFvXZtWhLRK3i0lDOuWLgVTP7GXC+rzjY9tKnBikPuu10I/rCOXcdgJmdAlwMXEn0teMT59zH4I1gAEVE5+8D8Eb6gOvw/mhBlLbFzLoB1wO7gNXAQ8DVRGdbmgELgRzgdTP7AO/3E3VtAf7pnOsHYGapwDpgMNHZllpRnx1x7YAY6rfVZ0dkW9Rn17ItMZWgBxFse+nTgpSHnXNuj9/DBLxPzqf5HkdTO5yZJQOTgcuB/0sU/j78/AfwJHCj73G0tmUw3gYxpwETgRHA677noq0tl+KNSFwEtMMbxSgXVW0pT4p8rgSyiN5/Yw0RdW2OlT4bYq7fVp/9fXmkUJ9duTyomJuDHkCw7aVrte10OJlZF6AlXicTre0oBrYA3wBPAacHqBPx7TCzScAKoAwvrjgC//+J+LY45/7onOsNtAceBE4GfhKgasS3Be/rePB+Hyt8Py8PUC8a2uLvWrw9IKL1/31DRG2bY6TPhhjot9VnH1ceKdRnVy4PqimMoAfbXrogSHlE8H0N9B/Alc65PDM76HsqqtrhnDsG/NX3VekHwBLgh0RZO/BGLYb5PT4Vbx4dRF9bAHDOHQYeMbPewFrgaaKvLYW++7Z4X5kC9PfdR1tbADCzFsAZzrm/m9l5vuKobEs9qc8Osxjpt9VnVy6PFOqzK5cHFasj6Mb3n0yCbS8d6dtOP473VclTZvYmsJ0oa4eZXVW+mgHeaMxuvI4+qtrhcyfehWBD8C5o24jXqURdW8xsuJk94utUAL7GmwcXdW0B3sQb7esFpPrK7iE621LuKmCp7+do7b/qSn12hIihflt9dgS2BfXZtW5LTI2gm1kc3vymM4BEM7sJ7yuHQNtL51rdt51uFGY2EW8kxt8UoqwdwA+AX5jZG3ijYpc6574KEm8kt6PSXDMzKwIOR2tbgBTgJuB6M/tfYL4Lvg17RLfFObfdL74M4NfAy3h/AKKqLX7GA/8JEK2/l9pSnx1Z7fCJiX5bfXZktkV9du3bYr6lYkREREREJALE6hQXEREREZGopARdRERERCSCKEEXEREREYkgStBFRERERCKIEnQRERERkQiiBF1EREREJIIoQRcRERERiSBK0EVEREREIogSdBEJKTNLMrNXzSzPzB4LdzwiIk2R+uLopp1ERSSkzOwC4FvgJ8AdzrlW4Y1IRKTpUV8c3eLDHYCIxBbn3N8AzGwv8L9hDkdEpElSXxzdNIIuIvVmZvOATlWKnwC+A7LwptGd45z7prFjExFpKtQXxx4l6CJyQphZT2ATcJNzbkG44xERaYrUF0cnXSQq1TKzM8zsdTObbWbPmtk0M7sn3HEFYmZvmdkmM0v0Kzunlse2MLM5ZubM7A+RFl+Acz3si/UOMxtoZkVmFhGftn3xpDrnPsf7o/APX3kfM1voi/u0cMYo0tRFat8e5X1xy0jqj9UXRzcl6BKUr/NaBixwzt3rnLsdaA6cGd7IgvoB0A1IAu8KduBPtTnQOXcE+POJCw1oQHzVeNY59w9gbwPPE0o/A943s58CM5xzXwI45z7F+/ckImEUyX17lPfFhyKsP1ZfHMV0kahUpw/QFRjE953XDOD+sEVUvXOBZOfcETNLAF7E64hr60SPeDQ0vuO47+eohX20ppxz7uZwxyAi1Yr0vj2a+2KIkP5YfXF00wi6VGc/Xkdzt++r0F7OuWPAXDM7r/xrPDO708w+MbN/mVnX8oPN7Fwz+5vvK9R3zWygr/zXvuN2mtlpZvaa7/FqM0v31f+fqsGY2RhfvR1m9piZ7Tez983sFDMbAuwCvvJVvwkY7Tsu28xuCtRAM7vczD43s7d9x/g/18fM/mRmT5vZK2aW6isP2nYzSzWzZWb2OzPbY2a/9B1TU3zbfecrNrMzzewk39ewD9f1l2Zm/c3skO9r4va+srt979Vs3+t8a2bPBjj2dDP7zlfnl77f20Ezu9HMpvu+Ft5qZh3rGpeIRIyI6tt9x8ZcX+w7Z6X+WH2x1JpzTjfdgt6AR/A6cgeUAL/FG3kA2OkrH4jXwTngcd9zKcDXwJ98j+cB+4CWQALe2qxFvsc/9B273Vf3DiA9QCzmq1cI9MYbBXHAMt/zX+EbyPA93uH/OMD5WgOHgbd9j+/1ne8PQCqwG5jqe+4A8IDfsQHbDkwB8oHOQC/gPr9jgsaH923WLt95TveVzQTigsT+cDXn6l8l1st85/0/vsdbgD1AsyDnfs9Xf4zvPXLAEaA7cLXv8cx6/nv6ie/408L9b1s33ZryLcL69pjpiwOcr6I/Vl+sW11uGkGXajnnHgBuweuA44E7gYW+p4/56vwDyPWVtfHdjwY64nV0+O7bAKOccyV4cwwTgR8D5wAfAF3N7FwgzTlXfj7/WMq/NtzrnPvMObcOr1O/3MySqfvXiiOBdGCD7/F6v+euxFuyariZ/QL4G94fn3LB2n4M7w/YOl/7X/BvQrBAnHOlfL9O7S3mzRE95Jwrq0uDzLta/7fOuUf8ivv67g/47vcD7fE6/ECO+e7XO+cOlh/rnNsG5Pget6tLXCISWSKpbycG+2II2B+rL5ZaU4IuQZlZXzNr65z7PdADmIXXsY00s9ODHFZ+XUMX332R7768oylfp7X8IqAbgDLgj77H9+GNiNRWru81g3Vw1TmlSoz+yr/OXe+ce9w5N8o592AN54sHXsYb9WgH/BUYUod4/h9QDEwCrgH+Uodjyy0B/s0qT+lZBhwFTjUzA07FG6naX4/zl4trwLEiEkYR2LfHYl8Mx/fH6oul1pSgS3Va4XVOOOdynXP3AeVz5WrqFMpHV8pHOhL9y51zH+B9zTgMrxN9FSgFRlG3q8uT8b6ePRDguZpG1Pf57lsGeK78fBeUF5hZ1U0gAmmP97Xu03jv0QPV1K0Un3PuW2ARcBIwyXlX2tfVBLz340kza+0772d4X5FeDDyI9zu8qh7nFpHYEGl9eyz2xVClP1ZfLHWhBF1qMsPMuvs9zgHW4s2dMwDfSID5ni+/X4I3F/FU3+POwEFfebn5wBfOuX845/YBy/FGEwJ9BeqvpZl1NLNTgLbAW865oirxgDf/EDNrbmZnBDjPX/A6z0vMu5L/Ml95mt9zo8zsJjMbCpznd2ywtg8DLnHO/Tve7m0HgxwTLL7yr1bfrOE9qKr8nP8EnsR7X570nf8kvFGsKc65X/lGoQpqcS7zi7X8vlmVxyISnSKpb4+lvrji/FTpj9UXS10oQZeanARsMLP/Z2ZP4nXKP8LrIMvnvo3GGxUAOMfM2jvnjuJdxNLDzH6D9zXqlc5b47bcK3gdebk/8f0cyOocxbuI6E1gDfAzMxuON2IC3txHgF/jTYFZiPcHpRLn3B68rzDb4M19LMT743QEb+7ijXgX8Py3L/ZFAGY2OFjbgQ7AIjN7AW806FbfMbWKzzn3HrDR997UipmdD5zsezgKWOz7eZKZ3Qmc5mvjp74VAfLNbL0F2JjDzHrjXVAF3sjOj3w/t/e1u9LvubYxikjEiZi+PVb6Yt/5g/bHeBfJqi+WWrHvr7sTiXzm7c72lXPutHDHEkpmFod3AdVhM7vPOfdEDfUfBh5yztU4emJm8XhfNSfh/eFLAc4Cljjnrmlw8HVgZj/Bu1irq3NuR2O+tohITdQXS6TQRkUSNfy+4otF04GpZnYZUOv5jmYWV4vVBX4K5DjnRvoddz/hufo/oeYqIiJho75YIoISdIkm1/ru25nZCOdcdlijCa1/AS3w1u69pxb1N+KtPzzWzN6sYd7+X4ARZjYLyMO7ICoF7w9RozGzPnhfo//RF4eISKRRXywRQVNcREREREQiiC4SFRERERGJIErQRUREREQiiBJ0EREREZEIogRdRERERCSCKEEXEREREYkgStBFRERERCKIEnQRERERkQiiBF1EREREJIL8f3G7XTSeoaQGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Paper figure 4 (Top) with caption\n",
    "# 'Evaluation of the SMP density model parametrization of Proksch et al. (2015) (P15; Top) \n",
    "#  and recalibrated coefficients for snow on sea ice (K20b; Bottom). \n",
    "#  In both cases the model is evaluated against manual density cutter measurements of snow density.'\n",
    "\n",
    "axis_value_size = 12\n",
    "axis_label_size = 14\n",
    "\n",
    "line_start = 100\n",
    "line_end = 700\n",
    "point_size = 8\n",
    "\n",
    "f, (ax1, ax2) = plt.subplots(1, 2, sharey=False, figsize=(12,5))\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "\n",
    "ax1.scatter(result['RHO'], result['mean_samp'], \n",
    "            s = point_size, color ='black', zorder = 1000)\n",
    "ax1.plot([line_start, line_end], [line_start, line_end], \n",
    "         'k-', color = 'gray' ,alpha= 0.8, zorder = 500)\n",
    "\n",
    "ax1.set_xlim(line_start,line_end)\n",
    "ax1.set_ylim(line_start,line_end)\n",
    "\n",
    "rho_bin_size = 20 #in kg m-3\n",
    "common_bin = np.arange(line_start,line_end, rho_bin_size)\n",
    "hist_kws = dict(density=True, bins=common_bin, histtype=\"stepfilled\", linewidth=1.25)\n",
    "ax2.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "\n",
    "ax2.hist(result['RHO'], alpha = 1, edgecolor=\"black\", \n",
    "         color = 'grey', label = 'Pit', **hist_kws)\n",
    "ax2.hist(result['mean_samp'], alpha = 0.6, edgecolor=\"black\", \n",
    "         color = 'deepskyblue', label = 'SMP', **hist_kws)\n",
    "\n",
    "ax1.text(550, 150,'N: %i \\nRMSE: %i \\nR$^2$: %0.2f'%(p2015_n, p2015_rmse, p2015_r2),  fontsize=12)\n",
    "ax1.set_ylabel('SMP P15 density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax1.set_xlabel('Snow pit density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax2.set_ylabel('PDF', fontsize=axis_label_size)\n",
    "ax2.set_xlabel('Snow density [kg m$\\mathregular{^{-3}}$]', fontsize=axis_label_size)\n",
    "ax2.legend(edgecolor = 'black',  fontsize=axis_value_size)\n",
    "ax2.set_xlim(line_start,line_end);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.41"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Error as a % of mean density\n",
    "np.round(rmse(result.error)/ result['RHO'].mean(),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "TYPE\n",
       "F    134.870551\n",
       "H     95.622354\n",
       "R    165.192983\n",
       "Name: error, dtype: float64"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# RMSE by layer type\n",
    "result.groupby('TYPE')['error'].apply(rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "campaign\n",
       "a     97.933150\n",
       "e    134.927592\n",
       "Name: error, dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Errors by campaign (a = alert, e = eureka)\n",
    "campaign_tag = []\n",
    "ice_tag = []\n",
    "lat_tag = []\n",
    "lon_tag = []\n",
    "for x in result['ID']:\n",
    "    ice_tag.append(np.asarray(pit_summary[pit_summary['ID']== x]['ICE'].values))\n",
    "    campaign_tag.append(np.asarray(pit_summary[pit_summary['ID']== x]['IOP'].values))\n",
    "    lat_tag.append(np.asarray(pit_summary[pit_summary['ID']== x]['LAT'].values))\n",
    "    lon_tag.append(np.asarray(pit_summary[pit_summary['ID']== x]['LON'].values))\n",
    "\n",
    "result['campaign'] = np.hstack(campaign_tag)\n",
    "result['ice_type'] = np.hstack(ice_tag)\n",
    "result['lat'] = np.hstack(lat_tag)\n",
    "result['long'] = np.hstack(lon_tag)\n",
    "\n",
    "result.groupby('campaign')['error'].apply(rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ice_type\n",
       "f    153.150085\n",
       "m     86.272081\n",
       "Name: error, dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Errors by ice type (f = fyi, m = myi)\n",
    "result.groupby('ice_type')['error'].apply(rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the dataset\n",
    "result.to_pickle(\"./output/smp_pit_filtered.pkl\")"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}