output_likelyhood.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9adfb973-5596-47f1-8642-76e5fe52cdc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import numpy  as np\n",
    "import pandas as pd\n",
    "import h5py\n",
    "\n",
    "import ehtplot\n",
    "from matplotlib import pyplot as plt, cm\n",
    "\n",
    "from common import dalt\n",
    "from common import hallmark as hm\n",
    "from common import viz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb2e6bf3-70dc-4abf-8e09-289286768b2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mag ['M' 'S']\n",
      "aspin [-0.94 -0.5   0.    0.5   0.94]\n",
      "Rhigh [  1  10  40 160]\n",
      "inc [ 10  30  50  70  90 110 130 150 170]\n"
     ]
    }
   ],
   "source": [
    "indices = ['mag', 'aspin', 'Rhigh', 'inc']\n",
    "\n",
    "def get(m):\n",
    "    def read(f, t, s, d):\n",
    "        return pd.read_csv(m+f'/stat_{f}_{t}.tsv', sep='\\t', usecols=indices+[s]).rename(columns={s:d}).set_index(indices)\n",
    "\n",
    "    return pd.concat([\n",
    "        read('86GHz',  'lc',    'med', 'Fnu_86GHz'   ),\n",
    "        read('86GHz',  'major', 'med', 'major_86GHz' ),\n",
    "        read('230GHz', 'mi3',   'med', 'mi3_230GHz'  ),\n",
    "        read('230GHz', 'major', 'med', 'major_230GHz'),\n",
    "        read('230GHz', 'minor', 'med', 'minor_230GHz'),\n",
    "        read('NIR',    'sed',   'med', 'Fnu_NIR'     ),\n",
    "        read('xray',   'sed',   'med', 'nuLnu_xray'  ),\n",
    "        # mring radius\n",
    "        # mring width\n",
    "        # mring asymmetry\n",
    "        # null locations\n",
    "    ], axis=1)\n",
    "\n",
    "df = get('cache/Illinois_thermal_w5')\n",
    "\n",
    "for k in indices:\n",
    "    globals()[k] = df.reset_index()[k].unique()\n",
    "    print(k, globals()[k][:16])\n",
    "    \n",
    "def edges(i):\n",
    "    m = (i[1:] + i[:-1]) / 2\n",
    "    l = [2*i[ 0] - m[ 0]]\n",
    "    r = [2*i[-1] - m[-1]]\n",
    "    return np.concatenate([l, m, r])\n",
    "\n",
    "Rhigh_edges = 10**edges(np.log10(Rhigh))\n",
    "inc_edges = edges(inc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "50d9320e-77b9-40d8-a9bf-5a872c535532",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "360\n",
      "32\n",
      "346\n",
      "18\n",
      "346\n",
      "317\n"
     ]
    }
   ],
   "source": [
    "lower = 34.380696128049614\n",
    "upper = 93.61106914382127\n",
    "\n",
    "print(len(df))\n",
    "print(np.sum(upper <= df.major_230GHz))\n",
    "print(np.sum(df.minor_230GHz <= upper))\n",
    "print(np.sum((df.minor_230GHz <= upper) & (upper <= df.major_230GHz)))\n",
    "print(np.sum((lower <= df.major_230GHz) & (df.minor_230GHz <= upper)))\n",
    "print(np.sum((lower <= df.minor_230GHz) & (df.major_230GHz <= upper)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "004ac95a-f544-402c-ba9f-864c223276a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Fnu_86GHz</th>\n",
       "      <th>major_86GHz</th>\n",
       "      <th>mi3_230GHz</th>\n",
       "      <th>major_230GHz</th>\n",
       "      <th>minor_230GHz</th>\n",
       "      <th>Fnu_NIR</th>\n",
       "      <th>nuLnu_xray</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mag</th>\n",
       "      <th>aspin</th>\n",
       "      <th>Rhigh</th>\n",
       "      <th>inc</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"6\" valign=\"top\">M</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">-0.50</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">160</th>\n",
       "      <th>30</th>\n",
       "      <td>2.278044</td>\n",
       "      <td>122.292752</td>\n",
       "      <td>0.098462</td>\n",
       "      <td>80.433761</td>\n",
       "      <td>71.208373</td>\n",
       "      <td>0.000143</td>\n",
       "      <td>5.555904e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>2.294351</td>\n",
       "      <td>124.047533</td>\n",
       "      <td>0.096340</td>\n",
       "      <td>80.520277</td>\n",
       "      <td>73.687309</td>\n",
       "      <td>0.000142</td>\n",
       "      <td>5.465083e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">0.00</th>\n",
       "      <th rowspan=\"4\" valign=\"top\">160</th>\n",
       "      <th>70</th>\n",
       "      <td>2.172176</td>\n",
       "      <td>140.504948</td>\n",
       "      <td>0.083838</td>\n",
       "      <td>79.035324</td>\n",
       "      <td>70.549255</td>\n",
       "      <td>0.000264</td>\n",
       "      <td>6.765862e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90</th>\n",
       "      <td>2.078028</td>\n",
       "      <td>143.003141</td>\n",
       "      <td>0.087216</td>\n",
       "      <td>76.296526</td>\n",
       "      <td>70.618663</td>\n",
       "      <td>0.000433</td>\n",
       "      <td>1.019600e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>2.104689</td>\n",
       "      <td>139.743133</td>\n",
       "      <td>0.092699</td>\n",
       "      <td>77.346543</td>\n",
       "      <td>71.366767</td>\n",
       "      <td>0.000252</td>\n",
       "      <td>7.046053e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>130</th>\n",
       "      <td>2.264433</td>\n",
       "      <td>134.096365</td>\n",
       "      <td>0.090839</td>\n",
       "      <td>80.733935</td>\n",
       "      <td>73.258668</td>\n",
       "      <td>0.000157</td>\n",
       "      <td>5.094390e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"16\" valign=\"top\">S</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">-0.50</th>\n",
       "      <th rowspan=\"5\" valign=\"top\">160</th>\n",
       "      <th>10</th>\n",
       "      <td>2.155446</td>\n",
       "      <td>98.970126</td>\n",
       "      <td>0.033885</td>\n",
       "      <td>74.656965</td>\n",
       "      <td>68.836112</td>\n",
       "      <td>0.000058</td>\n",
       "      <td>2.661791e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>2.148525</td>\n",
       "      <td>148.362848</td>\n",
       "      <td>0.050359</td>\n",
       "      <td>83.898539</td>\n",
       "      <td>60.813320</td>\n",
       "      <td>0.000079</td>\n",
       "      <td>4.205715e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>130</th>\n",
       "      <td>1.825985</td>\n",
       "      <td>151.620671</td>\n",
       "      <td>0.067492</td>\n",
       "      <td>79.273655</td>\n",
       "      <td>60.451933</td>\n",
       "      <td>0.000055</td>\n",
       "      <td>2.966531e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>1.793494</td>\n",
       "      <td>122.399787</td>\n",
       "      <td>0.060634</td>\n",
       "      <td>79.418426</td>\n",
       "      <td>66.419026</td>\n",
       "      <td>0.000025</td>\n",
       "      <td>1.330990e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>1.989254</td>\n",
       "      <td>96.365638</td>\n",
       "      <td>0.058050</td>\n",
       "      <td>79.132777</td>\n",
       "      <td>74.226127</td>\n",
       "      <td>0.000015</td>\n",
       "      <td>8.088204e+30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">0.00</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">10</th>\n",
       "      <th>130</th>\n",
       "      <td>1.708416</td>\n",
       "      <td>117.225247</td>\n",
       "      <td>0.088727</td>\n",
       "      <td>65.143758</td>\n",
       "      <td>55.550885</td>\n",
       "      <td>0.000032</td>\n",
       "      <td>6.555290e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>2.133763</td>\n",
       "      <td>130.702378</td>\n",
       "      <td>0.091887</td>\n",
       "      <td>78.968215</td>\n",
       "      <td>68.323343</td>\n",
       "      <td>0.000025</td>\n",
       "      <td>6.902032e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">160</th>\n",
       "      <th>130</th>\n",
       "      <td>1.914034</td>\n",
       "      <td>141.651613</td>\n",
       "      <td>0.064056</td>\n",
       "      <td>85.827636</td>\n",
       "      <td>56.365358</td>\n",
       "      <td>0.000025</td>\n",
       "      <td>1.779201e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>1.858771</td>\n",
       "      <td>104.861772</td>\n",
       "      <td>0.058441</td>\n",
       "      <td>76.669045</td>\n",
       "      <td>67.448547</td>\n",
       "      <td>0.000018</td>\n",
       "      <td>1.402985e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>1.819288</td>\n",
       "      <td>89.874705</td>\n",
       "      <td>0.050522</td>\n",
       "      <td>74.192265</td>\n",
       "      <td>70.727504</td>\n",
       "      <td>0.000021</td>\n",
       "      <td>1.425617e+31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">0.50</th>\n",
       "      <th rowspan=\"3\" valign=\"top\">10</th>\n",
       "      <th>30</th>\n",
       "      <td>2.038321</td>\n",
       "      <td>116.898198</td>\n",
       "      <td>0.072710</td>\n",
       "      <td>69.110696</td>\n",
       "      <td>60.472928</td>\n",
       "      <td>0.000036</td>\n",
       "      <td>3.044007e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>130</th>\n",
       "      <td>1.718483</td>\n",
       "      <td>104.790142</td>\n",
       "      <td>0.063397</td>\n",
       "      <td>56.477391</td>\n",
       "      <td>48.680879</td>\n",
       "      <td>0.000047</td>\n",
       "      <td>3.319026e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>2.174197</td>\n",
       "      <td>117.399384</td>\n",
       "      <td>0.070028</td>\n",
       "      <td>69.468492</td>\n",
       "      <td>61.231248</td>\n",
       "      <td>0.000045</td>\n",
       "      <td>3.524095e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">0.94</th>\n",
       "      <th rowspan=\"2\" valign=\"top\">40</th>\n",
       "      <th>10</th>\n",
       "      <td>1.848819</td>\n",
       "      <td>120.732230</td>\n",
       "      <td>0.075169</td>\n",
       "      <td>69.524496</td>\n",
       "      <td>63.202579</td>\n",
       "      <td>0.000130</td>\n",
       "      <td>3.412939e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>170</th>\n",
       "      <td>1.775968</td>\n",
       "      <td>120.583969</td>\n",
       "      <td>0.069612</td>\n",
       "      <td>70.607786</td>\n",
       "      <td>64.471256</td>\n",
       "      <td>0.000106</td>\n",
       "      <td>3.080134e+32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>160</th>\n",
       "      <th>10</th>\n",
       "      <td>2.293247</td>\n",
       "      <td>123.587652</td>\n",
       "      <td>0.073061</td>\n",
       "      <td>78.057396</td>\n",
       "      <td>67.928015</td>\n",
       "      <td>0.000023</td>\n",
       "      <td>4.523845e+31</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     Fnu_86GHz  major_86GHz  mi3_230GHz  major_230GHz  \\\n",
       "mag aspin Rhigh inc                                                     \n",
       "M   -0.50 160   30    2.278044   122.292752    0.098462     80.433761   \n",
       "                150   2.294351   124.047533    0.096340     80.520277   \n",
       "     0.00 160   70    2.172176   140.504948    0.083838     79.035324   \n",
       "                90    2.078028   143.003141    0.087216     76.296526   \n",
       "                110   2.104689   139.743133    0.092699     77.346543   \n",
       "                130   2.264433   134.096365    0.090839     80.733935   \n",
       "S   -0.50 160   10    2.155446    98.970126    0.033885     74.656965   \n",
       "                30    2.148525   148.362848    0.050359     83.898539   \n",
       "                130   1.825985   151.620671    0.067492     79.273655   \n",
       "                150   1.793494   122.399787    0.060634     79.418426   \n",
       "                170   1.989254    96.365638    0.058050     79.132777   \n",
       "     0.00 10    130   1.708416   117.225247    0.088727     65.143758   \n",
       "                150   2.133763   130.702378    0.091887     78.968215   \n",
       "          160   130   1.914034   141.651613    0.064056     85.827636   \n",
       "                150   1.858771   104.861772    0.058441     76.669045   \n",
       "                170   1.819288    89.874705    0.050522     74.192265   \n",
       "     0.50 10    30    2.038321   116.898198    0.072710     69.110696   \n",
       "                130   1.718483   104.790142    0.063397     56.477391   \n",
       "                150   2.174197   117.399384    0.070028     69.468492   \n",
       "     0.94 40    10    1.848819   120.732230    0.075169     69.524496   \n",
       "                170   1.775968   120.583969    0.069612     70.607786   \n",
       "          160   10    2.293247   123.587652    0.073061     78.057396   \n",
       "\n",
       "                     minor_230GHz   Fnu_NIR    nuLnu_xray  \n",
       "mag aspin Rhigh inc                                        \n",
       "M   -0.50 160   30      71.208373  0.000143  5.555904e+31  \n",
       "                150     73.687309  0.000142  5.465083e+31  \n",
       "     0.00 160   70      70.549255  0.000264  6.765862e+31  \n",
       "                90      70.618663  0.000433  1.019600e+32  \n",
       "                110     71.366767  0.000252  7.046053e+31  \n",
       "                130     73.258668  0.000157  5.094390e+31  \n",
       "S   -0.50 160   10      68.836112  0.000058  2.661791e+31  \n",
       "                30      60.813320  0.000079  4.205715e+31  \n",
       "                130     60.451933  0.000055  2.966531e+31  \n",
       "                150     66.419026  0.000025  1.330990e+31  \n",
       "                170     74.226127  0.000015  8.088204e+30  \n",
       "     0.00 10    130     55.550885  0.000032  6.555290e+32  \n",
       "                150     68.323343  0.000025  6.902032e+32  \n",
       "          160   130     56.365358  0.000025  1.779201e+31  \n",
       "                150     67.448547  0.000018  1.402985e+31  \n",
       "                170     70.727504  0.000021  1.425617e+31  \n",
       "     0.50 10    30      60.472928  0.000036  3.044007e+32  \n",
       "                130     48.680879  0.000047  3.319026e+32  \n",
       "                150     61.231248  0.000045  3.524095e+32  \n",
       "     0.94 40    10      63.202579  0.000130  3.412939e+32  \n",
       "                170     64.471256  0.000106  3.080134e+32  \n",
       "          160   10      67.928015  0.000023  4.523845e+31  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accept = pd.concat([\n",
    "    (1.7 <= df.Fnu_86GHz  ) & (df.Fnu_86GHz   <= 2.3),\n",
    "    (86  <= df.major_86GHz) & (df.major_86GHz <= 154),\n",
    "    (0   <= df.mi3_230GHz ) & (df.mi3_230GHz  <= 0.1),\n",
    "    # size_230GHz,\n",
    "                              (df.Fnu_NIR     <= 1.1e-3),\n",
    "                              (df.nuLnu_xray  <= 1e33),\n",
    "    # mring radius?\n",
    "    # mring width?\n",
    "    # mring asymmetry\n",
    "    # null locations\n",
    "], axis=1)\n",
    "\n",
    "df[accept.sum(axis=1) == 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d46934e0-52a5-47ff-ad53-127a77d4b7f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.31622776601683794"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rhigh_edges[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e72481a9-8978-4d65-b81d-5e3b5132ca4a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDAAAAGdCAYAAADkLoOBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5rElEQVR4nO3de7hkZ10n+u8v3UmTNEkjuZhwM4AJmJCBYCbMiKPxgiaRHtAxkiCMaKTHo5kZHPEIR4/X4RkdcfR4AOcJkgFFQOBRTDMBRBQZGeUSuSSBSczEOPSAJkHThNygO+/5Y1eFPk33rrV3V+21Vu3P53nq2buqa731e2vvb6+uX6/1rmqtBQAAAGDIjuq7AAAAAIBZNDAAAACAwdPAAAAAAAZPAwMAAAAYPA0MAAAAYPA0MAAAAIDB08AAAAAABk8DAwAAABg8DQx6U1UPr6rfr6q7q+pvquq5qzz3a6rqj6tqb1XdXFXfeZjnnVFV91XV6xdXOfRrjdl57yQTn5/cbtzIWmGR1pKF9TwflsU6smLfwdJbRy5mfh7xWWTxNDDo0yuTfCHJVyb53iS/UVVnH/ykqtqa5A+SvD3Jw5PsSvL6qjrzMGN+aGEVwzB0ys4BrmitPXRye8KGVAgbY61ZWOvzYVms53ffvoNl1zkXa/g84rPIgmlgjEBVnVJVV1fV31XV56pqd1Wd0HddR6Kqtif5F0n+79ba51trf5bk6iTPP8TTn5jkEUl+tbW2v7X2x0nef/Bzq+rSJHcmec8ia2c8ZAeO3FBztNYsyA6LtixZgXlaolzM/Dzis8jG0MAYhxOS/L9JHpPk9CQnJflXi3ihqnp7Vd15mNvb5/hSZybZ31q76YDHPpbkUF3POsxjTzqg7hOS/HySH5tjjYzfZs/O1H+oqjuq6v1VdcEca2FzGGqO1pqF9WQH1mJZsjJl38E8LEsuVv084rPIxtnadwHM1lq7OcnNk7v3V9W7k3xF1+2r6ilJ/mlr7Tc6vNYz11Xk2j00yd6DHtub5PhDPPd/JLktyY9X1a8m+aYk35jkTw54zi8keU1r7VNVh/r7hc1IdpIkP5HkE1k5RPLSJLur6imttf+5uBJZJgPO0VqzsNbnw5osUVYS+w7mZIlyMevziM8iG8QRGCNQVZdMut+3VdWdSV6S5KYZmz2otfbRLqE/UrWy4FM7zO3PDnr657PSkT3QCUnuOnjc1toXkzw7yXck+dusdDbfnGTP5HWfkuRbk/zqPOfD+G327CRJa+0DrbW7Wmv3t9Zel5XDHS+e/yxYVgPO0ZqysI7nw5osUVbsO5ibZcnFap9HfBbZWI7AGLiq+uYkv5TkOUk+Mnn41iQfrap/kuT/SXJfkj9N8sdZCVPLyjlal7fWPjY57O+ZWVl05keT7E/y+CTf21q7/qDXe0eSf3aYcv5ba+2iw9XaWrtgDVO7KcnWqjqjtfZXk8eenOSGw4z98ax0Oad1/vckr5vcvSArh6T9r0nH86FJtlTVWa21p66hJpaI7BxWy6EPg4QvM/AcrTULR5odOKwly8qh2HewZsuWi1U+j1wQn0U2TmvNbcC3rAT1T7NyONPJSX4rK4fzHZOVQ5UunjzvqKyE58+ysoP5miS7J392QZKXT76+Z/LYRUn+U89ze1OSNybZnuTpWTls6+zDPPcfJXlIkuOSvDjJXyfZNvmz45KcesDt5UnemuTkvn9+bv3dZKclycOSfPskO1uzssL23Ume0PfPx20ct6HnaC37kfU8382t622ZsmLf4Tav2zLlYvL8Q34eic8iG3pzCsnw/U6So7NyqNLbk/xVkk+01r6Q5BVJLq6q30ly4eT5H2krPpmVAB3so5Ovn8oazj9bkB9OcmxWzid7Y5L/o7X2YNezqt5RVf/X5O7zk3xm8txvSfKM1tr9SdJau6e19rfTW1YOCbuvtXb7Bs6F4ZGdlfn/+yS3J7kjyb9O8uzW2o0bXzIjNfQcrWU/MvP5cASWKSv2HczLMuUiOcznEZ9FNpZTSAautXZbkq876OFfmHz9XGvtiqo6Jsm1WdnBPKVWjl06Myt/WXzZkAd83+uhgK21v8/KuWSH+/OLDvj+x5P8eMdxf/ZIa2P8ZOdB/3jhBbG0hp6jNWZh5vNhvZYpK5MPXfYdHLFlysXkfqfPIz6LLJYGxrj9q6r6rqwc9vTayWN7k+zOStfyhT3VBUMnO3Dk5Ai6kRX4cnLBulRrbfazGIXpIjettRf3XAqMiuzAkZMj6EZW4MvJBV1ZAwMAAAAYPEdgAAAAAIM3+jUwHvrQh7YnPvGJfZcxd7fffntOPvnkvsuYq2WcU9L/vK699to7WmszC1jGrPT93i/KMs5rCHPqkpVlzEkyjPd/3pZxTkn/89rM+5Sk//d/EZZxTkn/8+qalapqxx133EaUtGH27duXrVtH/zHyyyzjvIYwp3vuuae11uZ25sfoj8DYvn17u/vuu/suY+7OO++8fPjDH+67jLlaxjkl/c+rqq5trZ0363nLmJW+3/tFWcZ5DWFOXbKyjDlJhvH+z9syzinpf16beZ+S9P/+L8Iyzinpf15ds7Jly5a2f//+jShpw/T93i/KMs5rCHOqqntaa9vnNd5oW0xVtTPJzm3btvVdCgyarMBscgLdyAp0M83KylVBgXkZ7SKerbXdrbVdW7Zs6bsUGDRZgdnkBLqRFehmmhUNDJiv0TYwpk466aS+S1iIXbt29V3C3C3jnJLxzGsZszKW936tlnFeY5nTMuYkGc/7vxbLOKdkPPOSlfFYxjkl45lX3+sPLMJY3vu1WsZ5LeOcRr8Gxnnnndf6Pq8H+tT1HExZYbPrkhU5YbOzT4FuNvt6MdDVvNfAGP0RGAAAAMDy08AAAAAABk8DAwAAABi80a8q8z/+9q78k1/8077LmOljv/uWvkvo7MnPuWQh43oP+nXdp/bmsT/6X+c65t++9x1zHS9JTr3gormPyWIs4uef9Ps7MJZ9yqKM6e/pMVnGfcoisrLZf/9++7U/uZBxn/+Cly1k3LHoO3/33PuFHHvuFb3W0JdF/U6PyWbP3yI4AgMAAAAYPA0MAAAAYPBGewpJVe1MsnPbiY/suxQYtGlWtu44re9SYLDsU6AbWYFupllJ+f9imKfRJqq1tru1tmvrQ+Z2SVlYStOsHLVNVuBw7FOgG1mBbqZZGfHHLRgkiQIAAAAGTwMDAAAAGDwNDAAAAGDwNDAAAACAwdPAAAAAAAZPAwMAAAAYPA0MAAAAYPA0MAAAAIDB08AAAAAABk8DAwAAABg8DQwAAABg8DQwAAAAgMEbdAOjqp5dVa+uqj+oqm/rux4AAACgHxvewKiqq6rqtqq6/qDHL6yqG6vq5qp6SZK01t7WWnthkhckec5G1woAAAAMQx9HYLw2yYUHPlBVW5K8MslFSc5KcllVnXXAU35q8ucAAADAJrR1o1+wtfa+qjr9oIfPT3Jza+2WJKmqNyV5VlV9MskvJnlHa+0vDzXeF+/em+tfsevB+6ecvzOnnL9zIbXDUNz2wd257YO7p3dP6rLN/nv35tNveNGD948/58Icf86Fh98AlsBd170zd133zundmVmxT2EzWs8+RVbYjNaTleSB3H/jmx+8t+XEs7P1pLPnXhsMyb47bsj+z94wvTvXnkO11uY5XrcXXWlgvL219qTJ/e9OcmFr7Qcn95+f5GlJbkryfUk+lOSjrbX/fPBYjz/7ye2X3vDOgx8enGeefdpCxn37DZ9ZyLhj4X1NLnnKI65trZ0363nbvvKM9ojn/toGVHRkLv6Gx/VdQmfXvO+WuY+52ee/SLf+2jNnZuWhj3pCe9IVV25USUfkxRedOfcxX/6Om+Y+ZjKuWje7D7z0gk77lLFk5WO/+5aFjPvk51yykHHHYlHv6yIs6mfVNSvbjj2u/c5f/M+5vvYi/v5bxN/TLMbY9n8feOkF97TWts9rvA0/AuMw6hCPtdbaryf59Y0uBgAAABiWoTQw9iR59AH3H5Xk06ttUFU7k+w89dGnL7AsGL9pVrbuWMzRKrAMpjnZduIj+y4FBk1WoJsH//119DF9lwJLZSiXUf1QkjOq6rFVdUySS5NcvdoGrbXdrbVdxx1/woYUCGM1zcpR2+Z25BYsnWlOtj5ETmA1sgLdPPjvry1b+i4Flkofl1F9Y5I/T/KEqtpTVZe31vYluSLJu5J8MsmbW2s3rDYOAAAAsHn0cRWSyw7z+DVJruk6jlNIoBunkMBsDouHbmQFunEKCSzGUE4hWTOnkEA3TiGB2RwWD93ICnTjFBJYjNE2MAAAAIDNYyhXIVkzp5BAN04hgdkcFg/dyAp04xQSWIzRHoHhFBLoxikkMJvD4qEbWYFunEICizHaBgYAAACweWhgAAAAAINnDQxYctbAgNmc1w/dyAp0Yw0MWIzRHoFhDQzoxhoYMJvz+qEbWYFurIEBizHaBgYAAACweWhgAAAAAINnDQxYctbAgNmc1w/dyAp0Yw0MWIzRHoFhDQzoxhoYMJvz+qEbWYFurIEBizHaBgYAAACweWhgAAAAAIOngQEAAAAMngYGAAAAMHiuQgJLzlVIYDZXVoBuZAW6cRUSWIzRHoHhKiTQjauQwGyurADdyAp04yoksBijbWAAAAAAm4cGBgAAADB4GhgAAADA4GlgAAAAAIOngQEAAAAMnsuowpJzGVWYzaUhoRtZgW5cRhUWY7RHYLiMKnTjMqowm0tDQjeyAt24jCosxmgbGAAAAMDmoYEBAAAADJ4GBgAAADB4o13Ec+phDzk6zzx78y5OuIi5v/2Gz8x9zGQxtS7KmGrt6nGnPDRX/cjT+y5jpnNPf1jfJXT2vHPm/3uyqPl/5NY75z7mIua/SF/3a/287osvOrOfF16H9774Gxcy7iL2K4uodVH7v0V4+Ttu6ruEuVtIVi76yfmPyaje176zsvWomvu/K8f07//NbhE/q7F9Tjn2pfMdzxEYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4o21gVNXOqrryzr17+y4FBm2albs+JytwONOc7Lvv7r5LgUGTFejmwazs3993KbBURtvAaK3tbq3tetiOHX2XAoM2zcrxJ8gKHM40J1sfsr3vUmDQZAW6eTArW7b0XQoMRlU9dbVblzFGfxlVAAAAYPA+nOSGJLdP7tcBf9aSfPOsATQwAAAAgEX7sST/Ism9Sd6U5Pdba59fywCjPYUEAAAAGIfW2q+21r4+yRVJHp3kPVX15qp6StcxNDAAAACADdFa++skf5DkD5Ocn+TMrts6hQQAAABYqKp6XJJLkzwryaeychrJy1pr93UdQwMDAAAAWLSbk3w8K0dffC7JY5L8cNXKWp6ttf80awANDAAAAGDRfj4rVxtJkoeuZ4B1NTCq6qeSPGny4q219tz1jAMAAAAsv9bazx7pGOtdxPOU1tqlrbXLNC8AAACArqrqeQd+7WrNDYyq+o4kO6rqsqq6uKouXusYAAAAwKb17w762sl6TiE5Kcl7khw9+R4AAABgrWotT15zA6O19rqqekdWGhhJcldVfa619ra1jgUAAADQxXrXwPhAkm9P8m1JPpbkW+ZW0QGq6nFV9ZqqeusixgcAAADGYb0NjCckOS3JI5KcmeSzXTesqquq6raquv6gxy+sqhur6uaqekmStNZuaa1dvs4aAQAAgCWxrsuoJvm5JC864Pvb1rDta5O8IslvTR+oqi1JXpnkGUn2JPlQVV3dWvvEOusDAAAAhummydcb17LRmhsYVfX0JDuS/Mnkoce31jq/aGvtfVV1+kEPn5/k5tbaLZPXeFOSZyWZ2cC4/Y7b8/Snnffg/R/4wV25/IW7upYDo/SaV1+Zq37zyundTovp/sPf35Ef+M5vevD+P3/O9+XZl75g/sXBgLztTa/N1b/7uundmVn54t17c/0rvrQPOeX8nTnl/J0Lqg6G4bYP7s5tH9w9vdtpnyIrbEbrysq+fT6rsOkc9FnlkD2H1tqlB37taj1HYDx8cpunRyb51AH39yR5WlWdmORlSc6tqpe21v7DwRuefNLJef8HPjzncmDYLn/hl3Z+xx5dd3TZ5iseflKu+v0/mf1EWCLPvvQFDzbqvu6Mr5iZlaO378iTrrhy1tNgqRzYfPjASy/otE+RFTajdWVl61afVdh0Dvqssu/gP6+qpx5is71J/qa19mXPP9B6rkKye/az1uxQl05prbXPJvmhBbweAAAAsPFeleSpST6elV7Akybfn1hVP9Ra+8PDbbiuRTyr6ier6o2T2xvWM8ZB9iR59AH3H5Xk0zNq2FlVV965d+8cXh6W1zQrd31OVuBwpjnZd9/dfZcCgyYr0M2DWdm/v+9SYIhuTXJua+281trXJjk3yfVJvjXJf1xtw/VeheQrW2uXTW7PXecYB/pQkjOq6rFVdUySS5NcvdoGrbXdrbVdD9uxYw4vD8trmpXjT5AVOJxpTrY+ZHvfpcCgyQp082BWtmzpuxQYoie21m6Y3plcvOPc6ZqYq1nPIp7fkWRHVV2WlfNU0lq7Zg3bvzHJBUlOqqo9SX6mtfaaqroiybuSbEly1YETAgAAAJbCjVX1G0neNLn/nCQ3VdW2JF9cbcP1LOJ5UpI/TnL05PtDrV9xWK21yw7z+DVJ1tII2Zlk5+Me/9VreXnYdKZZeeRjHtt3KTBY05xsO/GRfZcCgyYr0M00K8ds29Z3KTBEL0jyw0lelJV+wp8leXFWmhffdNit0uEUkqp6flXdXlV7qur7Wmuvy8q1Ws9M8m8n9zecU0igG6eQwGwOi4duZAW6cQoJHF5r7d7W2q+01r6ztfbs1trLW2v3tNYeaK19frVtu6yB8dNJLk7ylCSPrap3J3lLVo7AeNER1g4AAABsElV1RlW9tao+UVW3TG9dtu1yCsnnW2sfmrzQzyX5uyRnttbuXH/JR84pJNCNU0hgNofFQzeyAt04hQRW9V+S/EySX83KKSPfn45LU3Q5AuPUqtpVVd+Y5CuT7Om7eZE4hQS6cgoJzOaweOhGVqAbp5DAqo5trb0nSbXW/qa19rNJvrnLhl2OwPiZJP8oyfcmOSfJ8VX1R0k+kuQjrbU3rK9mAAAAYJO5r6qOSvJXk6uR/u8kp3TZcGYDo7V25YH3q+pRWWlonJPkoiQaGAAAAEAXL0pyXJJ/k+QXsnIayb/ssuGaL6PaWtuTZE/WcMnTRbAGBnRjDQyYzXn90I2sQDfWwIBVtSS/neSrsnJxkCR5dVYOlFhVlzUwBskaGNCNNTBgNuf1QzeyAt1YAwNW9TtZWcjzXyR55uS2s8uGaz4CAwAAAGCdbm+tXb2eDau1Nu9iNtTXfu157f0f+HDfZUBvjj26rm2tnTfreWPJykduvXPuY557+sPmPibj0yUrY8kJLMqy7VNYjLff8Jm+S+jdJU95RKesbN++vX12790bURIM0rFH1z2ttf/fYXtV9S1JLkvyniT3Tx9vrf3erPEcgQEAAABslO9P8sSsrH/xwOSxlmR5GxgW8YRuZAVmkxPoRlagG4t4wqqe3Fo7Zz0bWsQTlpyswGxyAt3ICnRjEU9Y1V9U1Vnr2XC0R2AAAAAAo/P1Sb6vqv46K2tgVJLWWpt5GVUNDAAAAGCjXLjeDTUwAAAAgA3RWvub9W472jUwqmpnVV155969fZcCgyYrMJucQDeyAt1Ms7Jv//6+S4GlMtoGhkWkoBtZgdnkBLqRFejGIp6wGKNtYAAAAACbhwYGAAAAMHgaGAAAAMDgaWAAAAAAg6eBAQAAAAzeaBsYLuMF3cgKzCYn0I2sQDcuowqLMdoGhst4QTeyArPJCXQjK9CNy6jCYoy2gQEAAABsHhoYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIOngQEAAAAMngYGAAAAMHijbWBU1c6quvLOvXv7LgUGTVZgNjmBbmQFuplmZd/+/X2XAktltA2M1tru1tquh+3Y0XcpMGiyArPJCXQjK9DNNCtbt2zpuxRYKqNtYAAAAACbhwYGAAAAMHgaGAAAAMDgaWAAAAAAg6eBAQAAAAyeBgYAAAAweBoYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIO3te8CVlNV25O8KskXkry3tfY7PZcEAAAA9GDDj8Coqquq6raquv6gxy+sqhur6uaqesnk4e9K8tbW2guT/PONrhUAAAAYhj5OIXltkgsPfKCqtiR5ZZKLkpyV5LKqOivJo5J8avK0/RtYIwAAADAgG34KSWvtfVV1+kEPn5/k5tbaLUlSVW9K8qwke7LSxPhoDtNsuf2O2/P0p5334P0f+MFdufyFu+ZfOAzIa159Za76zSund0/qso2ssBmtNStywmZknwLdvPutr88f/d7rp3c7ZeWL+/bJCpvOQfuVufYcqrU2z/G6vehKA+PtrbUnTe5/d5ILW2s/OLn//CRPS/ITSV6R5L4kf3aoNTC+9mvPa+//wIc3qnQYnGOPrmtba+fNet5YsvKRW++c+5jnnv6wuY/J+HTJylhyAouybPsUFuPtN3ym7xJ6d8lTHtEpK9u3b2+f3Xv3RpQEg3Ts0XVPa237vMYbyiKedYjHWmvt7iTff8gNqnYm2fm4x3/1QguDsZMVmE1OoBtZgW6mWTlm27a+S4GlMpTLqO5J8ugD7j8qyadX26C1tru1tuthO3YstDAYO1mB2eQEupEV6Gaala1btvRdCiyVoTQwPpTkjKp6bFUdk+TSJFf3XBMAAAAwEH1cRvWNSf48yROqak9VXd5a25fkiiTvSvLJJG9urd0wY5ydVXXlnXv3Lr5oGDFZgdnkBLqRFehmmpV9+11IEeZpwxsYrbXLWmuntdaObq09qrX2msnj17TWzmytPb619rIO4ziEETqQFZhNTqAbWYFunEICizGUU0gAAAAADksDAwAAABi80TYwnIMJ3cgKzCYn0I2sQDfWwIDFGG0DwzmY0I2swGxyAt3ICnRjDQxYjNE2MAAAAIDNQwMDAAAAGLzRNjCcgwndyArMJifQjaxAN9bAgMUYbQPDOZjQjazAbHIC3cgKdGMNDFiM0TYwAAAAgM1DAwMAAAAYvNE2MJyDCd3ICswmJ9CNrEA31sCAxRhtA8M5mNCNrMBscgLdyAp0Yw0MWIzRNjAAAACAzUMDAwAAABg8DQwAAABg8DQwAAAAgMEbbQPDKtjQjazAbHIC3cgKdOMqJLAYo21gWAUbupEVmE1OoBtZgW5chQQWY7QNjKnb77i97xIW4jWvvrLvEuZuGeeUjGdey5iVt73ptX2XsBBj+Z1ai7HMaRlzkozn/V+LZZxTMp55ycp4LOOckuTdb3193yV08sV9+/ouYe6W9XdqGee1jHMafQPjs3fc0XcJC3HVby7fL9syzikZz7yWMStX/+7r+i5hIcbyO7UWY5nTMuYkGc/7vxbLOKdkPPOSlfFYxjklyR/93jgaGPuXsIGxrL9TyzivZZzT6BsYAAAAwPKr1lrfNRyRqroryY1917EAJyVZtv/eWMY5Jf3P66taayfPetKSZqXv935RlnFeQ5jTzKwsaU6SYbz/87aMc0r6n9dm3qck/b//i7CMc0r6n1fXrDyQ5N4NqGcjbU2yfIeWLOe8hjCnY1trcztwYvQNDAAAAGD5OYUEAAAAGDwNDAAAAGDwNDAAAACAwdPAAAAAAAZPAwMAAAAYPA0MAAAAYPA0MAAAAIDB29p3AUfqoQ99aHviE5/Ydxlzd/vtt+fkk0/uu4y5WsY5Jf3P69prr72jtTazgGXMSt/v/aIs47yGMKcuWVnGnCTDeP/nbRnnlPQ/r828T0n6f/8XYRnnlPQ/r65Zqap23HHHbURJG2bfvn3ZunX0HyO/zDLOawhzuueee1prbW4HTlRrbV5j9WL79u3t7rvv7ruMuTvvvPPy4Q9/uO8y5moZ55T0P6+qura1dt6s5y1jVvp+7xdlGec1hDl1ycoy5iQZxvs/b8s4p6T/eW3mfUrS//u/CMs4p6T/eXXNypYtW9r+/fs3oqQN0/d7vyjLOK8hzKmq7mmtbZ/XeKNtMVXVziQ7t23b1ncpMGiyArPJCXQjK9DNNCtV1XcpsFRGuwZGa213a23Xli1b+i4FBk1WYDY5gW5kBbqZZkUDA+ZrtA2MqZNOOqnvEhZi165dfZcwd8s4p2Q881rGrIzlvV+rZZzXWOa0jDlJxvP+r8UyzikZz7xkZTyWcU7JeObV9/oDizCW936tlnFeyzin0a+Bcd5557W+z+uBPnU9B1NW2Oy6ZEVO2OzsU6Cbzb5eDHQ17zUwRn8EBgAAALD8NDAAAACAwdPAAAAAAAZv9KvK/OUn/leOPfeKuY75Az/9I3Mdb5Gu+vlX9l3CUlrE70DfP6tFZGURTr3gor5LWEp/+9539F3CKPzDvV/MWz/2mbmO+fJ33DTX8Vicj/3uW/ouYTSu+9TePPZH/2vfZcx08Tc8ru8SWINr3ndL3yXM3T33fmHu//5axL+VFpGVX9n5NXMfc2x+bPcn5z7mMuZkLRyBAQAAAAyeBgYAAAAweKM9haSqdibZWcec0HcpMGiyArNNc3Lqo0/vuxQYtGlWtu44re9SYNCmWUn5/2KYp9EmqrW2u7W2K1u29V0KDJqswGzTnBx3vEYfrGaalaO2be+7FBi0B//9Nd6PWzBIEgUAAAAMngYGAAAAMHgaGAAAAMDgaWAAAAAAg6eBAQAAAAyeBgYAAAAweBoYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIM36AZGVT27ql5dVX9QVd/Wdz0AAABAPza8gVFVV1XVbVV1/UGPX1hVN1bVzVX1kiRprb2ttfbCJC9I8pyNrhUAAAAYhj6OwHhtkgsPfKCqtiR5ZZKLkpyV5LKqOuuAp/zU5M8BAACATWjrRr9ga+19VXX6QQ+fn+Tm1totSVJVb0ryrKr6ZJJfTPKO1tpfHnK8fffm/hvf/OD9LSeena0nnb2Q2mEo9t1xQ/Z/9obp3ZO6bCMrbEZrzcrn/uGz+YnnfqnH/q3f9bw847uft6DqYBjWs0/Zf+/efPoNL3rw/vHnXJjjz7nw8BvAErjrunfmruveOb3bKSvJA/79xaZzUFbm2nPY8AbGYTwyyacOuL8nydOS/Osk35pkR1V9dWvtPx+8YW09Ntue8D0bU+UR+JWdX7OgcV+xkHE3ux/b/cm+S1jV1pO+tPO776OvvKPLNic/5qtyyS+/ZaF1MVzP+5F/P/cxX3/dZ+Y+5iK96rvOmpmVE77ixPzSG94562lr8uKLzpzreFPPPPu0uY/59hvG9TOdt2e+eDH71DG9r5c85RGd9ikPP+WUTb1Ped4588/fuac/bO5jJslHbr1z7mMu6u//i7/hcQsZd26+4YeT/HCSbvuUFUfN/bPKIt6nRX1W2ezG9L5e9fPzPeHhqEd8+8o3t928b57jDqWBUYd4rLXWfj3Jr290MQAAAMCwDKWBsSfJow+4/6gkn15tg6ramWRnHXPCIuuC0Ztm5YRTH9N3KTBY05yc+ujT+y4FBs0+BbqZZiU16Is+wugMJVEfSnJGVT22qo5JcmmSq1fboLW2u7W2K1u2bUiBMFbTrGzbfnzfpcBgTXNy3PGa4rAa+xTo5sHPKoP5uAXLoY/LqL4xyZ8neUJV7amqy1tr+5JckeRdST6Z5M2ttRtWGwcAAADYPPq4Csllh3n8miTXdB3HKSTQjcN9YTankEA39inQjVNIYDFGmyinkEA3DveF2ZxCAt3Yp0A3TiGBxZAoAAAAYPCGchWSNXMKCXTjcF+YzSkk0I19CnTjFBJYjNEmyikk0I3DfWE2p5BAN/Yp0I1TSGAxJAoAAAAYPA0MAAAAYPCsgQFLzvnKMJs1MKAb+xToxhoYsBijTZQ1MKAb5yvDbNbAgG7sU6Aba2DAYkgUAAAAMHgaGAAAAMDgWQMDlpzzlWE2a2BAN/Yp0I01MGAxRpsoa2BAN85XhtmsgQHd2KdAN9bAgMWQKAAAAGDwNDAAAACAwdPAAAAAAAZPAwMAAAAYPFchgSVnxXiYzVVIoBv7FOjGVUhgMUabKFchgW6sGA+zuQoJdGOfAt24CgkshkQBAAAAg6eBAQAAAAyeBgYAAAAweBoYAAAAwOBpYAAAAACD5zKqsORc8g5mcxlV6MY+BbpxGVVYjNEmymVUoRuXvIPZXEYVurFPgW5cRhUWQ6IAAACAwdPAAAAAAAZPAwMAAAAYvNEu4jl17tc8Ju//wCvmOuaP7f7kXMdLkq/4x1fMfcwk+YGf/pGFjDsWV/38Kxcy7iLe13/40Hx/T6eOPbrbe/DoHQ/Jr+z8mrm+9kduvXOu4yXJ66/7zNzHfN45p819zLFZxPu6KPP+PZ161UJGne3l77ipp1deuzHVugjPf8HL+i6hsyc/55JeX3/vXffnmvfdMtcx//a975jreEly6gUXzX3MJHOf+9gs4me1KIv6HeiqtmzpvYYuxvJvurEZ079BF/V7eutH5/t5zREYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4o21gVNXOqrryzr17+y4FBk1WYLZpTu6563N9lwKDNs3KA/ff3XcpMGjTrLT2QN+lwFIZbQOjtba7tbbrYTt29F0KDJqswGzTnBx3/Al9lwKDNs3KUdu2910KDNo0K1Wj/bgFgyRRAAAAwIaqqqOr6tyqOqXrNhoYAAAAwEJV1X+uqrMn3+9I8rEkv5XkI1V1WZcxNDAAAACARftnrbUbJt9/f5KbWmvnJPnaJP9nlwE0MAAAAIBF+8IB3z8jyduSpLX2t10H0MAAAAAAFu3OqnpmVZ2b5OlJ3pkkVbU1ybFdBti6wOIAAAAAkuRfJfn1JKcmedEBR158S5L/2mUADQwAAABgoVprNyW58BCPvyvJu7qMsa5TSKrqp6rqTVX1xqp6w3rGAAAAADafqnregV+7Wu8RGKe01i5d57YAAADA5vXvkrz+gK+drLmBUVXfkWTH5Dqte5OktXbNWscBAAAANrVay5PXcwTGSUnek+ToyfcAAAAAC7XmBkZr7XVV9Y6sNDCS5K6q+lxr7W1zrQwAAABgYl2LeCb5QJJvT/JtST6WlcuezF1VPa6qXlNVb13E+AAAAMA4rLeB8YQkpyV5RJIzk3y264ZVdVVV3VZV1x/0+IVVdWNV3VxVL0mS1totrbXL11kjAAAAMDw3Tb7euJaN1nsVkp9L8qIDvr9tDdu+NskrkvzW9IGq2pLklUmekWRPkg9V1dWttU+ssz4AAABggKZXNV3r1U3XcxWSpyfZkeRPJg89vrXWuWvSWntfVZ1+0MPnJ7m5tXbL5DXelORZSWY2MG6/4/Y8/WnnPXj/B35wVy5/4a6u5cAovebVV+aq37xyerfTYrqywma01qx87h8+m5947oUP3v/W73penvHda7o8OYzObR/cnds+uHt6t9M+Zf+9e/PpN7zowfvHn3Nhjj/nwsNvAEvgruvembuue+f0bqestAf2ywqbzkFZ+bKeQ1U99RCb7U3yN621fauNvZ4jMB4+uc3TI5N86oD7e5I8rapOTPKyJOdW1Utba//h4A1PPunkvP8DH55zOTBsl7/wS82HY4+uO7psIytsRmvNyglfcWJ+6Q3vnPU0WCqnnL8zp5y/M0nygZde0GmfsuXYHXnEc39tkWXB4BzYfLj1157ZKSt11BZZYdM5KCuHaki8KslTk3w8K5dRfdLk+xOr6odaa394uLHXcxWS3bOftWaHuvZra619NskPLeD1AAAAgI13a5LLW2s3JElVnZXkx5P8QpLfSzK/BsbkBX4yK12SZKXR8Nz1jHOAPUkefcD9RyX59IwadibZ+bjHf/URvjQsN1mB2aY5OfXRp/ddCgzaNCtbd5zWdykwaNOsZMvRfZcCQ/TEafMiSVprn6iqc1trt1Qd6tiGL1nvVUi+srV22eR2pM2LJPlQkjOq6rFVdUySS5NcvdoGrbXdrbVdD9uxYw4vD8tLVmC2aU6OO/6EvkuBQZtm5aht2/suBQZtmpWq9X7cgqV2Y1X9RlV94+T2qiQ3VdW2JF9cbcM1J6qqviPJjqq6rKourqqL17j9G5P8eZInVNWeqrp8slDHFUneleSTSd58YEcGAAAAWAovSHJzVq5s+qNJbpk89sUk37Tahus5heSkJH+c5OjJ96sf43GQ1tplh3n8miTXdB3HYfHQjazAbE4hgW6cQgLdOIUEDq+1dm+SX5ncDvb51badeQRGVT2/qm6fHC3xfa211yW5McmZSf7t5P6Gc1g8dCMrMJtTSKAbp5BAN04hgcOrqjOq6q1V9YmqumV667Jtl0T9dJKLkzwlyWOr6t1J3pKVIzBetM6aAQAAgM3nvyT5jST7snLKyG8l+e0uG3Y5heTzrbUPJUlV/VySv0tyZmvtznWVOicOi4duZAVmcwoJdOMUEujGKSSwqmNba++pqmqt/U2Sn62q/5bkZ2Zt2OUIjFOraldVfWOSr0yyp+/mReKweOhKVmA2p5BAN04hgW6cQgKruq9WwvFXVXVFVX1nklO6bNjlCIyfSfKPknxvknOSHF9Vf5TkI0k+0lp7wzqLBgAAADaXFyU5Lsm/SfILWTmN5F922XBmA6O1duWB96vqUVlpaJyT5KIkGhgAAABAFy0ra158VVbW1kySV2elz7CqNV9GtbW2J8merOGSp4vgvH7oRlZgNmtgQDfWwIBurIEBq/qdJD+e5LokD6xlw9GelOW8fuhGVmA2a2BAN9bAgG6sgQGrur21dnVr7a9ba38zvXXZcM1HYAAAAACs089U1W8meU+S+6cPttZ+b9aGo29gfGrvffmx3Z/su4zeXPXzr5z7mD/w0z8y9zGTxdQ6Jn3/no4lK9e875ZRjEly8Tc8biHj9vl7+nd778vL33FTb6+/FmOpc0ye/JxL+i6hsxdfdOZCxr3kpQsZtpNTL7iovxdnTfys+jWWf9f8ys6vWci4H7n1zrmPee7pD5v7mEn///4fsO9P8sSsrH8xPYWkJVn+BgYAAAAwGk9urZ2zng1He1JWVe2sqivvv/uuvkuBQZMVmG2ak3333d13KTBo06w8cL+swGqmWWltTesTwmbxF1V11no2HG0DY7owzrbtx/ddCgyarMBs05xsfYiFCWE1FvGEbiziCav6+iQfraobq+rjVXVdVX28y4ZOIQEAAAA2yoXr3VADAwAAANgQXS+ZeiiOaQIAAAAGb7RHYFTVziQ7Tzj1MX2XAoMmKzDbNCfbTnxk36XAoE2zsnXHaX2XAoM2zUq2HN13KbBURnsEhoUJoRtZgdks4gndWMQTurGIJyyGRAEAAACDp4EBAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIPnMqqw5GQFZnMZVejGZVShG5dRhcUY7REYLg0J3cgKzOYyqtCNy6hCNy6jCoshUQAAAMDgaWAAAAAAg6eBAQAAAAyeBgYAAAAweBoYAAAAwOBpYAAAAACDp4EBAAAADJ4GBgAAADB4W/suYL2qameSnSec+pi+S4FBkxWYbZqTbSc+su9SYNCmWdm647S+S4FBm2YlW47uuxRYKqM9AqO1tru1tmvb9uP7LgUGTVZgtmlOtj5ke9+lwKBNs3LUNlmB1UyzUjXaj1swSBIFAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIOngQEAAAAMngYGAAAAMHgaGAAAAMDgaWAAAAAAg6eBAQAAAAyeBgYAAAAweBoYAAAAwOBpYAAAAACDp4EBAAAADN7WvgtYTVVtT/KqJF9I8t7W2u/0XBIAAADQgw0/AqOqrqqq26rq+oMev7Cqbqyqm6vqJZOHvyvJW1trL0zyzze6VgAAAGAY+jiF5LVJLjzwgarakuSVSS5KclaSy6rqrCSPSvKpydP2b2CNAAAAwIBs+CkkrbX3VdXpBz18fpKbW2u3JElVvSnJs5LsyUoT46M5TLPl3s/9fd7y45c8eP+sZ1ySs7/te+ZfOAzIDX/45nzi3W+Z3j2pyzaywma01qx88e69uf4Vux68f8r5O3PK+TsXVB0Mw7vf+vr80e+9fnq30z5l/7178+k3vOjB+8efc2GOP+fCw28AS+Cu696Zu6575/Rup6y0B/bLCpvOQVmZa89hKGtgPDJfOtIiWWlcPC3Jryd5RVV9R5Ldh9rw2BMenkt++S2H+iNYWmd/2/c82Hx41XeddUeXbWSFzWitWTl6+4486YorF14XDMkzvvt5ecZ3Py9JcslTHtFpn7Ll2B15xHN/bZFlweAc2Hy49dee2SkrddQWWWHTOSgr++Y59lAaGHWIx1pr7e4k33/IDap2Jtl5wqmPWWhhMHayArNNc7LtxEf2XQoM2jQrW3ec1ncpMGjTrGTL0X2XAktlKJdR3ZPk0Qfcf1SST6+2QWttd2tt17btxy+0MBg7WYHZpjnZ+pDtfZcCgzbNylHbZAVWM81K1VA+bsFyGEqiPpTkjKp6bFUdk+TSJFf3XBMAAAAwEH1cRvWNSf48yROqak9VXd5a25fkiiTvSvLJJG9urd0wY5ydVXXl/XfftfiiYcRkBWab5mTffXf3XQoM2jQrD9wvK7CaaVZae6DvUmCpbHgDo7V2WWvttNba0a21R7XWXjN5/JrW2pmttce31l7WYRyHxUMHsgKzOYUEunEKCXTjFBJYDIkCAAAABk8DAwAAABi8oVxGdc1cGhK6kRWYzWVUoRuXUYVuXEYVFmO0R2A4rx+6kRWYzRoY0I01MKAba2DAYkgUAAAAMHhOIYElJyswm1NIoBunkEA3TiGBxRjtERgOi4duZAVmcwoJdOMUEujGKSSwGBIFAAAADJ4GBgAAADB4GhgAAADA4FnEE5acrMBsFvGEbiziCd1YxBMWY7RHYFiYELqRFZjNIp7QjUU8oRuLeMJiSBQAAAAweBoYAAAAwOBpYAAAAACDZxFPWHKyArNZxBO6sYgndGMRT1iM0R6BYWFC6EZWYDaLeEI3FvGEbiziCYshUQAAAMDgjb6Bce/n/r7vEhZi3x039F3C3C3jnJLkhj98c98ldLKMWbnrunf2XcJCLOO8xpKTL969t+8SFuK2D+7uu4S5W8Y5JeOZ1/57lzMry/j37zLOKRnPvNoD+/suYe7Gsk9fq9e8+sq+S5i7ZfxZjb6Bcd/n/qHvEhZi/2eX78P+Ms4pST7x7rf0XUIny5iVsfzjZa2WcV5jycm+u+/su4SFGMuH4rVYxjkl45nXA/d+ru8SFmIZ//5dxjklI5rXEjYwxrJPX6urfnP5GhjL+LMafQMDAAAAWH7VWuu7hiNSVXclubHvOhbgpCR39F3EnC3jnJL+5/VVrbWTZz1pSbPS93u/KMs4ryHMaWZWljQnyTDe/3lbxjkl/c9rM+9Tkv7f/0VYxjkl/c+ra1YeSHLvBtSzkbYm2dd3EQuwjPMawpyOba3N7cCJ0TcwAAAAgOXnFBIAAABg8DQwAAAAgMHTwAAAAAAGb7QNjKq6sKpurKqbq+olfddzJKrq1qq6rqo+WlUfnjz28Kp6d1X91eTrV/Rd5yxVdVVV3VZV1x/w2GHnUVUvnfz8bqyqb++n6tUdZk4/W1X/e/Lz+mhVXXzAnw1uTrIyLMuYk0RWhmQZcpLIygF/Nqh5LUtOElkZyu/UoYw9J4msDM0y5iRZjqysWWttdLckW5L8zySPS3JMko8lOavvuo5gPrcmOemgx/5jkpdMvn9Jkl/qu84O8/iGJE9Ncv2seSQ5a/Jz25bksZOf55a+59BxTj+b5MWHeO7g5iQrw7stY05WmZes9DOX0edkld8pWelxXsuUk8l8ZGWgWRlzTiY1ycrAbsuYk1XmNZqsrOc21iMwzk9yc2vtltbaF5K8Kcmzeq5p3p6V5HWT71+X5Nn9ldJNa+19Sf7+oIcPN49nJXlTa+3+1tpfJ7k5Kz/XQTnMnA5niHOSlYFZxpwksjICo8pJIisTQ5vXsuckkZVBGHlOElkZnGXMSbIUWVmzsTYwHpnkUwfc3zN5bKxakj+sqmuratfksa9srX0mSSZfT+mtuiNzuHmM/Wd4RVV9fHLY1vRwsyHOaYg1HYllzcqy5iSRlT4sa04SWel7XkOr50jJyvh+hmPISTLMmo7EsmZlWXOSjCcrazbWBkYd4rG24VXMz9Nba09NclGSH6mqb+i7oA0w5p/hbyR5fJKnJPlMkl+ZPD7EOQ2xpiOx2bIy9p+frPRjs+UkGf/PbyxZGVo9R0pWVozlZziWnCTDrOlIbLasjP3nN6asrNlYGxh7kjz6gPuPSvLpnmo5Yq21T0++3pbk97NyKM/fVdVpSTL5elt/FR6Rw81jtD/D1trftdb2t9YeSPLqfOnQqyHOaYg1rdsSZ2XpcpLISl+WOCeJrPQ9r6HVc0RkJcmIfoYjykkyzJrWbYmzsnQ5SUaXlTUbawPjQ0nOqKrHVtUxSS5NcnXPNa1LVW2vquOn3yf5tiTXZ2U+3zd52vcl+YN+Kjxih5vH1UkuraptVfXYJGck+WAP9a3Z9C+6ie/Mys8rGeacZGUcli4niaz0YclzkshK3/NaipwksjKg36nORpSTRFbGYulykowuK2vX9yqi670luTjJTVlZPfUn+67nCObxuKysBvuxJDdM55LkxCTvSfJXk68P77vWDnN5Y1YOU/piVjp8l682jyQ/Ofn53Zjkor7rX8OcfjvJdUk+npW/CE4b8pxkZVi3ZczJKvOSlY2fw1LkZJXfKVnpv/7R52QyD1kZyO/UGuY0mpxMapKVAd2WMSerzGtUWVnrrSYTAQAAABissZ5CAgAAAGwiGhgAAADA4GlgAAAAAIOngQEAAAAMngYGAAAAMHgaGAAAAMDgaWAAAAAAg6eBsWSq6vNHul1V/fd1jvGwqvrhgx5b11iHGf/rqurn5jUem5usQDeyArPJCXQjKxypaq31XQNzVFWfb609dKO2O2iM05O8vbX2pCMZBzaCrEA3sgKzyQl0IyscKUdgLKGqOr2qPllVr66qG6rqD6vq2AP+/F9W1cer6mNV9duH2P7zs8apqrdV1bWTx3dNNv3FJI+vqo9W1S8fONbk+39XVddPbi/qUutBdb2lqr5+Tm8TyAp0JCswm5xAN7LCEWmtuS3RLcnnk5yeZF+Sp0wee3OS502+PzvJjUlOmtx/+HS7A8eYfF1tnOl2xya5PsmJk+dff3A9k69fm+S6JNuTPDTJDUnOXe01DjG3TybZ0fd77LYcN1lxc+t2kxU3t9k3OXFz63aTFbcjvTkCY3n9dWvto5Pvr81K+JLkm5O8tbV2R5K01v5+neP8m6r6WJK/SPLoJGfMGOfrk/x+a+3u1trnk/xekn824zUeVFUPSXJ0a23vjNeBtZIV6EZWYDY5gW5khXXZ2ncBLMz9B3y/PyvdxySpJGtZ+OTLxqmqC5J8a5J/2lq7p6rem+QhM8apddR6oLOTfGJWsbAOsgLdyArMJifQjaywLo7A2Hzek+R7qurEJKmqh69jjB1J/mHyF8ITk/yTyeN3JTn+MNu8L8mzq+q4qtqe5DuT/Lc1vOY5ST6+jlphvWQFupEVmE1OoBtZYVUaGJtMa+2GJC9L8qeTw6r+0zqGeWeSrVX18SS/kJVDs9Ja+2yS908Wvvnlg173L5O8NskHk3wgyW+21j6yhtf0lwIbSlagG1mB2eQEupEVZnEZVQAAAGDwHIEBAAAADJ4GBgAAADB4GhgAAADA4GlgAAAAAIOngQEAAAAMngYGAAAAMHgaGAAAAMDg/X9bfNKQ1tArJwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x432 with 12 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "magmap = {'M':'MAD', 'S':'SANE'}\n",
    "\n",
    "def fill(ax, a, b, off, color):\n",
    "    L = inc_edges[a]\n",
    "    R = inc_edges[a+1]\n",
    "    B = Rhigh_edges[b]\n",
    "    T = Rhigh_edges[b+1]\n",
    "    \n",
    "    f = 0    + 0.25 * off / 5\n",
    "    ax.plot([f*L+(1-f)*R, f*L+(1-f)*R], [B, T], color=color)\n",
    "    f = 0.25 + 0.25 * off / 5\n",
    "    ax.plot([f*L+(1-f)*R, f*L+(1-f)*R], [B, T], color=color)\n",
    "    f = 0.5  + 0.25 * off / 5\n",
    "    ax.plot([f*L+(1-f)*R, f*L+(1-f)*R], [B, T], color=color)\n",
    "    f = 0.75 + 0.25 * off / 5\n",
    "    ax.plot([f*L+(1-f)*R, f*L+(1-f)*R], [B, T], color=color)\n",
    "\n",
    "def plot_notused(ax, df):\n",
    "    for i, (_, r) in enumerate(df.iterrows()):\n",
    "        a = i %  len(inc)\n",
    "        b = i // len(inc)\n",
    "        \n",
    "        if not r.Fnu_86GHz:\n",
    "            fill(ax, a, b, 0, 'r')\n",
    "\n",
    "        if not r.major_86GHz:\n",
    "            fill(ax, a, b, 1, 'g')\n",
    "            \n",
    "        if not r.mi3_230GHz:\n",
    "            fill(ax, a, b, 2, 'b')\n",
    "            \n",
    "        if not r.Fnu_NIR:\n",
    "            fill(ax, a, b, 3, 'y')\n",
    "            \n",
    "        if not r.nuLnu_xray:\n",
    "            fill(ax, a, b, 4, 'k')\n",
    "            \n",
    "    ax.set_xlim(inc_edges[0],   inc_edges[-1])\n",
    "    ax.set_ylim(Rhigh_edges[0], Rhigh_edges[-1])\n",
    "    ax.set_yscale('log')\n",
    "    \n",
    "def plot(ax, df):\n",
    "    \n",
    "    #print(df.set_index(indices).sum(axis=1))\n",
    "    #print(df.set_index(indices).sum(axis=1).values.reshape((4,9)))\n",
    "\n",
    "    v = df.set_index(indices).sum(axis=1).values.reshape((4,9))\n",
    "    ax.pcolor(inc, Rhigh, v, shading='auto', cmap='Blues')\n",
    "    \n",
    "    ax.set_xlim(inc_edges[0],   inc_edges[-1])\n",
    "    ax.set_ylim(Rhigh_edges[0], Rhigh_edges[-1])\n",
    "    ax.set_yscale('log')\n",
    "\n",
    "fig = viz.grid(\n",
    "    accept.reset_index(), plot, aspin=aspin, mag=mag,\n",
    "    figsize=(15,6),\n",
    "    xspace=0.05, yspace=0.05,\n",
    "    xtitle=r'$a_\\mathrm{{spin}}={}$', ytitle=r'mag={}',\n",
    "    xlabel='Inclination $i$',         ylabel=r'$R_\\mathrm{high}$',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2d209e4-07e1-4cf4-9808-034f3afccc28",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd33a870-cf79-4e4f-a4a1-5c13a642b9f8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}