component-clustering/model-variance-exploration.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# How resilient are our models?\n",
        "If we take many (>30) classifications of a galaxy and randomly pick 30, what model do we get, and how does this vary depending on the classifications we pick?\n",
        "\n",
        "*This script has been copied to run on all galaxies in `model-variance-exploration.py`*"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, load up our environment!"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "%load_ext autoreload\n",
        "%autoreload 2"
      ],
      "outputs": [],
      "execution_count": 1,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import json\n",
        "import os\n",
        "import re\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "from sklearn.model_selection import ShuffleSplit\n",
        "from shapely.geometry import MultiPolygon\n",
        "from descartes import PolygonPatch\n",
        "import lib.galaxy_utilities as gu\n",
        "import lib.python_model_renderer.parse_annotation as pa\n",
        "import lib.average_shape_helpers as ash\n",
        "import get_average_shape as gas\n",
        "from gzbuilderspirals import get_drawn_arms\n",
        "from gzbuilderspirals.oo import Pipeline"
      ],
      "outputs": [],
      "execution_count": 2,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now grab all classifications for our galaxy (need a galaxy from the validation subset in order to have the required number of classifications for the shuffle split)"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "dr8ids, ss_ids, validation_ids = np.load('lib/duplicate_galaxies.npy').T\n",
        "ss_id, val_id = ss_ids[94], validation_ids[42]\n",
        "print('Original id:', ss_id)\n",
        "print('Validation id:', val_id)"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Original id: 20902040\n",
            "Validation id: 21686557\n"
          ]
        }
      ],
      "execution_count": 3,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "all_cls = gu.classifications.query(\n",
        "    '(subject_ids == {}) or (subject_ids == {})'.format(\n",
        "        ss_id, val_id\n",
        "    )\n",
        ")\n",
        "all_models = all_cls['annotations']\\\n",
        "    .apply(json.loads)\\\n",
        "    .apply(ash.remove_scaling)\\\n",
        "    .apply(pa.parse_annotation)\n",
        "all_geoms = pd.DataFrame(\n",
        "    all_models.apply(gas.get_geoms).values.tolist(),\n",
        "    columns=('disk', 'bulge', 'bar')\n",
        ")"
      ],
      "outputs": [],
      "execution_count": 4,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For the purposes of demonstration, use 10 splits. `model-variance-exploration.py` uses 20 splits but takes ~4 hours to run."
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "N_SPLITS = 5"
      ],
      "outputs": [],
      "execution_count": 5,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now calculate the models for each of the splits, and save the results as a DataFrame"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ss = ShuffleSplit(n_splits=N_SPLITS, test_size=0.5, random_state=0)\n",
        "split_models = []\n",
        "pas = []\n",
        "\n",
        "gal, angle = gu.get_galaxy_and_angle(ss_id)\n",
        "\n",
        "for i, (train_index, _) in enumerate(ss.split(all_geoms)):\n",
        "    print('Split {}/{}'.format(i+1, N_SPLITS))\n",
        "    models = all_models.iloc[train_index]\n",
        "    drawn_arms = get_drawn_arms((ss_id, val_id), all_cls.iloc[train_index])\n",
        "    if len(drawn_arms) > 0:\n",
        "        p = Pipeline(drawn_arms, phi=angle, ba=gal['PETRO_BA90'],\n",
        "                     image_size=512, parallel=True)\n",
        "        pas.append(p.get_pitch_angle(p.get_arms()))\n",
        "    else:\n",
        "        pas.append((np.nan, np.nan))\n",
        "    geoms = all_geoms.iloc[train_index]\n",
        "    labels = list(map(np.array, gas.cluster_components(geoms)))\n",
        "    aggregate_disk, aggregate_bulge, aggregate_bar = gas.get_aggregate_components(\n",
        "        geoms, models, labels\n",
        "    )\n",
        "    split_models.append({\n",
        "        'disk': aggregate_disk if aggregate_disk else None,\n",
        "        'bulge': aggregate_bulge if aggregate_bulge else None,\n",
        "        'bar': aggregate_bar if aggregate_bar else None,\n",
        "    })"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Split 1/5\n",
            "Split 2/5\n",
            "Split 3/5\n",
            "Split 4/5\n",
            "Split 5/5\n"
          ]
        }
      ],
      "execution_count": 6,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "splits_df = []\n",
        "for model in split_models:\n",
        "    model_comps = {}\n",
        "    for key in ('disk', 'bulge', 'bar'):\n",
        "        if model.get(key, None) is None:\n",
        "            model[key] = {}\n",
        "        mu = model[key].get('mu', (np.nan, np.nan))\n",
        "        model_comps['{}-mux'.format(key)] = mu[0]\n",
        "        model_comps['{}-muy'.format(key)] = mu[1]\n",
        "        for param in ('roll', 'rEff', 'axRatio', 'i0', 'n', 'c'):\n",
        "            model_comps['{}-{}'.format(key, param)] = (\n",
        "                model[key].get(param, np.nan)\n",
        "            )\n",
        "    splits_df.append(model_comps)\n",
        "splits_df = pd.DataFrame(splits_df)"
      ],
      "outputs": [],
      "execution_count": 7,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pas = pd.DataFrame(\n",
        "    pas,\n",
        "    columns=('pa', 'sigma_pa'),\n",
        "    index=pd.Series(range(len(pas)), name='split_index')\n",
        ")\n",
        "pas.describe().loc[['mean', 'std']]"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 8,
          "data": {
            "text/plain": [
              "             pa  sigma_pa\n",
              "mean  19.428581  1.216343\n",
              "std    4.195982  0.406637"
            ],
            "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>pa</th>\n",
              "      <th>sigma_pa</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>mean</th>\n",
              "      <td>19.428581</td>\n",
              "      <td>1.216343</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>std</th>\n",
              "      <td>4.195982</td>\n",
              "      <td>0.406637</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 8,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "split_geoms = list(map(gas.get_geoms, split_models))"
      ],
      "outputs": [],
      "execution_count": 9,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "disks = [i[0] for i in split_geoms]\n",
        "MultiPolygon(disks)"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 10,
          "data": {
            "text/plain": [
              "<shapely.geometry.multipolygon.MultiPolygon at 0x1c2b9a69b0>"
            ],
            "image/svg+xml": [
              "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"195.57628355911993\" height=\"148.0398775044082\" viewBox=\"159.36833005740618 180.28027323046103 195.57628355911993 148.0398775044082\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,508.60042396533026)\"><g><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M 166.70421754328453,246.63742411232397 L 166.61189611515135,253.01206494291466 L 167.39156844792666,259.39786209279237 L 169.0357258708206,265.7333168445701 L 171.52853424906664,271.9574153056802 L 174.84598647528202,278.01021600550325 L 178.95613367101566,283.83342716450045 L 183.81939287193967,289.37096807592087 L 189.38892823342331,294.5695091936625 L 195.61110208531613,299.3789857249427 L 202.4259914920434,303.753079781643 L 209.7679653432088,307.64966644687456 L 217.56631641710703,311.03121946096024 L 225.74594232987357,313.8651726198249 L 234.22806881253888,316.12423340536094 L 242.9310083502741,317.78664582726435 L 251.77094687780775,318.8363999451537 L 260.66275095472906,319.2633860530177 L 269.52078764701923,319.0634920412259 L 278.25974921903617,318.23864299837555 L 286.79547469363456,316.7967826716683 L 295.0457603684057,314.7517969642448 L 302.9311514822182,312.12338020637253 L 310.3757074080652,308.93684548827054 L 317.3077330027678,305.2228808812098 L 323.66046907042465,301.01725389461154 L 329.37273528993524,296.3604670153926 L 334.3895194149694,291.2973676468564 L 338.6625070719485,285.87671620369235 L 342.1505470538814,280.1507165225289 L 344.8200476289892,274.17451311044704 L 346.6453000474404,268.0056600732269 L 347.60872613064777,261.70356683785485 L 347.70104755878094,255.32892600726422 L 346.92137522600564,248.94312885738645 L 345.2772178031116,242.6076741056088 L 342.78440942486594,236.383575644495 L 339.4669571986539,230.3307749446759 L 335.3568100029206,224.50756378567507 L 330.49355080199655,218.97002287425465 L 324.9240154405089,213.77148175651632 L 318.70184158861616,208.9620052252361 L 311.8869521818928,204.58791116853246 L 304.5449783307274,200.69132450330096 L 296.74662725682526,197.30977148921858 L 288.56700134405907,194.4758183303503 L 280.0848748613934,192.21675754481794 L 271.38193532366216,190.55434512291112 L 262.54199679612486,189.5045910050215 L 253.65019271920718,189.0776048971578 L 244.79215602691306,189.27749890895296 L 236.05319445489647,190.1023479517997 L 227.51746898029774,191.54420827851044 L 219.26718330553416,193.589193985931 L 211.3817921917181,196.21761074380294 L 203.93723626586717,199.40414546190823 L 197.00521067116443,203.11811006896903 L 190.65247460351156,207.32373705556398 L 184.94020838400104,211.98052393478292 L 179.92342425896652,217.04362330332276 L 175.6504366019841,222.46427474648283 L 172.1623966200512,228.19027442764633 L 169.4928960449431,234.16647783973178 L 167.6676436264919,240.33533087695196 L 166.70421754328453,246.63742411232397 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M 167.4039872928804,238.05276414873128 L 166.7922843103247,243.86047804289456 L 167.05151263847256,249.76993449900868 L 168.17917576681117,255.72422222096128 L 170.1644136829699,261.6659981630452 L 172.98810746058427,267.5380397748362 L 176.6230633849794,273.28379608562636 L 181.03427484353637,278.8479323211926 L 186.17925945849908,284.17686280792117 L 192.0084682154881,289.21926703213313 L 198.4657626476057,293.92658388469374 L 205.48895547954288,298.2534793310292 L 213.01040952508117,302.15828300266975 L 220.95768907009665,305.6033895056779 L 229.2542574681599,308.5556205811386 L 237.82021423023053,310.9865446299434 L 246.5730645100528,312.8727505246244 L 255.42851357464278,314.19607307131014 L 264.3012786086288,314.94376795047634 L 273.105910034386,315.10863445170276 L 281.75761443818254,314.68908482045515 L 290.1730711771079,313.68915954903855 L 298.2712348024062,312.118488464424 L 305.974115571468,309.9921979877777 L 313.2075305316525,307.3307654587365 L 319.9018179427162,304.1598219274415 L 325.9925081574141,300.5099053135132 L 331.42094449947194,296.41616630919145 L 336.1348481594325,291.91802985895407 L 340.0888216681635,287.058815475748 L 343.2447860993134,281.88532005039866 L 345.57234779019757,276.4473671719745 L 347.04909104933927,270.7973272993806 L 347.6607940318982,264.98961340521976 L 347.40156570374677,259.08015694910495 L 346.27390257541174,253.12586922715303 L 344.2886646592498,247.18409328506664 L 341.46497088163505,241.31205167327744 L 337.8300149572403,235.5662953624855 L 333.41880349868296,230.00215912692101 L 328.2738188837202,224.67322864019246 L 322.44461012673514,219.63082441597936 L 315.98731569461756,214.92350756341875 L 308.96412286267673,210.5966121170826 L 301.4426688171421,206.69180844544275 L 293.4953892721227,203.24670194243572 L 285.1988208740598,200.2944708669732 L 276.6328641119927,197.86354681816908 L 267.88001383217045,195.97734092348804 L 259.02456476758044,194.6540183768023 L 250.15179973359085,193.90632349763547 L 241.3471683078333,193.74145699641088 L 232.6954639040407,194.16100662765734 L 224.2800071651153,195.1609318990739 L 216.18184353981704,196.7316029836885 L 208.47896277075523,198.85789346033482 L 201.24554781056685,201.51932598937714 L 194.55126039950312,204.69026952067213 L 188.4605701848056,208.34018613459864 L 183.03213384274738,212.43392513892218 L 178.31823018278715,216.93206158915774 L 174.36425667405973,221.79127597236447 L 171.20829224290628,226.9647713977132 L 168.88073055202563,232.40272427613797 L 167.4039872928804,238.05276414873128 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M 326.41391185556034,252.50716951088853 L 326.0340141493831,249.116521417032 L 324.99596262348973,245.768076201365 L 323.30975428153056,242.49408122175876 L 320.99162823161805,239.32606683961373 L 318.0639092948217,236.2945427652444 L 314.5547930049628,233.4287042326859 L 310.4980740701904,230.7561508336347 L 305.9328209114726,228.30262071828798 L 300.902999412346,226.09174272289334 L 295.4570495034331,224.1448088111354 L 289.64741865941585,222.48056902088575 L 283.5300568011831,221.11505089109383 L 277.1638774675291,220.06140510781276 L 270.6101904455375,219.32977885592135 L 263.93211132383726,218.92721809616944 L 257.19395365498747,218.85759970873826 L 250.4606095807875,219.12159415675393 L 243.79692488549242,219.71665902937013 L 237.2670744954552,220.63706352656533 L 230.93394443951985,221.87394364989055 L 224.85852622215512,223.41538756760661 L 219.0993294419651,225.2465503321516 L 213.71181831223387,227.3497968450843 L 208.7478775102576,229.70487169273474 L 204.2553124995204,232.28909421690497 L 200.27738913696922,235.0775769420066 L 196.852416999134,238.0434652550593 L 194.01338043998197,241.1581960302877 L 191.78762093355888,244.3917727076418 L 190.1965737606315,247.7130541760644 L 189.25556157522323,251.0900546794184 L 188.97364683910297,254.490251856802 L 189.3535445452765,257.8808999506567 L 190.3915960711735,261.2293451663237 L 192.07780441312903,264.5033401459318 L 194.39593046304162,267.6713545280768 L 197.32364939983793,270.70287860244434 L 200.83276568969683,273.5687171350029 L 204.88948462446922,276.24127053405584 L 209.45473778319067,278.6948006494025 L 214.48455928231365,280.90567864479726 L 219.93050919122655,282.8526125565552 L 225.7401400352438,284.51685234680303 L 231.85750189347655,285.88237047659675 L 238.2236812271305,286.93601625987606 L 244.77736824912583,287.66764251176915 L 251.45544737082236,288.0702032715193 L 258.19360503967215,288.1398216589523 L 264.9269491138721,287.87582721093486 L 271.5906338091709,287.2807623383204 L 278.12048419920444,286.36035784112346 L 284.45361425513977,285.12347771780003 L 290.5290324725009,283.5820338000822 L 296.2882292526946,281.75087103553903 L 301.6757403824258,279.6476245226063 L 306.6396811844057,277.2925496749558 L 311.13224619513926,274.70832715078564 L 315.1101695576904,271.9198444256822 L 318.5351416955257,268.9539561126294 L 321.3741782546777,265.8392253374028 L 323.59993776110076,262.6056486600488 L 325.19098493403175,259.2843671916262 L 326.1319971194364,255.9073666882704 L 326.41391185556034,252.50716951088853 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M 169.50714011909122,244.04169183930827 L 169.1746582200389,250.56453859416393 L 169.69172424032791,257.12336192587713 L 171.05335855150997,263.65499678034706 L 173.24644787083366,270.0965399433701 L 176.2498715492914,276.3859558329781 L 180.03470497501468,282.46267393601283 L 184.56449813311326,288.2681721352704 L 189.79562663924264,293.74654030948255 L 195.67771186626322,298.84501877835066 L 202.15410611796045,303.5145064071083 L 209.1624381773281,307.71003347728896 L 216.63521397550437,311.3911947696963 L 224.50046659650383,314.5225386887596 L 232.68244935801232,317.07390868077016 L 241.10236529330504,319.0207336579837 L 249.67912600913027,320.3442646316322 L 258.3301326112622,321.031755274914 L 266.9720711770042,321.07658467712406 L 275.52171511379385,320.47832110660624 L 283.8967266768152,319.24272616858957 L 292.01644992647914,317.38169931777566 L 299.8026874891789,314.913163260072 L 307.18045364081183,311.86089134715917 L 314.07869646023187,308.2542786260982 L 320.4309820981289,304.12805874896765 L 326.1761345712973,299.52196946880565 L 331.25882492079813,294.4803699433263 L 335.63010406008436,289.0518135319782 L 339.24787418146195,283.28858020054975 L 342.07729418098904,277.24617303652934 L 344.0911151973512,270.98278372402353 L 345.26994303329354,264.5587321260256 L 345.6024249323423,258.0358853711695 L 345.08535891205685,251.4770620394567 L 343.7237246008712,244.94542718498639 L 341.5306352815475,238.50388402196336 L 338.5272116030938,232.21446813235218 L 334.74237817737054,226.13775002931743 L 330.2125850192683,220.33225183005948 L 324.98145651313894,214.8538836558473 L 319.09937128612205,209.7554051869796 L 312.6229770344248,205.08591755822198 L 305.61464497505665,200.8903904880449 L 298.14186917687687,197.20922919563708 L 290.2766165558778,194.07788527657027 L 282.09463379437295,191.5265152845601 L 273.67471785907657,189.57969030734608 L 265.09795714325094,188.25615933370125 L 256.4469505411194,187.5686686904158 L 247.80501197538103,187.5238392882062 L 239.2553680385878,188.12210285872357 L 230.88035647556603,189.35769779674382 L 222.76063322590932,191.21872464755864 L 214.97439566320634,193.6872607052582 L 207.59662951157338,196.7395326181711 L 200.69838669214977,200.34614533923167 L 194.34610105425276,204.47236521636222 L 188.6009485810879,209.0784544965246 L 183.51825823158708,214.12005402200396 L 179.14697909229682,219.54861043335524 L 175.52920897091968,225.31184376478012 L 172.69978897139612,231.35425092880092 L 170.68596795502992,237.6176402413099 L 169.50714011909122,244.04169183930827 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"2.0\" opacity=\"0.6\" d=\"M 327.6788274168497,250.6142926060495 L 327.2137576720351,247.32649691057506 L 326.08346858141,244.09678326490172 L 324.29884544680556,240.95625557697434 L 321.8770751645746,237.93515884521736 L 318.8414807064087,235.06258788265797 L 315.2212965063693,232.36620711796854 L 311.0513869172962,229.87198417190842 L 306.3719104480318,227.6039397749281 L 301.22793301503947,225.5839164344107 L 295.6689939329836,223.83136807938058 L 289.7486288240575,222.36317270853655 L 283.52385404064444,221.1934698458998 L 277.05461756673077,220.33352436946092 L 270.4032216859961,219.79161802425116 L 263.6337229767991,219.5729696646111 L 256.8113154122934,219.679684993771 L 250.00170250683553,220.11073628477277 L 243.27046455520676,220.86197227804956 L 236.68242705856417,221.92615816031662 L 230.30103641938172,223.29304523978453 L 224.1877489179326,224.94946964666298 L 218.401438854733,226.87947910842007 L 212.9978315587936,229.06448657888726 L 208.0289667222944,231.4834492416772 L 203.5426972298696,234.1130711640096 L 199.5822283091777,236.92802764927364 L 196.1857014409328,239.90120912770584 L 193.38582703560562,243.00398223635898 L 191.2095694142915,246.20646557405212 L 189.67788712754646,249.47781747558687 L 188.8055311130866,252.7865330338467 L 188.6009026362195,256.10074750928084 L 189.06597238103402,259.38854320475525 L 190.19626147165914,262.6182568504286 L 191.98088460626357,265.75878453835594 L 194.40265488849454,268.7798812701129 L 197.43824934666057,271.6524522326742 L 201.0584335466999,274.34883299736174 L 205.22834313577297,276.8430559434219 L 209.9078196050373,279.1111003404022 L 215.05179703802978,281.13112368092146 L 220.61073612008565,282.8836720359516 L 226.5311012290117,284.3518674067956 L 232.75587601242475,285.5215702694305 L 239.22511248633845,286.38151574587124 L 245.87650836706953,286.92342209108114 L 252.64600707627014,287.142070450721 L 259.4684146407722,287.03535512156134 L 266.2780275462337,286.60430383055933 L 273.00926549785885,285.8530678372827 L 279.59730299450507,284.78888195501554 L 285.9786936336875,283.42199487554575 L 292.091981135133,281.7655704686693 L 297.8782911983362,279.8355610069121 L 303.28189849427565,277.65055353644493 L 308.25076333077476,275.2315908736531 L 312.73703282319957,272.6019689513207 L 316.69750174389145,269.7870124660567 L 320.09402861213647,266.81383098762626 L 322.8939030174635,263.7110578789713 L 325.0701606387777,260.50857454128 L 326.60184292552276,257.23722263974526 L 327.47419893998256,253.92850708148362 L 327.6788274168497,250.6142926060495 z\" /></g></g></svg>"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 10,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "bulges = [i[1] for i in split_geoms]\n",
        "MultiPolygon(bulges)"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 11,
          "data": {
            "text/plain": [
              "<shapely.geometry.multipolygon.MultiPolygon at 0x1c2b992ef0>"
            ],
            "image/svg+xml": [
              "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"243.07172385765887 239.81520399256036 28.92464314060146 28.839475822004175\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,508.4698838071249)\"><g><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"0.5784928628120292\" opacity=\"0.6\" d=\"M 251.76025596195367,243.1473007924967 L 250.73737568073713,243.6333453032246 L 249.77107285609247,244.22864521247266 L 248.87065351245465,244.92746745668956 L 248.04478918032345,245.72308199564907 L 247.30143338467974,246.60782662641458 L 246.64774504819573,247.5731807745405 L 246.09001954689396,248.60984755186382 L 245.63362808223567,249.70784329061985 L 245.28296595351466,250.856593691618 L 245.04141022872142,252.04503566052716 L 244.9112872215341,253.26172385151654 L 244.89385008764472,254.49494089219917 L 244.98926675618327,255.7328102283272 L 245.19661831247282,256.9634105015057 L 245.5139078476767,258.17489035839003 L 245.9380796901217,259.3555825856953 L 246.46504883309004,260.4941164718408 L 247.0897402756638,261.5795273131177 L 247.80613789776018,262.6013620097773 L 248.60734239865573,263.5497797350997 L 249.48563774102334,264.41564670793315 L 250.43256546059106,265.19062415599643 L 251.4390061257772,265.8672486228203 L 252.49526716279792,266.4390038449079 L 253.59117620045694,266.90038350691805 L 254.71617903563617,267.2469442704892 L 255.85944127604412,267.47534856601925 L 257.00995268133477,267.583396735283 L 258.15663319774114,267.57004821534633 L 259.2884396650519,267.4354315597524 L 260.39447216827773,267.18084320048405 L 261.4640790097941,266.8087349626127 L 262.4869592910102,266.32269045188497 L 263.453262115655,265.72739054263735 L 264.3536814592927,265.02856829842005 L 265.1795457914241,264.23295375946 L 265.9229015870676,263.34820912869503 L 266.57658992355164,262.38285498056905 L 267.1343154248534,261.34618820324573 L 267.59070688951164,260.2481924644897 L 267.9413690182331,259.09944206349144 L 268.1829247430259,257.9110000945824 L 268.3130477502136,256.6943119035929 L 268.33048488410304,255.46109486291027 L 268.2350682155645,254.22322552678224 L 268.0277166592749,252.9926252536037 L 267.7104271240711,251.78114539671935 L 267.2862552816256,250.60045316941427 L 266.7592861386577,249.46191928326863 L 266.1345946960839,248.37650844199175 L 265.4181970739876,247.3546737453321 L 264.61699257309203,246.40625602000972 L 263.7386972307247,245.5403890471771 L 262.7917695111563,244.76541159911318 L 261.7853288459705,244.0887871322891 L 260.72906780894937,243.5170319102017 L 259.63315877129037,243.05565224819156 L 258.50815593611117,242.7090914846204 L 257.3648936957032,242.4806871890903 L 256.2143822904128,242.372639019826 L 255.06770177400617,242.38598753976328 L 253.93589530669598,242.52060419535746 L 252.82986280347004,242.7751925546254 L 251.76025596195367,243.1473007924967 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"0.5784928628120292\" opacity=\"0.6\" d=\"M 248.55829090762637,251.77912306128394 L 248.47088069418393,252.35358583590317 L 248.4878656490475,252.9481794515936 L 248.60908219781558,253.55717764679798 L 248.83336295886065,254.17471543598447 L 249.15854798584957,254.79484559266245 L 249.58150556920765,255.41159592442597 L 250.09816239619562,256.01902678842976 L 250.70354277913586,256.6112882933978 L 251.39181657400667,257.18267663727477 L 252.15635532791407,257.72768903795605 L 252.98979611471796,258.24107672808884 L 253.88411244403812,258.71789550357124 L 254.83069156073827,259.15355333893973 L 255.82041739047304,259.54385461108467 L 256.8437583324631,259.8850405053951 L 257.8908590540208,260.1738252151976 L 258.95163540279475,260.4074275858709 L 260.0158715226675,260.58359789888567 L 261.07331823803287,260.70063953782267 L 262.11379175896303,260.7574253277174 L 263.1272717566717,260.75340839037165 L 264.1039978647548,260.68862741108705 L 265.0345636768604,260.56370626610754 L 265.91000733551823,260.379848014348 L 266.72189783972857,260.1388233112801 L 267.46241624011213,259.84295335655156 L 268.12443093966425,259.4950875395654 L 268.7015663749346,259.09857599830275 L 269.18826441618694,258.6572373556603 L 269.5798378952277,258.17532194402736 L 269.8725157453987,257.6574708722598 L 270.0634793190077,257.1086713292646 L 270.15088953245,256.53420855464566 L 270.13390457758646,255.93961493895517 L 270.0126880288184,255.33061674375085 L 269.7884072677734,254.71307895456414 L 269.46322224078443,254.09294879788615 L 269.04026465742635,253.47619846612264 L 268.52360783043844,252.86876760211885 L 267.9182274474981,252.27650609715073 L 267.22995365262733,251.70511775327384 L 266.4654148987199,251.16010535259255 L 265.631974111916,250.6467176624597 L 264.7376577825959,250.16989888697714 L 263.7910786658957,249.73424105160882 L 262.80135283616096,249.34393977946388 L 261.77801189417096,249.00275388515325 L 260.7309111726132,248.71396917535077 L 259.6701348238393,248.4803668046775 L 258.6058987039665,248.30419649166288 L 257.5484519886011,248.18715485272588 L 256.507978467671,248.130369062831 L 255.4944984699627,248.13438600017702 L 254.5177723618792,248.19916697946155 L 253.5872065497736,248.32408812444106 L 252.7117628911158,248.50794637620055 L 251.89987238690543,248.74897107926853 L 251.15935398652192,249.04484103399705 L 250.49733928696975,249.39270685098313 L 249.9202038516994,249.78921839224586 L 249.4335058104471,250.2305570348883 L 249.04193233140634,250.7124724465212 L 248.7492544812353,251.2303235182888 L 248.55829090762637,251.77912306128394 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"0.5784928628120292\" opacity=\"0.6\" d=\"M 244.22646591424743,251.98133889774874 L 244.1430069369404,253.2230845135361 L 244.1885109805891,254.46819327340046 L 244.36253981637654,255.7046740993393 L 244.66341745148034,256.92061900506934 L 245.08824626979964,258.1043177764456 L 245.63293493757732,259.244370747224 L 246.2922378051864,260.3297985840793 L 247.0598054255999,261.35014802358836 L 247.92824570303813,262.29559254286767 L 248.88919508289197,263.15702699436474 L 249.9333990973237,263.92615529340395 L 251.05080149085651,264.5955703140165 L 252.2306410676043,265.1588252236149 L 253.46155532747855,265.6104955695131 L 254.73168989327,265.9462315193721 L 256.0288126747839,266.16279975246294 L 257.34043167055455,266.25811459831226 L 258.6539152726492,266.2312581228506 L 259.9566139159565,266.0824889686177 L 261.2359819004131,265.8132398638943 L 262.4796982129571,265.4261038247457 L 263.67578518561083,264.924809182858 L 264.8127238469832,264.3141836796695 L 265.8795648562497,263.6001079725771 L 266.8660339512974,262.7894590009931 L 267.76263089548047,261.8900437576607 L 268.56072097009803,260.91052410304263 L 269.25261813145147,259.86033334687335 L 269.83165903164604,258.7495854002272 L 270.29226719027747,257.58897737302686 L 270.6300066989854,256.38968655502765 L 270.84162494167185,255.163262772405 L 270.92508391897877,253.92151715661856 L 270.87957987533014,252.67640839675374 L 270.705551039543,251.43992757081537 L 270.4046734044389,250.22398266508486 L 269.9798445861195,249.04028389370853 L 269.4351559183419,247.90023092293006 L 268.77585305073325,246.8148030860749 L 268.0082854303193,245.79445364656584 L 267.1398451528811,244.84900912728654 L 266.1788957730273,243.98757467578895 L 265.134691758596,243.2184463767503 L 264.0172893650627,242.54903135613768 L 262.8374497883149,241.98577644653932 L 261.60653552844065,241.5341061006411 L 260.33640096264924,241.19837015078207 L 259.0392781811357,240.98180191769126 L 257.7276591853647,240.88648707184188 L 256.41417558327,240.9133435473036 L 255.11147693996276,241.06211270153642 L 253.83210895550616,241.3313618062599 L 252.58839264296304,241.7184978454086 L 251.39230567030847,242.21979248729565 L 250.255367008936,242.8304179904847 L 249.18852599966948,243.54449369757708 L 248.20205690462222,244.35514266916107 L 247.3054599604387,245.25455791249345 L 246.5073698858211,246.2340775671115 L 245.81547272446775,247.28426832328086 L 245.23643182427358,248.39501626992705 L 244.77582366564172,249.55562429712734 L 244.43808415693383,250.7549151151265 L 244.22646591424743,251.98133889774874 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"0.5784928628120292\" opacity=\"0.6\" d=\"M 248.00066795855128,251.2744853424008 L 247.89941735737193,251.85784388050382 L 247.90434296100335,252.4640225758127 L 248.0153973331893,253.08718359612067 L 248.23151095961748,253.72132556014765 L 248.55060254792556,254.36034133412255 L 248.96959907168537,254.99807684682133 L 249.48446536532532,255.62839035663387 L 250.09024298498414,256.2452115998972 L 250.78109796103217,256.8426002508492 L 251.55037698239084,257.4148031302159 L 252.39067147155265,257.95630961147043 L 253.29388893323724,258.4619046911816 L 254.25133088953748,258.9267192123449 L 255.25377665102445,259.34627675702785 L 256.29157211702255,259.7165367567222 L 257.35472274987694,260.0339334052259 L 258.43298982781624,260.29540999930714 L 259.5159890494393,260.4984483764275 L 260.59329054021373,260.6410931660231 L 261.6545192978708,260.72197062079226 L 262.6894551093486,260.7403018466342 L 263.688130977027,260.69591030382253 L 264.64092910637055,260.58922350718177 L 265.53867353054,260.42126890888426 L 266.3727184799743,260.19366400352504 L 267.13503164587405,259.9086007507642 L 267.8182715357304,259.56882446555716 L 268.41585817590556,259.1776073792715 L 268.9220364803815,258.7387171263086 L 269.33193167538053,258.2563804597264 L 269.6415962461038,257.7352425452966 L 269.8480479534546,257.18032222601755 L 269.9492985546339,256.59696368791475 L 269.9443729510025,255.99078499260588 L 269.83331857881655,255.3676239722979 L 269.61720495238836,254.73348200827093 L 269.29811336408034,254.0944662342958 L 268.87911684032053,253.45673072159724 L 268.3642505466806,252.8264172117847 L 267.75847292702167,252.2095959685214 L 267.0676179509737,251.61220731756936 L 266.29833892961506,251.0400044382027 L 265.4580444404533,250.49849795694792 L 264.55482697876863,249.992902877237 L 263.5973850224684,249.52808835607368 L 262.5949392609814,249.10853081139072 L 261.5571437949833,248.7382708116964 L 260.4939931621289,248.42087416319265 L 259.41572608418966,248.15939756911143 L 258.3327268625666,247.95635919199106 L 257.2554253717921,247.81371440239548 L 256.1941966141351,247.7328369476263 L 255.1592608026582,247.71450572178435 L 254.16058493497889,247.758897264596 L 253.2077868056353,247.8655840612368 L 252.31004238146582,248.03353865953426 L 251.4759974320316,248.26114356489353 L 250.7136842661318,248.54620681765437 L 250.03044437627557,248.8859831028612 L 249.4328577361003,249.2772001891471 L 248.92667943162434,249.71609044211 L 248.51678423662537,250.19842710869187 L 248.20711966590216,250.71956502312173 L 248.00066795855128,251.2744853424008 z\" /><path fill-rule=\"evenodd\" fill=\"#ff3333\" stroke=\"#555555\" stroke-width=\"0.5784928628120292\" opacity=\"0.6\" d=\"M 248.95234076578362,252.06851960224537 L 248.89626158640044,252.60891661216493 L 248.93609629234822,253.16274411272718 L 249.0714612536313,253.72466844234856 L 249.30105283167575,254.28927796255613 L 249.62265993407743,254.8511351750234 L 250.03318530863118,255.40482908765148 L 250.52867537156442,255.9450273253821 L 251.10435828270897,256.4665274838838 L 251.75468990092995,256.9643072315488 L 252.4734071772365,257.4335726772933 L 253.25358847136636,257.8698045383478 L 254.08772021096814,258.26880166342215 L 254.9677692514047,258.6267214920877 L 255.88526023933426,258.9401170607355 L 256.83135723499174,259.2059701987202 L 257.7969488071183,259.42172059499273 L 258.77273578102916,259.5852904552957 L 259.7493207947464,259.69510451245833 L 260.71729880073167,259.7501051970795 L 261.66734764163664,259.7497628224997 L 262.59031782777527,259.69408068597386 L 263.477320651706,259.58359503691497 L 264.31981379135124,259.41936991252146 L 265.1096835772217,259.2029868905159 L 265.83932313149296,258.9365298576869 L 266.50170562639534,258.6225649409211 L 267.09045195640823,258.2641157939975 L 267.5998921725325,257.86463447814924 L 268.02512008699784,257.4279682168266 L 268.3620405225347,256.9583223448323 L 268.60740875117045,256.46021980864816 L 268.7588617427303,255.93845760798786 L 268.81494092211335,255.39806059806858 L 268.77510621616557,254.84423309750633 L 268.63974125488255,254.28230876788473 L 268.41014967683816,253.71769924767693 L 268.0885425744365,253.1558420352099 L 267.6780171998828,252.60214812258158 L 267.1825271369495,252.06194988485095 L 266.6068442258049,251.5404497263495 L 265.95651260758393,251.0426699786845 L 265.2377953312773,250.57340453294 L 264.45761403714755,250.13717267188548 L 263.6234822975457,249.73817554681114 L 262.7434332571092,249.38025571814558 L 261.8259422691797,249.0668601494976 L 260.8798452735225,248.80100701151315 L 259.914253701396,248.5852566152406 L 258.93846672748515,248.4216867549377 L 257.9618817137675,248.31187269777473 L 256.9939037077822,248.25687201315378 L 256.0438548668772,248.2572143877336 L 255.1208846807395,248.3128965242596 L 254.23388185680784,248.42338217331832 L 253.3913887171626,248.58760729771183 L 252.60151893129216,248.8039903197174 L 251.87187937702097,249.07044735254618 L 251.2094968821186,249.38441226931195 L 250.62075055210562,249.74286141623577 L 250.11131033598136,250.14234273208388 L 249.686082421516,250.57900899340643 L 249.3491619859792,251.04865486540098 L 249.10379375734337,251.54675740158513 L 248.95234076578362,252.06851960224537 z\" /></g></g></svg>"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 11,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "bars = [i[2] for i in split_geoms]\n",
        "MultiPolygon(bars)"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 12,
          "data": {
            "text/plain": [
              "<shapely.geometry.multipolygon.MultiPolygon at 0x1c2b972a58>"
            ],
            "image/svg+xml": [
              "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"100.0\" height=\"100.0\" viewBox=\"229.50764881725465 236.47328718318005 53.94461153192955 32.96172596470575\" preserveAspectRatio=\"xMinYMin meet\"><g transform=\"matrix(1,0,0,-1,0,505.90830033106585)\"><g><path fill-rule=\"evenodd\" fill=\"#66cc99\" stroke=\"#555555\" stroke-width=\"1.078892230638591\" opacity=\"0.6\" d=\"M 231.5055973925113,256.066328198162 L 235.897210366945,238.4712357584367 L 281.45431177392754,249.84197213290383 L 277.0626987994939,267.43706457262914 L 231.5055973925113,256.066328198162 z\" /></g></g></svg>"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 12,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# splits_df.describe().transpose()[['count', 'mean', 'std']]"
      ],
      "outputs": [],
      "execution_count": 24,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Big picture\n",
        "\n",
        "We've run this code on all galaxies in the validation subset, with 20 splits for each galaxy. What do the results show us?"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "files = [f for f in os.listdir('model-variances')\n",
        "         if re.match(r'[0-9]+_.*?\\.pickle', f)]\n",
        "\n",
        "pas = [os.path.join('model-variances', f) for f in files if '_pa.' in f]\n",
        "comps = [os.path.join('model-variances', f) for f in files if '_components.' in f]\n",
        "\n",
        "pa_ids = [re.search(r'([0-9]+)_.*?', f).group(1) for f in pas]\n",
        "comp_ids = [re.search(r'([0-9]+)_.*?', f).group(1) for f in pas]"
      ],
      "outputs": [],
      "execution_count": 14,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pa_dfs = list(map(pd.read_pickle, pas))\n",
        "comp_dfs = list(map(pd.read_pickle, comps))"
      ],
      "outputs": [],
      "execution_count": 15,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Pitch angle sample error"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We'll look at pitch angle first. The big question is whether our predicted error on pitch angle from the model fit ties in with our observed variation from the classification splits"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pa_details_df = pd.concat([\n",
        "    pd.DataFrame(\n",
        "        df.describe().loc[['mean', 'std']].values.T,\n",
        "        index=['{}-pa'.format(i), '{}-sigma_pa'.format(i)],\n",
        "        columns=['mean', 'std'],\n",
        "    )\n",
        "    for i, df in zip(pa_ids, pa_dfs)\n",
        "])\n",
        "pa_details_df.info()"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "<class 'pandas.core.frame.DataFrame'>\n",
            "Index: 196 entries, 587726031710912518-pa to 587741600410632336-sigma_pa\n",
            "Data columns (total 2 columns):\n",
            "mean    158 non-null float64\n",
            "std     150 non-null float64\n",
            "dtypes: float64(2)\n",
            "memory usage: 4.6+ KB\n"
          ]
        }
      ],
      "execution_count": 16,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pa_details_df.head()"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 17,
          "data": {
            "text/plain": [
              "                                  mean       std\n",
              "587726031710912518-pa        19.612688  2.068204\n",
              "587726031710912518-sigma_pa   1.153470  0.285118\n",
              "587725550661992687-pa        16.759515  0.609073\n",
              "587725550661992687-sigma_pa   0.743140  0.096737\n",
              "587725577498984562-pa              NaN       NaN"
            ],
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mean</th>\n",
              "      <th>std</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>587726031710912518-pa</th>\n",
              "      <td>19.612688</td>\n",
              "      <td>2.068204</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>587726031710912518-sigma_pa</th>\n",
              "      <td>1.153470</td>\n",
              "      <td>0.285118</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>587725550661992687-pa</th>\n",
              "      <td>16.759515</td>\n",
              "      <td>0.609073</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>587725550661992687-sigma_pa</th>\n",
              "      <td>0.743140</td>\n",
              "      <td>0.096737</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>587725577498984562-pa</th>\n",
              "      <td>NaN</td>\n",
              "      <td>NaN</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 17,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.hist(pa_details_df['std'].dropna().values[::2],\n",
        "         density=True, alpha=0.8, label='Sample error');\n",
        "plt.hist(pa_details_df['mean'].dropna().values[1::2],\n",
        "         density=True, alpha=0.8, label=r'Mean of errors on fit $\\sigma_{\\phi\\mathrm{\\ fit}}$');\n",
        "plt.legend();\n",
        "plt.xlabel('Degrees')\n",
        "plt.ylabel('Density')\n",
        "plt.savefig('duplicates_plots/spiral-error-comparison.pdf', bbox_inches='tight')"
      ],
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHl1JREFUeJzt3XmUVeWd7vHvQ4EgIBqpiiYyxgAChSCUyhWNpB0a0aCJzRUUh3ZAwSEdWzt2t0ttG1vTcWGiAYV4FZxIlBgvGtIqRg0oJhaFCoIDIEoZcxkMRgRF9Hf/OIedqqKGQ3F2Hap4PmvVqrPn34Za9dS7373frYjAzMwMoFWhCzAzs92HQ8HMzBIOBTMzSzgUzMws4VAwM7OEQ8HMzBIOBTMzSzgUzMws4VAwM7NE60IXsLOKi4ujR48ehS7DzKxZWbRo0fqIKGlovWYXCj169KC8vLzQZZiZNSuS3s1lPV8+MjOzhEPBzMwSDgUzM0s0uz4Fs93V559/TmVlJZ9++mmhS7E9WLt27ejSpQtt2rRp1PYOBbM8qaysZJ999qFHjx5IKnQ5tgeKCDZs2EBlZSU9e/Zs1D58+cgsTz799FM6d+7sQLCCkUTnzp13qbXqUDDLIweCFdqu/gw6FMzMLOE+BbOUfOeOBXnd3+OXH93gOjfddBMPPfQQRUVFtGrVimnTpnHkkUfmtY6qhg8fzq233kpZWVlqx7Cm5VAAmHbsjvMufr7p6zDbBQsXLuSJJ56goqKCtm3bsn79erZu3Vrosnbatm3baN26dZ3TuW5njeN/QbMW4oMPPqC4uJi2bdsCUFxcnCy78cYbefzxx9myZQtHHXUU06ZNQxLDhw/nsMMOY/78+XzyySfcd9993HzzzSxZsoQzzjiDSZMmsXr1akaMGMGQIUOoqKigf//+3HfffbRv377a8Z966imuv/56PvvsMw4++GDuvfdeOnbsWG2dlStXcumll7Ju3Trat2/Pz3/+cw455BDOO+882rVrx+LFixk2bBidOnVi5cqVrFq1im7dunHvvfcyYcIEysvLad26NZMnT+bb3/42M2bM4NFHH2XTpk188cUXPP+8/5jbVe5TMGshTjzxRNasWUPv3r2ZOHFitV+Ql112GS+//DJLly5ly5YtPPHEE8myvfbai/Lyci655BJOPfVUpkyZwtKlS5kxYwYbNmwA4M0332TixIksX76cTp06MXXq1GrHXr9+PZMmTWLevHlUVFRQVlbG5MmTd6hx/Pjx3HHHHSxatIhbb72ViRMnJssqKyt58cUXk+2WLVvGvHnzmDVrFlOmTEESS5YsYdasWZx77rnJHTYVFRXMnj3bgZAnDgWzFqJjx44sWrSI6dOnU1JSwhlnnMGMGTMAePbZZznyyCMZMGAAv/vd73j99deT7UaNGgXAgAED6N+/P1/72tdo27Yt3/jGN1izZg0AXbt2ZdiwYQCMGzeOBQuq95e89NJLLFu2jGHDhjFo0CBmzpzJu+9WH39t06ZNvPjii4wePZpBgwZx8cUX88EHHyTLR48eTVFRUbW69t57bwAWLFjAuHHjADjkkEPo3r07b731FgAnnHAC+++//y7/+1mGLx+ZtSBFRUUMHz6c4cOHM2DAAGbOnMmYMWOYOHEi5eXldO3alRtuuKHafezbLze1atUq+bx9etu2bcCOtznWnI4ITjjhBGbNmlVnbV9++SX77bcfr7zySq3LO3ToUO90XXJdz3KTWktB0j2S1kpaWsfysyS9JmmJpBclDUyrFrM9wZtvvsnbb7+dTL/yyit07949CYDi4mI2bdrE7Nmzd3rf7733HgsXLgTgoYce4uijq98JNXToUF544QVWrFgBwCeffJL8Jb9dp06d6NmzJ4888giQCZJXX301p+Mfc8wxPPjggwC89dZbvPfee/Tp02enz8MalmZLYQbwM+C+Opa/AxwbEX+RdBIwHUjv3jmzJpbLLaT5tGnTJi6//HI2btxI69at+eY3v8n06dPZb7/9uOiiiygtLeXAAw/k8MMP3+l99+nThylTpnD++efTr18/JkyYUG15SUkJM2bMYOzYsXz22WcATJo0id69e1db78EHH2TChAlMmjSJzz//nDFjxjBwYMN/D06cOJEJEyYwYMAAWrduzYwZM6q1aix/FBHp7VzqATwREaUNrPcVYGlEHNTQPsvKyiLvL9nxLamWB8uXL6dv376FLiPvVq9ezSmnnMLSpbU2+m03VNvPoqRFEdHgAyW7S0fzBcBv61ooabykcknl69ata8KyzMz2LAUPBUnfJhMKP6xrnYiYHhFlEVFWUtLgK0bNLI969OjhVsIepKB3H0k6FLgbOCkiNhSyFjMzK2BLQVI34FHg7Ih4q6H1zcwsfam1FCTNAoYDxZIqgeuBNgARcRdwHdAZmJq953lbLp0gZmaWntRCISLGNrD8QuDCtI5vZmY7r+AdzWZmtvtwKJi1MJKScYIgM6R0SUkJp5xySkHquf322+nbty9nnXVWQY6fpprndtRRRwGwcePGHQYNbC489pFZWmp7KHJX5PhAZYcOHZLRUPfee2+efvppDjqowedCUzN16lTmzZtHly5ddnrbiCAiaNWqVa3TuW6Xlprn9uKLLwJ/C4Wqo8A2F24pmLVAI0eO5De/+Q0As2bNYuzYv3XxPfDAAxxxxBHJSKVffPEFAKeddhpDhgyhf//+TJ8+Hcg8zdy3b18uuugi+vfvz4knnsiWLVtqPebkyZMpLS2ltLSUn/zkJwBccsklrFq1ipNOOonbbrtth21qq2X16tX06dOHc845h9LSUubPn19tes2aNbUeq+Z2a9as4eSTT2bgwIGUlpbyy1/+Mue6cznv2s5t+/sjrrnmGlauXMmgQYO4+uqrdzjmq6++yre+9S369etHq1atkMR1111Xa31NbnuiNpevIUOGRN7d9a0dv8x20rJly6rPqO3nale+ctShQ4d49dVX4/TTT48tW7bEwIED49lnn42TTz45li1bFqecckps3bo1IiImTJgQM2fOjIiIDRs2RETE5s2bo3///rF+/fp45513oqioKBYvXhwREaNHj477779/h2OWl5dHaWlpbNq0KT7++OPo169fVFRURERE9+7dY926dbX+e9VWyzvvvBOSYuHChRERO0zXdaya682ePTsuvPDC5HgbN27Mue5cz7vmuXXo0CGpuX///rX+/2zZsiX69OkTf/jDHyIi4tprr42rrroqvvzyy1rXb4wdfhYjAiiPHH7HuqVg1gIdeuihrF69mlmzZjFy5Mhk/jPPPMOiRYs4/PDDGTRoEM888wyrVq0CMtfHBw4cyNChQ1mzZk0y4mrPnj0ZNGgQAEOGDGH16tU7HG/BggV897vfpUOHDnTs2JHvfe97zJ8/v94a66ule/fuDB06NFm36nR9x6q63oABA3j66af54Q9/yPz589l33313qu5czrsx5s2bx+DBgzniiCOAzP/Vhx9+mAxHfuONNyb/DtutXLmSUaNGMXPmTObPn8+MGTOqvSgpn9ynYNZCjRo1iquuuornnnsueYNaRHDuuedy8803V1v3ueeeY968eSxcuJD27dszfPjwZMjtqqORFhUV1Xn5aGfVVcvq1avz8m6F3r17U1FRwdy5c7n22ms57rjjduoSTVrnvXTpUgYMGJBMV1RUMHjwYLZs2cKVV17J888/zx//+EceeughOnXqBMCdd97Jtm3beP/99+ncuTMLFixg8+bNAHm/gcAtBbMW6vzzz+f666+v9gvouOOOY/bs2axduxaADz/8kHfffZePPvqIr3zlK7Rv35433niDl156aaeOdcwxx/DYY4+xefNmPvnkE379619zzDHH1LtNXbXk61h/+tOfaN++PePGjePqq6+moqIiL3XnYp999uHjjz+udVnnzp157bXXgMy7IR599FHGjBnD1KlTOfvssznyyCM5/vjjk3Ug84t/5MiRfP3rXwfg6KOP5swzz0zljrI9q6WQ77tBzHZjXbp04Yorrqg2r1+/fkyaNIkTTzyRL7/8kjZt2jBlyhRGjBjBXXfdRd++fenTp0+1Sze5GDx4MOedd15ySeTCCy/ksMMOq3ebumo58MADG3Wsmpd3lixZwtVXX02rVq1o06YNd955Z6P3tbM6d+7MsGHDKC0t5aSTTuLHP/5xsmzs2LHMmTOH0tJSiouLmTVrFp07d2bx4sWMGjWKrl27snTpUs4555w695/mXVWpvk8hDbv0PoWdCQW/T8F2Ukt9n4I1jUceeYTbb7+dNm3acP7551d71uS5555j6dKldOzYkeLiYrp27cpNN93EmWeeyWmnnbbDvnblfQp7VkvBzGw3NXr0aBYtWsQPfvADDjjggGrLtr93u6qHH344lTocCmZmu4lbbrml0CW4o9nMzP7GoWBmZgmHgpmZJRwKZmaWcCiY5VFzu8XbWp5d/Rl0KJjlSbt27diwYYODwQomItiwYQPt2rVr9D58S+rOqO3hNz/kZlldunShsrKSdevWFboU24O1a9euUe+u2M6hYJYnbdq0oWfPnoUuw2yX+PKRmZklHApmZpZwKJiZWcKhYGZmidRCQdI9ktZKWlrHckm6XdIKSa9JGpxWLWZmlps0WwozgBH1LD8J6JX9Gg/s+AYMMzNrUqmFQkT8HviwnlVOBe6LjJeA/SR9La16zMysYYXsUzgIWFNlujI7bweSxksql1TuB4PMzNLTLDqaI2J6RJRFRFlJSUmhyzEza7EKGQrvA12rTHfJzjMzswIpZCjMAc7J3oU0FPgoIj4oYD1mZnu81MY+kjQLGA4US6oErgfaAETEXcBcYCSwAtgM/GNatZiZWW5SC4WIGNvA8gAuTev4Zma285pFR7OZmTUNh4KZmSUcCmZmlnAomJlZwqFgZmYJh4KZmSUcCmZmlkjtOYVmb9qxha7AzKzJuaVgZmYJh4KZmSUcCmZmlnAomJlZwqFgZmYJh4KZmSUcCmZmlnAomJlZwqFgZmYJh4KZmSUcCmZmlnAomJlZwqFgZmYJh4KZmSUcCmZmlnAomJlZItVQkDRC0puSVki6ppbl3SQ9K2mxpNckjUyzHjMzq19qoSCpCJgCnAT0A8ZK6ldjtWuBhyPiMGAMMDWteszMrGFpthSOAFZExKqI2Ar8Aji1xjoBdMp+3hf4U4r1mJlZA9IMhYOANVWmK7PzqroBGCepEpgLXF7bjiSNl1QuqXzdunVp1GpmZhS+o3ksMCMiugAjgfsl7VBTREyPiLKIKCspKWnyIs3M9hRphsL7QNcq012y86q6AHgYICIWAu2A4hRrMjOzeqQZCi8DvST1lLQXmY7kOTXWeQ84DkBSXzKh4OtDZmYFklooRMQ24DLgSWA5mbuMXpd0o6RR2dX+GbhI0qvALOC8iIi0ajIzs/q1TnPnETGXTAdy1XnXVfm8DBiWZg1mZpa7Qnc0m5nZbiSnUJD0qKSTa7szyMzMWo5cf8lPBc4E3pZ0i6Q+KdZkZmYFklMoRMS8iDgLGAysBuZJelHSP0pqk2aBZmbWdHK+HCSpM3AecCGwGPgpmZB4OpXKzMysyeV095GkXwN9gPuB70TEB9lFv5RUnlZxZmbWtHK9JfXn2dtLE5LaRsRnEVGWQl1mZlYAuV4+mlTLvIX5LMTMzAqv3paCpAPJjGy6t6TDAGUXdQLap1ybmZk1sYYuH/09mc7lLsDkKvM/Bv4tpZrMzKxA6g2FiJgJzJR0ekT8qolqMjOzAmno8tG4iHgA6CHpyprLI2JyLZuZmVkz1dDlow7Z7x3TLsTMzAqvoctH07Lf/6NpyjEzs0LKdUC8/5bUSVIbSc9IWidpXNrFmZlZ08r1OYUTI+KvwClkxj76JnB1WkWZmVlh5BoK2y8znQw8EhEfpVSPmZkVUK7DXDwh6Q1gCzBBUgnwaXplmZlZIeQ6dPY1wFFAWUR8DnwCnJpmYWZm1vR25h3Nh5B5XqHqNvfluR4zMyugXIfOvh84GHgF+CI7O3AomJm1KLm2FMqAfhERaRZjZmaFlevdR0uBA9MsxMzMCi/XlkIxsEzSH4HPts+MiFGpVGVmZgWRayjc0JidSxpB5l3ORcDdEXFLLev87+z+A3g1Is5szLHMzGzX5RQKEfG8pO5Ar4iYJ6k9mV/0dZJUBEwBTgAqgZclzYmIZVXW6QX8KzAsIv4i6auNPREzM9t1uY59dBEwG5iWnXUQ8FgDmx0BrIiIVRGxFfgFOz7bcBEwJSL+AhARa3Mt3MzM8i/XjuZLgWHAXwEi4m2gob/qDwLWVJmuzM6rqjfQW9ILkl7KXm4yM7MCybVP4bOI2CplXtGcfYAtH7entgZ6AcPJvPLz95IGRMTGqitJGg+MB+jWrVseDmtmZrXJtaXwvKR/A/aWdALwCPB4A9u8D3StMt0lO6+qSmBORHweEe8Ab5EJiWoiYnpElEVEWUlJSY4lm5nZzso1FK4B1gFLgIuBucC1DWzzMtBLUk9JewFjgDk11nmMTCsBScVkLietyrEmMzPLs1zvPvpS0mPAYxGxLsdttkm6DHiSzJ1K90TE65JuBMojYk522YmSlpEZPuPqiNjQqDMxM7NdVm8oKNOJcD1wGdlWhaQvgDsi4saGdh4Rc8m0KqrOu67K5wCuzH6ZmVmBNXT56Adk7jo6PCL2j4j9gSOBYZJ+kHp1ZmbWpBoKhbOBsdlOYAAiYhUwDjgnzcLMzKzpNRQKbSJifc2Z2X6FNumUZGZmhdJQKGxt5DIzM2uGGrr7aKCkv9YyX0C7FOoxM7MCqjcUIqLeQe/MzKxlyfXhNTMz2wM4FMzMLOFQMDOzhEPBzMwSDgUzM0s4FMzMLOFQMDOzhEPBzMwSDgUzM0s4FMzMLOFQMDOzhEPBzMwSDgUzM0s4FMzMLOFQMDOzhEPBzMwSDgUzM0s4FMzMLJFqKEgaIelNSSskXVPPeqdLCklladZjZmb1Sy0UJBUBU4CTgH7AWEn9allvH+D7wB/SqsXMzHKTZkvhCGBFRKyKiK3AL4BTa1nvP4EfAZ+mWIuZmeUgzVA4CFhTZboyOy8haTDQNSJ+k2IdZmaWo9aFOrCkVsBk4Lwc1h0PjAfo1q1buoXtrGnH1j7/4uebtg4zszxIs6XwPtC1ynSX7Lzt9gFKgeckrQaGAnNq62yOiOkRURYRZSUlJSmWbGa2Z0szFF4GeknqKWkvYAwwZ/vCiPgoIoojokdE9ABeAkZFRHmKNZmZWT1SC4WI2AZcBjwJLAcejojXJd0oaVRaxzUzs8ZLtU8hIuYCc2vMu66OdYenWYuZmTXMTzSbmVnCoWBmZgmHgpmZJRwKZmaWcCiYmVnCoWBmZgmHgpmZJRwKZmaWcCiYmVnCoWBmZgmHgpmZJRwKZmaWcCiYmVnCoWBmZgmHgpmZJRwKZmaWcCiYmVki1Tev7dGmHbvjvIufb/o6zMx2glsKZmaWcCiYmVnCoWBmZgmHgpmZJRwKZmaW2CPvPnp77aaCHLdXQY5qZpa7VFsKkkZIelPSCknX1LL8SknLJL0m6RlJ3dOsx8zM6pdaS0FSETAFOAGoBF6WNCcillVZbTFQFhGbJU0A/hs4I62aCu07dywoyHEfv/zoghzXzJqfNFsKRwArImJVRGwFfgGcWnWFiHg2IjZnJ18CuqRYj5mZNSDNUDgIWFNlujI7ry4XAL9NsR4zM2vAbtHRLGkcUAbUMjYESBoPjAfo1q1bE1ZmZrZnSbOl8D7Qtcp0l+y8aiQdD/w7MCoiPqttRxExPSLKIqKspKQklWLNzCzdUHgZ6CWpp6S9gDHAnKorSDoMmEYmENamWIuZmeUgtctHEbFN0mXAk0ARcE9EvC7pRqA8IuYAPwY6Ao9IAngvIkalVdOeqlB3PYHvfDJrblLtU4iIucDcGvOuq/L5+DSPb2ZmO2e36Gi2lsvPZpg1Lx77yMzMEg4FMzNLOBTMzCzhUDAzs4RDwczMEr77qAlN3vj9HeZdud9PC1CJmVnt3FIwM7OEQ8HMzBIOBTMzSzgUzMws4VAwM7OEQ8HMzBIOBTMzSzgUzMws4VAwM7OEQ8HMzBIe5qKZqG2IDPAwGWaWX24pmJlZwqFgZmYJXz6yFsnvhjZrHLcUzMws4VAwM7OEQ8HMzBKp9ilIGgH8FCgC7o6IW2osbwvcBwwBNgBnRMTqNGsya6ncj2L5kFpLQVIRMAU4CegHjJXUr8ZqFwB/iYhvArcBP0qrHjMza1iaLYUjgBURsQpA0i+AU4FlVdY5Fbgh+3k28DNJiohIsS6z1BTqr/VC2tNaKIX8P26Kc06zT+EgYE2V6crsvFrXiYhtwEdA5xRrMjOzejSL5xQkjQfGZyc/k7S0kPXk1zHbPxQD6+ta64mGt9+d1XtuzVhLPS9oBuemKxq12W5/XvWp55xzOa/uuRwjzVB4H+haZbpLdl5t61RKag3sS6bDuZqImA5MB5BUHhFlqVRcQC31vKDlnltLPS9ouefm82pYmpePXgZ6SeopaS9gDDCnxjpzgHOzn/8B+J37E8zMCie1lkJEbJN0GfAkmVtS74mI1yXdCJRHxBzg/wD3S1oBfEgmOMzMrEBS7VOIiLnA3Brzrqvy+VNg9E7udnoeStsdtdTzgpZ7bi31vKDlnpvPqwHy1RozM9vOw1yYmVmiWYWCpBGS3pS0QtI1ha4nHyR1lfSspGWSXpdU+yvWmilJRZIWS6r7rtpmSNJ+kmZLekPSckn/q9A15YOkH2R/DpdKmiWpXaFraixJ90haW/UWdkn7S3pa0tvZ718pZI2NUcd5/Tj7s/iapF9L2q+x+282oZDjsBnN0TbgnyOiHzAUuLSFnNd23weWF7qIFPwU+J+IOAQYSAs4R0kHAVcAZRFRSuYGkeZ888cMYESNedcAz0REL+CZ7HRzM4Mdz+tpoDQiDgXeAv61sTtvNqFAlWEzImIrsH3YjGYtIj6IiIrs54/J/HKp+eR3sySpC3AycHeha8knSfsC3yJz9xwRsTUiNha2qrxpDeydfW6oPfCnAtfTaBHxezJ3NVZ1KjAz+3kmcFqTFpUHtZ1XRDyVHRUC4CUyz4U1SnMKhVyGzWjWJPUADgP+UNhK8uYnwL8AXxa6kDzrCawD7s1eGrtbUodCF7WrIuJ94FbgPeAD4KOIeKqwVeXdARHxQfbzn4EDCllMSs4HftvYjZtTKLRokjoCvwL+KSL+Wuh6dpWkU4C1EbGo0LWkoDUwGLgzIg4DPqF5XoaoJnt9/VQyofd1oIOkcYWtKj3ZB2Vb1O2Xkv6dzCXpBxu7j+YUCrkMm9EsSWpDJhAejIhHC11PngwDRklaTeZS399JeqCwJeVNJVAZEdtbdLPJhERzdzzwTkSsi4jPgUeBowpcU779P0lfA8h+X1vgevJG0nnAKcBZuzIyRHMKhVyGzWh2JInMtenlETG50PXkS0T8a0R0iYgeZP6vfhcRLeKvzoj4M7BGUp/srOOoPiR8c/UeMFRS++zP5XG0gA70GqoOrXMu8H8LWEveZF9o9i/AqIjYvCv7ajahkO1E2T5sxnLg4Yh4vbBV5cUw4Gwyf0m/kv0aWeiirEGXAw9Keg0YBPxXgevZZdmWz2ygAlhC5vdDs30CWNIsYCHQR1KlpAuAW4ATJL1NpmV0S3372B3VcV4/A/YBns7+Drmr0fv3E81mZrZds2kpmJlZ+hwKZmaWcCiYmVnCoWBmZgmHgpmZJVJ9yY7Z7k7SF2Ruv2xD5knQ+4DbIqKlDc1hlhOHgu3ptkTEIABJXwUeAjoB1+/qjiUVRcQXu7ofs6bky0dmWRGxFhgPXKaMouw49S9nx6m/GEBSK0lTs+PXPy1prqR/yC5bLelHkiqA0ZIOlvQ/khZJmi/pkOx6JZJ+ld33y5KGZecfW+UhxsWS9inQP4ftodxSMKsiIlZl393xVTKDw30UEYdLagu8IOkpYAjQg8x7Pb5K5gn7e6rsZkNEDAaQ9AxwSUS8LelIYCrwd2Tex3BbRCyQ1I3Mk/p9gauASyPihewgiZ+mf9Zmf+NQMKvbicCh21sBwL5AL+Bo4JFsv8OfJT1bY7tfQjLy7VHAI5mhhABom/1+PNCvyvxO2fVfACZLehB4NCIq839aZnVzKJhVIekbwBdkRs8UcHlEPFljnYbGpvok+70VsHF7n0UNrYChEVGzJXCLpN8AI8m0TP4+It7Y2fMwayz3KZhlSSoB7gJ+lh16+ElgQnZocyT1zr5M5wXg9GzfwgHA8Nr2l30vxjuSRme3l6SB2cVPkRlUb/uxt3d2HxwRSyLiR2RGBj4khVM1q5NDwfZ0e2c7dV8H5pH5Zf0f2WV3kxkSuyL7kvRpZFrXvyLzToVlwANkRhX9qI79nwVcIOlV4HX+9grZK4CybAf2MuCS7Px/krQ0O/rq5+zCG7TMGsOjpJo1gqSOEbFJUmfgj8Cw7HsWzJo19ymYNc4TkvYD9gL+04FgLYVbCmZmlnCfgpmZJRwKZmaWcCiYmVnCoWBmZgmHgpmZJRwKZmaW+P8433drRQNlxQAAAABJRU5ErkJggg==\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 18,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's make a histogram of the mean fitted standard deviation minus the sample standard deviation of pitch angles\n",
        "\n",
        "$$\\overline{\\sigma_{\\phi\\mathrm{\\ fit}}} - \\mathrm{s}_\\phi$$"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(5, 3), dpi=100)\n",
        "plt.hist(\n",
        "    pd.Series(pa_details_df['mean'].values[1::2] - pa_details_df['std'].values[::2]).dropna().values,\n",
        "    bins=20, density=True,\n",
        ")\n",
        "plt.xlabel(r'$\\overline{\\sigma_{\\phi\\mathrm{\\ fit}}} - \\mathrm{s}_\\phi$')\n",
        "plt.ylabel('Density');"
      ],
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 500x300 with 1 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAdcAAAEsCAYAAACR0fkbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHdxJREFUeJzt3X20HXV97/H3B2qoIgEVSHhsLFpppVTAi419uFGKUula19parK2K3NpCbZWLT2QpIvEB1DbLrqK25doIXa1SlUIrFWisuqTSIGCAilgRoogkSkXCY7iS7/1j5uju5uyTc/aZk332Oe/XWrM4M/Ob3/5OSPLJb/ZvZlJVSJKk7uwy6gIkSVpoDFdJkjpmuEqS1DHDVZKkjhmukiR1zHCVJKljhqskSR0zXCVJ6pjhKklSxwxXSZI6Ni/CNcmrk2xK8lCSDUmOnqLti5Jck+T7Se5PsjHJy/raJMmaJHcmeTDJ+iRPnfszkSRpHoRrkhOAtcBZwJHA9cDlSfYdcMj3gHcCK4HDgXXAuiTP72nzRuA1wMnAs4D72z5/fE5OQpKkHhn1g/uTbAC+WFV/1K7vAtwO/HlVnTPNPq4DLq2qM5IE+Dbwp1X1J+3+PYEtwIlV9dFp9Bdgf+DeYc5JkrSg7AF8u2YQmD82h8XsUJIlwFHA2RPbqmp7kvU0I9MdHR/gucDTgDe1m58MLAfW9/R5TxviK4FHhWuS3YDdejbtB9w80/ORJC1YBwJ3TLfxSMMV2BvYlWZU2WsLcOigg9qR6B00gfgI8IdV9S/t7uU9ffT3uZzJrQbO7N94++23s3Tp0qnqlyQtYFu3buWggw6CGV7JHHW4Dute4BnA44FjgLVJbq2qzw7Z39k03/tO2AP41tKlSw1XSdKMjTpc76IZeS7r274M2DzooKraDtzSrm5M8tM0o8/P9hy3DLizr8+NA/rbBmybWG+uNkuSNJyRzhauqoeBa2lGn8APJzQdA1w1g6524Uffmd5GE7C9fS6lmTU8kz4lSRrKqEeu0FyOPT/JNcDVwKnA7jS32JDkAuCOqlrdrq8GrgG+ThOoLwBeBpwCUFWV5H3AW5J8jSZs304zg/jinXhekqRFauThWlUXJtkHWEMz4WgjcFxVTUxIOhjY3nPI7sAHaGZuPUgzq/d3q+rCnjbvadv9FbAXcGXb50NzeS6SJME8uM91PmovI99zzz33OKFJkhaxrVu3sueeewLsWVVbp3vcyJ/QJEnSQmO4SpLUsZF/5ypJGr0Vp1/aST+bzjm+k37GnSNXSZI6ZrhKktQxw1WSpI4ZrpIkdcxwlSSpY4arJEkdM1wlSeqY4SpJUscMV0mSOma4SpLUMcNVkqSOGa6SJHXMcJUkqWOGqyRJHTNcJUnqmOEqSVLHDFdJkjr2Y6MuQJK0cKw4/dJZ97HpnOM7qGS0HLlKktQxw1WSpI4ZrpIkdcxwlSSpY4arJEkdmxfhmuTVSTYleSjJhiRHT9H2VUk+n+Tudlnf3z7Jh5NU33LZ3J+JJEnzIFyTnACsBc4CjgSuBy5Psu+AQ1YBHwGeA6wEbgeuSHJAX7vLgP16lt/uvHhJkiYx8nAFTgPOq6p1VXUTcDLwAHDSZI2r6neq6gNVtbGqbgZ+j+Y8julruq2qNvcsd8/lSUiSNGGk4ZpkCXAUsH5iW1Vtb9dXTrObxwGPAb7Xt31Vku8k+WqSDyZ5Uhc1S5K0I6N+QtPewK7Alr7tW4BDp9nHu4Fv0xPQNJeELwJuAw4B3gV8KsnKqnqkv4MkuwG79WzaY5qfLUnSo4w6XGclyenAS4BVVfXQxPaq+mhPsxuT3AB8neb72k9P0tVq4Mw5LFWStIiM+jvXu4BHgGV925cBm6c6MMnrgdOB51XVDVO1rapb2896yoAmZwN79iwH7rBySZIGGGm4VtXDwLX0TEZKMjE56apBxyV5I3AGcFxVXbOjz0lyIPAk4M4BdWyrqq0TC3DvjE5EkqQe8+Gy8Frg/CTXAFcDpwK7A+sAklwA3FFVq9v1NwFrgJcCm5Isb/u5r6ruS/J4mku8n6AZ/R4CvAe4Bbh8p52VJGnRGnm4VtWFSfahCczlwEaaEenEJKeDge09h5wCLAE+3tfVWcDbaC4zHw68AtiLZrLTFcAZVbVtjk5DkqQfGnm4AlTVucC5A/at6ltfsYO+HgSe31VtkiTN1KgnNEmStODMi5GrJGk4K06/dNQlaBKOXCVJ6pjhKklSxwxXSZI6ZrhKktQxw1WSpI4ZrpIkdcxwlSSpY4arJEkdM1wlSeqY4SpJUscMV0mSOma4SpLUMcNVkqSOGa6SJHXMcJUkqWOGqyRJHTNcJUnq2I+NugBJmo4Vp1/aST+bzjm+k36kqThylSSpY4arJEkdM1wlSeqY4SpJUscMV0mSOma4SpLUMcNVkqSOzYtwTfLqJJuSPJRkQ5Kjp2j7qiSfT3J3u6zvb5/GmiR3JnmwbfPUuT8TSZLmQbgmOQFYC5wFHAlcD1yeZN8Bh6wCPgI8B1gJ3A5ckeSAnjZvBF4DnAw8C7i/7fPH5+IcJEnqNfJwBU4DzquqdVV1E00gPgCcNFnjqvqdqvpAVW2sqpuB36M5j2OgGbUCpwLvqKpLquoG4OXA/sAL5/50JEmL3UjDNckS4Chg/cS2qtrerq+cZjePAx4DfK9dfzKwvK/Pe4ANg/pMsluSpRMLsMcMT0WSpB8a9ch1b2BXYEvf9i00ATkd7wa+zY/CdOK4mfS5GrinZ/nWND9bkqRHGXW4zkqS04GXAL9eVQ/NoquzgT17lgM7KE+StEiN+q04dwGPAMv6ti8DNk91YJLXA6cDv9J+rzph4rhlwJ19fW6crK+q2gZs6+l7OrVLkjSpkY5cq+ph4FrayUgASSYmJ1016LgkbwTOAI6rqmv6dt9GE7C9fS6lmTU8sE9Jkroy6pErNLfhnJ/kGuBqmpm+uwPrAJJcANxRVavb9TcBa4CXApuSTHyPel9V3VdVleR9wFuSfI0mbN9O873sxTvxvCRJi9TIw7WqLkyyD01gLqe5dHtcVU1MSDoY2N5zyCnAEuDjfV2dBbyt/fk9NAH9V8BewJVtn7P5XlaSpGkZebgCVNW5wLkD9q3qW18xjf4KeGu7SFKnVpx+aSf9bDrn+E760fwz1rOFJUmajwxXSZI6ZrhKktQxw1WSpI4ZrpIkdcxwlSSpY0OFa5Kf7LoQSZIWimFHrrck+UyS3/UF5JIk/XfDhuuRwA00jy7cnOQvkxzdXVmSJI2vocK1qjZW1WuB/YGTgP2AK5P8R5LT2scZSpK0KM1qQlNV/aCqLgJeDLwJeArwJ8DtSS5Isl8HNUqSNFZmFa5JnpnkAzTvTT2NJlgPAY6lGdVeMusKJUkaM0M9uD/JacArgacB/wy8HPjnqpp4e81tSU4ENnVQoyRJY2XYt+KcAvw18OGqunNAm+8A/3vI/iVJGlvDhuuxwDd7RqoAJAlwUFV9s6oeBs6fbYGSJI2bYb9z/Tqw9yTbnwjcNnw5kiSNv2HDNQO2Px54aMg+JUlaEGZ0WTjJ2vbHAtYkeaBn967As4CNHdUmSZ1bcfqloy7hh+ZTLerWTL9zPaL9b4CfBR7u2fcwcD3N7TiSJC1aMwrXqnoOQJJ1wGurauucVCVJ0hgbarZwVb2y60IkSVooph2uSS4CTqyqre3PA1XVi2ZdmSRJY2omI9d7aCYyTfwsSZImMe1w7b0U7GVhSZIGG+o+1ySPTfK4nvWfSHJqkud1V5okSeNp2IdIXELzsH6S7AVcDbwOuCTJKR3VJknSWBo2XI8EPt/+/JvAZuAnaAL3NTPpKMmrk2xK8lCSDUmOnqLt05N8om1fSU6dpM3b2n29y80zqUmSpNkYNlwfB9zb/vw84KL2If7/ThOy05LkBGAtcBZNYF8PXJ5k3yk+91bgdJpAH+TLwH49yy9OtyZJkmZr2HC9BXhhkoOA5wNXtNv3BWbyYInTgPOqal1V3QScDDwAnDRZ46r6YlW9oao+Cmybot8fVNXmnuWuGdQkSdKsDBuua2gec7gJ2FBVV7Xbnwd8aTodJFkCHAWsn9jWjn7XAyuHrGvCU5N8O8mtSf42ycGz7E+SpGkb9glNH09yJc0l1+t7dn0a+IdpdrM3zcP+t/Rt3wIcOkxdrQ3AicBX2/rOBD6f5LCquneyA5LsBuzWs2mPWXy+JGmRG/Zl6VTVZvq+96yqq2dd0SxV1ad6Vm9IsgH4BvBbwIcGHLaaJoQlSZq1ocI1ye40k4qOofme9b9dXq6qn5xGN3cBjwDL+rYvY+rJSjNSVd9P8p/AU6ZodjbNxKoJewDf6qoGSdLiMuzI9f8C/xP4G+BOfvRYxGmrqoeTXEsT0BcDJNmlXT93yLoeJcnjgUPaWgfVso2eCVLJoHfBS5K0Y8OG668Cx1fVv83y89cC5ye5huZBFKcCuwPrAJJcANxRVavb9SXAz7THLgEOSPIM4L6quqVt8yfAP9FcCt6f5jafR4CPzLJWSZKmZdhwvRv43mw/vKouTLIPzezj5cBG4LiqmpjkdDCwveeQ/fnvs5Ff3y6fA1a12w6kCdInAd8FrgR+vqq+O9t6JUmajmHD9QxgTZJXVNUDsymgqs5lwGXgqlrVt74JmPKabVW9ZDb1SJI0W8OG6+tovsfckmQT8P96d1bVkbOsS5KksTVsuF7caRWSJC0gwz5E4qyuC5EkaaEY9vGHJNkrye8lOTvJE9ttRyY5oLvyJEkaP8M+ROJwmmcA3wOsAM6jmT38IpoZvi/vqD5JksbOsCPXtcCHq+qpwEM92/8Z+OVZVyVJ0hgbNlz/B/CXk2y/g+Z+VUmSFq1hw3UbsHSS7T9F8+AGSZIWrWHD9R+BtyZ5TLte7TtT3w18opPKJEkaU8OG6+uAx9OMUh9L8/jBW4B7gTd3U5okSeNp2Ptc7wGOTfILwM/RBO11VbW+y+IkSRpHMw7X9rVwJ9LcdrOC5nVztwGbk6SqZvz6OUmSFpIZXRZO86LTf6R5n+sBwI3Al4GfAD4M/EPH9UmSNHZmOnI9keY+1mOq6jO9O5I8F7g4ycur6oKO6pMkaezMdELTbwPv6g9WgKr6V+Ac4He6KEySpHE103A9HLhsiv2fopngJEnSojXTcH0isGWK/VuAJwxfjiRJ42+m4bor8IMp9j/C8O+IlSRpQZhpEAb4cJJtA/bvNst6JEkaezMN1/On0caZwpKkRW1G4VpVr5yrQiRJWiiGfbawJEkawHCVJKljhqskSR0zXCVJ6pjhKklSxwxXSZI6NvJwTfLqJJuSPJRkQ5Kjp2j79CSfaNtXklNn26ckSV0babgmOQFYC5wFHAlcD1yeZN8BhzwOuBU4HdjcUZ+SJHVq1CPX04DzqmpdVd0EnAw8AJw0WeOq+mJVvaGqPgoMegTjjPqUJKlrIwvXJEuAo4D1E9uqanu7vnJn9plktyRLJxZgj2E+X5IkGO3IdW+at+z0v8JuC7B8J/e5GrinZ/nWkJ8vSdLILwvPF2cDe/YsB462HEnSOBvlu1fvonn/67K+7csYMFlprvqsqm30fIebZMiPlyRphCPXqnoYuBY4ZmJbkl3a9avmS5+SJM3UKEeu0Nwyc36Sa4CrgVOB3YF1AEkuAO6oqtXt+hLgZ9pjlwAHJHkGcF9V3TKdPiVJmmsjDdequjDJPsAamglHG4HjqmpiQtLBwPaeQ/YHvtSz/vp2+Rywapp9SpI0p0Y9cqWqzgXOHbBvVd/6JmCHX4hO1ackSXPN2cKSJHXMcJUkqWOGqyRJHTNcJUnqmOEqSVLHDFdJkjpmuEqS1DHDVZKkjhmukiR1zHCVJKljhqskSR0zXCVJ6pjhKklSxwxXSZI6NvJXzkkLxYrTL+2kn03nHN9JP/NFV78u0jhx5CpJUscMV0mSOma4SpLUMcNVkqSOGa6SJHXMcJUkqWOGqyRJHTNcJUnqmOEqSVLHDFdJkjpmuEqS1LF5Ea5JXp1kU5KHkmxIcvQO2r84yc1t+xuTvKBv/4eTVN9y2dyehSRJjZGHa5ITgLXAWcCRwPXA5Un2HdD+2cBHgA8BRwAXAxcnOayv6WXAfj3Lb8/JCUiS1Gfk4QqcBpxXVeuq6ibgZOAB4KQB7V8LXFZV762qr1TVGcB1wB/1tdtWVZt7lrvn7AwkSeox0nBNsgQ4Clg/sa2qtrfrKwcctrK3fevySdqvSvKdJF9N8sEkT+qobEmSpjTq97nuDewKbOnbvgU4dMAxywe0X96zfhlwEXAbcAjwLuBTSVZW1SP9HSbZDditZ9Me0z0BSZL6jTpc50RVfbRn9cYkNwBfB1YBn57kkNXAmTuhNEnSIjDq71zvAh4BlvVtXwZsHnDM5hm2p6pubT/rKQOanA3s2bMcOGXVkiRNYaThWlUPA9cCx0xsS7JLu37VgMOu6m3fOnaK9iQ5EHgScOeAOrZV1daJBbh32ichSVKf+XBZeC1wfpJrgKuBU4HdgXUASS4A7qiq1W37PwM+l+R1wKXAS4BnAr/ftn88zSXeT9CMZg8B3gPcQjPxSZKkOTXycK2qC5PsA6yhmZS0ETiuqiYmLR0MbO9p/4UkLwXeQTNR6WvAC6vqP9omjwCHA68A9gK+DVwBnFFV23bCKUmSFrmRhytAVZ0LnDtg36pJtn0M+NiA9g8Cz++yPkmSZmLUE5okSVpwDFdJkjpmuEqS1DHDVZKkjhmukiR1zHCVJKljhqskSR0zXCVJ6pjhKklSxwxXSZI6ZrhKktQxw1WSpI4ZrpIkdWxevBVnIVtx+qWz7mPTOcd3UMnC5K/v5Lr4dZFGpavfv6P8s+3IVZKkjhmukiR1zHCVJKljhqskSR0zXCVJ6pjhKklSxwxXSZI6ZrhKktQxw1WSpI4ZrpIkdcxwlSSpY4arJEkdM1wlSerYvAjXJK9OsinJQ0k2JDl6B+1fnOTmtv2NSV7Qtz9J1iS5M8mDSdYneercnoUkSY2Rh2uSE4C1wFnAkcD1wOVJ9h3Q/tnAR4APAUcAFwMXJzmsp9kbgdcAJwPPAu5v+/zxuToPSZImjDxcgdOA86pqXVXdRBOIDwAnDWj/WuCyqnpvVX2lqs4ArgP+CJpRK3Aq8I6quqSqbgBeDuwPvHCOz0WSpNG+LD3JEuAo4OyJbVW1Pcl6YOWAw1bSjHR7Xc6PgvPJwHJgfU+f9yTZ0B770Unq2A3YrWfTHgBbt26dyelMavu2B2bdRxd1LFTz6de3i1pg/vy+k8ZdF3+Whu1jpOEK7A3sCmzp274FOHTAMcsHtF/es58dtOm3Gjizf+NBBx00oPnOtef7Rl3Bwjbffn3nWz3SuOr4z9IewLSTdtThOl+czaNHw08EvrcTa9gD+BZwIHDvTvzc+WAxnzss7vNfzOcOi/v8x+nc9wC+PZMDRh2udwGPAMv6ti8DNg84ZvMO2m/u2XZnX5uNk3VYVduAbX2bd+q12OarYgDurapFdR14MZ87LO7zX8znDov7/Mfs3Gdc30gnNFXVw8C1wDET25Ls0q5fNeCwq3rbt47taX8bTcD29rmUZtbwoD4lSerMqEeu0FyOPT/JNcDVNDN9dwfWASS5ALijqla37f8M+FyS1wGXAi8Bngn8PkBVVZL3AW9J8jWasH07zZD+4p12VpKkRWvk4VpVFybZB1hDM+FoI3BcVU1MSDoY2N7T/gtJXgq8A3gX8DXghVX1Hz3dvocmoP8K2Au4su3zobk+n1nYRnOvb//l6cVgMZ87LO7zX8znDov7/Bf0uaeqRl2DJEkLynx4iIQkSQuK4SpJUscMV0mSOma4SpLUMcN1Hkjy5iRfSPJAku8PaHNwkkvbNt9J8t4kI5/tPReS/FSSS5LclWRrkiuTPGfUde0sSY5vX734YJK7kyy6W8iS7JZkY5JK8oxR1zPXkqxI8qEkt7X/37+e5Kz2+esLzkxfMzqODNf5YQnwMeCDk+1MsivNPb1LgGcDrwBOpLl9aSH6JM1tYs+lebHD9cAnkwx6NvSCkeQ3gL+huc/754BfAP5upEWNxnuY4ePmxtyhNH8f/wHwdOD/0Lwh7F2jLGouzPQ1o+PKW3HmkSQnAu+rqr36tv8qTeDsP3H/b5KTgXcD+7RPuloQkuwNfBf45ar6fLtt4oHZx1bV+qmOH2ftlYhNwJlV9aERlzMy7e/3tcBvAF8GjqiqSR9dupAleQNwSlX95Khr6VL7hrIvVtXEa0J3AW4H/ryqzhlpcR1y5DoeVgI39jxYA5rX7C2l+VfuQvJfwFeBlyfZvQ2cPwC+Q/OozIXsSOAAYHuSLyW5M8mnkhw26sJ2liTLgPOAl9G813kx25Od+/KQOdfzmtHeV4Jub9cHvWZ0LBmu42HQa/Ym9i0Y1VxK+RXgCJo3ZTwEnEbzhK27R1nbTjAxQnkbzRPIfg24G/hskieOqqidJc2T3D8M/EVVXTPickYqyVOAPwb+ctS1dGyq14wuqL/LDNc5kuScdjLGVMugd9YuONP99Wj/gn0/zUj1l4CjaZ4J/U9J9hvlOQxrBr8XJv48vrOqPlFV1wKvBAp48chOYJZmcP5/TPNqr7NHXHJnhvl7IMkBwGXAx6rqvNFUrtlakLNN54k/pflX+FRunWZfm2lCpteynn3jYLq/Hs+lGbE9oec1VH+Y5FiaiVzj+J3MdM994h8PN01srKptSW6lecb2uJrJ//uVwLb86HVkANck+duqesXclDenZvT3QJL9gc8AX6B9GckCM8xrRseS4TpHquq7NBNzunAV8OYk+1bVd9ptx9JM8rlp8GHzx3R/PZI8rv1xe9+u7YzplZYZnPu1NA8xfxrNyyZI8hhgBfCNOSxxTs3g/F8DvKVn0/40cwtOADbMTXVzayZ/D7Qj1s/QzC14Zftd5IJSVQ+3v8+PoX1LWX70mtFzR1lb1wzXeSDJwcATaUYnu/bc13dLVd0HXEETon+T5I003028A3h/+6L3heQqmu8Zz0+yBngQeBXwZJrbkRasqtqa5C+As5LcThOob2h3f2x0le0cVfXN3vUk97U/fr2qvjWCknaaNlg/S/P//PXAPhOj96paUCM6dvCa0YXCcJ0f1tBc8pzwpfa/zwE+W1WPJPk1mvtgrwLuB84H3rpTq9wJququJMcB7wT+FXgMze0Y/6uqrh9pcTvHG4Af0Nzr+liaEdtzF8FkrsXuWOAp7dL/D4k8uvn4msZrRhcE73OVJKljY/kdliRJ85nhKklSxwxXSZI6ZrhKktQxw1WSpI4ZrpIkdcxwlSSpY4arJEkdM1wlSeqYjz+UxkCSkTxKraoW1KP3pJ3Fxx9KktQxLwtLktQxw1UaI0l+Psmnk/xXkupblo66vn5JfjPJjUkebGten2T3UdclzTXDVRoTSX6O5p2fXwJ+CTgO+B7waeCEqto6uuoeLcl+wEeAvwZ+GlgFXMQCe4WaNBm/c5XGRJLPAXdU1Ut7tp0LHFVVK0dX2eSSHAlcC6yoqm/soO0pwNKqevdOKU6aY45cpTGQZBnwi8AH+nbdD1Rf21OSvGmKvt6b5MtJ3pzkb9tteyX5/SmOOWeSy9D9y6F9h11PM6q+McnHkrwqyRMGfMSRwHWDPl8aN45cpTGQ5AXApTSju3t7tq8Hbqqq1/RsOw/4+6r6lwF9fRdYVlXbe7atAD5eVc8ccMw+wJN2UOatVfVw33EBng08D/h1YDnwrKq6rd1/GPDnwErgduBNVXXRDj5HmvccuUrjYeLP6g8nAyU5HPhl4O/a9cOSfAZ4GfCBJC/q7yTJPwBPAK5L8odJrml3vRP4mSQbk7y1/7iq+m5V3byD5eFJjquq+reqOhM4AniYJmRJ8ti29lOAbwEvBN4+1K+ONM/4EAlpPGwAHgTem+SdwCHA+4H3V9W/9wTVbwGfpAmqv6eZQPRDVfXrSe6qqme0o9WT2l1vBp42aOQ6jCTPAo4BrgC+AzwL2Af4StvkOOAL7c+3AF8D9ujq86VRMlylMVBV303yW8CfAjcA3wTOBda2TeZjUG2lGVmfCiwFvgG8rqo+1e4/HLgROKz97+E05yaNPcNVGhNV9UmaUelk5l1QVdVXaEJ/kHuBpwL7AjcBbwE+uBNKk+ac4SotDLMNqnvZ+SPdC4BLgKcDW2gucV+6k2uQ5oThKi0MswqqqvqvJNcluRH4WFWtmaM6ez/zLuAXkvwn8LOTTYiSxpW34kgLSBtUh41LUCXZC/hMVR0x6lqkLjlylRaINqjuH5dgBaiq79PcoiMtKI5cJUnqmA+RkCSpY4arJEkdM1wlSeqY4SpJUscMV0mSOma4SpLUMcNVkqSOGa6SJHXMcJUkqWOGqyRJHTNcJUnqmOEqSVLH/j+2gsvW7sbqGAAAAABJRU5ErkJggg==\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 19,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Looks like our fit is consistently under-reporting errors on the standard deviation, sometimes the sample error is 10° larger than the error on the fit."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Model component sample error\n",
        "\n",
        "Can we get a measure of the variance on our model parameters? Can we trust them at all?"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Q: How often do our splits result in the same model?\n",
        "\n",
        "*in the below plot, if volunteers always agreed on the components to use **and** the isophote of the components, we would see a bimodal distribution with points at either zero or one (either all models for a given galaxy have a component, or all do not)*"
      ],
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "normalised_counts = pd.DataFrame([\n",
        "    comp_dfs[i].count() / len(comp_dfs[i]) for i in range(len(comp_dfs))\n",
        "])[['disk-axRatio', 'bulge-axRatio', 'bar-axRatio']]\n",
        "fig, ax = plt.subplots(figsize=(10, 4), ncols=3, sharey=True)\n",
        "for i, k in enumerate(('disk-axRatio', 'bulge-axRatio', 'bar-axRatio')):\n",
        "    ax[i].hist(\n",
        "        normalised_counts[k], density=False, alpha=0.7,\n",
        "        label=k.split('-')[0], color='C{}'.format(i),\n",
        "    )\n",
        "    ax[i].set_title(k.split('-')[0].capitalize())\n",
        "    ax[i].set_xlabel('Fraction of models with {}'.format(k.split('-')[0].capitalize()))\n",
        "ax[0].set_ylabel('Frequency')\n",
        "plt.savefig('duplicates_plots/model-consistency-comparison.pdf', bbox_inches='tight');"
      ],
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 720x288 with 3 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAmEAAAEWCAYAAAAuOkCvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHxNJREFUeJzt3XmcZHV57/HPlxmVICDLTBDZBhVE4hJxcI0rLkQU0CDiCmYioqIx6o2Y5EaSm9zoNWpMcAmKDhpBFBfALIrsSRQY1mERQQQBt4nKGmR97h/ndKhpu5nqnq461dWf9+tVrz77eepUPdVP/c6vzklVIUmSpOHaoOsAJEmSFiKLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEXYGEvyiST/u4/lrknyvGHEJI2KJJXkkV3HIWnhsgibx9ri6fYktyS5Mcl/JjkkyQYAVXVIVf2fruOUBqUnB25N8ssk/5xku67jkkaR+TJ6LMLmv5dU1SbADsD7gHcDR3UbkjRUL6mqjYGtgZ8C/9BxPNIoW+98SbJ4zqNaoCzCxkRV3VRVJwKvAA5M8pgkK5P8FUCSJUm+3raY/SLJWRMtZr2SPDrJD5K8ctjPQVofVfUr4HhgV4Akpyf5g4n5SQ5K8u9TrZtkyyQnJbk5yblJ/qp32SS7JDm5zZ0rkuw/6OcjDdIU+bJXkgvaHLguyeETyyZZ1p6+X5Hkh8Cp3UQ9fqxmx0xVnZPkeuAZk2a9E7geWNqOPwVY655VSXYDvga8uaq+PuhYpbmUZCOaLyHfmcXqHwVuAx4KLAO+AVzbbvfBwMnAnwO/CzwWODnJJVV12fpHLg3fFPlyG/A64FLgMTTv8Qur6ms9qz0LeDRw7zBjHWcWYePpR8AWk6bdRdP8vENVXQWcNWn+M4AVwGuq6vSBRyjNna8luRt4MLAGeOFMVk6yCPg94DFV9d/AZUmOBp7dLvJi4Jqq+kw7fkGSLwMvB/5iDuKXhmnKfJn0uX9xkmNpiq7eIuzwqrptWIEuBJ6OHE/bAL+YNO0DwFXAN5NcneSwSfMPAf7TAkzz0L5VtRmwIXAocEaSh85g/aU0X0iv65nWO7wD8OT2VP6NSW4EXk3TaibNN1PmS5InJzktyZokN9H8T1gyad3rJm9M68cibMwk2Z2mCFur70tV3VJV76yqhwN7A+9IskfPIocA2yf58PCileZOVd1TVV8B7gF+h+b0ykY9i0xXNK0B7ga27ZnW+4ux64AzqmqznsfGVfWmOQxfGqop8uUY4ERgu6p6CPAJIJNXG26U488ibEwk2TTJi4EvAP9UVasnzX9xkkcmCXATTeL1nte/BdgTeGaS9w0rbmmupLEPsDlwOXAh8LIkG7XXA1sx1XpVdQ/wFeDwdtldaPrGTPg6sHOS1yZ5QPvYPcmjB/uMpMGZIl82AX5RVb9K8iTgVZ0GuEDYJ2z+O6k9v38vcBnwIZpvMJPtBBxBc+rll8DHquq03gWq6sYkzwdOS3JXVa3zQq/SCDgpyT0039KvBQ6sqkvbVt3daX6GfzHweWC6ixIfCqwEfgJcARwLLIemFTnJC2hy60M0X14vAt4xqCckDdB0+fJm4INJjgDOAL4IbNZhnAtCqmxdlKReSd4PPLSqDuw6Fknjy9ORkha89jpgj2tP0TyJ5tTlV7uOS9J483SkJDX9YY4FHkZz+vKDwAmdRiRp7Hk6UpIkqQOejpQkSerAvDgduWTJklq2bFnXYUgAnHfeef9VVUvXveTgmBMaJeaEtLZ+c2JeFGHLli1j1apVXYchAZDk2q5jMCc0SswJaW395oSnIyVJkjpgESZJktQBizBJkqQOWIRJkiR1wCJMkiSpAxZhkiRJHbAIkyRJ6oBFmCRJUgcswiRJkjowL66YL83GipXnzmq9ow7afY4jkSTp19kSJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDAy3CkvxRkkuTXJLk2CQbJtkxydlJrkpyXJIHDjIGSZKkUTSwIizJNsDbgOVV9RhgEXAA8H7gw1X1SOCXwIpBxSBJkjSqBn06cjHwG0kWAxsBPwaeCxzfzj8a2HfAMUiSJI2cgRVhVXUD8LfAD2mKr5uA84Abq+rudrHrgW2mWj/JwUlWJVm1Zs2aQYUpzRvmhLQ2c0Lz3SBPR24O7APsCDwMeDCwZ7/rV9WRVbW8qpYvXbp0QFFK84c5Ia3NnNB8N8jTkc8DflBVa6rqLuArwNOBzdrTkwDbAjcMMAZJkqSRNMgi7IfAU5JslCTAHsBlwGnAfu0yBwInDDAGSZKkkTTIPmFn03TAPx9Y3e7rSODdwDuSXAVsCRw1qBgkSZJG1eJ1LzJ7VfVe4L2TJl8NPGmQ+5UkSRp1XjFfkiSpAxZhkiRJHbAIkyRJ6oBFmCRJUgcswiRJkjow0F9HSpJGyDGvmN16rzpubuOQBNgSJkmS1AmLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdsAiTJEnqgEWYJElSByzCJEmSOmARJkmS1AGLMEmSpA5YhEmSJHXAIkySJKkDFmGSJEkdGGgRlmSzJMcn+W6Sy5M8NckWSU5OcmX7d/NBxiBJkjSK+irCkjx2ltv/CPBvVbUL8HjgcuAw4JSq2gk4pR2XJElaUPptCftYknOSvDnJQ/pZoV3umcBRAFV1Z1XdCOwDHN0udjSw7wxjliRJmvf6KsKq6hnAq4HtgPOSHJPk+etYbUdgDfCZJBck+VSSBwNbVdWP22V+Amw11cpJDk6yKsmqNWvW9PVkpHFmTkhrMyc03/XdJ6yqrgT+DHg38Czg79u+Xi+bZpXFwG7Ax6vqCcBtTDr1WFUF1DT7O7KqllfV8qVLl/YbpjS2zAlpbeaE5rt++4Q9LsmHafp0PRd4SVU9uh3+8DSrXQ9cX1Vnt+PH0xRlP02ydbvdrYGfrUf8kiRJ81K/LWH/AJwPPL6q3lJV5wNU1Y9oWsd+TVX9BLguyaPaSXsAlwEnAge20w4ETphl7JIkSfPW4j6X2wu4varuAUiyAbBhVf13VX3uftZ7K/D5JA8ErgZeT1P4fTHJCuBaYP9ZRy9JkjRP9VuEfQt4HnBrO74R8E3gafe3UlVdCCyfYtYe/QYoSZI0jvo9HblhVU0UYLTDGw0mJEmSpPHXbxF2W5LdJkaSPBG4fTAhSZIkjb9+T0e+HfhSkh8BAR4KvGJgUUmSJI25voqwqjo3yS7AxC8dr6iquwYXliRJ0njrtyUMYHdgWbvObkmoqs8OJCpJkqQx11cRluRzwCOAC4F72skFWIRJkiTNQr8tYcuBXdvbDEmSJGk99fvryEtoOuNLkiRpDvTbErYEuCzJOcAdExOrau+BRCVJkjTm+i3CDh9kEJIkSQtNv5eoOCPJDsBOVfWtJBsBiwYbmiRJ0vjqq09YkjcAxwP/2E7aBvjaoIKSJEkad/12zH8L8HTgZoCquhL4zUEFJUmSNO76LcLuqKo7J0aSLKa5TpgkSZJmod8i7IwkfwL8RpLnA18CThpcWJIkSeOt3yLsMGANsBp4I/AvwJ8NKihJkqRx1++vI+8FPtk+JEmStJ76vXfkD5iiD1hVPXzOI5IkSVoAZnLvyAkbAi8Htpj7cCRJkhaGvvqEVdXPex43VNXfAXsNODZJkqSx1e/pyN16RjegaRnrtxVNkiRJk/RbSH2wZ/hu4Bpg/zmPRpIkaYHo99eRzxl0IJIkSQtJv6cj33F/86vqQ3MTjiRJ0sIwk19H7g6c2I6/BDgHuHIQQUmSJI27fouwbYHdquoWgCSHA/9cVa8ZVGCSJEnjrN/bFm0F3Nkzfmc7TZIkSbPQb0vYZ4Fzkny1Hd8XOHowIUmSJI2/fn8d+ddJ/hV4Rjvp9VV1weDCkiRJGm/9no4E2Ai4uao+AlyfZMcBxSRJkjT2+irCkrwXeDfwnnbSA4B/GlRQkiRJ467flrCXAnsDtwFU1Y+ATQYVlCRJ0rjrtwi7s6oKKIAkDx5cSJIkSeOv3yLsi0n+EdgsyRuAbwGfHFxYkiRJ463fX0f+bZLnAzcDjwL+vKpO7mfdJIuAVcANVfXitkP/F4AtgfOA11bVnfe3DUmSpHGzziKsLaK+1d7Eu6/Ca5I/BC4HNm3H3w98uKq+kOQTwArg47PYriRJ63ToKYfOet0j9jhiDiOR1rbO05FVdQ9wb5KHzHTjSbYF9gI+1Y4HeC5wfLvI0TQXfpUkSVpQ+r1i/q3A6iQn0/5CEqCq3raO9f4O+GPu+yXllsCNVXV3O349sE3/4UqSJI2Hfouwr7SPviV5MfCzqjovybNnGliSg4GDAbbffvuZri6NHXNCWps5ofnufouwJNtX1Q+rajb3iXw6sHeSFwEb0vQJ+wjNLywXt61h2wI3TLVyVR0JHAmwfPnymsX+pbFiTkhrMyc0362rT9jXJgaSfHkmG66q91TVtlW1DDgAOLWqXg2cBuzXLnYgcMJMtitJkjQO1lWEpWf44XO0z3cD70hyFU0fsaPmaLuSJEnzxrr6hNU0wzNSVacDp7fDVwNPmu22JEmSxsG6irDHJ7mZpkXsN9ph2vGqqk2nX1WSJEnTud8irKoWDSsQSZKkhaTfe0dKkiRpDlmESZIkdcAiTJIkqQMWYZIkSR2wCJMkSeqARZgkSVIHLMIkSZI6YBEmSZLUAYswSZKkDliESZIkdcAiTJIkqQMWYZIkSR2wCJMkSeqARZgkSVIHLMIkSZI6YBEmSZLUAYswSZKkDliESZIkdcAiTJIkqQMWYZIkSR2wCJMkSeqARZgkSVIHLMIkSZI6YBEmSZLUAYswSZKkDliESZIkdcAiTJIkqQMWYZIkSR2wCJMkSeqARZgkSVIHLMIkSZI6YBEmSZLUAYswSZKkDgysCEuyXZLTklyW5NIkf9hO3yLJyUmubP9uPqgYJEmSRtUgW8LuBt5ZVbsCTwHekmRX4DDglKraCTilHZckSVpQBlaEVdWPq+r8dvgW4HJgG2Af4Oh2saOBfQcVgyRJ0qgaSp+wJMuAJwBnA1tV1Y/bWT8BtppmnYOTrEqyas2aNcMIUxpp5oS0NnNC893Ai7AkGwNfBt5eVTf3zquqAmqq9arqyKpaXlXLly5dOugwpZFnTkhrMyc03w20CEvyAJoC7PNV9ZV28k+TbN3O3xr42SBjkCRJGkWD/HVkgKOAy6vqQz2zTgQObIcPBE4YVAySJEmjavEAt/104LXA6iQXttP+BHgf8MUkK4Brgf0HGIMkSdJIGlgRVlX/DmSa2XsMar+SJEnzgVfMlyRJ6oBFmCRJUgcswiRJkjpgESZJktQBizBJkqQOWIRJkiR1wCJMkiSpAxZhkiRJHbAIkyRJ6sAgb1skSdKCdegph8563SP2OGIOI9GosiVMkiSpAxZhkiRJHbAIkyRJ6oBFmCRJUgfsmC9J88kxr+g6AmlK/hBh5mwJkyRJ6oAtYZIkacEZhZY7W8IkSZI6YBEmSZLUAYswSZKkDliESZIkdcAiTJIkqQMWYZIkSR2wCJMkSeqARZgkSVIHvFirJEnTWJ8LenZlFC5Cqv7YEiZJktQBW8IkSRIwP1v+5jNbwiRJkjpgESZJktQBizBJkqQO2CdMGjfHvGL2677quLmLQ5L6tFB/0WlLmCRJUgdsCZN0n/VpRVsfXbXAddVq2NVxlsbQfP5Fpy1hkiRJHeikCEuyZ5IrklyV5LAuYpAkSerS0E9HJlkEfBR4PnA9cG6SE6vqstlsb8XKc2cdy1EH7T7rdRearo7z+uxXC4Sn9jSG5vMpNvWvi5awJwFXVdXVVXUn8AVgnw7ikCRJ6kyqarg7TPYD9qyqP2jHXws8uaoOnbTcwcDB7eijgCum2eQS4L8GFO76GMW4jKk/64pph6paOqxgJvTkxBJgS6bPia6N4ms6wdhmZ9RzAubf/4lRjAlGM675GFNfOTGyRdgMtreqqpbPZYxzYRTjMqb+jGJMvYxv9oxtdkY5tn6MYvyjGBOMZlzjHFMXpyNvALbrGd+2nSZJkrRgdFGEnQvslGTHJA8EDgBO7CAOSZKkzgz915FVdXeSQ4FvAIuAT1fVpeuxySPnJrI5N4pxGVN/RjGmXsY3e8Y2O6McWz9GMf5RjAlGM66xjWnofcIkSZLkFfMlSZI6YREmSZLUgXlThPVzq6Mk+ye5LMmlSY7pOqYkH05yYfv4XpIbBx1Tn3Ftn+S0JBckuTjJi0Ygph2SnNLGc3qSbYcQ06eT/CzJJdPMT5K/b2O+OMluQ4hpXcfpQUmOa+efnWRZz7z3tNOvSPLCDmJ7R5t/F7ev5Q498+7pyYU5/yFOH7EdlGRNTwx/0DPvwCRXto8D5zq2PuOb9rNikMdufXJgGMetH6OYM6OYK6OaI6OWG0PPiaoa+QdNB/7vAw8HHghcBOw6aZmdgAuAzdvx3+w6pknLv5XmRwijcKyOBN7UDu8KXDMCMX0JOLAdfi7wuSEcq2cCuwGXTDP/RcC/AgGeApw9AsfpzcAn2uEDgON6XseLgAcBO7bbWTTk2J4DbNQOv2kitnb81o6P20HAEVOsuwVwdft383Z482HHN2n5tT4rBnzsZpUDwzhuc/jaDzVnRjFXRjVHRjE3hp0T86UlrJ9bHb0B+GhV/RKgqn42AjH1eiVw7IBj6jeuAjZthx8C/GgEYtoVOLUdPm2K+XOuqs4EfnE/i+wDfLYa3wE2S7L1AEPq5zjtAxzdDh8P7JEk7fQvVNUdVfUD4Kp2e0OLrapOq6r/bke/Q3MNwGFYn1uhvRA4uap+0X52nAzs2XF8w/qsWJ8cGMZx68co5swo5sqo5sjI5cawc2K+FGHbANf1jF/fTuu1M7Bzkv9I8p0kg/5A6CcmoDnVRvNN69Sp5ncQ1+HAa5JcD/wLzbeLrmO6CHhZO/xSYJMkWw44rnXp+zUe4v7+Z5mquhu4ieYWRoOOdabbX0HzbXHChklWtbm57xzGNZPYfq89fXB8kokLRg/jNV7fz4pBHrt1mS72YefGdEYxZ0YxV0Y1R+ZjbsxpTgz9OmEDtJjmlOSzab5VnJnksVU1lH5Y63AAcHxV3dN1IK1XAiur6oNJngp8LsljqureDmN6F3BEkoOAM2nuojAqx0szkOQ1wHLgWT2Td6iqG5I8HDg1yeqq+v4QwzoJOLaq7kjyRpqWkecOcf/9muqzoutjpwEZsVwZ9RwZy9yYLy1h/dzq6HrgxKq6q21a/h5NUdZlTBMOYEinF+gvrhXAFwGq6tvAhjQ3I+0spqr6UVW9rKqeAPxpO63rAnrYt9jqZ3//s0ySxTSnk3/e57qDjo0kz6N5/fauqjsmplfVDe3fq4HTgScMM7aq+nlPPJ8CntjvusOIr8evfVYM+Nity3Sxj8rt50YxZ0YxV0Y1R+ZjbsxtTvTTUa3rB00r19U0TZETnfd+a9IyewJHt8NLaJoFt+wypna5XYBraC+MOyLH6l+Bg9rhR9P0CRtYfH3GtATYoB3+a+Avh3S8ljF9B8y9WLsD5jkj8Nq9hbU7GX+xHf4t1u5kfDVz2zG/n9ieQNPJdqdJ0zcHHtTzOl/J/XS+HVBsW/cMvxT4Tju8BfCDNsbN2+Ethv26tsv92mfFoI/dbHNgGMdtDl/7oebMKObKqObIqObGMHNiqAmzngflRTStW98H/rSd9pc03yJoD8iHgMuA1cABXcfUjh8OvG/EjtWuwH+0b/gLgReMQEz7tUn0PZpvYQ8aQkzHAj8G7qJpSV0BHAIc0vOe+mgb82pg+Qgcpw1pfkl6FXAO8PCedf+0Xe8K4Hc7iO1bwE/b99SFNC3TAE9rj99F7d8VHcT2N8ClbQynAbv0rPv77fG8Cnh9F69rO/5rnxWDPnbrkwPDOG5z9NoPPWdGMVdGNUdGLTeGnRPetkiSJKkD86VPmCRJ0lixCJMkSeqARZgkSVIHLMIkSZI6YBEmSZLUgbEswibdWf3CJMvmYJv7Jtm1Z/wv24vtDUySY9tbSPzRAPdxUJIj1neZGe7zP9u/y5K8aqb7SXJ6kivaY/PdJEck2Wzy9u9n/VvXJ/75wByY0T7mcw5cmOTyJAf3sc7hSd61fpGPH3NlRvtYELkyTON026Jet1fVb083M8niau4hNhP7Al+nuQ4ZVfXn6xHfOiV5KLB7VT1ykPvpQlU9rR1cBrwKOGYWm3l1Va1K8kCa69ucQHvrj57tL2TmwAib4xzYAvh+kpXV3ARZM2OujLCucyXJohrgLQfHsiVsKm3VfGKSU4FTkmyc5JQk5ydZnWSfnmVf136juCjJ55I8Ddgb+EBbTT8iycok+7XL75HkgnY7n07yoHb6NUn+omcfu0wR14ZJPtPOvyDJc9pZ3wS2aff3jEnrrEzy8TQ3Lb06ybPb/V6eZGXPcq9st3tJkvf3TH99ku8lOQd4es/0pUm+nOTc9vF0Jkny8nZ7FyU5c4r5H02ydzv81SSfbod/P8lft8MTLVHvA57RPseJb28PS/JvSa5M8v+mfjXv0ybSHwPbJ3l87/aTbJ3kzHb7l0xxHJck+XaSvda1n3FgDoxnDgAbA7fR3mu1Z9sk2a/3ePRM3719fS9M8oEkl7TTF7Xj57bz39jH/seOubJgcuXjaW4AfmmSv+iJ65ok709yPvDyPrY7e3N55dtRebQHeOJKxF9tpx1Ec/XbLdrxxcCm7fASmivchuZWFt8DlrTzJpZfCezXs4+VNFd535DmFkk7t9M/C7y9Hb4GeGs7/GbgU1PE+k7g0+3wLsAP220uY/rbJqwEvtDGuw9wM/BYmqL6POC3gYe121raPtdTab6dbd0z/YE0V84/ot3uMcDvtMPbA5f3HLuJZVYD27TDm00R2wHAB9rhc7jv1hefAV7YDt/a/n028PWedQ+iuYXFQ9pjcC2w3RT7OJ1JV68Hvga8YtL238l9V2BeBGwyMR/YCjgbeH7X71dzwByYZQ5cAVwM3A68sWferT3D+wEr2+HDgXe1w5cAT22H3zdxnIGDgT9rhx8ErAJ27Pr9bK6YKwwmVyZer0Xtco/reS3+eBjvv3FtCbu9qn67fby0Z/rJVfWLdjjA/01yMc0tJLah+cf8XOBLVfVfAD3LT+dRwA+q6nvt+NHAM3vmf6X9ex5Nokz2O8A/tfv6Ls0baed1P0VOqubdshr4aVWtrqp7aW47sQzYHTi9qtZU05T++TauJ/dMvxM4rmebzwOOSHIhcCKwaZKNJ+33P4CVSd5A88ad7Cyabyu70jTF/zTJ1sBTgfvtq9U6papuqqpftevv0Mc60Lyek50LvD7J4cBjq+qWdvoDgFNokuzkPrc/35gD458Dr66qx9H8A3xXkr5yJU3/yU2q6tvtpN7TOy8AXtc+/7OBLYGd+tnuPGauLNxc2b9t7bqApqDetWed4xiCce0TNp3beoZfTVPdP7Gq7kpyDU01Pdcm7kp/D3N7vCe2e2/P8MT4Ypr7Xs3UBsBT2jf0/0juq2+q6pAkT6a5iel5SZ5YVT/vmX9D+yG/J3AmzU1N96f5NnML69b7XPo6ZkkW0Xyzu7x3elWdmeSZbawrk3yoqj4L3E3zIfdC4Iw+Yhon5sD9m3c5UFVr2n8kT6b5p9x7L7qZvp6haY35xgzXG0fmyv2b17mSZAPgXTR96X7ZnprtfU1vm2obc21cW8L68RDgZ21CPYf7KuhTgZcn2RIgTUc+gFuATabYzhXAsiQTHSJfy8z+sZ9Fk+Ak2ZmmUr9iJk9kGucAz0rT72kR8Mo2rrPb6VsmeQBrn+/+JvDWiZEkv9ZZNckjqursajqargG2m2Lf3wHeTpNUZ9G80c+aYrnpjmnf2ufwN8B1VXXxpHk70Hzr+yTNTcF3a2cVzY1Wd0ny7vXZ/zxnDoxHDmwEPIHmhsLQtCY8uv0n89LJy1fVjcAt7T9HaE4JTfgG8Kb2uJBk5yQPXp/4xoS5Mn65silNoXVTkq2A312fbc/WQmsJ6/V54KQkq2n6PXwXoKoubTsEnpHkHppmyoNozql/MsnbaM7t0y7/qySvB76UZDHNKbBPzCCOjwEfb+O4Gzioqu7o/TYxG1X14ySHAafRfLv956o6AZqfqgPfBm6k6Qcx4W3AR9sm98U0SXHIpE1/IMlO7TZPobmD/WRnAS+oqquSXEvz7WaqpLoYuCfJRTT9Fn45g6f4+SR30PRb+RZNX4fJng38ryR30fQDe93EjKq6J8krgROT3FJVH5vBvseFOTD/c+B2mhxYWVXntdMPo/ll3hqa13XyKSKAFTSv5b00/2xvaqd/iub01PlpXoA1NH2DFjpzZQxzJckFNK/ldTSnTocuzWliSdJCkWTjqpr4FfFhwNZV9YcdhyUtOAu5JUySFqq9kryH5n/AtTStN5KGzJYwSZKkDizkjvmSJEmdsQiTJEnqgEWYJElSByzCJEmSOmARJkmS1IH/D8eM1QFJFZxIAAAAAElFTkSuQmCC\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 20,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "component_parameter_variances = pd.DataFrame([\n",
        "    df['std']\n",
        "    for df in \n",
        "    (\n",
        "        comp_df.describe().transpose()\n",
        "        for comp_df in comp_dfs\n",
        "    )\n",
        "], index=comp_ids)"
      ],
      "outputs": [],
      "execution_count": 21,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.figure(figsize=(10, 13))\n",
        "component_parameter_variances.drop([\n",
        "    'disk-i0', 'disk-n', 'disk-c',\n",
        "    'bulge-i0', 'bulge-n', 'bulge-c',\n",
        "    'bar-i0', 'bar-n', 'bar-c',\n",
        "], axis=1).hist(ax=plt.gca());"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/Users/tlingard/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py:3267: UserWarning: To output multiple subplots, the figure containing the passed axes is being cleared\n",
            "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 720x936 with 16 Axes>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAlMAAALyCAYAAAD+AVGjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xu8XFV99/HP1wTkEuSqpyEBDwqiKCXYiFjUHrlUBAVsfVCKSBCb9ikoSPpAvLRQBRtaERF8qFGQUMNNLkIRrRhzXhQfpRKIBIjIpUESQiKYAEFEA7/nj71OmAznsufMzN57Zr7v12temX2d356zMvObtdZeSxGBmZmZmY3Py8oOwMzMzKyTOZkyMzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7MmOJkyMzMza0LXJFOSlkk6sOw42k3SOyTdV3YclumVcmfFcHmydqhauZK0u6TFkp6W9AlJm0v6D0lPSvp22fGNR9ckU1UlaVDS7yStk/S4pGslTW7g+JC069ByRPxXROzenmjNzMyaI2mGpOfT917tY8e0y6nAwojYKiK+AnwA6AO2j4j/VVrgTXAyVUOZdrwnJ0bEJGBXYBLwxTa8hnWoNpY760EuT9YOecuVpInp6U8iYlLd49G07dXAPTWHvRr4ZUSsb3HYhem2/3BvkXSvpDWSvilpM0nbSrpR0q/T+hslTR06INUcnSXpx8BvgdfUn1TSayX9SNITqXZpvqRtarb9RtKb0/KO6bUG6s8TEWuB7wDTas69j6SfSForaaWkCyRtmrbdknb7ecrqPyhpQNLymuPfkK5hraR7JB3WijfSGtKucjdD0o8lnZv+vg9J+tO0/hFJqyUdW3fOj9Udf2t6/qep7O6UlvdKcb2+nW+MjUsnlKevSjqn7vw3SPpkO94Qa4kiytUTwBmjBSHpR8C7gAvS99rlwD8CH0zLx7fyoovSbcnU0cC7gdcCrwM+S3aN3yTLfHcGngUuqDvuGGAmsBXw8DDnFfDPwI7AG4CdSAUmIh4ETgO+JWmL9FrzImLwJSeRtgf+AnigZvXzwCeBHYC3AQcAf5fO/c60z14pq7+y7nybAP8B/AB4FfBxYL4kNwMWq13lDuCtwF3A9sBlwBXAW8hqOT9M9oE0aawAI+L/AV8D5knaHPgW8A8R8YvcV2lFqXx5AuYBRynVVEjaATgwndOqqd3l6iGyprqzRgsiIvYH/ovUYhMRRwFfAK5Myxc1fmkVEBFd8QCWAX9bs3wI8OAw+00D1tQsDwKfa/C1jgDurFt3A7CE7IPq5XXn/y3wJBDAYmDnUc59MnBdzXIAu9YsDwDL0/N3AI8BL6vZfjlwRtl/j155tLPcATOA+2uW90zloa9m3RPAtJpzfqzu+FtrljcBFqVy+n1AZb9/fnR0eVoKHJSenwjcVPb750dp5epXw6xbD6yteTxYd97asnUG8K2y36dmHt1WM/VIzfOHgR0lbSHpa5IelvQUcAuwjaQJwx2n7G65oc5y96R1fZKukLQineNbZDVJtb4OvAk4PyKeq9v2iYjYGvhjYFugthr1dalq9bF07i8Mc+6R7Ag8EhEv1F33lJzHW2u0pdwlq2qePwsQEfXr8tQkEBF/AC4hK6fnRPoUs8rpiPJEVjv14fT8w8C/5zzOytHOclV77iE/jYhtah6vbenVVEy3JVM71TzfGXgUmAXsDrw1Il4BDDWdqWbfDV8qkd0tN9RZ7o1p9RfSPnumc3y49vhULf5l4CLgDEnbDRdcRCwBzgS+Kmno+AuBXwC7pXN/ui620TwK7KSNOwXuDKzIeby1RrvKXaOeAbaoWf6j2o2SpgCnk1XrnyPp5eN8HWuvjihPZD8qD5e0F1n3h++M83WsGO0sVz3/w6zbkqkTJE1NycxngCvJ2nmfBdam9aeP47xbAeuAJ9MX0v+p234ecHtEfAz4LvBvo5xrHlm78lBH8a2Ap4B1qTPw/67bfxXDdPpLbiNrQjxV0ibKOr2/j6wfhBWnXeWuUYuBv0i/NncFNnTkTMn7JWQJ//HASuDzBcRkjat8eQKIiOXAz8hqpK6JiGcLiMnGryrlqit1WzJ1GVln7IeAB8lqgb4MbA48DvyUrK9Io/4JeDNZv6fvAtcObZB0OHAwLyZBpwBvlnT0cCeKiN+TJV//kFb9PfBXwNNkTYVX1h1yBlmn4bWSjhzmXO8D3pOu7/8CHwl3Ki5au8pdo84Ffk+WgM8D5tds+wTZTQr/kJr3jgOOk/SOAuKyxnRCeRoyj6zvlZv4qq/ocvU2vXScqbe08PyVInebMDOz8ZD0TrLmvle7D571sm6rmTIzswKkoVlOAr7hRMp6nZMpMzNriKQ3kN3uPpmsqcisp+UZGn4nSQvTyKn3SDoprT8jDRWwOD0OaX+4ZmZWtohYGhFbRsSfRsRTZcdjVrYx+0wpm5R3ckTcIWkrskH/jgCOBNZFhOeZMzMzs541cawdImIl2W3URMTTkpbiQSHNzMzMgAbv5pPUTzZC6pvIhgCYQTZG0u3ArIhYM9rxO+ywQ/T392+07plnnmHLLbdsIOTO1E3XuWjRoscj4pVlx9HL5aleJ1+3y1PxuvG6hq7J5ak5VY+x6Pjylqcxa6aGpFG+rwFOjoinJF1INuhfpH/PAT46zHEzySZJpK+vjy9+ceNWwXXr1jFpUt7ZCzpXN13nu971rpEmuyxUf38/t99++0brBgcHGRgYKCegEnXydUtyeSpYN17X0DW5PDWn6jEWHV/e8pQrmUq3wF4DzI+Ia2Hj+ZwkfR24cbhjI2IuMBdg+vTpUf8mVP0P1yq9cp1mZma9Js/dfCKbgmJpRHypZv3kmt3eD9zd+vDMzMzMqi1PzdR+wDHAEkmL07pPA0dJmkbWzLcM+Ju2RGhmZmZWYXnu5ruVjWeQHnJT68MxMzMz6yy5O6CbmZl1O0nLyCaefx5YHxHTJW1HNgl9P1lLzJFj3b1uvaVjkqn+2d/Nve+yOYe2MRKrsiUrnmRGzrLicmJjcXnqWe+KiMdrlmcDCyJijqTZafm0Rk/q8tS9PDefmZnZ6A4H5qXn88hmATHbwMmUmZnZiwL4gaRFaZxEgL40GwjAY0BfOaFZVXVMM5+ZmVkB3h4RKyS9CrhZ0i9qN0ZESBp26pD6QaoHBwc32t63Oczac32uIOqPLcq6detKe+08qhqfkykzM7MkIlakf1dLug7YB1glaXJErExjLK4e4dhRB6k+f/71nLMk39fusqMHxtynHao+wHRV43Mzn5mZGSBpS0lbDT0H/pxsQOobgGPTbscC15cToVWVa6bMzMwyfcB12cQfTAQui4jvS/oZcJWk44GHgSNLjNEqyMmUmZkZEBEPAXsNs/4J4IDiI7JO4WY+MzMzsyY4mTIzMzNrgpMpMzMzsyY4mbJSSJog6U5JN6blXSTdJukBSVdK2rTsGM3MzPJwMmVlOQlYWrN8NnBuROwKrAGOLyUqMzOzBjmZssJJmgocCnwjLQvYH7g67eK5r8zMrGM4mbIyfBk4FXghLW8PrI2IoXkWlgNTygjMOpObjc2sTB5nygol6b3A6ohYJGlgHMd3/NxX7VDV+aoKNNRs/Iq0PNRsfIWkfyNrNr6wrODMrLs5mbKi7QccJukQYDOyL7/zgG0kTUy1U1OBFcMd3A1zX7VDVeerKkJNs/FZwCk1zcZ/lXaZB5yBkykzaxM381mhIuJTETE1IvqBDwE/ioijgYXAB9JunvvKGuFmYzMrlWumrCpOA66QdCZwJ3BRyfFYB3Cz8fh0Y7NwN16TdY4xkylJOwGXkk0AGcDciDhP0nbAlUA/sAw4MiLWtC9U6zYRMQgMpucPAfuUGY91JDcbj0M3Ngt34zVZ58jTzLcemBURewD7AidI2gOYDSyIiN2ABWnZzKwwbjY2syoYM5mKiJURcUd6/jTZHTNTgMPJOnaCxwUys2o5jawz+gNkfajcbGxmbdNQnylJ/cDewG1AX0SsTJseI2sGNDMrhZuNzawsuZMpSZOAa4CTI+Kp7O7jTESEpBjhuFE7eObtNJi3EyhUsyOoO0eamZl1p1zJlKRNyBKp+RFxbVq9StLkiFgpaTKwerhjx+rgmbfT4IzZ380TKlDNjqDuHGlmZtadxuwzlQbAuwhYGhFfqtl0A1nHTnAHTzMzM+tRee7m2w84Bthf0uL0OASYAxwk6X7gwLRsZmbW0TzXozVqzGa+iLgV0AibD2htOGZmZqXzXI/WEE8nY2ZmltTM9fiNtDw01+PVaRcPBWQv4WTKzMzsRZ7r0RrmufnMzMzwXI9Q/WF8qhqfkykzM7NMz8/1WPVhfKoan5v5zMzM8FyPNn5OpszMzEbnuR5tVKU38y1Z8WRDo5ubmZm1m+d6tEa4ZsrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mrFCSNpP035J+LukeSf+U1u8i6TZJD0i6UtKmZcdqZmaWh5MpK9pzwP4RsRcwDThY0r7A2cC5EbErsAY4vsQYzczMcnMyZYWKzLq0uEl6BLA/cHVaPw84ooTwrMO4ptPMqsDJlBVO0gRJi4HVwM3Ag8DaNCM7wHJgSlnxWUdxTaeZla70ufms90TE88A0SdsA1wGvz3uspJnATIC+vj4GBwc32t63Oczac/0wR75U/bGdbN26dV11PXlFRAAj1XT+VVo/DzgDuLDo+MysNziZstJExFpJC4G3AdtImphqp6YCK0Y4Zi4wF2D69OkxMDCw0fbz51/POUvyFetlRw+MuU+nGBwcpP696BWSJgCLgF2Br+KaTjMr2JjfOpIuBt4LrI6IN6V1ZwB/Dfw67fbpiLipXUFa95D0SuAPKZHaHDiIrElmIfAB4ArgWOD68qK0TuKazsZ1Y01mN16TdY48P+EvAS4ALq1bf25EfLHlEVm3mwzMS7UJLwOuiogbJd0LXCHpTOBO4KIyg7TO45rO/LqxJrMbr8k6x5ifEhFxi6T+9odivSAi7gL2Hmb9Q8A+xUdkncw1nWZWBc30mTpR0keA24FZEbFmuJ1aWY2eVxWrel0FbdYWruk0s9KNN5m6EPg82V0znwfOAT463I6trEbPq4rV7a6CNms913SaWRWMa5ypiFgVEc9HxAvA1/GHlpmZdTgPAmvjNa5kStLkmsX3A3e3JhwzM7PSeBBYG5cxkylJlwM/AXaXtFzS8cC/SFoi6S7gXcAn2xynmZlZW3m6KxuvPHfzHTXManfmNDOzruNBYG08PAK6mZlZ0uuDwFb9zvOqxudkyszMrE6vDgJb9TvPqxrfuDqgm5mZdRtJr0w1UtQMAruUFweBBQ8Ca8NwzZSZmVnGg8DauDiZMjMzw4PA2vi5mc/MzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMoKJWknSQsl3SvpHkknpfXbSbpZ0v3p323LjtXMzCwPJ1NWtPXArIjYA9gXOEHSHsBsYEFE7AYsSMtmZmaV52TKChURKyPijvT8aWApMAU4HJiXdpsHHFFOhNZJXNNpZlXgZMpKI6mfbIb224C+iFiZNj0G9JUUlnUW13SaWekmjrWDpIuB9wKrI+JNad12wJVAP7AMODIi1rQvTOs2kiYB1wAnR8RTkjZsi4iQFCMcNxOYCdDX18fg4OBG2/s2h1l7rs8VQ/2xnWzdunVddT15pQR8ZXr+tKTams6BtNs8YBA4rYQQzawHjJlMAZcAFwCX1qwb+tU3R9LstOwPKstF0iZkidT8iLg2rV4laXJErJQ0GVg93LERMReYCzB9+vQYGBjYaPv586/nnCV5ijUsO3pgzH06xeDgIPXvRa9xTaeZlWXMb52IuCV9SNXyrz4bF2VVUBcBSyPiSzWbbgCOBeakf68vITzrUK7pbEw31mR24zVZ58j3E/6l/KvPxms/4BhgiaTFad2nyZKoqyQdDzwMHFlSfNZhXNPZuG6syezGa7LOMd5kaoPRfvVBa3/55VXFXyf+1ZSJiFsBjbD5gCJjsc7nmk5rJUk7kXVp6QMCmBsR57mfsI1lvMlUrl990NpffnlV8ReifzWZtYVrOq2Vhu4OvUPSVsAiSTcDM3A/YRvFeLMY/+ozs9K5ptNayXeH2njlGRrhcrJCtIOk5cDp+FefmZl1sfHcHdoNNzRUvUtKVePLczffUSNs8q8+MzPrOuO9O7QbbmioepeUqsbnEdDNzMyS0e4OTdtH7SdsvcnJlJmZGbnuDgX3E7ZhtPY2OjMzs87lu0NtXJxMmZmZ4btDbfzczGdmZmbWBCdTZmZmZk1wMmVmZmbWBCdTZmZmZk1wMmVmZmbWBCdTZmZmZk1wMmVmZmbWBCdTZmZmZk3oykE7+2d/N9d+y+Yc2uZIzMzMrNu5ZsrMzMysCU6mzMzMzJrgZMrMzMysCU6mzMzMzJrgZMrMzMysCU6mrFCSLpa0WtLdNeu2k3SzpPvTv9uWGaOZmVkjnExZ0S4BDq5bNxtYEBG7AQvSslkuTtDNrGxNJVOSlklaImmxpNtbFZR1r4i4BfhN3erDgXnp+TzgiEKDsk53CU7QzaxEraiZeldETIuI6S04l/WmvohYmZ4/BvSVGYx1Fifo1kqu6bTx6MoR0K1zRURIipG2S5oJzATo6+tjcHBwo+19m8OsPdfneq36YzvZunXruup6WiBXgt6r5akby0sLr+kS4ALg0pp1QzWdcyTNTsunteLFrDs0m0wF8IP05fe1iJjbgpis96ySNDkiVkqaDKweacdUxuYCTJ8+PQYGBjbafv786zlnSb5ivezogTH36RSDg4PUvxeWGS1B79Xy1I3lpVXXFBG3SOqvW304MHTyecAgTqasRrPJ1NsjYoWkVwE3S/pFqnLfoJW//Frt/PnX5953zylbN/Va3fhLsIVuAI4F5qR/8/9hzIaXO0E3y8FdEWxUTSVTEbEi/bta0nXAPsAtdfu07JdfmZr91dmNvwTHQ9LlZL/wdpC0HDidLIm6StLxwMPAkeVFaF3CCbq1xWg1nd3QbFz1H/5VjW/cWYykLYGXRcTT6fmfA59rWWTWlSLiqBE2HVBoINY1nKBbAXLVdHZDs3HVf/hXNb5mqoT6gOskDZ3nsoj4fkuiMjPLyQm6FcA1nTaqcSdTEfEQsFcLYzEzMyuVazptPKrfWcnMzKwgrum08fB0MmZmZmZNcM2U9az+2d/Ntd+yOYe2ORLrBi5PZr3LNVNmZmZmTXAyZWZmZtYEJ1NmZmZmTXAyZWZmZtYEJ1NmZmZmTXAyZWZmZtYEJ1NmZmZmTXAyZWZmZtYEJ1NmZmZmTXAyZWZmZtYEJ1NmZmZmTfDcfGZj8JxrZmY2GtdMmZmZmTXByZSZmZlZE9zMZ2ZWIDcbm3UfJ1M5lfkB2OrXznu+Rs5pZmatU9bn/qw91zOQa8/OUNT3XVPNfJIOlnSfpAckzW7mXGYuT9ZKLk/WSi5PNppx10xJmgB8FTgIWA78TNINEXFvq4Kz3tEN5alXa/yqeN0uT52ritfdDeXJ2quZmql9gAci4qGI+D1wBXB4a8KyHuTyZK3k8mSt5PJko2ommZoCPFKzvDytMxsPlydrJZcnayWXJxtV2zugS5oJzEyL6yTdV7fLDsDj7Y6jKDp7xE1tv85RXrvV53x1618pn24pT234W3Xydbs8NWkc5akjrmssddc9dE09VZ5a/VnyCdjhEx+udNloW9lt5vOpmWRqBbBTzfLUtG4jETEXmDvSSSTdHhHTm4ijI/TKdTbB5akJvXrdo3B5GkU3Xlebr6lnylPVY6xqfM008/0M2E3SLpI2BT4E3NCasKwHuTxZK7k8WSu5PNmoxl0zFRHrJZ0I/CcwAbg4Iu5pWWTWU1yerJVcnqyVXJ5sLE31mYqIm4CbmoxhxCrRLtMr1zluLk9N6dXrHpHL06i68braek09VJ6qHmMl41NElB2DmZmZWcfyRMdmZmZmTSgsmRprKH5JL5d0Zdp+m6T+omJrpRzXeYqkeyXdJWmBpNJu4+02vTLdg6SLJa2WdHfNuu0k3Szp/vTvtmXG2A26oTx1Y1mRtJOkhelz9B5JJ6X1lb6uKpYnScskLZG0WNLtaV2p72MjZVaZr6T39C5Jby4y1lqFJFM1Q/G/B9gDOErSHnW7HQ+siYhdgXOBNoya1F45r/NOYHpE/DFwNfAvxUbZnXK+993iEuDgunWzgQURsRuwIC3bOHVRebqE7isr64FZEbEHsC9wQvrbVPa6Kl6e3hUR02qGGyj7fbyE/GX2PcBu6TETuLCgGF+iqJqpPEPxHw7MS8+vBg6QpILia5UxrzMiFkbEb9PiT8nGK7Hm9cx0DxFxC/CbutW1/3/mAUcUGlT36Yry1I1lJSJWRsQd6fnTwFKy0cirfF2dVJ5KfR8bLLOHA5dG5qfANpImFxPpxopKpvIMxb9hn4hYDzwJbF9IdK3T6JQDxwPfa2tEvaPXp3voi4iV6fljQF+ZwXSBbi5PXVNWUneQvYHbqPZ1VbU8BfADSYvS6O1QzfdxpJgq8762fToZG56kDwPTgT8rOxbrLhERknybro2pk8uKpEnANcDJEfFUbUNGJ19Xwd4eESskvQq4WdIvajdW8X2sYkxQXM1UnqH4N+wjaSKwNfBEIdG1Tq4pByQdCHwGOCwinisotm6X673vYquGqrfTv6tLjqfTdXN56viyImkTskRqfkRcm1ZX+boqWZ4iYkX6dzVwHVlzZBXfx5Fiqsz7WlQylWco/huAY9PzDwA/is4bBGvM65S0N/A1skSqCoW0W/T6dA+1/3+OBa4vMZZu0M3lqaPLSupLexGwNCK+VLOpytdVufIkaUtJWw09B/4cuJtqvo8jxXQD8JF0V9++wJM1zYHFiohCHsAhwC+BB4HPpHWfI0sqADYDvg08APw38JqiYiv4On8IrAIWp8cNZcfcLY/h3vtufACXAyuBP5D1ETierH/hAuD+VMa2KzvOTn90Q3nqxrICvJ2sr89dNZ+jh1T9uqpWnoDXAD9Pj3tqvq9KfR8bKbOAyO6SfBBYQnanfCnvp0dANzMzM2uCR0A3MzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7MmOJkyMzMza4KTKTMzM7Mm9HQyJWmZpAPLjgNAUr+kkDQxLQ9K+ljZcXWz8f79JV0i6cx2xFQlkr4n6diy4+g2LnfWSi5P1dDTyZRZt5A0IOkFSeskPS3pPknHNXD8GZK+VbsuIt4TEfNaH62ZWXdxMtUkZcZ8H4dqnMza6NGImAS8Avgk8HVJu5cck5lZ13MyBW+RdK+kNZK+KWkzSdtKulHSr9P6GyVNHTogNcGdJenHwG+B19SfVNIMST+WdK6kJ4AzJL1M0mclPSxptaRLJW1d4LXaSw33958h6dbanVIT7K7DnUDSqZJWSnpU0sdq95X0cklflPQrSask/ZukzUcKRtJsSQ+m2qV7Jb2/ZtuFkq6pWT5b0gJJqj1HZG4CfgP8cc3+50l6RNJTkhZJekdafzDwaeCDqWbr52n9hqZml92Wq1q5Wybp/0i6S9Izki6S1Jeaep+W9ENJ26Z9ByQtH+b4A9PzmySdU7PtCkkXj/+tshy6sjxJ+iNJv5W0fc22Nyv7bt6kmTes1ZxMwdHAu4HXAq8DPkv2vnwTeDWwM/AscEHdcccAM4GtgIdHOPdbgYeAPuAsYEZ6vIssAZs0zHmtWMP9/XNLicgpwIHArsBA3S5z0nmnpe1TgH8c5ZQPAu8Atgb+CfiWpMlp2yxgz/Qh+Q7geODYiIi6mF4m6TBgB+CBmk0/S3FsB1wGfFvSZhHxfeALwJURMSki9homrhm47LZS1codwF8CB6Xj3gd8jyzJfiXZZ+Incob3UeAYSftLOhrYBzgp57E2Pl1ZniLiMWAQOLJm9THAFRHxh7GOL1RE9OwDWAb8bc3yIcCDw+w3DVhTszwIfG6Mc88AflW3bgHwdzXLuwN/ACYC/UAAE2te42Nlv0fd/Bjp75/+drfW7RvArun5JcCZ6fnFwD/X7Lfr0L6AgGeA19ZsfxvwPw3EuBg4vGb5rWQ1Tg8DR9WsHwBeANYCzwHPAyePce41wF7p+RnAt+q2byiDo5Xdsv+OnfaoYrlLMR1ds3wNcGHN8seB79SUteXDHH9gzfJfAo8AjwNvL/s97+ZHt5cn4IPAj9PzCcBjwD5lv+/1D9dMZf/hhzwM7ChpC0lfS00aTwG3ANtImjDccZLekZpH1km6Z4RzA+zIxrVYD5MlUn0tuRIbj5f8/Rs8fse6c9Q+fyWwBbBI0lpJa4Hvp/VDd8sNlZuj07qPSFpcs/+byGqYAIiI28hqOwVcVRfLoxGxDVmfqa8A+9dulPT3kpZKejKde+vac+e4Tpfd1qlUuUtW1Tx/dpjlSQ3E9x9kX3z3RcStY+1sTevm8nQ9sIekXchqup6MiP/OeWxh3Ckadqp5vjPwKFlzyu7AWyPiMUnTgDvJvsCGbGhaiYj/YviCEXXLj5I1Hda+3nqyQjYVK8Nwf/9nyD48AJD0R6Mcv5KN/3a153uc7EPjjRGxov7AiHhP7bKkVwNfBw4AfhIRz0taTE25k3QC8PIU56nAPw9z3ucknQbcJ+mIiPhOahY8NZ37noh4QdKamnPXl9V6o5Vda1xlyt041Mc5gfTFWuMsYCmwi6SjIuLyJl/TRte15SkififpKuDDwOuBf2/y9drCNVNwgqSpkrYDPgNcSdYP6llgbVp/eote63Lgk5J2kTSJF/uprG/R+a1xw/39fw68UdI0SZuRNYGN5CrgOElvkLQF8A9DGyLiBbLk6FxJrwKQNEXSu0c415ZkSc2v077HkdVMkZZfB5xJ9qFyDHBqSvRfIiJ+D5zDi/0atiJLfn4NTJT0j2Q1WENWAf0a+c5Ul93WqlK5a9Qvgc0kHZo6AX+WLMEnvdY7geOAjwDHAudLmtKi17bhdW15Si7moPvaAAAgAElEQVQla7Y8DCdTlXUZ8AOyppMHyb6svgxsTpaR/5SsSrMVLiYrCLcA/wP8jqzt2Mrzkr9/RPwS+BzwQ+B+YMRmioj4HlmT2kKyzt4/TZueS/+eNrQ+NRn/kKzWc7hz3UuWAP2ELLnZE/gxbBha41vA2RHx84i4n6wz579Lqv/gGXIxsLOk9wH/SVaOf0nWDPA7Nq7K/3b69wlJd4xwLpfd1qlMuWtURDwJ/B3wDWAFWc3CcgBJryD74jsxIlakWvuLgG9KG991ai3VleWpZp8fk/UJvSMiRrrhq1RKnbrMrAUkvQG4G3i5a22sKC531kpVLE+SfgRcFhHfKDuW4bhmyqxJkt6vbByWbYGzgf+oygeQdS+Xu/GTdLGy8dLurlm3naSbJd2f/h0aB0mSviLpAWXjJr25vMjbp8rlSdJbgDeTNV9WkpMps+b9DbCarHr9eeB/lxuO9QiXu/G7BDi4bt1sYEFE7EY2FMjstP49wG7pMRO4sKAYi1bJ8iRpHlmz4skR8XTZ8YzEzXxmZtZzJPUDN0bEm9LyfcBARKxUNlDuYETsLulr6fnl9fuVFLpVkGumzMzMoK8mQXqMF8dQm8LGN2ssT+vMNnAyZWYdT9IESXdKujEt7yLpttTP5UpJm5Ydo3WOyJpsGm62kTRT0u3pMbMNoVlFFTpo5w477BD9/f0brXvmmWfYcsstiwyjVN1wvYsWLXo8IuoH6StcJ5anKsdXVmwtKk8nkQ0SOTR21tnAuRFxhaR/I5vHcNS+Lp1YnlqpW661ifK0StLkmma+1Wn9CjYexHJqWvcSETEXmAtZeZo+ffrXard3y3sMvXMtuctTkXPX/Mmf/EnUW7hw4UvWdbNuuF7g9qjAXEidWJ6qHF9ZsTVbnsi+3BaQTZ9zI9mo7o/z4jyXbwP+c6zzdGJ5aqVuuda85YlsPtS7a5b/FZidns8G/iU9P5Rskl4B+wL/nef83V6eeuVa8pYnTydjZp3uy2RT5WyVlrcH1saLt3W7j4ttRNLlZBPs7iBpOdksF3OAqyQdTzaw7ZFp95vIJg9+APgt2ejuZhtxMmWFStMa3EI2XcBE4OqIOF3ZJJZXkH0RLgKOiWxKFLMRSXovsDoiFkkaGMfxM8lud6evr4/BwcGNtq9bt+4l67pVL11rRBw1wqYDhtk3gBPaG5F1OidTVrTngP0jYl2ah+lWSd8DTqHBPi5mwH7AYZIOATYj6zN1HrCNpImpdipXH5fp06fHwMDARtsHBwepX9eteulazVrNd/NZoVIz9Lq0uEl6BFl/l6vT+nnAESWEZx0mIj4VEVMjoh/4EPCjiDiabI6xD6TdjgWuLylEM+sBpddMLVnxJDNmf3fM/ZbNObSAaKwIkiaQNeXtCnyVbMTdXH1cWtkss2TFk7n223PK1rn2y6PKTSlVjm0cTgOukHQmcCfZZLtt1Z/jcwz8WdbL8n7fgctJpyk9mbLeExHPA9MkbQNcB7y+gWNb1iyT+0Pt6Hzny6PKTSlVji2PiBgEBtPzh4B9yozHzHqHm/msNBGxlqw55m2kPi5p04h9XMzMzKrGyZQVStIrU40UkjYHDiIbbNF9XMzMrCO5mc+KNhmYl/pNvQy4KiJulHQvBfdxMTMzawUnU1aoiLgL2HuY9e7jYmZmHSl3M58nEjUzMzN7qUb6TA1NJDpkaCLRXYE1ZIMsmpmZmfWUXMmUpKlkkz1+Iy0LD7JoZmZmlrtmamgi0RfSsicSNTMzMyNHB/R2TyTatznM2nP9MEdurFtGZu6yUabNzMx6Xp67+do6kej586/nnCVjh9HKUajL1OmjTJuZmdnGxmzm80SiZmZmZiNrZgT004BTJD1A1ofKgyyamZlZz2lo0E5PJGpmZma2Mc/NZ2ZmZtYEJ1NmZmZmTXAyZWZmZtYEJ1NmZmaJpE9KukfS3ZIul7SZ56K1sTiZskJJ2knSQkn3pg+sk9L6MyStkLQ4PQ4pO1Yz6y2SpgCfAKZHxJuACWRDAnkuWhuVkykr2npgVkTsAewLnCBpj7Tt3IiYlh43lReimfWwicDmkiYCWwAr8Vy0NgYnU1aoiFgZEXek508DS/G8jmZWARGxAvgi8CuyJOpJYBGei9bG0NA4U2atJKkf2Bu4jWzaohMlfQS4naz2ak150ZlZr5G0LXA4sAuwFvg2cHADx7dkLlqo/ny03TTPbCuuxcmUlULSJOAa4OSIeErShcDngUj/ngN8dJjjRv2wauQ/RRkfalX+AKpybGYFORD4n4j4NYCka8l+6BU6Fy1Ufz7abppnthXX4mTKCidpE7JEan5EXAsQEatqtn8duHG4Y8f6sGrkP8WM2d/NtV8rP9Sq/AFU5djMCvIrYF9JWwDPAgeQ1ZQPzUV7BZ6L1obhPlNWKEkim8dxaUR8qWb95Jrd3g/cXXRsZtbbIuI2so7mdwBLyL4j5+K5aG0Mrpmyou0HHAMskbQ4rfs0cJSkaWTNfMuAvyknPDPrZRFxOnB63WrPRWujcjJlhYqIWwENs8lDIZiZWUdyM5+Zdaw0OvV/S/p5GgT2n9J6j1htZoVxMmVmnew5YP+I2AuYBhwsaV88YrWZFcjJlJl1rMisS4ubpEfgEavNrEBOpsyso0makG5mWA3cDDyIR6w2swK5A7qZdbSIeB6YJmkb4Drg9XmP7fRBYFvJg7aajZ+TKTPrChGxVtJC4G20aMTqqg8C20oetNVs/NzMZ2YdS9IrU40UkjYHDiKbPHtoxGrwiNVm1maumTKzTjYZmCdpAtmPw6si4kZJ9wJXSDoTuBOPWG1mbeRkysw6VkTcBew9zHqPWG1mhXEyZTaG/rx9YeYc2uZIzMysisbsM+URhs3MzMxGlqcDukcYNjMzMxvBmMmURxg2MzMzG1muoRE8wrCZmZnZ8HJ1QG/nCMN9m+cbObhbRubt9VGGJe0EXAr0kdVwzo2I8yRtB1wJ9APLgCMjYk1ZcZqZmeXV0N187Rhh+Pz513POkrHDqOqowY3yKMOsB2ZFxB2StgIWSboZmAEsiIg5kmYDs4HTSozTzMwslzx383mEYWuZiFgZEXek50+TlaUpwOFkfe/AffDMzKyD5KmZ8gjD1haS+skGXLwN6IuIlWnTY2TNgGZmZpU3ZjLlEYatHSRNAq4BTo6IpyRt2BYRISlGOG7UPniN9EnL01evEXlet8p95qocW1mWrHgy9wTGZta7PAK6FU7SJmSJ1PyIuDatXiVpckSslDSZ7M7RlxirD14jfdJa/SWZp19flfvMVTk2s6Kkbi3fAN5EdpPMR4H78A0yNopcQyOYtYqyKqiLgKUR8aWaTTeQ9b0D98Ezs/KcB3w/Il4P7EXWr3M22Q0yuwEL0rLZBk6mrGj7AccA+0tanB6HAHOAgyTdDxyYls3MCiNpa+CdpD7AEfH7iFiLb5CxMbiZzwoVEbcCGmHzAUXGYmZWZxfg18A3Je0FLAJOwjfI2BicTJmZmWUmAm8GPh4Rt0k6j7omvWZukMk7SDVUf6DqbrphpRXX4mTKzMwssxxYHhG3peWryZKpltwgk3eQaqj+QNXddMNKK67FfabMzMyAiHgMeETS7mnVAcC9+AYZG4NrpszMzF70cWC+pE2Bh4DjSANWSzoeeBg4ssT4rIKcTJmZmSURsRiYPswm3yBjI3Izn5mZmVkTXDNlZlag/pwj7y+bc2ibIzGzVnHNlJmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZWceStJOkhZLulXSPpJPS+u0k3Szp/vTvtmXHambdy8mUmXWy9cCsiNgD2Bc4QdIeZPOpLYiI3YAF1E1Wa2bWSk6mrFCSLpa0WtLdNevOkLRC0uL0OKTMGK1zRMTKiLgjPX8aWApMAQ4H5qXd5gFHlBOhmfUCJ1NWtEuAg4dZf25ETEuPmwqOybqApH5gb+A2oC8iVqZNjwF9JYVlZj3AI6BboSLilvSlZ9YykiYB1wAnR8RTkjZsi4iQFCMcNxOYCdDX18fg4OBG2/s2h1l7rm9X2KOqj6Xd1q1bV/hrmnULJ1NWFSdK+ghwO1kfmDVlB2SdQdImZInU/Ii4Nq1eJWlyRKyUNBlYPdyxETEXmAswffr0GBgY2Gj7+fOv55wl5XxMLjt6YMx9WmlwcJD66zezfJxMWRVcCHweiPTvOcBHh9txrJqERn5dt7rGIc/rVvnXf5VjG4myKqiLgKUR8aWaTTcAxwJz0r/XlxCemfWIMZMpSTsBl5L1OQhgbkScJ2k74EqgH1gGHOnaBBuPiFg19FzS14EbR9l31JqERn5dz8g54WxeeWoSqvzrv8qxjWI/4BhgiaTFad2nyZKoqyQdDzwMHFlSfGbWA/LUTA3denyHpK2ARZJuBmaQ3Xo8R9JssluPT2tfqNathppj0uL7gbtH299sSETcCmiEzQcUGYuZ9a4xk6n0JbcyPX9aUu2txwNpt3nAIE6mbAySLicrNztIWg6cDgxImkZW87kM+JvSAjQzM2tQQ32mfOuxNSsijhpm9UWFB2JmZtYiuZOpsm897rSOsSPpxE6+lk9/jj5Ys/Zcz4zZ32XZnEMLiMjMzIqQK5mqwq3HRd8m3C4d2snXzMzMRjDmCOg5bj0G33psZmZdQtIESXdKujEt7yLpNkkPSLpS0qZlx2jVkmc6maFbj/evmzttDnCQpPuBA9OymZlZpzuJbJ7HIWeTTXm1K7AGOL6UqKyy8tzN51uPzcysJ0iaChwKnAWcklpn9gf+Ku0yDziDbLBhM8ATHZuZmdX6MnAq8EJa3h5YGxFDd0otJxseyGwDTydjZmYGSHovsDoiFkkaGMfxLZs4u+p3fXfTnemtuBYnU2ZmZpn9gMNSv+DNgFcA5wHbSJqYaqemAiuGO7iVE2dX/Q72brozvRXX4mY+MzMzICI+FRFTI6If+BDwo4g4GlgIfCDt5rvX7SVcM2VWgjwDfAIe3NOsGk4DrpB0JnAnnrXB6jiZMjMzqxMRg2RzzhIRDwH7lBmPVZub+czMzMya4GTKzMzMrAlOpszMzMya4GTKCiXpYkmrJd1ds247STdLuj/9u22ZMZqZmTXCyZQV7RLg4Lp1s4EFEbEbsCAtm5mZdQQnU1aoiLgF+E3d6sPJ5rsi/XtEoUGZmZk1wcmUVUFfRKxMzx8D+soMxszMrBEeZ8oqJSJCUoy0fay5rxqZYynvHFmt1MjcXFDs/FzdNNeWmVmRnExZFaySNDkiVkqaDKweacex5r5qZI6lGTlHIW+lWXuuzz03FxQ7P1c3zbVlZlYkN/NZFdxANt8VeN4rMzPrME6mrFCSLgd+Auwuabmk44E5wEGS7gcOTMtmZmYdwc18VqiIOGqETQcUGoiZmVmLuGbKzMzMrAlOpszMzMya4GTKzDqapygys7I5mTKzTncJnqLIzErkDujWVZaseLKU8aOsPBFxi6T+utWHAwPp+TxgEDitsKDMrKeMWTPlKnQz60CeosjMCpOnZuoS4ALg0pp1Q1XocyTNTsv+1WdmlTPaFEVjTU/U6PQ/rVT01D6eTshs/MZMplyFbmYdKNcURWNNT3T+/Osbmv6nlYqcSgg8nZBZM8bbAd1V6GZWZZ6iyBomaSdJCyXdK+keSSel9e7aYqNq+ifXaFXo0Lpq9G6pfnZVullrpSmKBoAdJC0HTiebkuiqNF3Rw8CR5UVoHWQ9MCsi7pC0FbBI0s3ADNy1xUYx3mQqVxU6tK4avegq73ZxVbpZa3mKImuV1OKyMj1/WtJSYAru2mJjGG8zn6vQzcysa6W+wnsDt+GuLTaGMauEXIVuZma9RNIk4Brg5Ih4StKGbUXdHVr17iDd1GWlFdeS524+V6GblaQ/5wCky+Yc2uZIzHqDpE3IEqn5EXFtWl343aFV79rSTV1WWnEtnk7GKkPSMklLJC2WdHvZ8ZhZb1FWBXURsDQivlSzyV1bbFSeTsaq5l0R8XjZQZhZT9oPOAZYImlxWvdp3LXFxuBkyszMDIiIWwGNsNldW2xEbuazKgngB5IWpY6cZmZmleeaKauSt0fECkmvAm6W9IuIuKV2hyrPpZZHu+JrxV013XR3jplZkZxMWWVExIr072pJ1wH7ALfU7VPZudTymLXn+rbE14o7f7rp7hwzsyK5mc8qQdKWafoGJG0J/Dlwd7lRmZmZja26P+Gt1/QB16XB8SYCl0XE98sNyczMbGxOpqwSIuIhYK+y4zAzM2uUm/nMzMzMmuBkyszMzKwJTqbMzMzMmuBkyszMzKwJ7oBuZlZB/bO/m3vfZXMObWMkZjYW10yZmZmZNcHJlJmZmVkTnEyZmZmZNaGn+0x1Sp+EvHG634SZmVnxejqZMutFIyXns/Zcz4yabU7OzczycTOfmZmZWRNcM9ViYzXJ1f/6L0OnNG+amZl1AtdMmZmZmTXBNVNmXaCR2sYyX9s1nWbWjZpKpiQdDJwHTAC+ERFzWhKV9SSXJ2ulXipPrUimx9sFodUJclWT814qT9a4cSdTkiYAXwUOApYDP5N0Q0Tc26rgapX9H6zMX/69oOjyZN3N5clayeXJxtJMzdQ+wAMR8RCApCuAwwEXLhsPlydrJZcnayWXpw6VpyJkqFa2mYqYZjqgTwEeqVlentaZjYfLk7WSy5O1ksuTjartHdAlzQRmpsV1ku6r22UH4PGWvubZrTxba32iDdc7pB3XPcI5X936V8qnjPLUSu38+zerPjaXJ6DCf69WG2/ZLPPztpvLU5W/x5Ku+b8xVPabKU/NJFMrgJ1qlqemdRuJiLnA3JFOIun2iJjeRBwdpdeutwE9UZ6qHF+VYxuHnihPrdRL1zoOLk91fC0ba6aZ72fAbpJ2kbQp8CHghmaCsZ7m8mSt5PJkreTyZKMad81URKyXdCLwn2S3il4cEfe0LDLrKS5P1kouT9ZKLk82lqb6TEXETcBNTcYwYpVol+q1682tR8pTleOrcmwN65Hy1Eq9dK0Nc3l6CV9LDUVEKwIxMzMz60mem8/MzMysCaUmU5IOlnSfpAckzS4zlnaTdLGk1ZLuLjuWbjBW2ZH0cklXpu23SeqvWHzvlHSHpPWSPlCx2E6RdK+kuyQtkFTareZl6pXPJ0k7SVqY/ub3SDqp7Ji6USeXp5HKiKTtJN0s6f7077Zlx5qXpAmS7pR0Y1reJX1XPJC+OzZt5HylJVM1w/O/B9gDOErSHmXFU4BLgIPLDqIb5Cw7xwNrImJX4FygsFFbcsb3K2AGcFlRcTUQ253A9Ij4Y+Bq4F+KjLEKeuzzaT0wKyL2APYFTujiay1FF5SnkcrIbGBBROwGLEjLneIkYGnN8tnAuek7Yw3Zd0huZdZMbRiePyJ+DwwNz9+VIuIW4Ddlx9El8pSdw4F56fnVwAGSVJX4ImJZRNwFvFBQTI3EtjAifpsWf0o2pk6v6ZnPp4hYGRF3pOdPk33BeHTv1uro8jRKGan9nJ0HHFFOhI2RNBU4FPhGWhawP9l3BYzjWspMpjw8v41XnrKzYZ+IWA88CWxfSHTVLtuNxnY88L22RlRNVf4btk1qDt8buK3cSLpO15SnujLSFxEr06bHgL6SwmrUl4FTefHH7PbA2vRdAeP4+7gDupkNS9KHgenAv5Ydi7WfpEnANcDJEfFU2fFY9YxWRiIbGqDywwNIei+wOiIWtfK8bZ+bbxS5huc3G0aesjO0z3JJE4GtgSeKCa/SZTtXbJIOBD4D/FlEPFdQbFVS5b9hy0nahOxLcn5EXFt2PF2o48vTCGVklaTJEbFS0mRgdXkR5rYfcJikQ4DNgFcA5wHbSJqYaqca/vuUWTPl4fltvPKUnRuAY9PzDwA/iuIGVaty2R4zNkl7A18DDouITvhwbIcq/w1bKvUXuQhYGhFfKjueLtXR5WmUMlL7OXsscH3RsTUqIj4VEVMjop/s7/CjiDgaWEj2XQHjuJbSkqmU/Q0Nz78UuKqbh+eXdDnwE2B3ScslNXSngL1opLIj6XOSDku7XQRsL+kB4BQKvMskT3yS3iJpOfC/gK9JKqTs53zv/hWYBHxb0mJJHfOh3yo99vm0H3AMsH/6ey9Ov9qtRbqgPI1URuYAB0m6HzgwLXeq04BT0nfG9mTfIbl5BHQzMzOzJrgDupmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZmZmZNcHJlJmZmVkTnEyZmZmZNcHJlJmZmVkTuj6ZkrRM0oHjOO4SSWe2I6bxkHSmpMclPZaW3y/pEUnrJO1ddny9pBvKVG0skgbSpMvWZkPvu6R3SLovx/5nSPpWEbG1Wvpsek3ZcVhvlbuydH0y1SnSF/Sz6QNo6HFB2rYzMAvYIyL+KB3yReDEiJgUEXeWFbeZNS4i/isidi87jiHpy/MP6XNnraT/J+ltDRw/KOljtevSZ9NDrY/Wxqtq5a6bOJkqmaSJNYvvSx9AQ48T0/qdgSciYnXNvq8G7iksUOsYdWXKLK8rI2ISsAOwEPh2yfGYdYxeSabeIuleSWskfVPSZpJmSLq1didJIWnX4U4g6VRJKyU9KuljtftKermkL0r6laRVkv5N0uYjBZNqoU6TdBfwzGhffqk56WZgx/Sr8XJJ64AJwM8lPTiO98OaV/kyJekNqcZgraR7JB3W0nfAxiRpb0l3SHpa0pXAZmn9Rk2r6W+3Iu13n6QDhjnXJun//zWSNh3h9b4t6TFJT0q6RdIb0/pNJS2W9PG0PEHSjyX9Y/05ImI9MB+YIumVaf9tJd0o6depzN8oaWradhbwDuCCuhr12vK8taRL0/EPS/qspF75/ilckeVOUn/6Wx+nrOvJGkl/K+ktku5Knz8X1Oy/URNizfETJW0nabmk96VtkyQ9IOkjLX6LWq5XCvPRwLuB1wKvAz7byMGSDgZOAQ4EdgUG6naZk847LW2fArzkQ6rOUcChwDbpw2tYEfFD4D3Ao6m26qj06xFgr4h4bSPXYi1T6TIFCPgP4AfAq4CPA/MluYq/IOmL5zvAvwPbkdX0/OUw++0OnAi8JSK2IitXy+r22Tyd6zngyIj4/Qgv+z1gN7K/+R1kSRFp/w8Dn5P0BmA22Q+ys0aI+yPAE8CatPplwDfJasR3Bp4FLkjn/gzwX7zY7eDE+nMC5wNbA68B/iyd/7gRrsGaUFK5A3grWdn7IPBl4DNkn29vBI6U9GdjxR4RvwE+Cnxd0quAc4HFEXHpWMeWrVeSqQsi4pH0hzqL7EunEUcC34yIeyLit8AZQxskCZgJfDIifhMRTwNfAD40xjm/kmJ6tmbdd1IWP/T46wbjtOJUvUztC0wC5kTE7yPiR8CN44jTxm9fYBPgyxHxh4i4GvjZMPs9D7wc2EPSJhGxLCJqa5xfAXwfeBA4LiKeH+kFI+LiiHg6Ip4jK1N7Sdo6bbsbOJPsy/HvgWPqznWkpLVkidJfAx8Y+qEXEU9ExDUR8dtUHs8iS4rGJGkCWdn9VIptGXAOcEye461hhZe75PMR8buI+AHwDHB5RKyOiBVkyXauG6XS8d8GFgCHAH+T57iy9Uoy9UjN84eBHRs8fse6c9Q+fyWwBbBoKAkiK4BD1ePf04sdyo8e4RxDjoiIbWoeX28wTitO1cvUjsAjEfFCXZxTGozTxm9HYEVERM26h+t3iogHgJPJkp/Vkq6QVFue9gX+mCwx3nAubXyzys6p6W6OpAclPcWLtQw71JxrHlnt0k0RcX9dKFdFxDZAH3A38Cc1r7WFpK+lJrqngFuAbVKiNJYdyL7ca6/dZbF9Ci13Nfuvqnn+7DDLk8hvLvAm4JKIeKKB40rTK8nUTjXPdwYeJcuctxhaKemP6g+qsRKYOsL5HicrKG+sSYK2HmqKi4j31HQon19zXG1Bt85T9TL1KLBTXb+UnYEVOa7NWmMlWb8j1azbebgdI+KyiHg7WaITwNk1m38A/DOwQFJfzTG1N6v8Cvgr4HCyppWtgf60a+3r/1+yGsp3S3r7CLE8TlYzeoakyWn1LGB34K0R8QrgnXXnHu3z7HHgD+nahrgstk/R5a5RG31OAht9TqYEfS5wKfB3GqHPadX0SjJ1gqSpkrYja8e9Evg58EZJ0yRtRk0zyzCuAo5T1qF3C+AfhjakX/5fB85NbbxImiLp3W26FquGqpep24DfAqemDqQDwPuAKxo4hzXnJ8B64BPpb/AXwD71O0naXdL+kl4O/I4ska6tUSQi/gW4jOyLbYf6cyRbkfVteYLsy+oLda9zDFlt0wzgE8A8ScPWFkTEfcB/AqfWnPtZYG0q86fXHbKKrD/UcOd6nqy8nyVpK0mvJusv6HGM2qPocteoxcA7U23q1sCn6rZ/miyx+yjwr8ClOWtAS9UrydRlZFn2Q2Ttv2dGxC+BzwE/BO4Hbh3p4Ij4HvCV/9/e/UdLXtd3nn++FlAJqICQOx1Amz0yOkRGyHYILolzA5IhYoRkidEwBBLOdpITDa49J7TuOaOJ2d1mV0TG5CTpCEM7g4pBHAi6STqEu45nEyIg4VfHAZk2dk9DxwjKZTOaJu/9o76Nl+v9Ufd+69a36tbzcc49t74/qupdVe+qetfn+/l8vvSGCz8K/EWz6VvN/ysPrm+awP+U3q+4lfrDeU2on17FbWg4Rjqnmo6iP0Fv8MLX6LVI/FxV/XW/t6F2mtfgp+gVL1+n1zH3lgV2fSG9AQdfAx6n13l8/hcMVfV+ev2d/rQpaOb7KL3DOXuBh/lOTh2cq+5D9HJgtqo+BtxNr4PvYv4vYHNT0H8IOLyJ8S/oHXae61rgovRGcv3bBW7rHfRaJB6j9774GHD9EvetVeog71Ya3056Pz7vB+6h11IKQJL/gV6h/XNNEX4VvcJqa9v7XWt5/mFV9aMZDfMg8MJaYiSe1C9zSpLG126/oKsAACAASURBVKS0TLWW3ulbXpjkaHrV8h/6pac2zClJWh8spvr3i8B+eod0ngV+udtwtA6YU5K0DniYT5IkqQVbpiRJklqwmJIkSWphqGeXP/bYY2vjxo3PW/fMM89wxBFHDDOMgZnU2O+5556vVdVxAw5pxdZbPrU1ro99lPNpXIzra3/QIOM3n1bHHFpYv/k01GJq48aN3H333c9bNzMzw/T09DDDGJhJjT3Jd52aoAvrLZ/aGtfHPsr5NC7G9bU/aJDxm0+rYw4trN988jCfJElSCxZTkiRJLVhMSZIktWAxJUmS1ILFlCRJUgtDHc23kAf2foPLtn5m2f12bzt/CNForSV5EfA5emcsPxS4uarem+Qk4BPAy+idSfyS5uznK9JvPoE5peVt7DOXwHzSYPWbe+bdaLBlSsP2LeDsqnotcBpwXpIz6Z3o95qqeiXwJHB5hzFKktQ3iykNVfXMNouHNX8FnA3c3KzfAVzYQXiSJK2YxZSGLskhSe4D9gM7gS8DT1XVgWaXPcDxXcUnSdJKdN5nSpOnqp4FTktyFPBp4NX9XjfJZmAzwNTUFDMzM8/bPnU4bDn1wALX/G7zrzvuZmdn191jkqRxYDGlzlTVU0nuBF4HHJXk0KZ16gRg7yLX2Q5sB9i0aVPNP33Ah2+8lasf6C+td188vew+42TcTwchSePKw3waqiTHNS1SJDkcOBfYBdwJXNTsdilwazcRSpp0TVeELya5vVk+KcldSR5NclOSF3Qdo0ZLX8VUkv8lyUNJHkzy8SQvMrm0ShuAO5PcD3wB2FlVtwNXAu9K8ii96RGu6zBGSZPtCno/8g5ytLGWtGwxleR44FeBTVX1GuAQ4K2YXFqFqrq/qk6vqn9eVa+pqt9o1j9WVWdU1Sur6qer6ltdxypp8iQ5ATgf+EizHBxtrGX0e5jvUODwJIcC3wPsw+SSJK0/HwJ+DfjHZvllONpYy1i2p25V7U3yAeBvgL8H/oTeDNUmlyRp3UjyJmB/Vd2TZHoV119ytPFKDHtU8riPBu46/mWLqSRHAxcAJwFPAX8AnNfvHQxqKPsovshdv3htjHPskrRGzgLenOSNwIuAlwDXMqDRxivR92mxBjQqedxHA3cdfz9jyN8A/Jeq+luAJLfQS7ihDmUfxWHsXb94bYxz7JK0Fqrq3cC7AZqWqX9dVRcn+QN6o40/gaONtYB++kz9DXBmku9pOuKdAzyMQ9klSZPB0cZa0rLFVFXdRa+j+b3AA811tmNySepYkhOT3Jnk4Wb6liua9cck2Znkkeb/0V3HqvFSVTNV9abmsqONtaS+poquqvcC7523+jHgjIFHJEn9OwBsqap7k7wYuCfJTuAy4I6q2pZkK7CV3g9ASRo4Z0CXNLaqal9V3dtcfpreRIvH0xs0s6PZzalbJK0piylJ60KSjcDpwF3AVFXtazY9Dkx1FJakCeCJjiWNvSRHAp8C3llV3+yNlempqkpSi1xvyalb+p3rB/qfvuWBvd/oa79Tj39p3/c97lOdjHv8ksWUpLGW5DB6hdSNVXVLs/qJJBuqal+SDcD+ha673NQt/c71A/1P37IW8weN+1Qn4x6/5GE+SWOrma7lOmBXVX1wzqbb6E3ZAk7dImmN2TIlaZydBVwCPJDkvmbde4BtwCeTXA58BXhLR/FJmgAWU5LGVlV9Hsgim88ZZiySJpeH+TRUS0yy+L4ke5Pc1/y9setYJUnqhy1TGrbFJlkEuKaqPtBhbJIkrZjFlIaqmftnX3P56SQHJ1mUJGkseZhPnZk3ySLA25Pcn+R6z6UmSRoXtkypEwtMsvg7wPuBav5fDfzCAtdbcpLFqcP7n2hxvU0S6MSHktQNiykN3UKTLFbVE3O2/z5w+0LXXW6SxQ/feCtXP9BfWq9kUsRx4MSHktQND/NpqBabZLGZpfqgnwQeHHZskiSthi1TGrbFJll8W5LT6B3m2w38YjfhSZK0MhZTGqolJln87LBjkSRpEDzMJ0mS1ILFlCRJUgsWU5IkSS1YTEmSJLVgB3RJGoCNWz/TdQiSOmLLlCRJUgsWU5IkAUlelOQvk/xVkoeS/Hqz/qQkdyV5NMlNSV7QdawaLX0VU0mOSnJzkr9OsivJ65Ick2Rnkkea/56YVpI0zr4FnF1VrwVOA85LciZwFXBNVb0SeBK4vMMYNYL6bZm6Fvijqno18FpgF7AVuKOqTgbuaJYlSRpL1TPbLB7W/BVwNnBzs34HcGEH4WmELdsBPclLgdcDlwFU1beBbye5AJhudtsBzABXrkWQkiQNQ5JDgHuAVwK/DXwZeKqqDjS77AGOX+S6m4HNAFNTU8zMzKw6ji2nHlh+J2h1H3PNzs4O7La60HX8/YzmOwn4W+DfJXktvSS7Apiqqn3NPo8DU2sToiRJw1FVzwKnJTkK+DTw6hVcdzuwHWDTpk01PT296jgu63N06O6LV38fc83MzNAm3q51HX8/xdShwA8A76iqu5Jcy7xDelVVSWqhKy9XqU8d3l8FPooVc9eVcBvjHLskrbWqeirJncDrgKOSHNq0Tp0A7O02Oo2afoqpPcCeqrqrWb6ZXjH1RJINVbUvyQZg/0JXXq5S//CNt3L1A8uHMajqe5C6roTbGOfYpYOSXA+8CdhfVa9p1r0P+J/ptagDvKeqPJG2lpXkOOAfmkLqcOBcep3P7wQuAj4BXArc2l2UGkXLdkCvqseBryZ5VbPqHOBh4DZ6SQUml/qU5MQkdyZ5uBl6fEWz3tGhWo0bgPMWWH9NVZ3W/FlIqV8bgDuT3A98AdhZVbfT6w/8riSPAi8DruswRo2gfmdAfwdwYzO3xmPAz9MrxD6Z5HLgK8Bb1iZErTMHgC1VdW+SFwP3JNlJb4DDHVW1LclWeq2fDmjQkqrqc0k2dh2H1oequh84fYH1jwFnDD8ijYu+iqmqug/YtMCmcwYbjta7ZtDCvuby00l20RsZ4+hQDdLbk/wccDe94v3JrgOStH55bj51pmlROB24C0eHanB+B3g/vfmB3g9cDfzCQjsuN0Cm3+Hpa2ElA0TGfUDJuMcvWUypE0mOBD4FvLOqvpnkuW3DGB0KozlCtA2/kHqq6omDl5P8PnD7EvsuOUCm3+Hpa2Elg27GfUDJuMcvWUxp6JIcRq+QurGqbmlWD3V0KIzmCNE2/ELqOZhHzeJPAg92GY+k9c9iSkOVXhPUdcCuqvrgnE0HR4duw9Gh6lOSj9Pra3dskj3Ae4HpJKfRO8y3G/jFzgKUNBEspjRsZwGXAA8kua9Z9x56RZSjQ7UiVfW2BVY7bF3SUFlMaaiq6vNAFtns6FBJ0thZdtJOSZIkLc5iSpIkqQWLKUmSpBYspiRJklqwA7om1sY+J2Tcve38NY5EkjTObJmSJElqwWJKkiSpBYspSZKkFuwzJQ1Iv32wwH5YkrSe2DIlSZLUgsWUJElSCxZTkiRJLVhMSZIktWAxpaFKcn2S/UkenLPufUn2Jrmv+XtjlzFKkrQSFlMathuA8xZYf01Vndb8fXbIMUmStGoWUxqqqvoc8PWu45AkaVD6nmcqySHA3cDeqnpTkpOATwAvA+4BLqmqb69NmJoAb0/yc/RybEtVPdl1QJImS5ITgY8CU0AB26vq2iTHADcBG4HdwFtW8xm1krnoNF5WMmnnFcAu4CXN8lX0Ds18IsnvApcDvzPg+DQZfgd4P70Pr/cDVwO/sNCOSTYDmwGmpqaYmZl53vapw2HLqQcGGtz8+1jMSu6339tcidnZ2TW5XWmCHKD3Y+7eJC8G7kmyE7gMuKOqtiXZCmwFruwwTo2YvoqpJCcA5wP/G/CuJAHOBn622WUH8D4sprQKVfXEwctJfh+4fYl9twPbATZt2lTT09PP2/7hG2/l6gcGO7H/7ounl90H4LKVzIDe522uxMzMDPOfD42vlbRi3HDeEWsYyeSoqn3Avuby00l2AccDFwDTzW47gBkspjRHv986HwJ+DXhxs/wy4KmqOvhTfA+9hJNWLMmG5kMM4CeBB5faX5LWWpKNwOnAXcDUnM+ox+kdBlzoOku2nA+61RwG18o97i3bXce/bDGV5E3A/qq6J8n0Su9gUIdlRvFF7vrFa6Or2JN8nN4vvGOT7AHeC0wnOY3eYb7dwC8OPTCNrSTXAwc/p17TrBtIHxdNpiRHAp8C3llV3+wdjOmpqkpSC11vuZbzlbRe92tQrdzj3rLddfz9tEydBby5mfvnRfT6TF0LHJXk0KZ16gRg70JXHtRhmbU4LNJW1y9eG13FXlVvW2D1dUMPROvJDcBv0es4fNBW7OOiVUhyGL1C6saquqVZ/cTBFvQkG4D93UWoUbTs1AhV9e6qOqGqNgJvBf6sqi4G7gQuana7FLh1zaKUpEUsMt3GBfT6ttD8v3CoQWksNf2BrwN2VdUH52y6jd73HPh9pwW0mWfqSnqd0R+l14fK1gVJo6KvPi7SPGcBlwBnzzsjwzbg3CSPAG9olqXnrGjYU1XN0BvFQFU9Bpwx+JAkaXCW6uPSRYfhtTDoPpAP7P1G3/ueevxLW9/fqPQ/rarPA1lk8znDjEXjZbBjyCVpNPTVx6WLDsNr4YbzjhhoH8hhT/Mxzv1PJbCYkrQ+Hezjsg37uGgd63c+st3bzl/jSCab5+aTNNaa6Tb+HHhVkj1JLsc+LpKGyJYpSWNtkek2wD4ukobElilJkqQWLKYkSZJasJiSJElqwWJKkiSpBYspSZKkFiymJEmSWrCY0lAluT7J/iQPzll3TJKdSR5p/h/dZYySJK2ExZSG7QbgvHnrtgJ3VNXJwB3NsiRJY8FiSkNVVZ8Dvj5v9QXAjubyDuDCoQYlSVILFlMaBVNVta+5/Dgw1WUwkiSthKeT0UipqkpSi21PshnYDDA1NcXMzMzztk8dDltOPTDQmObfx2JWcr/93uZKzM7OrsntSpKWZjGlUfBEkg1VtS/JBmD/YjtW1XZgO8CmTZtqenr6eds/fOOtXP3AYNN698XTy+4DcFmfZ29fyW2uxMzMDPOfD0nS2vMwn0bBbcClzeVLgVs7jEWSpBWxmNJQJfk48OfAq5LsSXI5sA04N8kjwBuaZUmSxoKH+TRUVfW2RTadM9RAJEkaEFumJEmSWrCYkiRJasFiSpIkqQX7TEnSmHtg7zdWNDWHFpfkeuBNwP6qek2z7hjgJmAjsBt4S1U92VWMGj3LtkwlOTHJnUkeTvJQkiua9Z6cVpK03tyA5w/VCvVzmO8AsKWqTgHOBH4lySmYXJKkdcbzh2o1lj3M15wzbV9z+ekku4Dj6SXXdLPbDmAGuHJNopSkFUqyG3gaeBY4UFWbuo1IY8zzh2pJK+ozlWQjcDpwFyaXpNH3o1X1ta6D0Pqx1PlDlzt36KDPG7oSy523c9zP7dl1/H0XU0mOBD4FvLOqvpnkuW1tkqvfE9OO4ovc9YvXxjjHLklD1tf5Q5c7d2iXgwSWOx/ouJ/bs+v4+yqmkhxGr5C6sapuaVYPJLn6PTHtWpwYtq2uX7w2xjl2qU8F/EnzQ+/3ms+i5xnlloSV6PdH6VoYxI+yMfhxd/D8odvw/KFawLJVTHpNUNcBu6rqg3M2mVySRtkPV9XeJN8L7Ezy103n4ueMckvCSmw59UBfP0rXwiB+6I7Sj7vm/KHTwLFJ9gDvpfc998nmXKJfAd7SXYQaRf28+84CLgEeSHJfs+49mFySRlhV7W3+70/yaeAM4HNLX0uTbr2eP3TjMj8Mtpx64LkfD7u3nT+MkNaVfkbzfR7IIpvHOrkkrU9JjgD+u2YE8hHAjwG/0XFYktYpZ0DXyBjVoezL/aIblfue+8tyKRPyq3MK+HQzUOZQ4GNV9UfdhiRpvbKY0qhxKLtaq6rHgNd2HYekyeCJjiVJklqwmNIoOTiU/Z5myLokSSPPw3waJcsOZR/UJLBd+/CN/c0ksuXU/m9znCfAlaRxZjGlkdHPUPZBTQK7HvU719AoToArSePMw3waCUmOSPLig5fpDWV/sNuoJEla3mT+hNcocii7JGksWUxpJDiUXRpP/c6FNiHzm2lCeZhPkiSpBYspSZKkFiymJEmSWrCYkiRJasEO6NKEGfSJm+1YLGnS2TIlSZLUgsWUJElSCxZTkiRJLdhnSpIkrZgTtn6HLVOSJEktWExJkiS1MDaH+VYynHsSmhQlaZws9Rm+5dQDXLaKKTv8rF8bg54+ZRLYMiVJktRCq5apJOcB1wKHAB+pqm0DiaolO8WNp1HNJy1tVFuNzScNkvmkpay6mEpyCPDbwLnAHuALSW6rqocHFZwmh/mkQTKfNEjmUztd/uAa1n23Ocx3BvBoVT1WVd8GPgFc0OL2NNnMJw2S+aRBMp+0pDaH+Y4HvjpneQ/wQ+3CGX/jcIhx49bP9N3hc4hxmk8aJPNJg2Q+aUlrPpovyWZgc7M4m+RL83Y5FvjaWsexkFzV+iZWHfsA7ruVX+0z9kXifMWg4+nXKOdT1/p9Tbs0hvk0FsbhtV/KauM3nwZnVHKoxXdj6/jb5FObYmovcOKc5ROadc9TVduB7YvdSJK7q2pTizg6Y+wDNfH51NYkP/YFDCSfxsW4v/ZjEP+6z6cxeA2W1HX8bfpMfQE4OclJSV4AvBW4bTBhaQKZTxok80mDZD5pSatumaqqA0neDvwxvaGi11fVQwOLTBPFfNIgmU8aJPNJy2nVZ6qqPgt8tmUMY9kk2jD2ATKfWpvkx/5dBpRP42LcX/uRj38C8mnkX4NldBp/qqrL+5ckSRprnk5GkiSphU6LqSTnJflSkkeTbO0yluUkOTHJnUkeTvJQkiua9cck2Znkkeb/0V3HupgkhyT5YpLbm+WTktzVPP83NR0rx9Y45dMgJbk+yf4kD3YdiwZruZxO8sLmvfto817e2KzfmOTvk9zX/P3usGNv4lgu/tcnuTfJgSQXzdt2afO5+kiSS4cX9eRY7Htt3Mz/butCZ8XUnOn5fxw4BXhbklO6iqcPB4AtVXUKcCbwK028W4E7qupk4I5meVRdAeyas3wVcE1VvRJ4Eri8k6gGYAzzaZBuAM7rOggNVp85fTnwZPMevobee/qgL1fVac3fLw0l6Dn6jP9vgMuAj8277jHAe+lNjHkG8N5R/qE6xhb7Xhs387/bhq7Llqmxmp6/qvZV1b3N5afpvXDH04t5R7PbDuDCbiJcWpITgPOBjzTLAc4Gbm52GdnY+zRW+TRIVfU54Otdx6GB6yen537+3Ayc07y3R8Gy8VfV7qq6H/jHedf9l8DOqvp6VT0J7MQfDAO3xPfa2Jj/3daVLouphabnH4sXsWlKPx24C5iqqn3NpseBqY7CWs6HgF/jOx9aLwOeqqoDzfLYPP+LGNt8khbRT04/t0/zXv4Gvfc2wEnNoY//J8mPrHWwC2jznvT9PGTzvtfGyfzvtk7YAX2FkhwJfAp4Z1V9c+626g2NHLnhkUneBOyvqnu6jkXSUOwDXl5VpwPvAj6W5CUdx6QRtdT32igbpe+2LoupvqbnHyVJDqOXcDdW1S3N6ieSbGi2bwD2dxXfEs4C3pxkN72m9rOBa4Gjkhyca2zkn/9ljF0+ScvoJ6ef26d5L78U+Luq+lZV/R1A80XzZeCfrnnEi8TWWMl70vfzkCzyvTYuvuu7Lcl/6CKQLoupsZqev+mHcB2wq6o+OGfTbcDBkSaXArcOO7blVNW7q+qEqtpI73n+s6q6GLgTODiCZiRjX4GxyiepD/3k9NzPn4vovbcryXFNB3CS/PfAycBjQ4r7oDbvyT8GfizJ0U3H8x9r1mmAlvheGwuLfLf9qy5i6ayYao7vH5yefxfwyRGfnv8s4BJ6le/B4cZvBLYB5yZ5BHhDszwurgTeleRRev0srus4nlUbw3wamCQfB/4ceFWSPUnGdlSmvmOxnE7yG0ne3Ox2HfCy5j38Lr4zmvj1wP1J7qPXMf2XqmqogxT6iT/JDybZA/w08HtJHmqu+3Xg/fQKsi8AvzHs+CfEYt9rWiFnQJckSWrBDuiSJEktWExJkiS1YDElSZLUgsWUJElSCxZTkiRJLVhMSZIktWAxJUmS1ILFlCRJUgsWU5IkSS1YTEmSJLVgMSVJktSCxZQkSVILE1tMJbkhyW8m+ZEkX+pj//cl+Q/DiE3jbdRyK8nhSf4wyTeS/EGz7jeTfC3J42t1vxqcEcyp3UneMIz7ksbBoV0H0LWq+k/Aq7qOQ+vPsHIryUbgvwDPzNt0eVXdBFwETAEvq6oDSV4ObAFeUVX71zo+DY6fV9JomvhiShpnSea+h4+qqgML7PYK4D/P2fZy4O8spLSQJIcukkeSFjExh/mSnJ7k3iRPJ7kJeFGzfjrJnjn7XZlkb7Pfl5Kcs8BtHZbk40k+leQFC2zfmKSS/HySryZ5MskvJfnBJPcneSrJb83Z/3nN5HOuf2iSn05yz7zbf1eSWwf01KiljnLr8iR/A/zZMrH9OvBvgJ9JMpvkF4GdwPc1yze0evBaE6OQU0nenOSh5vNqJsk/W7tHLI23iSimmg+Q/wj8e+AY4A+A/2mB/V4FvB34wap6MfAvgd3z9jm8ua1vAW+pqm8vcdc/BJwM/AzwIeB/Bd4AfD/wliT/oo/wbwNOmvdBdgnw0T6uqzXWYW79C+CfNbezqKp6L/C/AzdV1ZFV9XvAjwP/tVm+rI+HqSEahZxK8k+BjwPvBI4DPgv84ULFmKQJKaaAM4HDgA9V1T9U1c3AFxbY71nghcApSQ6rqt1V9eU5218C/BHwZeDnq+rZZe73/VX136rqT+j1Z/l4Ve2vqr3AfwJOXy7wqvoWcBPwrwCSfD+wEbh9uetqKLrKrfdV1TNV9fdz1n2taUU4+GdLwngahZz6GeAzVbWzqv4B+ABwOPA/tnto0vo0KcXU9wF7q6rmrPvK/J2q6lF6v8TeB+xP8okk3zdnlzOBfw5sm3tbzeGSg38vn7P/E3Mu//0Cy0f2Gf8O4GeThF6r1CebIkvd6yq3vrpALMdW1VFz/na1eFzqzijk1PfNvc+q+sdm+/Grf1jS+jUpxdQ+4PimGDno5QvtWFUfq6ofptdpt4Cr5mz+E+D/AO5IMjXnOkfO+fubVcT3DPA9c5b/ybyY/gL4NvAjwM/Sa/7XaOgqt+Z+0Wp9GYWc+q/NbQLQxHIisHeVj0la1yalmPpz4ADwq01nzJ8Czpi/U5JXJTk7yQuB/0av9egf5+5TVf8n8DF6H1DHDii++4DXJ3l5kpcC715gn48CvwX8Q1V9fkD3q/ZGPbc0fkYhpz4JnJ/knCSH0ZtK41vA/7uqRyStcxNRTDWdLn8KuAz4Or3+ALcssOsLgW3A14DHge9lgcKmqt5Pr1PnnyY5ZgDx7aTXL+p+4B4W7g/174HXAE6ON0JGLLeemncI510rvL5GwCjkVFV9iV4/zQ83t/8TwE8s04Fdmlh5/mF5japmVM5+4Aeq6pGu45EkST0T0TK1Tvwy8AULKUmSRoszoI+BJLuBABd2HIokSZrHw3ySJEkteJhPkiSpBYspSZKkFobaZ+rYY4+tjRs3DvMuV+2ZZ57hiCOO6DqMvgw71nvuuedrVXXc0O5wEccee2wdd9xxY/M6DcM45e1Bo5RPcz+fxvG5HKZRfX5GJZ80WYZaTG3cuJG77757mHe5ajMzM0xPT3cdRl+GHWuS7zq1RRc2btzIBz7wgbF5nYZhnPL2oFHKp7mfT+P4XA7TqD4/o5JPmiwe5pMkSWqh72IqySFJvpjk9mb5pCR3JXk0yU1JXrB2YUqSJI2mlbRMXQHMPQv9VcA1VfVK4Eng8kEGJkmSNA76KqaSnACcD3ykWQ5wNnBzs8sOnFBSkiRNoH47oH8I+DXgxc3yy4CnqupAs7wHOH6hKybZDGwGmJqaYmZmZlWBPrD3G33ve+rxL13Vfcw1Ozu76liHbZxiHSUbt36mr/12bzt/jSPRqOo3R8A8kSbZssVUkjcB+6vqniTTK72DqtoObAfYtGlTrXb0x2Ur+VC7eHX3MdeojlRZyDjFKknSetNPy9RZwJuTvBF4EfAS4FrgqCSHNq1TJwB71y5MSZKk0bRsn6mqendVnVBVG4G3An9WVRcDdwIXNbtdCty6ZlFKkiSNqDbzTF0JvCvJo/T6UF03mJAkSZLGx4pmQK+qGWCmufwYcMbgQ5IkSRofzoAuSZLUgsWUJElSCxZTkiRJLVhMSZIktWAxJUmS1ILFlCRJUgsWU5IkSS1YTGmokrwoyV8m+askDyX59Wb9SUnuSvJokpuSvKDrWDUekhyV5OYkf51kV5LXJTkmyc4kjzT/j+46Tknrl8WUhu1bwNlV9VrgNOC8JGcCVwHXVNUrgSeByzuMUePlWuCPqurVwGuBXcBW4I6qOhm4o1mWpDVhMaWhqp7ZZvGw5q+As4Gbm/U7gAs7CE9jJslLgdfTnM6qqr5dVU8BF9DLIzCfJK0xiykNXZJDktwH7Ad2Al8GnqqqA80ue4Dju4pPY+Uk4G+Bf5fki0k+kuQIYKqq9jX7PA5MdRahpHVvRefmkwahqp4FTktyFPBp4NX9XjfJZmAzwNTUFLOzs8zMzKwqji2nHlh+J1j17XehzfMxpg4FfgB4R1XdleRa5h3Sq6pKUgtdeX4+zX3uZmdn2XLqs30HMmHP+yTmmrQoiyl1pqqeSnIn8DrgqCSHNq1TJwB7F7nOdmA7wKZNm+rII49kenp6Vfd/2dbP9LXf7otXd/tdmJmZWfXzMab2AHuq6q5m+WZ6xdQTSTZU1b4kG+i1gn6X+fk097mbmZnh6s8/03cg45QngzCBuSYtysN8GqokxzUtUiQ5HDiXXofhO4GLmt0uBW7tJkKNk6p6HPhqklc1q84BHgZuo5dHYD5JWmO2TGnYNgA7khxCr5j/ZFXdnuRh4BNJfhP4Ik2HYqkP7wBubKbTp4YTdwAADOxJREFUeAz4eZrcSnI58BXgLR3GJ2mds5jSUFXV/cDpC6x/DDij7e1v7PPQndaPqroP2LTApnOGHYukyeRhPkmSpBYspiRJklqwmJIkSWrBYkqSJKkFiylJkqQWLKYkSZJa6HxqBIeyS5KkcWbLlCRJUgsWU5IkSS0sW0wleVGSv0zyV0keSvLrzfqTktyV5NEkNzWncpAkSZoo/bRMfQs4u6peC5wGnJfkTOAq4JqqeiXwJHD52oUpSZI0mpYtpqpntlk8rPkr4Gzg5mb9DuDCNYlQkiRphPXVZyrJIUnuA/YDO4EvA09V1YFmlz3A8WsToiRJ0ujqa2qEqnoWOC3JUcCngVf3ewdJNgObAaamppiZmXne9i2nHljgWu3Mv4/VmJ2dHcjtDMM4xSpJ0nqzonmmquqpJHcCrwOOSnJo0zp1ArB3ketsB7YDbNq0qaanp5+3/bI1mGdq98XTy+6znJmZGebHOqrGKVZJktabfkbzHde0SJHkcOBcYBdwJ3BRs9ulwK1rFaQkSdKo6qdlagOwI8kh9IqvT1bV7UkeBj6R5DeBLwLXrWGckiRJI2nZYqqq7gdOX2D9Y8AZaxGUJEnSuHAGdEmSpBYspiRJklqwmJIkSWrBYkqSJKkFiylJkqQWLKYkSZJasJiSJElqwWJKkiSpBYspSZKkFiymJEmSWrCYkiRJasFiSpIkqQWLKUmSpBYspjRUSU5McmeSh5M8lOSKZv0xSXYmeaT5f3TXsWp8JDkkyReT3N4sn5TkriSPJrkpyQu6jlHS+mUxpWE7AGypqlOAM4FfSXIKsBW4o6pOBu5olqV+XQHsmrN8FXBNVb0SeBK4vJOoJE0EiykNVVXtq6p7m8tP0/sCPB64ANjR7LYDuLCbCDVukpwAnA98pFkOcDZwc7OL+SRpTVlMqTNJNgKnA3cBU1W1r9n0ODDVUVgaPx8Cfg34x2b5ZcBTVXWgWd5Dr2CXpDVxaNcBaDIlORL4FPDOqvpmrzGhp6oqSS1yvc3AZoCpqSlmZ2eZmZl5bvuWUw8sdLVW5t7+qJv/fKx3Sd4E7K+qe5JMr+L6z8unuc/d7OwsW059tu/bmqTnHSYv16SlWExp6JIcRq+QurGqbmlWP5FkQ1XtS7IB2L/QdatqO7AdYNOmTXXkkUcyPT393PbLtn5m4PHuvnh62X1GxczMzPOejwlwFvDmJG8EXgS8BLgWOCrJoU3r1AnA3oWuPD+f5j53MzMzXP35Z/oOZJzyZBAmMNekRXmYT0PV9Ge5DthVVR+cs+k24NLm8qXArcOOTeOnqt5dVSdU1UbgrcCfVdXFwJ3ARc1u5pOkNWUxpWE7C7gEODvJfc3fG4FtwLlJHgHe0CxLq3Ul8K4kj9LrQ3Vdx/FIWsc8zKehqqrPA1lk8znDjEXrS1XNADPN5ceAM7qMR9LksGVKkiSpBYspSZKkFiymJEmSWli2mPJcapIkSYvrp2XKc6lJkiQtYtliynOpSZIkLW5FfaY8l5okSdLz9T3P1KDOpTb/XE6jei61cTrv1DjFKknSetNXMTXIc6nNP5fTqJ5LbZzOOzVOsUqStN70M5rPc6lJkiQtop+WqYPnUnsgyX3NuvfQO3faJ5NcDnwFeMvahChJkjS6li2mxvFcahv7PHS4e9v5axyJJEla75wBXZIkqQWLKUmSpBYspiRJklroe54paVLZB0+StBRbpiRJklqwmJIkSWrBYkqSJKkFiylJkqQWLKYkSZJasJiSJElqwWJKkiSpBYspSZKkFiymJEmSWrCYkiRJasFiSpIkqQWLKUmSpBYspiRJklqwmJIkSWrBYkqSJKkFiykNVZLrk+xP8uCcdcck2Znkkeb/0V3GqPGR5MQkdyZ5OMlDSa5o1ptTkobGYkrDdgNw3rx1W4E7qupk4I5mWerHAWBLVZ0CnAn8SpJTMKckDZHFlIaqqj4HfH3e6guAHc3lHcCFQw1KY6uq9lXVvc3lp4FdwPGYU5KGyGJKo2CqqvY1lx8HproMRuMpyUbgdOAuzClJQ3Ro1wFIc1VVJanFtifZDGwGmJqaYnZ2lpmZmee2bzn1wJrHuJi5cXRl/vMxKZIcCXwKeGdVfTPJc9uWyqn5+TT3uZudnWXLqc/2HcOkPe+TmmvSQiymNAqeSLKhqvYl2QDsX2zHqtoObAfYtGlTHXnkkUxPTz+3/bKtn1nrWBe1++LpZfdZazMzM897PiZBksPoFVI3VtUtzeq+cmp+Ps197mZmZrj688/0HccovP7DNIm5Ji3Gw3waBbcBlzaXLwVu7TAWjZH0mqCuA3ZV1QfnbDKnJA3NssWUQ9k1SEk+Dvw58Koke5JcDmwDzk3yCPCGZlnqx1nAJcDZSe5r/t6IOSVpiPo5zHcD8FvAR+esOzjseFuSrc3ylYMPT+tNVb1tkU3nDDWQNbBxBYcYd287fw0jmRxV9Xkgi2we+5ySNB6WbZlyKLskSdLiVtsBve9hx0uNloHRHX01TiNVxilWSZLWm9aj+ZYbyr7UaBkY3dFX4zRSZZxilSRpvVntaL4nmuHGLDeUXZIkaT1bbTHlsGNJkiT6mxrBoeySJEmLWLbP1Hoeyi5JktSWM6BLkiS1YDElSZLUgsWUJElSCxZTkiRJLVhMSZIktWAxJUmS1ILFlCRJUgsWU5IkSS1YTEmSJLVgMSVJktSCxZQkSVILy56bT5I02jZu/Uxf++3edv4aRyJNJlumJEmSWrCYkiRJasFiSpIkqQWLKUmSpBYspiRJklpwNJ/UgX5HX63UllMPcNkyt+2Irm458k5afya6mFrqQ23+l5IfbJIkaSETXUxJ0qhaq9ZLSYNnnylJkqQWLKYkSZJa8DBfn+w0KkmSFmIxJUlac/4g1XrmYT5JkqQWWrVMJTkPuBY4BPhIVW0bSFRalbUY/TPMX4nm02gxn1ZmUkffDfpxr+T2bMXSqFh1y1SSQ4DfBn4cOAV4W5JTBhWYJov5pEEynyQNU5uWqTOAR6vqMYAknwAuAB4eRGCaOOaTBsl8WoD9lqS10abP1PHAV+cs72nWSathPmmQzCdJQ7Pmo/mSbAY2N4uzSb601vc5CL8KxwJfW+n1ctUaBLO8VcXaj0UezyvW4r76MT+ffvRHf/TvWKPHPo76yduOcnSp+x6ZfJr3+bRm76tR12eOdP78jFo+aXK1Kab2AifOWT6hWfc8VbUd2N7ifjqR5O6q2tR1HP0Yp1iXsKp8WiePfWB8Pp7T+vPJ53JpPj/Sd7Q5zPcF4OQkJyV5AfBW4LbBhKUJZD5pkMwnSUOz6papqjqQ5O3AH9Mbenx9VT00sMg0UcwnDZL5JGmYWvWZqqrPAp8dUCyjZpwOTY5TrItaZT6ti8c+QD4fjQF8PvlcLs3nR2qkqrqOQZIkaWx5OhlJkqQWJraYSnJikjuTPJzkoSRXLLDPdJJvJLmv+fs3XcQ6J57dSR5oYrl7ge1J8m+TPJrk/iQ/0EWcw5LkvCRfah7v1q7jGabF8jfJMUl2Jnmk+X9017GOo0nOreUkuT7J/iQPdh2LNCom9jBfkg3Ahqq6N8mLgXuAC6vq4Tn7TAP/uqre1FGYz5NkN7Cpqhac2yXJG4F3AG8Efgi4tqp+aHgRDk9zupD/DJxLb0LGLwBvm/v6rWeL5S9wGfD1qtrWFAFHV9WVHYY6diY9t5aT5PXALPDRqnpN1/FIo2BiW6aqal9V3dtcfhrYxfjPkHwBvQ+4qqq/AI5qvnTXo+dOF1JV3wYOni5kIiyRvxcAO5rddtArsLQyE51by6mqzwFf7zoOaZRMbDE1V5KNwOnAXQtsfl2Sv0ryfyf5/qEG9t0K+JMk9zQzN883SafQmKTHuqR5+TtVVfuaTY8DUx2FNc7MLUkrsuankxl1SY4EPgW8s6q+OW/zvcArqmq2OYT2H4GThx3jHD9cVXuTfC+wM8lfN78SNaHm52+S57ZVVSWZzOP4kjREE90yleQwel9EN1bVLfO3V9U3q2q2ufxZ4LAkxw45zLnx7G3+7wc+Te9wxFx9nUJjnZikx7qgRfL3iYOHdpv/+7uKb4xNfG5JWpmJLabS+wl/HbCrqj64yD7/pNmPJGfQe77+bnhRPi+WI5qOxiQ5AvgxYP5omtuAn2tG9Z0JfGPOIZ/1ZqJPF7JE/t4GXNpcvhS4ddixrQMTnVuSVm6SD/OdBVwCPJDkvmbde4CXA1TV7wIXAb+c5ADw98Bbq7vhj1PAp5va7lDgY1X1R0l+aU68n6U3ku9R4P8Dfr6jWNecpwtZNH+3AZ9McjnwFeAtHcU3tsytpSX5ODANHJtkD/Deqrqu26ikbk3s1AiSJEmDMLGH+SRJkgbBYkqSJKkFiylJkqQWLKYkSZJasJiSJElqwWJKkiSpBYspSZKkFiymJEmSWvj/AU/oATB23z2vAAAAAElFTkSuQmCC\n"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "execution_count": 22,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    }
  ],
  "metadata": {
    "kernel_info": {
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.7",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3"
    },
    "nteract": {
      "version": "0.14.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}