"
],
"text/plain": [
" NAME DATE POS_7DAYAVG DTH_NEW DTH_NEW_2WK\n",
"0 Adams 2020-03-16 0.0 0.0 0.0\n",
"1 Ashland 2020-03-16 0.0 0.0 0.0\n",
"2 Barron 2020-03-16 0.0 0.0 0.0\n",
"3 Bayfield 2020-03-16 0.0 0.0 0.0\n",
"4 Brown 2020-03-16 0.0 0.0 0.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.read_csv(\"wi-covid.csv\")\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For clustering, we'll want each county on its own row. Historical data should be spread across columns, so that we can group rows based on year-long trends. Pandas's `pivot` function can automatially take three columns and put the values of one to the index, the values of another to the columns, and the values of the third to the cells of the final table:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
DATE
\n",
"
2020-03-16
\n",
"
2020-03-17
\n",
"
2020-03-18
\n",
"
2020-03-19
\n",
"
2020-03-20
\n",
"
2020-03-21
\n",
"
2020-03-22
\n",
"
2020-03-23
\n",
"
2020-03-24
\n",
"
2020-03-25
\n",
"
...
\n",
"
2021-02-21
\n",
"
2021-02-22
\n",
"
2021-02-23
\n",
"
2021-02-24
\n",
"
2021-02-25
\n",
"
2021-02-26
\n",
"
2021-02-27
\n",
"
2021-02-28
\n",
"
2021-03-01
\n",
"
2021-03-02
\n",
"
\n",
"
\n",
"
NAME
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Adams
\n",
"
0.0
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
...
\n",
"
3.00
\n",
"
2.71
\n",
"
2.57
\n",
"
3.14
\n",
"
2.14
\n",
"
1.86
\n",
"
1.29
\n",
"
1.29
\n",
"
1.57
\n",
"
1.29
\n",
"
\n",
"
\n",
"
Ashland
\n",
"
0.0
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
...
\n",
"
0.43
\n",
"
0.43
\n",
"
0.71
\n",
"
0.57
\n",
"
0.71
\n",
"
0.71
\n",
"
0.57
\n",
"
0.57
\n",
"
0.43
\n",
"
0.29
\n",
"
\n",
"
\n",
"
Barron
\n",
"
0.0
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.0
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
0.00
\n",
"
...
\n",
"
7.71
\n",
"
7.29
\n",
"
7.57
\n",
"
7.71
\n",
"
8.14
\n",
"
7.71
\n",
"
8.29
\n",
"
9.00
\n",
"
8.86
\n",
"
8.57
\n",
"
\n",
"
\n",
"
Bayfield
\n",
"
0.0
\n",
"
0.0
\n",
"
0.00
\n",
"
0.25
\n",
"
0.2
\n",
"
0.17
\n",
"
0.14
\n",
"
0.14
\n",
"
0.14
\n",
"
0.14
\n",
"
...
\n",
"
0.14
\n",
"
0.14
\n",
"
0.29
\n",
"
0.29
\n",
"
0.14
\n",
"
0.29
\n",
"
0.29
\n",
"
0.29
\n",
"
0.43
\n",
"
0.43
\n",
"
\n",
"
\n",
"
Brown
\n",
"
0.0
\n",
"
0.0
\n",
"
0.33
\n",
"
0.50
\n",
"
0.4
\n",
"
0.50
\n",
"
0.43
\n",
"
0.43
\n",
"
0.43
\n",
"
0.29
\n",
"
...
\n",
"
23.29
\n",
"
23.57
\n",
"
18.86
\n",
"
18.71
\n",
"
21.57
\n",
"
20.86
\n",
"
20.14
\n",
"
19.57
\n",
"
19.71
\n",
"
18.00
\n",
"
\n",
" \n",
"
\n",
"
5 rows × 352 columns
\n",
"
"
],
"text/plain": [
"DATE 2020-03-16 2020-03-17 2020-03-18 2020-03-19 2020-03-20 \\\n",
"NAME \n",
"Adams 0.0 0.0 0.00 0.00 0.0 \n",
"Ashland 0.0 0.0 0.00 0.00 0.0 \n",
"Barron 0.0 0.0 0.00 0.00 0.0 \n",
"Bayfield 0.0 0.0 0.00 0.25 0.2 \n",
"Brown 0.0 0.0 0.33 0.50 0.4 \n",
"\n",
"DATE 2020-03-21 2020-03-22 2020-03-23 2020-03-24 2020-03-25 ... \\\n",
"NAME ... \n",
"Adams 0.00 0.00 0.00 0.00 0.00 ... \n",
"Ashland 0.00 0.00 0.00 0.00 0.00 ... \n",
"Barron 0.00 0.00 0.00 0.00 0.00 ... \n",
"Bayfield 0.17 0.14 0.14 0.14 0.14 ... \n",
"Brown 0.50 0.43 0.43 0.43 0.29 ... \n",
"\n",
"DATE 2021-02-21 2021-02-22 2021-02-23 2021-02-24 2021-02-25 \\\n",
"NAME \n",
"Adams 3.00 2.71 2.57 3.14 2.14 \n",
"Ashland 0.43 0.43 0.71 0.57 0.71 \n",
"Barron 7.71 7.29 7.57 7.71 8.14 \n",
"Bayfield 0.14 0.14 0.29 0.29 0.14 \n",
"Brown 23.29 23.57 18.86 18.71 21.57 \n",
"\n",
"DATE 2021-02-26 2021-02-27 2021-02-28 2021-03-01 2021-03-02 \n",
"NAME \n",
"Adams 1.86 1.29 1.29 1.57 1.29 \n",
"Ashland 0.71 0.57 0.57 0.43 0.29 \n",
"Barron 7.71 8.29 9.00 8.86 8.57 \n",
"Bayfield 0.29 0.29 0.29 0.43 0.43 \n",
"Brown 20.86 20.14 19.57 19.71 18.00 \n",
"\n",
"[5 rows x 352 columns]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df = pd.pivot(df, index=\"NAME\", columns=\"DATE\", values=\"POS_7DAYAVG\")\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Clustering: KMeans and Agglomerative"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's use KMeans clustering to group similar rows. Note that for this example, we're just clustering based on total cases (for many use cases, it may make sense to scale by population or other otherwise standardize the data)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,\n",
" n_clusters=8, n_init=10, n_jobs=None, precompute_distances='auto',\n",
" random_state=None, tol=0.0001, verbose=0)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"km = KMeans(n_clusters=8)\n",
"km.fit(df)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 1, 4, 1, 7, 1, 1, 4, 4, 1, 4, 1, 3, 5, 1, 1, 4, 5, 1, 5, 1, 4,\n",
" 1, 1, 1, 1, 1, 4, 1, 5, 1, 5, 1, 1, 1, 4, 5, 1, 1, 1, 2, 4, 1, 1,\n",
" 0, 4, 1, 1, 4, 4, 1, 5, 1, 5, 1, 4, 1, 4, 5, 4, 1, 1, 1, 1, 4, 1,\n",
" 5, 6, 4, 1, 0, 4], dtype=int32)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"clusters = km.predict(df)\n",
"clusters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's put those cluster IDs back in the GeoDataFrame so we can plot it. The cluster IDs are in an order corresponding to the DataFrame we fit to (`df`). This is probably the same as the order in the `wi` DataFrame, but we should explicitly assign the index to each cluster ID before adding it to a different DataFrame, just to be safe."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"wi.plot(column=wi[\"cluster\"], cmap=\"tab10\", figsize=(8,8))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compactly, we could have done all the fitting and plotting in just three lines (note that clusters will probably be similar, but the colors for each cluster might change each time we run this):"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAHSCAYAAABclZoLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3wcZ53/38/M7qpbsqqbXOIa20mc2InTG0kINYELoQQCHFwu1IM7OAIHP45ywB1HOEJNgEAgIQmEhCTgFCchxXHce5dcJVm2JKu3LTPP749dyZLVVrs7O8+snvfrZUuanfLd2dnn85RvEVJKNBqNRqPRjA/DbQM0Go1Go/EiWkA1Go1Go0kALaAajUaj0SSAFlCNRqPRaBJAC6hGo9FoNAmgBVSj0Wg0mgTwuW1AqiktLZWzZ8922wyNRqPRZAibN29uklKWnbk94wR09uzZbNq0yW0zNBqNRpMhCCGODrddT+FqNBqNRpMAWkA1Go1Go0kALaAajUaj0SSAFlCNRqPRaBJAC6hGo9FoNAmgBVSj0Wg0mgTQAqrRaDQaTQJoAdVoNBqNJgG0gGo0Go1GkwBaQDUajUajSQAtoBqNRqPRJIAWUI1Go9FoEkALqEaj0Wg0CaAFVKPRaDSaBNACqtFoNBpNAmgB1Wg0Go0mAbSAajQajUaTAFpANZoxCFs2TZ1Bt83QaDSKoQVUoxkDKeH+NYdH3WfdoVOjvt4btvjy4zt4771v8PlHt6XSPI1G4xI+tw3QaFQn4DP44MWzBm2TUiKE6P+7ND8w6jl++MIBHt5QA8DS6ZNSb6RGo0k7egSq0cTBtKIcAE609bL/RAcv7G0Y9Pq88oJRj3/tQFP/79efPSX1Bmo0mrSjBVSjGYHesEVNc/egbXc9voNXDjQwZVI2z+ysH/TaY5trRzxXW0+4//d3nDc1tYZqNBpX0AKq0YzAd1ft5QO/WkdzV6h/m2VL/uGCGRTm+JlTljdo/111bVQ3dAx7rtyA2f97yLKdMVij0aQVLaAazQjYEmqae3h+9wkguu55/4cvpDgvwMySXBZNmdS/HWBSto/PP7p9yHksW/Z78d62cmb/cRqNxttoAdVMODYeaeZ4a8+Y+9128Uw+culsblo2HYB7XqzmnG88x+ajLf373PzT13lkY9Q56PxZk1lWWTTkPE9tr8M0DD5y6Wy+ffPSFL0LjUbjNqKv95wprFixQm7atMltMzQK8uC6o6w92MT2mja+dfMSrl1UAUSnXrP9xpiOQCu/8wI3LpnC19+xBMOIeuA+ua2OlXNKmFKYjW1LhGCQdy5AbUs3UyZl4zN1f1Wj8SJCiM1SyhVnbtdhLJoJw09equZEey9FuX6uXlDev33hlAL+vq+hX0Ab2nt5taqJW5bPGHT885+7ikk5PoQQvLTvJMdbeweFt/SJ6pnMmJzrwLvRaDRuo7vEmgnD45+8lGsWlnHejCKCkdOOPJ29Ec6eenpdsjMY4aoFZUB0ffNkey8Ahbl+hBA8u6uef/3jdsoLsoZcozds8ejGY+w+3ubwu9FoNG6jR6CaCcO0ohx+89GLBm3rDVtc9f2/0x2yuGZROb/44HLOKsvvf/3/Xqjit2uPsKyyiF9/eAWmIXho/TF++oELuGxeKQAHGzt5fEst22paqW7o5GR7kJK8AOu/8iY9bavRZDBaQDUZS0tXkMl5p0eJf9laR7bf4E1nV+A3DX73xhHufeUQ7b0RAFbvOUl3KEKO38RnGnT0hrn/9cN09EZ45UAjv3n9CP905Vn87h8vQgjBa1WNvLi3gUc2HqM3PDg05evvOFuLp0aT4WgB1WQsA8WzJ2Tx7b/toakzxPSiHK5ZVMajG2sIW4Od6P772X2s3nOS/3jbYt553jT++cqz+N/nDwDwX6v2cuPSKVQW53KyvZdPPriFjmBk0PGTsn3cuqKS6xfrZAkaTaaju8iajGdXXRufeXgLTZ3RhAh1rT08uO7YEPEEeHDdMU62B3lo3VEA7rxqLudMLyQ3YFJWkEVNSzQz0S9eOUhHMEKWz+CK+aXMmJzDLz64nLedO42vvn0xOQMSJ2g0msxEj0A1GUtPyOKux3fw5Lbj4z72YGMnD284xq0rKvnTnZfgNw3MAV62t62cyeGmLm5ZPoO3LJ2KEQtfESKatq8wxz/ofJYtBx2v0Wi8jxZQTcay8UhzQuIJ0NQZ4suP76SjN8wdV87t336kqQu/Lxoz+tszHJIA3rxkaKL42pZu6tt6uXB2cUK2aDQaNdECqvEkT26rY/msyUNiLJ/ddYKQZVOaF+A//rIzqWucV1nEpOzTI8mntx/ny4/v5LvvPofpseoso7HhcDOPbqzhZHsvP37/+UnZotFo1EMLqMZTdAYjbDzcTEG2r188719zmH+8fA4Av193hNerRy9uHS83L5vG+y6aySsHGllWWcTbz53K9YsryPZH1zdtW/LnLbWcPXUSS6cXDjq26mQHX39qN6X5AX7z0Qvxa49cjSbj0AKq8QzNXSHe+qPXkEhe/9K1AOytb+fhDceYV55PS3eI3cfbU3KtkrwANy+bztqDTXzkNxu470MrBoknwA9W7+enfz9IaX6AF/71KopyTxfVnl9RwDP/ckVKbNFoNGqiBVTjCaSUfOnPOzgRywp054Ob+djlZ/GhX68nYktuv39DSq/3qWvmMTkvwE9eqmZhRQHXLiof9PqGw83c+8ohILpe+qFfb+Bnt11AZbFO26fRTBS0gGqUpSdksf9kB0dPdbGmqonVe072v/bC3gZe2Nvg2LWnFWUTithUN3Ry/0cuHORB+7cd9dz15x1EbMm5MwpZWFHA7NI8phZmO2aPRqNRDy2gGuXoDVs8sPYIP/l7NR29kbEPSDHTi3K47uwKLClZ9S9XUJofTciws7aN/31+P28cOsVlc0u4/ZLZXL2wbEj1FY1GMzHQAqpRimDE4kO/Xs/GIy1j7+wAPkNw3+3L8ZkGPiAr30RKyf97cje/jyVXWHvXtUyLwwtXo9FkNmO6BgohFgohtg341y6E+JwQolgIsVoIURX7OTm2vxBC3COEqBZC7BBCXDDgXB+O7V8lhPjwgO3LhRA7Y8fcI2Jd+pGuocks2nrCQDTZwHdX7XNNPAHuef/5LJk22KN29/H2fvHsw7Yzq46uRqMZP+MqqC2EMIE6YCXwKaBZSvk9IcRdwGQp5ZeEEG8FPgO8Nbbfj6SUK4UQxcAmYAUggc3AcillixBiA/BZYD2wCrhHSvmMEOJ/hrvGaDbqgtreoralmw/9egMCaOgI0hlM/5QtwMKKAs6rLOS//+FchBA0dPTyyv5G/ue5/TR2BAftaxqCX92+gmvOcCzSaDSZSaoKar8JOCilPCqEuAm4Orb9AeBl4EvATcDvZFSZ1wkhioQQU2P7rpZSNscMWg3cKIR4GZgkpVwX2/474Gbgmdi5hruGxsPYtmR7bSvP7T7Jfa8exO3BXI7f5M+fvJT8rNNfh39/bAcv72/sf/1nt11AbWsPFQVZXLmgbFA4i0ajmZiMV0DfBzwc+71CSlkf+/0EUBH7fTpQM+CY2ti20bbXDrN9tGsMQghxB3AHwMyZM8f3jjKY+rYemrtCzCnNIzeQ/uXu5q4QAZ9BfpYP25bUtfbw/57cxSsHGl0XzYEU5wXwnZGn9u3nTuPl/Y2U5gf41k1L9WhTo9EMIe5WVQgRAN4JfPnM16SUUgjhaJM42jWklPcB90F0CtdJO7yAbUu+9+w+fr3mMJYtuXxeKfd/5EICPuez4di25NndJ7j3lYPsrGvDllCWn0V3KEJXyHL8+uPlkrNK+OltFwwZUa492MRZZXn86L3nc86MwhGO1mg0E5nxDEveAmyRUvYF450UQkyVUtbHpmj7gvLqgMoBx82Ibavj9HRs3/aXY9tnDLP/aNfQjEJnKMJ9rx7q/3tNdRP/9qft/Oi9yzAcqgjSG7Z4Ye9JfvHKQXbVDc4G1NIdIqLSkHMAPlNQnBcYsv3uW5e5YI1Go/ES4xHQ93N6+hbgKeDDwPdiP58csP3TQohHiDoRtcUE8DngOwM8aW8AviylbI559l5M1InoduDHY1xDcwahiM13Vu2lsSPIvPL8Ia8/vf04ZflZfO3tZ6csbrGjN8wzu06wpqqJv+9rGFJcug9VxRPg9eomGjp6KS84nQThcFMXpfkBCrL9oxyp0WgmOnEJqBAiD7ge+OcBm78H/FEI8THgKHBrbPsqoh641UA38FGAmFB+C9gY2++bfQ5FwCeB3wI5RJ2HnhnjGpoz+MoTO3lsc+2o+9z/+mHCls21Z5dz5fyyhOpThiI2r1U18ujGGv6+v2HYotQD8ZtizH3cxJbwp021fOqaeUA0ZeDUwmx21LYxvzyfycOMTjUajQbGGcbiBSZqGMu+E+28tK+B/3l2f1z737RsGv9zy7lk+eLzJg1FbJ7cVse9rx6iuqEzbrvyAyadCq59DqQ0P8CaL107ZB304Ilm5k7RNTw1molOqsJYNIoytyyfV2JhF/Hw5LbjVJ3s5GtvX8yFsyfT3B1iTVUTW4+1srWmhbaeMJWTcynI9nGqM8Se+na6ExBC01Q/zV1TZ4i/7qjnluUzBm2fXVZE1cl25pYVOLZ2rNFovIsWUI/RG7ZoaA8ys2Rw1Y9vPL2bB9cdG9e59tS38/5frsMQDBtWUtPck4ypAJiG+nUwpxflcPXCsiHbTdNgblkBG480c+HsYi2iGo1mEOq3bppBdIcsvvCn7eyqa+ONg6foDVtIKdle05bwOZ308TEVT7R+9tRJ/P5jF/UnjD8TwxBcNKeYrTWtWAo7Q2k0mvSjR6AeI2LZbDjSzDt/soYcv4khBAXZPo639bpt2rCoOmjLz/LxX+9ayiVzSwZ54A7kVGeQn/79IK8caGBZ5WQeWneU77z7HJ2FSKPRAFpAPUc4NgqyJf2JCUYKH1EBv6nmJMfl80q5adn0Ufcpyc/ihiUV3P/6YQ42dgFQ29rDL29fQWGODnHRaCY6arZumhEJhtX2aD2TLL+aj9hwsbJnsvVYCx/77cZB2zYcbua9975BQ4eaI36NRpM+1GzdNCMSjNhumzAuVB2BXjqvZMTXwpbNr147xHvvWzds+sF9Jzq46SevjyucR6PRZB56CtdjeE1Az0zSrgKXnFXCpXNLh2yvbujghy9Usba6iZbu8KjnqG/r5T2/WMu9H1rBRXN0rKhGMxHRAuoxvDaFayjmhbuwooDvv+fc/r/Dlk1Nczc/fqmav2yrYzx5RVq6w3zwV+v5/nvOHXM9VaPRZB5aQD2G10agiuknM0tymV6UA0Rjar//3H5+veZwwucLWTb/8sg2dta28e83LkpLxRuNRqMGWkA9htcEVLUR6At7T3LTT9ZwxYIyXtzbwL4THSk576/WHGZNdRNfeevZ1LX28MiGY9ywZAqfuGquTsCg0WQoWkA9RjCip3CTQUqobenhp38/mPJz7zvRwQ9X72drLKnF9to2Xj3QyHfffQ5nlY3t9avRaLyFFlCPEQx7awSqmH4CTqcXHPyG1x9u5s3/9yo3L5tOtt/kzUumcPn8oQ5MGo3Ge2gB9Ri9HhuBqiigTs6o2gz1Qgpbkj/FSs1tPtrCXz9zef+07trqJho7g9oJSaPxIFpAPUav10agbhswDE5WiBnLi3dPfTu/X3eURVMKeGlfA5uOtrD5aAvBsM2tF1Y6ZpdGo0k9WkA9htfWQFUcgjqa4D6OMJj/WrWX0BnOYH3pGIMRK+4arRqNxl20z73HOLPh1Ywft51ih/sMq0528LHfbuRfH93ugkUajSYR9AjUY3gtjMVWsARYR6+DyfcTFOdHNtYAUJjj575XD9LaHWZyboAPrJxJXpb+mmo0KqK/mR7Da164qoWxAGOm6UuGZHP/tvWE+c6qff1/F+b49dqoRqMoegrXY3huDVRBHM1vn+IB94Prj6b2hBqNJmVoAfUYnvPCVW8AinDQKJliBd1R28bdqw+k9JwajSY1aAH1GF6LA1VvBdRZL1wn3u/L+xuobkhNykGNRpM6tIB6DK9VY1FwAIrPwTjQcZVziZMdtW28vL8R6cC5NRpN4mgB9Rg9HhNQFTEdjGORDnUZfvHKQbYca3Xk3BqNJjG0gHqMnpAW0GQxhXOPvVOjxKbOEK/sb3Dk3BqNJjG0gHoMr8WBqjjp6HIiooSpbe1x8OwajWa8aAH1GDoTUfI4KuoOnvz53Sc50dbr3AU0Gs240ALqMbw2AlVxCOp2Kr9E6QxG+PWaQ26bodFoYmgB9RjeS6SgnoI6mR3JdthT9ldrDtPW41wmJY1GEz9aQD2G16ZwVbTW0QGow6NbKeHBdTo7kUajAlpAPYaewk0ew8GnPh2hmj9cfYDukIMJ8TUaTVxoAfUYXhuBqoiTYSxOT+ECTJ+cQ45f1wzVaNxGC6iHsGxJRMHyYKORDkEZLw7qZ1pGoDXN3dzyizf48YtV7KhtVbJknEYzEdDlzDyEF0efCuonPgfncNPRYbAlbD7awuajLfxg9QFK8wNcMb+MKxeUcsX8Mkrzsxy3QaPRaAH1FN7zwFVzBOpkKj/LhdFgU2eIJ7bW8cTWOgDOnjqJKxeUctWCMpbPmkyWT0/3ajROoAXUQ3jOgQglfYhwMpd8twKpFvfWt7O3vp17XzlEbsDk5mXTWVCRz1ULy5ldkutoOTeNZiKhBdRDBD1WCxTUnMI1HByBqtbJ6Q5Z7Drexh82HIOn91BZnMMV88u4Yl4pl84tpTDX77aJGo1n0QLqIbw4hatiCS4nx18hxT+jmuYe/rD+GH9YfwxDwLkzirhyfilXLyrnvBlFjk5vazSZhhZQD6Ha6CYeVHQQdXIGM2wp+BmN8BnYErbVtLKtppV7XqqmMMfPZfNKuGJ+GZfPK6WyODe9dmo0HkMLqIfwooC64VQzFlI6p6ARS733G++Qu60nzKqdJ1i18wQAs0tyuXx+1LP30rklFGTr6V6NZiBaQD2EF6dwLRWncIVzNs0syUOIaL5dQdTj1xACISA3YLLxSItj1041R051c+TUMR5cdwyfIVhWWRQdnc4v5bwZhfhMHUaumdhoAfUQXowDtWz1bBYOroIebuoa8bWL5hQ7dt3RSMU6dMSWbDrawqajLfzwhQMUZPu4dG50uveqBWV6ulczIdEC6iH0FG5qmGiOMk50GDp6Izy3+yTP7T4JwJzSPK6YX8pl80q5ZG4Jk/R0r2YCoAXUQ3hxBKrimqCT5cyUJA1v93BTF4ebuvjdG0cxDcF5Mwq5fH4ZV8wvZVllEX493avJQLSAeghPCqiCI9B0CIpCl007li3ZcqyVLcdauefFKvKzfFx8Vkn/CHVuWZ5O5qDJCLSAeggh4KLZk7EBaYONxLKj/6SMhlCELZuQZROxJJLo+pdlS8KWJGLbBCN2WpMbqDiFO9Gabrffb2cwwgt7T/LC3uh075RJ2THv3mgyh7ICnbtX4020gHqI9p4wGxL04pxTmtfv4OI3BFk+A5/PwGcITCEwzai3qCnAMAwCptHvTWoYYCAQBpjitFcpgBBDV9gkgJTYRMv9zCzJ7R89S9kn5lFBF9FdidiSsGVjy6jo2rbEMAR+UxC2JEJEp4ODEYuwgtPCKqPa3TrR3stjm2t5bHMtEM3de8X8Ui6fV8pFc4rJ1qXaNB5BC6iHSMaJaOD0b9iWhEMWpCFva67foNuBFIR+Q1CQ7Scv2yTLZ+IzYqEjQuAbEDoixAAnGhHtCOQG3GmgI7Zk0ZSC2Bqs7LdNEu007DvR4cyFVVPQM+jL3Xvfq4cI+AxWzimOCWpZ9H5NMKcvjXfQAuohepMQIrfWT50aLYZtSXN3iObu8R+7Ytbk1BsUB5uPjjx7MK0oO42WqEsoYvNaVROvVTUB+yjND3DJ3FIum1vC5fNLmTFZh8to1EELqIfoCSc+YuwNR1JoSfyomUjBbQuGYjpplILvN16aOkM8vf04T28/jmkIPrhyJp+/fgFFuQG3TdNo0L7lHqI3GQF1YQQqhKK5cN02YDi0fo6JZUseeOMoV//vyzyw9oiaeYc1EwotoB6iJ4k1Szccb9SN/csUSYmTDHu7rd1hvv7Ubm78v1d5Yc9JJSv+aCYGqrZwmjNo6gxS396b0LF+JytIj4K6GX/Ua3CdnMIVDibPd5ODjV18/Heb+NCvN7DfKQcsjWYU9Bqo4vSGLe58cDNrqpoSTkrgMwzCVvoT0fuUFVD17OoMOrdG7WTyfBVYU93EW370Kh9YOZPPX7eAknwdV6pJD3oEqjgv72/g5f2NSWX08bk0AlU2ZZ6CZjkpoEq+4RRjS3hw3TGu/v7LPLD2sCcrF2m8hxZQxfnrjvqkz+HWVKrWz/hJJkRpLFT9HJygIxjhF68c5Lq7X2HVznq9PqpxFC2gihCxbJ7bfYLXqhr7t0kpaWgPJn1ut6ZSA4o6EU20JnWivV8hBDXNPXzyoS3ceu8b7KhtddskTYai10AV4d8f28HjW+sAuHHJFN6/ciY9IYsNR5qTPrdbI1BVvXAn0IAMmIjv9/Q73nikhXf+5HXuef/5vOPcqQghONLURV6WT+fg1SSNFlAFqG3pZtWu01O1z+4+wbO7T6Ts/O4JqJpNt6ojMlOATvObPMM97p99eCu/ePkgn7l2Hj96sYpv3rSU16ub+OuO40wpzObr71iCLSVZvsFpHm1b0hmKMCnbj5QSW8Lv3zhCWUE2bzt3anrekEZZtIC6jJSSX756yNE1MLemUlUNY1HTKufEc6Jp8kil0vbUt/OJh7YA8Kk/bKGxIxjbH9YfaqalO8THrziLO6+a23/M+sPN3PvqQS45q4QfvVjF5NwAda09FOcFuG5x+RDB1Uws1Jxjy2D+uLGG6+5+hf95dh8n2no51RXi9+uOOnpNwxDk+vVH3YeKfiVOrlOr2mFwkz7xhOjzUNXQSVNniMe31NLWE+5/bdXOel7e38h3n9lHd8iirrUHgOauEH/aVMvr1U287743iJyRFak75E7qTE160SPQNLL1WAtfeWInEVtS3dDJz185SFGO3/F0d4cau5hXnk91Q6ezFzqDYMSmND9AU2cordcdC6ngmCzgM4ikoTrOxCDxzzdsSb78+A7Or5zMyfZe1h5sGnHfr/5lV//v5/zn83zmTfO4/ZLZ7Dnezuce2conrp7LbStn6WoyGYwW0DQRjFhsOtIyKJ5TSmjpDo9yVOpw4ytc29LD+ZVFygloR2+Ei+YUj3hP+j8hKaPze6d/4DcEYdsGohskMlrTdOCxMvpTAHZsuNvn2NIn3n3TjFJGt/hNI1Yxp+9MsXJsQ2waHJYiY2Yaou9a0Sv17SOBLkdjTN1h5ZxiOoMRpIy+bxmrP4uUHD2VQImeGFk+g1U7T7Bq5/h8EHrCFt9/bj8PrTvWP0r92pO7WVPdxDdvWkrFJF1tJxPRAuowrxxo5MW9J3liSx0dLjZkbiU1UHG6NJm6m+dOL2RHXVsKrUkNK2ZNZsuxiROu0Rux2H28PeXnTSZhiZT0i2cfz+0+iWkIZpXk8dFLZ9PWE+aZXScIRWy+8OaFQLRzrddSvYkWUAcIRiy+/uRuNhxu5lBTl9vmAO5lI4rIzKqYEVK0AoiKHRUn8dKKft9o9o8bazjVFZ2NKcr109YTZmpRNk9sqeP8mUV8793n6ulej6EFNIVIKfnZywf5/nP73TZlCG6NQCMZFpcRdKkwueYMPJheqU88IVpRZqDzYFVDJy/sbeCGxRV87e2LycvSTbMX0J9SCvnrjnolxROGj41LB4rmUsg4PKgnSeHU23Uz9V9zV4hHNtbwWlUTH1g5k2WVRVw2r9Q1ezRjo5u3FFB1soNbf/EGn3l4q9umjIh7a6ATrGV3iwl2m0eK9Uz6vI6cdXzUtfbw/ef2U9OcuDOUJj3EJaBCiCIhxGNCiH1CiL1CiEuEEMVCiNVCiKrYz8mxfYUQ4h4hRLUQYocQ4oIB5/lwbP8qIcSHB2xfLoTYGTvmHhH7dox0DZXoCkb42AObUpJyz0ncGqHYGbY4p2hyJSUa/nTi1POs0gS9Hn2qT7xTuD8CnpVS3iKECAC5wFeAF6WU3xNC3AXcBXwJeAswP/ZvJfBzYKUQohj4OrCCqHf9ZiHEU1LKltg+/wSsB1YBNwLPxM453DWU4Z6XqjjmgZ6iUz32sbAyTEDduo9jIYTAENGZBsOIhrL4DIHPNDANMWgK3xCiX4D6wmsMEU24EQ2dOb2t79yC6OtSytOzGbEwG0OIWDhNNKSHvv2FGHRtIcA0jNi0vuBQY+eIIU5zSvMozPHTG7YIRWwitux/f0I4V4BcpeotW2taqSzOddsMzSiMKaBCiELgSuAjAFLKEBASQtwEXB3b7QHgZaLidhPwOxl9EtfFRq9TY/uullI2x867GrhRCPEyMElKuS62/XfAzUQFdKRrKMNftydfbiwduLUGqlKDlAoCPjVXPaIxp7ERfywUI5prR93kDLNLckcU0JK8AJuOtox4bHFewBGbVHKyvuvPO/jvZ/bx509cypRCHUeqIvG0BnOARuA3QoitQohfCSHygAopZZ96nAAqYr9PB2oGHF8b2zba9tphtjPKNZSguqFzSNyXqrg2AlWoQcpkVMyulAyW0+m5RkCl+9gdsjivspDSfGc6C5rkiUdAfcAFwM+llOcDXUSnUvuJjTYdffJGu4YQ4g4hxCYhxKbGxsbhdkk5q3bW89Z7XkvLtVKBWxOPlp1ZCtobVnNEJzy4CmoagoIsH5XFOUwvyqEkL0B2LGfzWAkNnOoQqjRhsmhKAd+++Rx82pVdWeL5ZGqBWinl+tjfjxEV1JOxqVliPxtir9cBlQOOnxHbNtr2GcNsZ5RrDEJKeZ+UcoWUckVZWVkcbyl5Xj3QGEu9phkd7zXsXsSLzloHG7voCEaoae6hrrWHU10hesM25QVZNHclX0g+EZLJRJRKKiZlcfetyxybqtakhjEFVEp5AqgRQiyMbXoTsAd4CujzpP0w8GTs96eA22PeuBcDbbFp2OeAG4QQk2PetDcAz8VeaxdCXBzzvr39jHMNdw3XuX6xUrPJypJpa6CKtK8ZTUNHkLrW3lH3capbZivwAd96S8cAACAASURBVBdk+/jl7StYPG2S26ZoxiBeL9zPAA/FPHAPAR8lKr5/FEJ8DDgK3BrbdxXwVqAa6I7ti5SyWQjxLWBjbL9v9jkUAZ8EfgvkEHUeeia2/XsjXMNVukMRHt5QM/aOmoxD1Q6BAu1+RqDCSP59F1Zy7owit83QxEFcAiql3EY0/ORM3jTMvhL41AjnuR+4f5jtm4Clw2w/Ndw13ObXrx3mhb0n3TbDE2Raw65CAzscKoycMoGwAqknH9lQQ2l+Fu9ZUamncBVHp/JLgLUHT7ltgmcIazfctKCornsOFToiHcEI331mH7aET1w9t397TXM3W461sPVYK/Mr8rlt5SwXrdSAFtBxU9vSzU4Fy1mNhVsNbDCiptdqxiHcb/jTiVPvNqJQT+T7z+1jTXUjbztnGhHb5p4Xq/rjZhdWFGgBVQAtoOOgprmbj/52I50eLFDsVnxbTyizBFTVAbU6zb63cSv+dDhsCa9Xn+L16lMsmlIwKOlEWUEWW4+1cP5M5bKbTii0gMbJ2uom/v3PO6ht8UbiBFVQJSwgVaga1+rFOFAViSjaQ/KdkUpsTXUTGw43854VM+gN2zR1BqkszuHC2cUsqCjg7KnagzcdaAEdg5rmbv7fk7t4rarJ02Lg1syUQjNiKUHZXLhuG5BuHHqulP2KD/MBhyybh9YfG7TtwXXHyPGbrPvKmyjM8afJuImLFtBR2HSkmVt+8UZKz3lbaRNFdduwAnnY2QVEsvIJ+3MJ+bIIGVmEhY8wPsIYWAhsBEHboNcWhG0I2xCyJWELQpYkaNkEI/aY3oOuCag7l3UMt3IKj4Vb5ercQqWUe+lgPJ9uT9hi/aFT3LBkimP2aKJoAR2FPfXtKT9ndk8zxvGq/gwWWSk6r1Vayc8K3j7i626EXwiReSPQUETNNzTB9NNRfAaolmRsvFP0n3hoC19929m8/6KZZPtNh6zSaAEdhf0nOlJ+TiMyfPWJpMnKG/VlN8qKmUIo5dWYyUw0AXXysVo6vTBWni32D4mUUQcjKaPfJSkllh39Zw/YVpafxcn2IBJJxJbYsfJvRqz2m884LYWG0TdzIDAMCJgGWT4TQxArAxct3WYIQW3L+EomWrbkG0/v4e7nD3DHlWfxmTfNT/Vt0qAFdFQ2OlAk2wg744QkzdHXO9yIb/OZwtPrxsMx0aYOVcXJGZVtNYmHqeVn+WjsdCeP73B0BCO8WtXIP14+h7ws3dynGp3mfwQilk1Nc+rFzgh2pfycANIcPWOJGzrmMzLv8VLUCXfCoWo3RkUns41HWrj+7ldo7FBH2DOFzGvhUsSa6iZ6HChdJYLjm4qJl7FGoBEXWn5TVY+bJLAVbbonWhiLqisDpoICCnDx3BKlYlwzBT2mH4H6I0e4bKqfoC3otQXdEegM2/RGJD1hK+GcmSLsTC/Q9gdgFL1348tzZuxaJmApGic40aaWVU3qLxTLCLVoSgE/fv/5zK8ocNuUjEQL6AgUVa/hgrUvjfi6LQzIzoesXAhkY/uzwRdA+rKwAjnYviwsM4uwP5uIkUXY9BPBh3AoDaCt4BRuJo5AQ4oK6ERDUf1EtYjcqYXZzCvPd9uMjEUL6AhY4fCorxvShp726D8Gz4W7Eb5carXzvoo2bAykMLARRISPiDCJYOL3dVM5zUdYCsK2ICQFIRtCFvRakqAliViSsCUJWXZKnH+MDBTQoKJhLBMNN7zK40G1R/7v+xv54QtV/Ov1C9w2JSPRAjoCZbPPYv8br7ltRtwURtopW/eHlJ1PAvizkb4A+AJg+pA+P5j+/m22L4A0oz9tXxZSmNimD8vwYxl+RCCbg8XT+88nOD1y6JtyFDH3fgZsP71PFDsWLiDpm4qW/a/bMuqRKW2JJaMhBXYstKDvXDK2jx0LPYhY0RCDiD12AoqBZGJcq1dR0VkH1Ow03r/mMG9eUsGSaYVum5JxaAEdgcYjh9w2YVykukERAOFeRLg34XPklVbwh4J3p8wmpxBE12tPx+VFtxlCYJoCUxiYhsBnQH27mp6ME82JSFU/XBUzQnUGI3z+0W38v7cv4fL5pW6bk1FoAR2BnEneSsasZI/cI2EsEgjbEuy+v4bHr+Dooo+J50TktgXDo6KAAhw42ckjG49pAU0x3mjhXMAXSFWSvTSh4PdWeERA48VnKniTJyhupKaMB4X7WDqlnwPoEegIXPGBDzNl7gI6TjWy7/VXaThykOkLFzN3+UVUbXiD4wf2um3iGaj3zRVGZn1hVfYqVlNOnENR/VQ6peK/3aAdiVKNFtARMAyThZdcDsCKt79r0GtLr72BZ37yAw5t2eiGacOi5BSuyKwRqIoOIn2oa5kzqDoCVW0t+qLZxfSELfaf7GBy7uihbprxk1ktXJrIzsvn5i9+jfOuf4vbpqhNpk3hqthJmaCoKqCK6SfzKvL5052XsPsbb9ZTuA6QWS1cGhGGwQVvvcltMwagXoOSaVO4PlN/XVRB1XwWiuknf1h/jH/63SZlnZu8jp7CTYKW+uNum9CPkqnNTJ+Kup4wCs/gkhfwsXJO8ZDtI93+gTG5fX+cue/AR2rwa6fja/tPJMSQi/l9BhsOp76iEYClaFZ/1aZwAV6rauLIqS7mlumMRKlGC2gSVMyZy6SyctobG9w2RU2nCtMHEbeNSB0qr4F2hSJsPNLithmDOHe6c4H7Sj7voN4QNEZrt0N1iCc4ek4qCfKLS3jH5+5SwvVO2qmvHJMsMsOmcJVttFFz5ONkuj1l10AVZZNinatMQQtokkyZt4CFl1zhthnYCk5pyQzzwlVymrwf9WxzsgKQrsw1Pv60udZtEzKSzGrhXOL8N7/dbROQCgooGTYC1Y32+Oh1oJ5uH7aiH4Z68wBRLp+nMxA5gRbQFDB90WIufvd7XbXBttRbbLSNzFpiD0bUmyY/jXpNd3uvc8+kmvKpLrdfMsttEzISLaAp4rL3foiP/ODnZOXmuXJ9O6LeCDTTBFTJZBUq4+CUtxsF4uNBRatuWjaNOaXutEuZTma1cC5TMqOSmeecR9X6tWm/tpIjUJFZU7imygIq1Gu6I7akrCCLvIBJtt8k228ghOgvL9df446o1vaVnxOc7qz0OQuFLZuwZffXrA34DYIRGyEgx29iGgIhYhV0DEFrd5i2ntFr+jqBSuvkhoDPX7eAT10zb0J1/vadaMdnCOaVFzh+LS2gKaanvd2V69oKeuHaGeZEFPCp3AipZ5stobEjSKPD1+kYZqp4xezJrnieqiOf8JW3ns3HrzjLbTMcoTsUYd+JDi6YORmIdlz6Ogk7atr49t/2sPGr15Hlc7YTrwU0xXS1OhM4PhZWRMERaIatEKiciUg9+YSIm45tbimZIgqa4zf52OVz3DbDMf62o543Dp7iVGeITUeb+ev2egCe+NSlzC7N4wtvXkgwYmsB9Rqm352Ezbal3gjUyrA1UKneMrPSuLlO6VacqCozpT1hix21bZxXWeS2KY6wvbaVqoZOVj+6jY7g6cHDP/9+Mz+/bTkXDZOVywnU7VJ7DCkl2577G6dqjrlyfRXXQCNk1hqoksO8GIoMfAYRUdTRx0kUWgLl3lcPum2CY2w83MLOurZB4glQ39pLxaT01XLWApoidr28mhfv/znSpWFKJJR+h4mxsDJsDVRh/VQSN8XEtUsr9JDcceXcMffpCkbYVdeWBmtSy9zyoV7FQsDHr5iTVoepzGrhXKK9qZGNT/7ZVRussHq5Lm2VWpMMR9U77ZZdqiZaSBfvOn86CypGTx5f39bDHb/fxLt+9jp/21GfJstSw20XzRyy7aoFZfzjZeld982sRSoXaD5ex6P/+SW621pds8EwTTXXQDOsf6ZySShV5cJnCsJW+q1z8oqmITDE6VAbU0T/Ng2B3xAU5vgRRNdhLVsSHiDmoTTFa9+yfAa5gZGb91DE5s4Ht7C9Jtpufe7RrdhS8o7zpqXFvmR5Zudxbji7jHkVBfzs5UMALJ46Ke0FH7SAJknHqUZXxROijku21eOqDcNhKTgumjIpmymTsjGMaIMnEP2OH/3Vuc78SwokkiyfQVlBdH1FSrD7mulYXGNfaGNfLKAVa0D78BkGflNgGAJTRK/b1wgPvFOy/7/BSKIlx/pGV4Y4bXvAVO9eg3sl4Np7wiyZVoBlRz8H245+PoYR++wkRCwbS0oiVuyftLEtyZLphWw+Otib3o7Fp1oy6hwV7a4O/ZDeODS6F77PEI6vDb/jvGnUtY7cHnQFI/zfCwf6xRMgbEme2n6ct50zVemqQ310RyQfu2wWSyuLOdLUTW1Ltytxv1pAk6SwfIrbJuAPBAj3qiig6o1AcwMm22rd7fA4Qbq8DsdLdNSe/hHokVPdCR9rS4lTA8V0eAfvrW/n1QONtHSFmFmcy32vHaKlK4Qtox3ILL/Ba1VNQ45bveckTV1ByguyHbcxWe6+dRkdvWF+8HwVP/7ABbjzlGkBTZotzzzptgkYPjU/xoiCAuqF3nUiTPQ1P6+Qjo+puqETgB+/VE0wYg2aQj/WPHLH4qI5xZTkpc+DNRlqmrupLM7lojnF7K1vZ6mDtWdHQ70WziNIKXn90d+z9Zmn3TYFw1QzXETFKdwM1U9Clg5S1QymMxgZ1/rzZ6+dj+mRL8i9rx7EsiU3Lp1CbYt7s29qDl0UR0rJuscfYd3jj7ptCgCGqebHGJHq9c9ULDydCtLlnDIRyMwnZHSuWVjG5fPVKnn2930NLJk+ifKCbO5efYD2njBZfoMv3rCQb998Tv9+y2cWEg6H8fv9abdRzZZXYXq7Onns21/j5KEqt03pxzDUEypQcwSqsCNtUqhancSbOPeQGEK9urJTJmVz963L3DZjCM/vOcEvXztEwGewq66dps4gAPPLC7hl+Qwg+tyXTcpxzUYtoOPENH00H1erurtQdAo3oqCAemSGaty4lbpuLFQO/RkZ5+6lzzSUmy249cJKJue5k4J0NPbWd7CtZrDDX5bPYFbxacF0e8pZC+g48Wdnc/ZlV7HjxWfdNqUfVRup82Ud0/MasHzZ2IaJNHzRn7H2yTJ9WMKHjYkUIBHYGEghouEgCGwEEiMavoHAim2zEdgyui0iBRZg2QKbaKhBXy+/73c7FmJSlKXmvUoWVZ8BRc0aHQdt9hkC1VKevGWp+5EEAznZ3ktZfla/M9RALp1bwoVzSgZVX3ETLaDjJNzby6Ftm9w2YzCKTuEGGg+RW1vjthmD6Dn7Ktay2G0zUo7bPfGRUKGRGy9OrpOrNk+wdPokZpXkum3GIN52z2tML8oZ9lPYWdfOwcZOzlKkQLgW0PEioPPU0BgqN1G1kVIxOxKGmtPdyaLoI6CsXa6hmIKunFMyasaidLP2YBNNnSGaOgeP0/OzfHQGI0wriiZCUaXNU3PoojD+rGxMxeIuhaIjUBlRT0BtI/2eeulA0SVQZe3SRJlXPnq+3HSz5ejQIujvPn86f//C1ZxXWUTVyU4iCoVsqaUEHmHKvAXU7dvjthn9qNIbOxNLwRJrUtHORrKo6kTkRZy8kyrNtC+aUsD7h0nK7ibNXUPT8TV2BikryOLT18zDtm0CDhfJHg9aQMdJfdV+DMNk+qIl/fNT0TRSEmL5UJE2tm0jbYm0LSzLIhIKYoXDREIhIqEgkVDqXAlUDWNJ5XtMFZk6As3U+FY3cPJOBnwGXSE1ZmacXPu0bMmrBxq5ZlH5uI6bP0wFmcVTJwFw/eIKQK2sW1pAx8nWZ5+mZs/OlJzLME1MfwCf34/h82H6fNGfhg9hGgjDwDBMhCEQwoj9FBBLgC5jX/Xs3Fymn70UZEy0pUTaNrYVwbZtbMtC2jZWJIxtWViRCHYkQiQcdrQMmm2r0VAMRCo6Wk+WDH1bGUfAp0ZnVwi4euH4xC1e7nv1IF3B6Hc/XgHtE9xTsVjPgZwplyql49QCOk5Mf+ripWzLwrZ6kk4EP33REur27U7iDALT748KuGFgBgJR4TZNDLNPxA1M0wDDjIZMGAY5efkEe3qiot7fgkfjRaSUnKjen9T7coQMK/Ldh6peuJrBFOb46Q1H1/AMARFL0hWKpD25wsIKZ6ZvpZT88rXDNHYE+cils0fdtzsUITfgo+pkB//xxC42HBm+ko1beW7jQQvoOLn2I3dwcPN6etrVqeKe/OhDYoVDp0ejXfEdNePspRzfr85acDxILaBpRVGzRmX94dMNeUGWSUcwdTMpB04OjW2EaDm6c2cUsWkYJxonOHCyg65ghLys1EqAEIIr5pXyl211VEwaXNUlbNn4DEFdaw/7T3Rw2bxo6sDHNteOKJ4Ak3PVXXbJzNbEQdoaTiglnlHcaqXUWYuIFztDBdSDOqUZQMiS/XVk04Et4YW9Jx05993vXcZzn7uSHL8x6D21doewbMna6lP8xxO7eN9969h6rIWphcOXT7vkrBJK8wPDll5TBT0CHSf+7ByEYSBtdVyp3cN7zXamjkBVRdUMSfHivS5i/Hz+0W2094T50CWzU37u+RUFQ0JkSvOzEEKw7tApTrT3cqK9ly8/vnPISLWPb928hOqGLj750Gbeed40JadydWsyTgrLKyitnOW2GYOQLn3N09ljThUZK6AeFypV8eAjHje2hK89uZtPPLiZ2pbEC5CPhBBiUIhd3+/FA/Lu7jvRwSsHGoc9fu3BU9y4dApPf+ZyJcUTtIAmRNmsOW6bMBiXvuRSem8UnqntoapfZK+H16Szk+jWs/nMrhN88U870nKttdVNca/z/uD5A3zj6d2cPWWSw1YljqrfO6U564KL3DZhEG4JmVQoHiteZIY+8uoOQL33jLiFmx/h9tpWWrudj9t+ekf9kAorI9HWE+Y3rx9h13HVfE5Ok5mtiYOcqq1hzSMPuG3GIE7V1kTjQNOMMAwMRUupjUSmxoEqi8fvd1r7iC7eq+6QxWce3kowxek3z0y719E7NNPQcBgCAmZUnopy1HXVUdcyBanbv5cnvvefBLvjjPNIE8GuTnDBqen4/j0UT6+kuU6tiiujYdrqpRcci7xAtAScjYyWZRswrShjZdrOdNYRYuj6Xd82IQaPdiR92bSi5+k7f9+alZQSQ4j+9l2I6MSsEGDEztT3Wt81+/7O8Rnk+M3+8xO7Vl8wvCGi1zREdLq3Kxzpj5McL7NLcsnymf3XHnhPjFio8mnbT7+Hvt8ZcF/6/kYOTpXZd68GErEkm48lH35yqjPIRXOK+88vgQ2HRw7vSDWvVTWx7BurKc4L8O4LpvO2c6eysKIgqVShPnPwGM1vxjdm+8KbF7KwooBHNtaQrVCy+zNR1zLF2Lf2VZ75yd3YCuZ3BVyb/5EKZhsajdyeU8Bkt82IG78p+GRDHMk7mnu5hpxhX5JIl9ci408+svucAlbVnEroKj5TsP9kR0LHJsPyWal5no6c6ubIqdPOPAXZ6W+ee8IWda09/Pilau595RAP/dNKLpxdnJJzN3YE2RpnR+Opbcf5+QeX8+2bl1JeMLyXrgroKdw4aDlxXG3xxL2E8paCFVdGwza81WdMhQ+Llxx5krHUrZAZpxyN3J79/v3HLiIr2IaVonns9973xqAOwnAU5wV4/a5ruXphGX/bcXzEEBdV8FZr4hJ2JKK0eEZxqfHwWDys5wTUbQPSTDI9erf0xksdlHg5d0YhK88qAUqGvNbQ3kt2wGRSdnwZgqSMJonoCo7dhlZMymZaYTa3XTyTXL/631U9Ao2DkhkzufTW29w2Q0lUTBg/GtJjAjrRMJIQI/fK+jnTzXHTy/1d508f8bWygizMcdzrp7Yf57+f3R/XSDZgRmNHKyfnUZKfFfc13EK3JnFyyT+8n/zJJTx/7z1umzIsbiU18NwINFMTKWQIyXw6ruXdzbwBKEumjZy4QAgRdw7dQ42dfOXxnXGXcBtv+TO30a3JOFh6zfXkFKgb1OsGtuWxEagWUKVJTkDd8qRz5rpuhlmnKpzlxb0N46p/2pdg3ivoEeg4aDh8kJ6OdrfNGBY9As1MvJguMRm86ETk1Ag0Xued4UKWkqUrBRVottW08tu1R+Lad2phNv92w0LOryxK+rrpRAvoODi8dZPbJoyCOw2t7TEBLYu0cH7lvOgf4vQPgQARDffoj3ccEOM4MIbxzFjAQa/JwbmJ5YAdRCzWse98fftPMg2m9UR3sPvqqcb2MRDQMnwJrExkOhFunAYCGwOJQEZ/ytM/BTaGtBGxG9uXHOMoea7Y7JRs90bi+25JCVMLs6hvG1qMOlFS0XF7YO0R6lrjq3X8X+9ayrWLKpK+ZrrRAjoOfIHUFdNOOW7lw/WYgAYivWyNM5VYuriqeBIVh+LL0JLpzOxsxd7y84SObbjs0ym2xjuk2oGqIE4P25F49UAjT2yti2vfmcW5XD6vLKnruUVcAiqEOAJ0ABYQkVKuEEIUA48Cs4EjwK1SyhYR/SR/BLwV6AY+IqXcEjvPh4Gvxk77bSnlA7Hty4HfAjnAKuBfpJRypGsk9Y6TIBJyPldk4rhVkcVbAqpiKr8sL1addgx9LxKhrTvM+TOLsGyJZUuyfEY0C5CMzodIGZ0StmX09ZqWHtp6hnbacnyCt507jdmluUnZ88yu+rj3/bcbFhDwedM3YTwj0GuklAMrm94FvCil/J4Q4q7Y318C3gLMj/1bCfwcWBkTw68DK4i29puFEE/FBPHnwD8B64kK6I3AM6NcwxUWX3ktm/72BMEutVL5gXtll7w2AlUxmbxfQVF3j2TuxcS9j10hi63H4p9ZKc6LjjCnFWazeEoelaUFFGT5CIfDfPyKuRQnkf3HsiUv7m2Ie39VS5XFQzJTuDcBV8d+fwB4mai43QT8TkYn0dcJIYqEEFNj+66WUjYDCCFWAzcKIV4GJkkp18W2/w64maiAjnQNV5hUVs7Ci69gx4vPumXCiLhVE9RrXrgqMp6YukwnmQ6OmHBpJ8ZPWUEWNy+bxoWzi5hfmsPs8qKUT/9uONxMQ0f867Gt3d5dvohXQCXwvBBCAvdKKe8DKqSUfeP0E0DfCvB0YGB28drYttG21w6znVGu4RqTp05z24RhcSfoemCKcG8gUG/ErEegp0kmq4/hsWcxncwrz+cz187jxqVTyPI5W0Hp12sOjWv/Lz62nUfvuISyAvUTJ5xJvAJ6uZSyTghRDqwWQuwb+GJsvdLRp3e0awgh7gDuAJg5c6aTZtDZkr7qCOPBjbVI0+fDinir96hiE2tO4KnHM5FJxFS6NQJV8ZkaiN8U/Oh9y0ZNjpAq1lY38cI4pm8B8rN8dPSGM1dApZR1sZ8NQogngIuAk0KIqVLK+tgUbd9dqwMqBxw+I7atjtPTsX3bX45tnzHM/oxyjTPtuw+4D2DFihWOPs+1e3c5efqEcSNXrzDUW08cCxUTKZhaP/tJZjBeIHuZUViAKQSmEb2vhogmZxCAaYBPgCkkpgATiSmi4TFHe3wcaUkwDERxBf389QvSIp4AW8ZZ1u3fb1zIJ66a62IaxuQYU0CFEHmAIaXsiP1+A/BN4Cngw8D3Yj+fjB3yFPBpIcQjRJ2I2mIC+BzwHSFEX+2fG4AvSymbhRDtQoiLiToR3Q78eMC5hruGKxxY/zoNh8c3PZEu3FiL9FoxbVDTiUjx9tczVK69f1DPfTwUX/qPHCGxEZBb/gfxcuOSKWm5jpQy7sQJANedXc4nr57nnEFpIJ7WpAJYI4TYDmwA/ialfJaoqF0vhKgCrov9DVEv2kNANfBL4JMAMeehbwEbY/++2edQFNvnV7FjDhJ1IGKUa6SdSCjE03d/V9mwDSvswlSqBzuNKoaxRCZYtqHRcOtOGEmsjav+6W0+mp7IPyEEX3zzQvxxTqmcaO8lYqnZnsbLmCNQKeUh4Lxhtp8C3jTMdgl8aoRz3Q/cP8z2TcDSeK/hBr5AgKnzFlJfvd9tU4bFDQE1fckFW7uBiiPQsPJNcBpxKK/s2NdN4jNQ/OO7e/UBrl5Y7tgaY3VDB+sPN/MPF8zgvRfOZOORFh7bXDvmcbvq2nl5fyPXLXbdNzRh1GtNFGb2sgvcNmFErEj610A9KaAKjkBDbmYN1wDJOSDZiitoTyiCmcKWfmCav9eqGnnffetZveck+050APFPTL3tnKnMLc9PnWEuoFP5jYM3OgtoX34zJjY+bExp48PCkDamHcGUEQw7+k9EwggrDFYYaUXAiiAjEWwrjIxEsMIhIqEgdoqEzw1vWC8u/Kvd1GncWhdIJgRG5Rl40xB846alFOelbvT5enUT04pyALjz95vpClm8vL+RU50h/vjPl4xZx9NnCD562Wy+8OaFjofUOI0W0HGwXU7h1eZxfuCC6F32wXA+CgY2OaYk35QUmDbl/gjFZohCESKfXnLsXgKRXsxwL0YkCLYFVgRsC2lFsCNhZCRMsKMtBe9wfHjSC1fBhVsdB+o+QibuhKdyxZxrFpZx07KRi2PHQ31bD90hi4pJ2Zxs62FWSR4/e7mao6e6B5Uq21nXxp+31I462p1Tmscvb1/OvPKCpGxSBS2g4yDHn3rBsDHosqDLgpNAdU8AGCMPpRH7N0CUP9tyb8ptGwvDNCmfsziWb1MSLSMiYz1yiZQSaVtI28K2LaRtR3+3IlhWBGlFt9uWjR2JKOug5TR+rZ+ncWkN1BinCPalLzYNgSkg22dgiGilHdOIpoMwYr8bsXJjQgyetZFS0hmy6R5Hvczx8r4LE4+Ll1Ky7lAzH3tgI1MLs3n4jov5xENbOHBy5OpAD60/xjvOnTri65+/fkHGiCdoAY2bU51Bntt90m0zhsXAdkV8hOGnvfXGBA4EfCB8YBL9B9EvrGFIDNPGMCTClBhCggGGkAjDRhiSvnwaQshYCTIZ3U40pCDaRsX2E3ZsX0DC7oLJ0KpW8gdDwVGxW7g1Q1Cy8298KrsAIiFEJAiWBdIG26Zv4l9KEFZoqIXVcEWCv5cTLgAAIABJREFU16279KM8Xp943tmRyAuYfOyKs7h2UXnC5/jqX3bx0PpjACyZWkBDe3BU8QTYW9/OnNJcPn7FHB7dWENH7+AlqkvOKknYHhXRAhonv15z2G0TRiTbUHcKaTwIIZBSYEUMnOqTd5XmAU1j7pdOkklfp0kNorsd0d0++j4OXNeJDkNxXoA/3XkJc8sSd9BZU9XUL54Aq/c28uK+xriOXbXzBEW5fj559Vye33OyP8n9dWc75wnsFlpA42TdoVNumzAiOS4tRXrRiUjFSWKfLmd2GrfCWFzCCa/wqxeUcVZp4sXFgxGLzzy8ZdC2nvD4urSt3WH++9n9LJpSwMziHK5cUMa3bhoSqeh5vOcF4gI1zd3sqhu9d+om2aY7siAUTIs3Jgq2z/7MmEDQJIDtwAP5gZUzk+rcPrqxhpYUVUjZd6KD3ICPz1+3wJMd7rHwYAuYXoIRiw//ZgMhhTNm5Lg0hSsM77mgq+gwqXPhTlxSLaBzSvNYMbs44eO31bTy9ad2p9AiWDilYMzQFq+iBXQM1lQ1cahRvQLaA8ky9Ag0bhQUK+1ENHGxU9wEv/WcxPPe2rbk/jWHU97JtDI4UYgHW8D0UtUwuteZCrgVBmGYHlxCV/C7rOXzNCrG6TpJKkegpiF4z/JE0+nDT/5ezVPbj6fMnj7uvGpuys+pClpAx+CW5TOYVTJGXKbLuFUNQhjeE1DVK2doJha+FPqb37xsOrMTdB6qae7mxy9VpcyWPlbMmszS6ekppeYGWkDHoDQ/i1WfvYIv3LAg7ioD6catD9GLAqoiaj5VLjHB+jeTQqnJIPbZa+fxv+85N+Hjv/H0bsJW6m/+jMk5KT+nSmgBjYO8LB+fvnY+j9xxMefPLHLbnCG4petCeNCJyG0DNKMiJ1gYS2538jHJ584o5PPXJ+7l+uyuE7ywtyFpO4YjXYW83UIL6DhYPquYx+68lLePkqrKDXxCe+HGi4q+1FrUT6Oil7ST+DuTF9APrpyVVIjI95/bl7QNI3HDEu+WKosHLaDjxDQE//feZbz/osRzTKaagHDLC9d7Aupc1tHEmVhjLs1AjJb6pI73GYL3rJiR1DkWVDiTm/aqBWXMKkk8oYMX0AKaAD7T4Fs3LeEd501z2xQAslwSUDwYxqJi7cYMjC/XxIlob0wqE1W2P/lOrC+VxUIH8Nk3zXfkvCqhvUASxGcafP+Wc9lb3061y6EupmmSW1yGbYWxw2GscAgr7HzCdE+OQNXTTwUl3T2kPbF6EwIozPFxqiux7+v88rykpm9rW7p5cW9iRTJ8hiA3YJLtN/H7DPyGwGca+AzB8lmTWT5rcsJ2eQUtoEmQ7Tf51e0reOdP1tDem5rC2IlQnzONZwtvGbxRSvyGJMeU5Bo22YYky5RkC5ssAX7Dxo8kIGz8wo7+xMbExkTix8KHjRErHG5iY9oWAhtD2hjSwspOPObMLVRcA7Un2sLfKNgqfkAOklUwiVnFWcMKaLbfYFKOPypOPgO/aWD6jGgJNdPAMAULpiSeML6jN8xdf97BzOJcQGLLWGdORsO9LBsitk3EkkQsSdi2iURsQpZN2JJEbEl7b2TYtu8TV2du7OdAtIAmyezSPJ769OX816q9rN7jTrkzY7geqBCEpSAcgXacGSl+Lr8cPx2OnNspVJzCbTeDLJrdwumxaKz2GoK+Wi19f4GITZ33feYDSm0hzlhQFQP+DURi+kysSISh49++uq59f9tI7Jh3jx0rHj3w775SenZ/0UshfAjDhxB+ED6E8NFXxFbKAFJGn0choqXChJBR71sBPn870xYujr7ngR0LEf2vryzd6W0M/kMM3DL4zg1+j6L/leh1RnguZF+JPIO6falNcQcgllzAjkWFXDCjEMsUBE3oNOGUKWkRktYxRpefPKeSxlCYsoB/3Nd+cN0x1lSnvkhGjt/k8nmlKT+vimgBTQGzS/P45e0r+OKftvOnzbVpv/6wApoGTA+6v6gnn9But3N062/Ses3pixZTt29PWq8ZD0VTptJ6IjnHGifIK0o8v+xo+I8fpfzsLtZOymPo0zn692uyz+CG0kIaQ+Of/rVtydqDzpT1e8+KGRmb+/ZMvOcFojBvSSIPZTK45YRieHDqUUWT3fn41Oz8SEXncJ0qWN9x7DDlkd6Ejl2aH82QFu/oM2JL3mjt5G+NrXxt82Feq0q9gAZ8BndceVbKz6sqegSaQqrGqNbuGG4JKGquKY6GilO4mtNIFXs4DpNoPtyKrPinbTsjFv+8+ygvNkfLMi4Ip77REAL+fOelzJisdurTVKIFNIV8YOVMzqss4j+f2s2hxq60lUBzq5qH8GDWGBULQyhokmsoOwJ18twJfo1mZgewpRx1Cedwd5CH60/x18Y2DvUEE7RwbAwB37/lPM6ZkdmZh85EC2gKKcj2c/FZJfzy9hUUZPt42z1rqGvtcf7CLk7hqpiYYDQsFUc4rnx+Ct4H1JxiBxw1LNFEYpP9JtVdvczPyx42lKUlHOGrVXX9o86BpLqbcl5lEf+wPLmEDl5EC6gDVBZHpzD+5U3z+fc/73D8emO1vwsq8snP8kXd1GXUXd2WEsuO/ovYNpHY73bMdT1sScKWjWVHfw43cnOpjndSRFRsoRU0yT0UvRlOPjcJdqByTJO5edns6ephSf7QadNvHzw+rHhC6u/yu8+fnuIzegMtoA5y6bwShIAbFlfwWlUT3aH0jtf8pqC8IIsDKVqbNQX4Y4HSPp/Bb7pamVRpEjAEAQz8JgSEQZYhKMzyEwnZ+ITAFNEHzRACP+DHwCS6hmoSC7SQp0MwJBKBQMQaLYHAQCKkwCB6LmlJhCT2b0DYRd/f0QiJaMMX+yltCEvnE0yMl4lWA3M0VJnCFcJAGAJhGBiGien3k51fgDCM6GhPCIQQCGH0r9v2vXbmOq5txeIkpURKibRtrEgEKxLGjkTIRjLVb5JvCLqlpC7OdmJ+TgBTCGZkBQZt77Fs/vfICR6qbx7xWCvFnodXLihL6fm8ghZQB5kyKZvVn7+KeeX5NHT0srb6FG09YR5afzRlogYj9yYXVBSw+/jwPdBEsCRYkVgDF7JoBepG2HdmcS7HmrtTdu2BLKjI50BTYvdvAYkHnmucx3ZJQHMKi2iZsQDLtqOJLdLo2r5442ssjv3esfJKHsoePWRmapafHy2ayYrCaJ7ZQv/gZjwsJc3h0RO7WCksQHH1wszPeTsSWkAdxGcazCuPNtjlBdncHJvmeP9FM/nxS1VcvbCM8oJs9ta3E/AZRCzJ+sOneGj9sXGNVkfyXMz2uRel5FB6TSAaw5bwsQrOEAo3pi0VvA8A0nZnVV3aNuE+8XYxOXE8oWElfh8L8rIxjeG/ZAEhePHU6B3nVN7liZDzdiS0gLpAwGfwbzcs7P+7b80U4LrFFbz3wko++KsNNHYGseJo8Uf6zrnZRjqZ3CGZfLbx3M9048YUrnp3IYptuSWgirjDxSGguzp7uHlrFauWL6DYP7QJf+B4Ew2h0UegqRrnF+b4Ob9SvRrJ6UILqILMKy/gmkXlPLzhGDl+g5yArz9RcyCWCzMQW4sUBtQMM1Wa5TMcm0KNBycFNLkRqHrSIV0Z8ah3HyC69ugGdppCzsZCxJmw4UhPiOs37udrc6fxtrIi/LGKLkd7gnz30NiZnFL1bqcV5SSVzN7raAFVlFuWT+fhDcfoCdv0hEPjPr6yONfVKjFmEiWaxiKZUBQVR6CuaJmCt8FNVHFeMqz4i1LUBcPcuecoxf5aLi3KxxSC55ra6I3jGY+k6Ov51qXuZF9TBS2ginLO9OSmRWpb3Bt9gsNTuEk0dgoOQN1ZA1WUaIL79GMrMoVrJDCF3Ry2+Gtj27iOCaXomVs4xZli3F5B58JVmPys/9/em8fJdZV33t/n3KrqfdWuliy1rLY2IyRbXoR3g5Exi81iMGGLcXBgnAy8JJPAJJ+QmYTPDGHeyTJDIAw7bzIYmAAelhiCDWHzvtuSbO2LtUutpdeqe8/7x73dakmt7lruVlXPV59S3zp17zmnTlWd3z3nPOd5yr+/WTmvPcSalM557BtCoZLZtlSOQBMgra2QmCu/lNxZSUxrwCMh3d92NJUeBaaW0BFoSsllDL9zTS9/868vlXV90kIRpXvBSvq67pYcs9sbEPH3l0qw/XRs/94Tu/pDq2exmATkbKD/GAtWrGIsdJoNgn/5WPZsDD90VzEkZcyT1PaZs4lLQPNYMlL5fcPyucneqCeNCmiKueeGpTy87Si/2VZazL5cxrBpf3j7P8shSluQSmaHX9h3/na5bHFX+RlXRPwCeuLgfk4c3D/pa42tyU3L2aRu/Kwdj2eaKHHdQIhviDhSKP/GYfGMZtqb6ltCdAo3xWQdw5fvvIy3XFKam6yu5iwjhdodgUaVsyTkESiJEehUJBkRJaqwYcWQVFzdiUgJRkSVkqtwn/gNy2fXtQUuqICmnsasw39/+xr+5JYVRV+zqDt5ryBR/q4i694TCwuXLgFNy3pg3JgILceLRWLcTpOt0NvJqvn1FXllMlRAq4S7ru7lQ9dfWNRd44ET5QXoDZNI94FG1sEnIxxps8JNy3pg3DgJ7UGdiMS4BpyJ0l1YnaAtWCUYI/zxzct55D++mr9622pes2I2DZOI6ZqFHexM0IHCGJGOQGtsCJo2Z/Jp2dIRN1HuXS6aGKdwM05l77dvtvqVru8V4CqksznH29ct5O3rFjJa8Nh2+BSP7jjGp360iVMjBXKOk3QVIycq/UxB95kKvEJ9Cmga1kCJcR9sJe+3p7OJ1XUWPHsyVECrmFzGsHxuO8vntnPN0pn84Nl9/PzFg0lXK3IKEU0xWqC90Y+b6noWdyz8VBA/NWwD0awjZIwhQ+mepqIkca88EsL+irKKTYGAxtj2lYxAL+5pT0d7JYwKaI2weGYL99ywlPdf1cvLx4f46L1P8fSe0ryThEmUv618Bab3U/HI9vPHTxzDCL4PYhl7Lhg5bb/rhx6146+Jv8kUrCXvBkHLrcW1BEHLXRYM7o7k/ZRPwhbcxiTiVL7e5MCpQEAX12n4srNRAa0xmnIOF85q5X+9dx33/NMTPLrjWCL1iHJLyGhEAloMnoXRM8LBVC42Dl6o4aXCwDiZ04Gg4yahkU0alkDjHNWZCtyFnRhO6LuRMtSIqEaZ3d7IvXev5967r6y5tYpKwpmlkWIjcMSJSXAt3Zhkyk7FGmiMI+9KRqD1HMJsIiqgNYwxwhVLZvDN313PZ37rkpoQ0kwahgkhk0w4s6kxCW5xkBr8jIsmRiOicq2OF81oLtm5S62iAloHNGYdXr96Hv/0gSu5YELw7iiJapCYrdB7ShqRFDouSCouZ90To4CWe6OyblG37iEN0FaoI1obMvzFbRdX9Siumut+PoxN2wpovY4Ck3/PtpCPraxyp6xvWD4r5JpUL2pEVGdcd9EsfvTha3j9//hltMY4EY2qalFAJX0D0ISdqtfeZ1w0o+FtaWp3DFe4QzjWxfG84OEfi/Vo4wQr5ozieAUcr4Bx8xh3FKcwgnHziDuKFEaRwgiSH3sMs2H5TaHVsdpRAa1D+ua0sWHVXP7v0y9HVkZUmpAKbzEhkzZXfpCcY31IcvSb/OcgTU2h5TXHwJrP/EVo+Y1h3QKQCz3fakSncOuUC2dFu48rqmW9WhRQm4KOO00kJd5pWIq2IyOhvXs3omY8undPNBlXISqgdUj/4Cj/69+2RVpGVA7fa9F4IY1GRImS0PRxGm5kxC3QEdJ3vGCjacctjz0USb7VSO31Rsq0NOUcrlsWrSFAVDElnRRu+aiUVBoRJdjOlWzwr4S0uKYL6zselT3v0z/5EaPDQxHlXl2ogNYhDRmHv3/XpfzeDUsjKyNsv7FjpKSPC5fkBz7nUpMNPQ0p+RzCavl8RDexw6dO8tT9P4gk72pDBbSOWRLhOmhUs5I1uASKkD5PRJLQKDBRauy7lY/whmDjLx6MbJapmqjDX4kyxq1reuhoykaStxvRjyst02xhkkZXfqmwqImbGnvLhQjfz+HdO9n30qboCqgSdBtLHbPn2CDHh8LfuD2zNcfhkyOh5wswkvfo6Wxib38NrcHUWMddrYgRjDHjD8dxznguImcciwiO44wfj93cTVzDdUYGwU4c3Nrg87bjT8eOJTh63dFdjBq/axbr7900hQLGc/2H69KUy5E/cRzxXP/huuB5GLeA5PNIfgQzMhxpe2197GHmX7Qi0jLSjgpoHXNqJHwzg5mtOXq6mnh6dzSh1Pb2D3H54q6aEtCfz7yWXcvXTEgpRlEFEYvBn0YSsQh+J+zhr0EXrJCfYInZmM1warSABNc4BnJiyYrFMZYskBFLTiwZ8cgstmSxOGPPsWTwyJigY8disDhYMrgYLGKtH70NsMEOVzthp6sVM64hNpi89jCn1SNIlzHVCc51ETwEa/1r/NLOfK1g/b9ucM5YuQVPuL55F4V8HhvEdXWtxfU88q6H63m+1figLzheWDE5raVt0+MlXzaXn0x7zvyLlvPyi8mOAE8eOZxo+WlABbSOWTW/g3esW8i9j1Uej7K9McNFc9voHxyNTDxPU1vTuEMmx56RaH+KTVnD/Bnt7NzVX3FeRqIzEouKtZkCw0NTjMiiWBqIdLkh+d9AQ4vGBNU10DrnxhWzK85jycwWhvMuj+04xpaDAyHUamrSsF8vTKLuCtubMqxe0MmTIYgnVJ94QnLO8aMyxkqDLcCcJX1JVyFxVEDrnBuXz+bP3rCyIg8/M1tzZwWZjpYk3cxFQdSdYVPGYWA03On6avsEktpbapyIZhYS/gCcTIYla9clW4kUoAJa52Qdw/uv7uUrd17GohnlhTrbdjj6UWctE/V2gGzG0JwNtyOvNqf+SW3LcTJRBQdPtv0vff1tNHdoUG0VUAWAa/pm8eXfvoxsGVHq57Q3RlCjKaiuvntaot4xcvDEMBv3HWduR4ifU5V9BkkJaC2OQBeuWs2r3v6u5CqQIlRAlXGWzGrl0297JS250u6aG7MOHU1qj1Y2EXeGC7ubWTyzlUNhbi2qsnXQpNZAjRPNCDQp/Vx62Xre/Md/hpOJZv94taG9nnIGt63t4eq+mfzjQ7v4zINbGHWnN+l/fOcxls9tY1G34Zm9UVvgVtZ5LO5qYGnzKNZO2FoxYZvFhN0U47vzvAmvAcEGisnrNTE02cQuewTDiYLDyTwcH/E4Plyg4J3e/xc2PZ2NdDTlaGvMsOvoAFsPhTvNHlWwgKhIzIgoIgGNW0JnLerlqne8m9616zAmqvdUfaiAKucws7WBD7+mj8t7u/nz+55n+5GBaYNvb9p/EoC+2a28dPBUHNUsievmO7SQZ8kjn0Py0W4wLwYL0NCCbW7n+NzreYJw15NyGYcX9p0AYNX8dvYdj8axRbWQ3BRuVGITzw2Mk81y3bvfz5rXvr4+3TtOgwqocl7WXziD+/+faxnOu/zpd5/j249PHwewrTFdX6nGjOG98iQNv/p10lU5AwEYGUBGBnCGT0HIAto5YUr9+ZdPhJo3BFscq2gQmtS2j6gENA4/tA0tLbz5jz5Bz/KVkZdVregthTItjVmHj7ymuD1fQ/noQ3OV0nW8btYQDS+mSzzPJdzOffncNmzEU3xp2IdYEontA41KQCPJ9gw2fPDDKp7ToAKqFMXM1gb6ZrdOe96uI4PMbW+IoUbFMSy5pKswPSFq0aUXdLFp/0k2B1PqUVF1kThqzIiIqAMQiNC7Rvd5TkfR3yoRcUTkSRH5fvC8V0QeFpEtInKviN9TiUhD8HxL8PriCXl8PEjfLCIbJqTfHKRtEZGPTUiftAwlfhqzDv/ykWu5/yPXMrvt/AI5MOqyoLu8/aTFUoremBSGCjubsEaLzVmDE2xDimMmoKpIaMQclRGRjdgd1KyFi8jktLudjlJuyz4MbJzw/FPAX1trlwLHgLuC9LuAY0H6XwfnISIrgTuAVcDNwN8HouwAnwFeB6wE3hmcO1UZSgI4Rlg2t40/fcPU0zqnhsN3Ur96ThuvXdjNunkdJbnyq4aJxrDquKqng0e2Hw0pt6mpOm9QSQloRIY3XsQj0EJ+NNL8a4WiPl0RWQC8HvhC8FyAG4FvB6d8FbgtOL41eE7w+quD828FvmGtHbHWbge2AJcHjy3W2m3W2lHgG8Ct05ShJMiGVXO45ILJjV6Wz23jpQOVTx/OaD599+sYYX0hxyufHeKGjaPMdx2asg5ZR2oywHa5xDmrWm1LoEndSkW2BhpWxJjzMDyg3sWKoViTyb8B/ghoC57PAPqttWNDjT1AT3DcA+wGsNYWROR4cH4P8NCEPCdes/us9CumKeMMRORu4G6ACy64oMi3pJRLQ8bh63ddwXee3Muf3/f8+H5GgOG8S7luca/u6WSOydBgYd6mAf5pXgP7To7w+p4ZtD97+gfd98wAfZw5vWSxWCN+qBAD1jGQgczuPA0L7hk/L5d9kIPbXyivghERlu5V4s+4VKpPQJMhqhFo1AJ66S23Rpp/rTCtgIrIG4CD1trHReT66KtUOtbazwOfB1i3bl2VWTdUJy0NGd595SKuu2gW7/7iw+w8Msi6RV08vvNYWfnNbm3g8h0FnIHT+xXfeTCDbcwh26bftykIMhYIEwB/DdBiGOb0mm1DVwrt5kIycDlwIr79rab69rEkU2xEI1DPi26Ne17fMi679a2R5V9LFPPLvQp4k4jswJ9evRH4W6BTRMYEeAGwNzjeCywECF7vAI5MTD/rmvOlH5miDCUlLOxu5mM3L0fE71TL6VI7m7K8XZpxBs7sFJwRj8zxwjnplWCjtl4sg7CMiA6fis9Zgg5AiyOqEahbiE5AL7v1beptqEim/XSttR+31i6w1i7GNwJ6wFr7LuBB4G3Bae8Dvhcc3xc8J3j9AevbvN8H3BFY6fYCfcAjwKNAX2BxmwvKuC+45nxlKCnida+Yxw9+/xpuX7eAd152AcvntnFBCZa4fZ3NNO2OafSUxu0XIY2OTo24XL64O5S8pkP3gRZJVEZEhfAN9cAX/MWvWBtJ3rVIJW5j/hj4hoj8JfAk8MUg/YvA10VkC3AUXxCx1j4vIt8EXgAKwD3WWhdARH4PuB9wgC9Za5+fpgwlZayc387K+e3cvu70ZMIj24/y0W8+xf7jw2esk55NZ2Qhn84ljcG4w6zRk7uP0dmcpX8wH2Ku51Jl8pncpx6RcLv5aD7fuUsvItsYc3SlKqYkAbXW/gz4WXC8Dd+C9uxzhoHbz3P9J4FPTpL+Q+CHk6RPWoZSHVze283P/8MNWGv5xUuH+fufbeHRHeeukT5x5BQXdzaS6Y+20wdSOQK1IY7meme2sCUOX8TVpqC1to3FjWYEumDFxZHkW6uk0KJCqSUcI2Qcww3LZ/OPv3MlC7qazjnnyOAo+xbU86bt8Dr3LQdPsbqnkzkp8gaVBqJ2bXheIhJuN6Ip3LbuGZHkW6uky/O3UtPkMoZPX7uM3/xkB1YEa8ALtl5kjkbTIZxDGtfuQqyTZ+GpPf0I0NPVxN5jQ6HlPZFia5wxgmMEYwSDb2gm4q+hmuCvgP+6gCP+uY4RTHDOxGOQwGDt9LXjeQWvwZlNKgiOCd+hflFENAKNytFBS1c8a+i1ggqoEivr1/dw+NFD7HohHo85Z5NUXMi4sUBPZ3QC2pR1GA5C3I1vaLH+GrPrnV5zLHh2yvXv6Vg8ozmUWKZX9NXWCDQKI6KWrm4Wv/KS0POtZeqjN1FSxYYPXMzKq+YlUnYqrUcjWpeN8p0WPEve9R+jwd+8Zyl44RrsuCG1TVJTuGGub59N2I7qr33XneQaz11iUc6PjkCV2Mk1Zbj+Xcs51T/CrudjHonWyQgUYPvhCN2xxXQjEta23eRMxwQnm8U4DsZkEMdgTMZ/7mQQ4z8XI4hxcDINiHEQcRgZOMaxfTvPm7OTzeG54c0wXLBqdWh51QsqoEoiiBFee9cqvvRHv8QrxNe9pXMKN5r3f+jkCI1Zw3A+fOcRmZjcBlb7CNTr6CPbeuWZaQQOs6aZhZ21cMfUAppxCMtuPdfUREtnV0i51Q9p7E2UOqGhOcucRe2xlplOAY2GJbNaIhFPiG8XS6Fcx8pnkdgUbiU3R3bq76pxwhv/dM6ZH9mWm1pGW0xJlJvuWkVTWza+AlMpoOF37qsXdDCztYGsE41wxNXXhjUC9RIT0EpuYKZuZCdTuYCKGLKNjcxZclHFedUjOoWrJEpbdyOvfPVCHvrutphKTJ8jhSiw1vJwhLFBTUxroGHNFJc7AhURf5uMMTjGQYwENl82eN1gjCBixs/1t9P45TkVdLHG6aBnxZWIZBBxQDIgDoLx/5osXT0WiwPWARys52CtwbMO1gpYwVrBeoLnOXiuwS0Y3IKDmzd4nl/XRWumjvGrTI4KqJI4l968mJaOBn761Y3Tn1whNsIoFmkiamvjuAJqhxW160h+LqsKHb74Wd/w2XoCnsW6grUW6wnW8w2XPM8/r9L3ebySi00PR/a3TX9eBYx9TXKN6jy+HNI4n6XUIcuunEtbd/Q+ONMYjSUKntlzfFKvT2ERxwA0TDulUddh6LBh+Ihh+Khh5Jhh9LgwetKQHxQKQwZ3RPDygnUFsRLbTcL5sBXsny2V/dsqkvq6RQVUSQUiwuJXzoy8HOvF5PEoBczvjHBPXwza0pLLcHQwHI871Tjv4IZkQFUMe1/sj62sWkIFVEkNV9/ex9J1syMtw4Y1Jxgq0XSUrpvG91o8JsQhaDWufNsYBTSN/kWqARVQJTUYI7zmzpWRWuV6dbIGCvD4rn5e0dNOLhP+zzyOmfAwp3Bj2rYaKl6MU7j5keq+2UoKFVAlVTiOYcHy6Bxap3MNNLre/dm9J1jY1YQTsoLE0bWHWecURrGbFi/GEWj/gcFY11xrBRVQJXXMXNAaWd7pgSGlAAAgAElEQVSptMKNeP5s66EBls4Kt03jMSIKr5AUfurT4sW43OAWPPZt1XXQUlEBVVJHb4TGRGkU0Di85IjA7LbwYoR6MQzpuppz9IRkSexW4SponCNQgCfu3xVrebWACqiSOjrnNNM+M5otLTaFc3meRL8Hb9P+kywMcVtLHOtzh06NcPD4cCh+d2PWolCI00c0wM7njnDk5VOxllntqIAqqUNEuGTDokjyTucaaDw8vqufy3vDWV+O4z7k6MAoec+SCcEdYTWOQN1C/N/Vx3+0U9dCS0AFVEklq67pYX5fZ/gZp3AEamO0EH1k+1EuX1x51I04W9EJYS00jinnsCnk419u8FwbX6SAGkBd+SmpJD/qcmDHidDzTeMUro35PvbxncdYPreNTftPxlpuMaxe0EFjxjmjE39qd+XGLdU4hZvEZMnsxW3pDDqfUlRAlVSyb0s/bgShuFLpSCHmDsu1lW8Rieo+ZCTv8cye8N3KFVJ441QM4oCNcSA6f2kEsz41jE7hKqnkxKGhSPJN5RpoAp370YGR2MsshqimWsMKixY3xom3i+4/OBhredWOjkCVVPLySJ7mC9sQz/oLbhbEWvDwBccF63rjkTX8k4LoGtYGUTX8F/0IHP5f8ZrINjZix8/zfO9EZXSwbQuXMHr0ICMDJVguipBrbCLT3ILT2IxkG9jX0FJy2ZVy+NQoTVmHoTLX2SoKFD0FUQ3Gq9UuxnEk1j2sO589wvIr58VYYnWjAqqkjoMnhvmDX7/IyeEKHL8bzjO/8kZoniTZWhyx4+s/1trxIMxj2mrEkjV+to3GcixvmDHfozXjYYCMWGxwvQQXGeP/HXYN/QXD8bw5RyWucLqB6GJ3TkbetVzc08aTu+pj87xXhVa4EK4/4GIYHaqfYAthoAKqpI6thwYqE89yEMFFJpiXTui4gkMPYcxl6JDnpx8pOBwplLCPc5L+MKmufbRQ/tgmqhnRqAxYqnUEKjFP4RYisDuoZXQNVEkdp0bq6y44KZvHF14+SUdTdI77yyGqtqjWNVAn5jjXS9bMirfAKkcFVEkdz+3V4L5xYIFlc9vKujaqtcow/d9OpEr1E2Pi66Kb2rKsunZ+bOXVAiqgSuo4FlIQ5Wohyb796d3HuLinveTrotoNFNWSX1RGT1FjIghFdz5WXdtDJhvzkLfKUQFVUkdDjJ1GGkhy2/pIwfLc3hOs7ulIsBan0TXQMzEhuDEshkzW8IrrFsRSVi1RXz2VUhV88LoLk65C3dGQLa0riHFmMRSq1bmOxNTOq29cSHN7Lp7Cagi1wlVSR71N4XbLELfPHUCCTa0STDhK8J9Yi7Eenphgm4xvLez70A222ojgXynBsRnbGXsGgr+VVhjbXisIllYOsXjuMNnCCBkvj1MYQQIXOBbBGgcrBrEWsR6CB20W8Vz/YV0QMx6azX8vnv9+rH8sngt4iPWnVMXzEDfvP6zr79t9yeUat4C4BX/h0nP948IoQ33r+eJAH+BP9TZmHZqyhgZHcAQc4wBC1giOCI4RMuK/Nq8hx6xFZny+3G9GGRfWs5+PtZA/IrZBOwFntaoN3vFYO8K56639BwYZOpkv9Wvhv88YtrE0tma5ZMMFkZdTi6iAKqnjiTrZmzhG8+AR5v7ma8nWobOL4ZMn8dz0WkA3Pv8A95h/Q9zRSae95y67k/6DXUy+qpznUMT1Ox9dc5vLFtA4/NJe8cZeGprTZY1dLVTZRIxSD/TOjN8zT5LE7Ux+Mgb7j9E+K91bGIz1MOcRzzQT1zpmObR2NbDiVWp5Wy7J/3IV5SwuW9zNX9x2cdLViI04w5lNRUNzK3Mv7Eu6GmUjJp1OAKSCadio10DX3bIYp8T1b+U0OoWrpJL3XLmIQyeG+bsHtiRdlcixKRlTHdj2EgDts+fS1tXN3s0vJFyj0hBJp6ltZdOw0X03mttzLF+vfm8rQW89lNTylkvqw6zeS9nPsDAyzLF9e5OuRs1QySgyyjXQJWtn4dTZlrGw0RGoklr29kcT0ixt2JTtsTCOQ1Ob71zBdV06585j4NhRTh09knDNpkZS6iyhIhGMcFTd1t0YWd71ggqoklr+2483J12FM2jMGlzP4llwS9iZL+K7qAt2SozjR1s700eOiMFkHEQMYgwiwuhQvDEaTx09coZYHtjqT+12zJnL6NAQTjZLfmiIkcGBWOs1LWmdwjWCkzUYIxhHEONvlxETfCccwRg//fTr/jm5xgyzFrWNPx/Lz9/BZLHuIfyNSS7YAtbmsd4onjeKdUfwvDzWLeB5Bax1sV4ez83jFUYZHbgKWJRUs9QEKqBKatl7LD0j0NaGzKRO7gV/T6LI6f1/YyFMx7B2amfm3385A70fOu/rv7/jc6lw5nr8wP7x4/ZZc1InoCLpNCLCWty8F3pcTydTYODQZ8q+PpO9MsTa1Cc6Aa6kljuv6k26CuPY8wiYBVwLBc//69rJdyFWWHjYOYZAGuuUTgGNah3TdSvrvgv58vamKqdRAVVSy29dfgFXL52ZdDWAtHbNSZKudVtI8Qg0Kgf5XmUZu4X0Os2oFlRAldTS0Zzl9nXpsMQtZc2zPkhfe6TViCgqRATjlL8Kl3ajsGpABVRJNSOFdIwq8m6C9UiZlS6cf0o7UVI6Ao1yK4pxyu/Cl6xdF2JN6hMVUCXVzGlPh6m9tZCURzZJ4XSpTeWIPJ1TklG2VNkjUBEWvXJtuJWpQ1RAlVTz6y2Hk67COI0abHgCaRTQNNYp2tVikylPQPsuW0+usSnk2tQfKqBKanlsx1G+/KsdSVdjnDhCS1ULqZzCtemcwo0SJ1N6FBUxhive8o4IalN/qIAqqcUCo0muPZ5NQqJhUzqySh8p+q7EhJQR2fzi61/DnF4NWh8GKqBKarlscTdvWduTdDXGSWzQlcbRXgqxobsqSD+lroG2dHZxzbvujKg29Yd6IlJSzSff/Ar29g9x8OQIrmex2HE9MYH3Hw/wPBu87k8vCjJ+rmf91wqexbP2tAu94LVi7WFUxk4TR6DnkrHpNCKKEmNK68JvvuejNLW2RVSb+kMFVEk1TTmHm1bO4S9/sDHScoxAU9Yh4wiDoy4F19LRnCVjhJacAwgHTg5HWodqIo0Cam39jUBLWQNdetl6Fq9Wy9swUQFVUs/8zuitBT0LA6NndsD9g76rszE74JyGfjpNCgWUlApolDPwxileQC+55U3RVaRO0R5BST1OSqxfvST2PqZ2/TMdn8lE6nINNJMr6rxZFyxmwYqLI65N/aECqqSe7z+zL+kqJIaT1hBd6dNP8NInoE5GIgmzJuIxY94hCiPFhbq7YPXaVE67VzsqoEqqeXTHUb7/zMtJVwOAnhimks8mfZ2esGDFxZw4dDDpipyDk+tOugrnMLu3g30vHQ8930zOZe8LX+fInq1FnT/rgsWh10HRNVAl5aya305rLsPJSWJxxo2XgB1uGhwWzOtbhue6ZBoaOXFwP3s2Ppd0lSbF2pakqzAJEX1+trQbq9auGdHUo85RAVVSTXMuww8/fA2PbD/KX92/iQMnRpKrTAJalpT/3THm9S1n30ubkq1EsaRtsA6RfWe8En1GtM5QAY0CFVAl9SzsbmZhdzO3vGIeX/n1Dv6/h3ayt38o9nokMRZMcgbXyeY4tm9PchUoleQH67HhlRhMu7GlNaKa1De6BqpUDU05hw9dfyEP/OF1vO3S+OOEeglMp5oERWHu0j6GT51KrgIlUz8K6pUQTDuTzdHS2RVhbeoXFVCl6mjIOHz6bau554Z4/XkmshyZ0Ai0Z8XF7N34fDKFl0k9BdQWkaL94C69fH3EtalfVECVqkRE+MhrLmLtBZ2xlZmEgCYlCu5ogmvNZZNCZ/JROlEoUkAvvv6m6CpR5+gaqFK1ZB3DJ297Bbf83S9iKc9NQEGzCd3iOjmH5a/uBfFF3IogY/sZgzT/9tsixiLGA+MiAmJcEBtcZ/1zACsWEcvWfzWcOnostLqKGHLNTZgUeoqK8htT1AhUhHkXLYuwFvWNCqhS1ayc385rV87hxy8ciLwst1TTxxBI6ge6d+Mm1ly1G5zK1kBlkuMlN1yNN9rjq4uAOK4vuk4eMf4DpwB4IAUQD8QXZawN0lwwwyAjWDOMCIzuvB12VFTd8In0pmv6+f2uufM1cHaEqIAqVc+li7r46aaDuBG72iu4CRgRJTSoWrh6CTjROPA3M38Z6trRmIyISZ8noijvuYoZgXb3xG9sV0+kb85DUUrkd6+7kDeunhd5OfkEfOFKAtPG81csZsblP4u93EoRk7yzjbOxEX5njONMe05+OP7tXvWEjkCVqufgyWF+9uKhyMspuPFP4cbpR7+1u4veaxsws3+UTqcE05DGEWiUnqScIoJp7928kcHj/TR3xGdsV0+ogCpVzcBIgfd/5dHx0GO1R7Qj0MbWFpbf5jsktw3PIpmTkZYXLenbxhLpCDQ7fSgzN5/noe/cy42//buR1aOeUQFVqpaBkQJ3ffVRntt7IpbykohmFrUz+SVXzYSWf/HLirSkGEhh5JooQ+A5RcYCPXk4+tmZemVaARWRRuDfgIbg/G9baz8hIr3AN4AZwOPAe6y1oyLSAHwNuBQ4ArzDWrsjyOvjwF2AC/x7a+39QfrNwN8CDvAFa+1/DdInLSOk965UMceH8nzga4/xyPajsZV5PiOl2a1Z1rWPcMBtZNAVRjzY1T9K4Tzn5zKGRd3NNDc4HB/Mc2qkgLVwZCDer/aitUvILvpBrGVGiqRvH6hXguGZMYI4gZOE4G7mjGPjn5PJgJMBm5uJbXGQrAOOwTqCNQIGrAjWAEbY3zbIidETtOfaw3+DdU4xI9AR4EZr7SkRyQK/FJEfAR8F/tpa+w0R+Ry+MH42+HvMWrtURO4APgW8Q0RWAncAq4D5wL+KyEVBGZ8BbgL2AI+KyH3W2heCaycrQ6lzPvj1x2MVz8l4xewGrtn6f3C27wZgol8k29LFyKJXMtLYwVCunW0NCxAj9A/m2XlkgJcOnrk9pDlr6G7J0dWcZeuhgfH0ckaF2cZGCqOj2ElMQHNNzTiNTcj8i+he/zQ2fcuGZdPYvoWZPRbruWA9LBbw/GPrYT0P6+XxvALWy/vrk+Ovuf7Duljrt5uIwVrrX+MWgutcrJs/45y3Xp6ho/8JxLo43jDijWBsHk+yeJLDXtCIK414phFLBsRivGEcbwjHHcDxBjHeUMl+j+9YdR3Pn9xe1LlmUy8fWP2B0gpQpmVaAbX+KvjYrz0bPCxwI/BbQfpXgT/HF7dbg2OAbwP/U/x5qFuBb1hrR4DtIrIFuDw4b4u1dhuAiHwDuFVENk5RhlLH7D46yOO7wtuIXyozmrOs6LCseepLyODk08cycIzGF35GI9CVyfCthVN3XoN5j8H8KEcHRlm7sIMXD5xiYNQtyRNRQ2s7u5Zt4JtHZjA757K0aRSLZUlmgNnuMQZMM18/MZ9+14FB+FLmGX8uqEaQzMvsee6Z2MttRmge2XFOumMDC9hi2riMOyVTguI+czj+dqkHiloDFREHfwp1Kf5ocSvQb60dsxvfA/QExz3AbgBrbUFEjuNPwfYAD03IduI1u89KvyK45nxlKHXMf/nRRkYL8U3XrZjZQIOxjHjC+uGNNL/4S2RkYPoLA4p1uTbGk7uPs2xOK0cGRpF8ccZRTjbLz5a8hceONoHAwXyGg3n/5/0bWoDZk1xV9aueZ2BtMrvyTEJ3IaV8ejctUnd+UVCUgFprXWCNiHQC3wGWR1qrEhGRu4G7AS644IKEa6NEydCoy4+fj97r0BgNGcOrn/gM4pZv5Vus0++JbD5wiqwjrOjpJLezmYZFy9nfsoB2d5DGvc+RP3WCbGs7A4f242Sz/HjZe3jhZEOpNSu5XqkmoSVQSaMP3gm0Zlu5YeENSVejJinJCtda2y8iDwLrgU4RyQQjxAXA3uC0vcBCYI+IZIAOfGOisfQxJl4zWfqRKco4u16fBz4PsG7duvSZ4imh8f/+ePN5jXPC5qIZDbz60AMViSf462TlkHctv9g1wN7ed7JtMAsjgeC19kErYC1/OG8jucPbeWGgVPFUwiLtUWA2LN5AW64t6WrUJNP+skVkVjDyRESa8I19NgIPAm8LTnsf8L3g+L7gOcHrDwTrqPcBd4hIQ2Bd2wc8AjwK9IlIr4jk8A2N7guuOV8ZSp1xdGCU93/lUb7wy+KMJipldmuW1278Kpmdla8dlTMCnci2odzkkbVF+G8nVvC97ldXlH+tYEuIkRluwckIaLGl7h/cH2k96pliftnzgAdF5Bl8sfuJtfb7wB8DHw2MgWYAXwzO/yIwI0j/KPAxAGvt88A3gReAfwHusda6wejy94D78YX5m8G5TFGGUmc8vO0ID2w6GEtZWUe4vuUoMhCOoVKlAjp15sLzA43R5V9FJCWgSY1Aiw3w/qu9v+IXe+KJWFRvFGOF+wywdpL0bZy2op2YPgzcfp68Pgl8cpL0HwI/LLYMpf5obYze58eslizLOyyv2P5jnIc2h5dxmVO4SmkkNgKtAv7g53/At974LRa1L0q6KjWFeiJSqoLOplxkeYvAXU2baXrugcjKUOIg3WuRYVOKn+ShwhD3bb2P31/7+9FVqA7RW2OlKugfis5Lz5vnDNP0fITimdAaWb0hmXSvRYaNU2L3/eDuB8lXaBCnnImOQJWq4MRQ+KGq1sxu4Krt38X8ZkfoeU8kyogcymmMSaqdE9p/WqLrooWtC8kY7fLDRFtTqQoOnhwu+1oj0N2c5dSIyxtmnmTmvmcww6dwHt4SYg2nwKZzn6CttX2gCc2n2YSaMVPC57e4fTGfeNUnIg9OUG+ogCpVQXdL+Wugd8w4Qvej34RMA/LSSIi1Kg4df8ZFUi2djCg5JYjhPWvuobuxO8La1CcqoEpVMLO1PEcBCztydD91n9/FFeIXT2BSp+5K+EhSQ8GEyJQw5L665+oIa1K/qBGRUhUsn1u6J5V1c3K8adMXkHz507+hoGugsZBUM0tCBTsljHx3ntgZYU3qFx2BKlXBjNYGFnY3sfvo0LTnvma+0Lf7FzgPPVtrq3zKFOh9yvl5eP/DrJq5Kulq1Bw6AlWqhmv7Zk35umOE18wXlj/6ZTK7VDynI+0+XKsFLyFHGaV8ep9/5vNsPLIxsrrUKzoCVaqGa/pm8uiWA7yqsAWxHqeaZ/KLgS6ODhW4vfMgXU98F9mawn1uavkYC0k1c1LWzF4JEjqQH+C3/+W3+fR1n+baBddGWKv6QgVUqRrWdMNNT/49xvWdKrTi+4z0xGBSulUEai5oWHpJTECrYyJvsDDIH/78D/nBm3/ArOapZ3OU4qiOT15RgK7ODpqbzrXGTbN4QsTO5CuixqS9zkag5ZQ6VBjiu1u+G3pd6pW0/rIV5Rwampu5+Iabkq5G6aR0CrfWVkCTW9NN5vM1ZZa748SOcCtSx6iAKlXF6NBg0lWoGdIp6xUgCYUVS6gbLcUT0UTac+0h16R+UQFVqoqTRw4nXYXSSen+inTWqvpIagq3lH2gE8l7KTS0q1JUQJWqomd59e1lS68z+Rr7+dfZGmi2zKUBDa4dHjX2C1JqHa8QflSWyEmtgNbWJG5S7yapKdxsmV8r3cYSHiqgSlXRMWdu0lUoGZtSK2Fba75jE1oDTaodyx2BOsYJuSb1i+4DVaqKldfcQFNbOw/98728vPmFpKtTFGkdgNqy7ThTSkIC6iZlRFTm293avzXcitQxOgJVqo7eNZfyjj//L1z2prcmXZXiSKuC1hiSlBWuTcoKtzwe2/8YR4ePhlqXekUFVKlKjHG49l13cunrb026KtOT0n2gtbYGmtgINKEp3HInYgu2wE93/TTUutQrOoWrVDWZXHlxQuMktfpZawKaEHEKqHUaINMAJkeLbWJ+Uw9ZkyUrObKSJUMGRzI4ZDAYZPzf2DODYHhs07O89cK3YhwdQ1WCCqhS1TS1dSRdharlqy+8m+2H3zC+H9TzoOAZ8h64nu/Xx3qCa09v1XDEN4hyDDjGkjVgsONabMQi4s9a28A3kOtCwQrWwkjB8PkNf4nnDWIkB+Ig4iA4gEHEAbIIGfydjg7+WCsDNoeQBZvBWgdsBqwD1mBxyB+fyYKVBnCx1sV6nu9GUXzR8NMKWK+A5+ax1vMfnn8+uIjI+AMRRIx/DIgRIHjdjImmcKAwTG7W6/z3awXPgkfw1xo8C3lXKHh+G7sWvEB0W7pWMDzYhuc5uNah4DmM5LOMFrIMjeYYzWfIFzKM5h1GCxnOuOnZDm+q4PN/cdYBlq+fV0EOigqoUtXkmpuSrsK0pHUJ9GQ+x94TpY7gKx+xPP29T+MMuBXnEyazFuxl97P3lnXt3grKXbx2Jft3tFaQQ/n8/H9vZubCVmYuKD1YveKj43elqnHzVbAvNK0Kmli1UtoeCWATbIvCqMcPP/ssp46NJFaHakdHoEpVk971xYlYFnafHikbEQQQEYycTjMiiDCeJhPenLWWgmeDaUGLFxz71xIc2+Bc8ADrWf+vtX6atThG/EeiDVcVH1osJP39PXlkmG9/6jFuvvti5i7R5ZBSUQFVqpqWzu6kqzAt1sLuo0NJV+McultyyRScynmvhEaCKZidGOgf4Z8//Thv/oNLmLe0M+nqVBWp/CorSrHM61uWdBWmxXop9USUVLlJD7smJanWSMd3o/eVs5jdq1FaSkUFVKlqWjq7aGhuSboa05D8KCNVpFE/E3NEnw4uf2Mvjm5pKRltMaXquXDdFUlXYUrSG40lIdI4Ak3sI0p+BCoCXfPSfhOaTlRAlaqnc26697JJGgWD5AaCNpW+zJMagia/ncda35hIKR0VUKXqueK2t9PY2pbOkQ2QzjnL5ITdOulrD5uQqtuUBLceHkhHPaoNFVCl6jGOwwf/4etc/Y73JF2V85DOKdyk7jesUQEdw3PTsQezpSP9LjHTiG5jUWoCJ5NhzYbXkx8ZZvOvf0H/gX1JV2kC6RMMSHIKN33tUc8C2tSWpbkjoS1NVY6OQJWaoaG5havveC/v+dTfMmtRb9LVOU369AJIcAo3hb2O9RJqC5v81Omlr1uMSeGsQDWQwq+yolRGrqmZDR/8MDMXLkq6KqkmsRFoCjtrz0tmBOrmkzXeyTU66lC+AnQKV6lJ5ixZynv/6n/w1E9+yLbHH6GhpZWXHv4Vnhu/1WNqrXATqpaXwvawSQlooTQBdbIGJ2NwMoKTMRjHjwzjR4gBY05HixE/eMzp758wvhwvwX9zFrfT0KQyUC7ackrNIsawdsMbWLvhDQDs27KZh7/zLXY++ySFkTjXntInGJCgJ6IUznsV8slUyptGQOf0tnNw14nx3S5u3sPNh7d3dNmVOvqsBBVQpW6Yt3QZt/2HP8UtFNj57JNse/xR9r20mYM7t6XCJ2ncJDeFm1DBU2C9ZCpl7XRiaCPdKjq/T33fVoIKqFJ3OJkMS9ZexpK1lwEwMjjAyy9uYt9Lm9j13NPs3bwxVEEVk0LFIElfuAkVPAVJ3T957nRGRNE1lgh0zmmOLP96QAVUqXsamlvoXXMpvWsu5VW3v4vBE8fZ/fyzPPfgj9nx9BMV559WAU0sAEkK10C9pEag3tTxbKNsqub2nFrfVkhKf9mKkhzN7R0sW381b/n4fwrFz25ajYiSIo2T5V4hma7QdUenfF0iFDjPS+MnUV3oCFRRzoOIsPKaG9j62MNJV6WmSOMUrpvQPlAvf5LuxS1Y1w+Sbj1wXQ/PtXgFj/3bjkdW9tDJPCePDtPW3RhZGbWOCqiiTMHSy9fT0NLCyMBA2XmkNRqLroGexnMTElA3z9G95X+3KmXPpqOseNX8xMqvdnQKV1GmwBiHOUv6KsojtWugCZHG2wlJqCt0C8l6IkrpvV3VoL9sRZmGNTfdUtH1Ivozm4ib1uZIYK3azecTnaEYGZjaiEmZmrR+lRUlNfRd8So2fOgjSVejZrApdSwhCdXLZJILqm1S6Ni/mlABVZQiWHXtjTjZbJlXp3OeLDFHCiltj6QaJJNJrj3UErcyVEAVpQjEGNa94c3lXesk42c1raTRFy6ASWitWkxyI9DRYZ3CrQQVUEUpkqve8R6yjU0lX2eMCuhEPEnnqCcpYy8nwSncU0eSjQZT7aiAKkqRiAjz+paVfp2Tzt1iyUVjSabc6TAJfU5ikruh2LP5GFanccsmnb9sRUkpXXPns+vZp0q6RlI6Ak1sH2hC5U5HUgJqipzCNUEIM8cxmIxgHME4Jghh5nstMo74z8fCmgUhzeBMj1h+un/+8ECeprZc+G+sDpC0bvIul3Xr1tnHHnss6WooNUp+ZBjreVhr8TwPr1DA81ysZ7Gei+d5wbGH5xbwPI8ChkOZLlzP4suH37EZEc721Db2cxTxrUJFwDGCCTo/iz3DWlQE8q7HaMF/DBc8hvMuIwWPodECrsd4GSJ+/mO/+MasYTjv4QWFWnt6VCoIFou14FnfS07Bs+RdS8H1/PzzLsN5l7zr4Xrgeh5573QBnrW4nsWzvjMJz1pcC5e3NrNg7yjWw/e441m/zaz1n7un+yRr/TpM3LBog2YcM0Y6LRRyOuRlUOcxDz9jXn6sZ88ocyLifonRwVNkG5vINTWTa2wk29hIrqmZbEMjmVwOJ5vFcTKYjIMYgzEOxvEffhxOg3Ey42nGGMQY/7pMdvy5cRxMJouTyZBp6MHJNGCcsVifBicrOBmHTM6QCWKARunWT5kaEXncWrvu7HQdgSpKCWQbynN71hNyPZTKsWPCav2bEjHXYdTgSykBFVBFUeoSMYKjozqlAtSISFEURVHKQAVUURRFUcpABVRRFEVRykAFVFEURVHKQAVUURRFUcpABVRRFEVRykAFVFEURVHKQAVUURRFUcpABVRRFEVRykAFVFEURVHKQAVUURRFUcpABURrg9oAAAaBSURBVFRRFEVRykAFVFEURVHKQAVUURRFUcpgWgEVkYUi8qCIvCAiz4vIh4P0bhH5iYi8FPztCtJFRP5ORLaIyDMicsmEvN4XnP+SiLxvQvqlIvJscM3fSRA6/XxlKIqiKErSFDMCLQB/YK1dCVwJ3CMiK4GPAT+11vYBPw2eA7wO6AsedwOfBV8MgU8AVwCXA5+YIIifBT4w4bqbg/TzlaEoiqIoiTKtgFpr91lrnwiOTwIbgR7gVuCrwWlfBW4Ljm8FvmZ9HgI6RWQesAH4ibX2qLX2GPAT4ObgtXZr7UPWWgt87ay8JitDURRFURKlpDVQEVkMrAUeBuZYa/cFL+0H5gTHPcDuCZftCdKmSt8zSTpTlKEoiqIoiZIp9kQRaQX+D/ARa+2JYJkSAGutFREbQf2KKkNE7safLgY4JSKbI6jCTOBwBPnWOtpu5aHtVjraZuWh7TY9iyZLLEpARSSLL57/aK395yD5gIjMs9buC6ZhDwbpe4GFEy5fEKTtBa4/K/1nQfqCSc6fqowzsNZ+Hvh8Me+lXETkMWvtuijLqEW03cpD2610tM3KQ9utfIqxwhXgi8BGa+1/n/DSfcCYJe37gO9NSH9vYI17JXA8mIa9H3itiHQFxkOvBe4PXjshIlcGZb33rLwmK0NRFEVREqWYEehVwHuAZ0XkqSDtPwL/FfimiNwF7ATeHrz2Q+AWYAswCNwJYK09KiJ/ATwanPefrbVHg+N/B3wFaAJ+FDyYogxFURRFSRTxDV+V6RCRu4OpYqUEtN3KQ9utdLTNykPbrXxUQBVFURSlDNSVn6IoiqKUQc0LqIg0isgjIvJ04IrwPwXpvSLycOA+8F4RyQXpDcHzLcHriyfk9fEgfbOIbJiQfnOQtkVEPjYhfdIyqgkRcUTkSRH5fvBc220aRGRH4JryKRF5LEhT15dTICKdIvJtEdkkIhtFZL222dSIyLLgOzb2OCEiH9F2ixFrbU0/AAFag+MsvhOIK4FvAncE6Z8DPhQc/zvgc8HxHcC9wfFK4GmgAegFtgJO8NgKLAFywTkrg2smLaOaHsBHgX8Cvj/Ve9J2O6PNdgAzz0r7K+BjwfHHgE8Fx7fgG81J8L18OEjvBrYFf7uC467gtUeCcyW49nVTlVEND3xPY78THOeATm2zktrPwXc2s0jbLcZ2T7oCMX/JmoEn8P3xHgYyQfp6/C014G+3WR8cZ4LzBPg48PEJed0fXDd+bZD+8eAh5yujWh74e3J/CtwIfH+q96Ttdka77eBcAd0MzAuO5wGbg+N/AN559nnAO4F/mJD+D0HaPGDThPTx885XRtofQAewncAmQ9usrDZ8LfArbbd4HzU/hQvj05BP4Tti+An+yKffWlsITpnoPnDc5WDw+nFgBqW7KJwxRRnVwt8AfwR4wfOp3pO222ks8GMReVx8L1mgri+nohc4BHxZ/OWCL4hIC9pmpXAH8L+DY223mKgLAbXWutbaNfgjqsuB5QlXKfWIyBuAg9bax5OuSxVytbX2EvzIRPeIyLUTX7T+bXvkri+jLiNEMsAlwGettWuBAc6KvKRtdn4CG4E3Ad86+zVtt2ipCwEdw1rbDzyIPy3YKSJjjiQmug8cd0UYvN4BHGFqF4WTpR+Zooxq4CrgTSKyA/gG/jTu36LtNi3W2r3B34PAd/Bv2g6I744SKd715fnSp3R9OUkZaWcPsMda+3Dw/Nv4gqptVhyvA56w1h4Inmu7xUTNC6iIzBKRzuC4CbgJPyTbg8DbgtPOdkU4ZoX2NuCB4A7rPuAO8a1Ne/Hjlj6C71mpT3zL0Rz+VMp9wTXnKyP1WGs/bq1dYK1djP+eHrDWvgtttykRkRYRaRs7xl+beg51fXlerLX7gd0isixIejXwAtpmxfJOTk/fgrZbfCS9CBv1A1gNPAk8g9+R/VmQvgS/I9+CP/XREKQ3Bs+3BK8vmZDXn+Cvn24msEYL0m8BXgxe+5MJ6ZOWUW0P/CAAY1a42m5Tt9USfIvip4Hnx94X/truT4GXgH8FuoN0AT4TtMGzwLoJeb0/aIMtwJ0T0tcF3+WtwP/ktEOUScuohgewBngs+J1+F98aVNts+nZrwZ+16ZiQpu0W00M9ESmKoihKGdT8FK6iKIqiRIEKqKIoiqKUgQqooiiKopSBCqiiKIqilIEKqKIoiqKUgQqooiiKopSBCqiiKIqilIEKqKIoiqKUwf8P+VeoL0i5ZKIAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
],
"text/plain": [
"NAME Adams Ashland Barron Bayfield Brown Buffalo Burnett \\\n",
"DATE \n",
"2020-03-16 0.0 0.0 0.0 0.00 0.00 0.0 0.0 \n",
"2020-03-17 0.0 0.0 0.0 0.00 0.00 0.0 0.0 \n",
"2020-03-18 0.0 0.0 0.0 0.00 0.33 0.0 0.0 \n",
"2020-03-19 0.0 0.0 0.0 0.25 0.50 0.0 0.0 \n",
"2020-03-20 0.0 0.0 0.0 0.20 0.40 0.0 0.0 \n",
"\n",
"NAME Calumet Chippewa Clark ... Vernon Vilas Walworth Washburn \\\n",
"DATE ... \n",
"2020-03-16 0.00 0.0 0.0 ... 0.0 0.0 0.0 0.0 \n",
"2020-03-17 0.00 0.0 0.0 ... 0.0 0.0 0.0 0.0 \n",
"2020-03-18 0.00 0.0 0.0 ... 0.0 0.0 0.0 0.0 \n",
"2020-03-19 0.25 0.0 0.0 ... 0.0 0.0 0.5 0.0 \n",
"2020-03-20 0.20 0.2 0.0 ... 0.0 0.0 0.6 0.0 \n",
"\n",
"NAME Washington Waukesha Waupaca Waushara Winnebago Wood \n",
"DATE \n",
"2020-03-16 0.00 0.00 0.0 0.0 2.00 1.00 \n",
"2020-03-17 0.00 0.50 0.0 0.0 1.00 0.50 \n",
"2020-03-18 0.67 0.67 0.0 0.0 0.67 0.33 \n",
"2020-03-19 0.50 2.25 0.0 0.0 1.00 0.25 \n",
"2020-03-20 0.60 2.40 0.0 0.0 0.80 0.20 \n",
"\n",
"[5 rows x 72 columns]"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dft = df.T\n",
"print(\"shape:\", dft.shape)\n",
"print(\"variance:\", dft.values.var())\n",
"dft.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have 352 rows, but of course there will be similarities across them. Can we approximate the table by finding a few key rows (called components), such that every row in approximately some weighted combination of those component rows? In particular, how much of the 3230.9 variance can be explained by these few rows?\n",
"\n",
"A principal component analysis (PCA) can help us find such components and tell us how much of the variance those components can explain.\n",
"\n",
"Passing 0.95 tells the PCA that we want it to find however many component rows are necessary to explain 95% of the variance in our original table (alternatively, we could pass an integer to require an exact number of components):"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PCA(copy=True, iterated_power='auto', n_components=0.95, random_state=None,\n",
" svd_solver='auto', tol=0.0, whiten=False)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca = PCA(0.95)\n",
"pca.fit(dft)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`explained_variance_ratio_` tells us how many components are necessary, and how much variance they each explain. In this case, one component row can explain 90.6% of the variance, and a second component row can explain an additional 4.5%."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.90638618, 0.04510856])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pca.explained_variance_ratio_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see what those two rows look like in a table. `components_` gives us a numpy array, but we can use the same column names from our original DataFrame."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
NAME
\n",
"
Adams
\n",
"
Ashland
\n",
"
Barron
\n",
"
Bayfield
\n",
"
Brown
\n",
"
Buffalo
\n",
"
Burnett
\n",
"
Calumet
\n",
"
Chippewa
\n",
"
Clark
\n",
"
...
\n",
"
Vernon
\n",
"
Vilas
\n",
"
Walworth
\n",
"
Washburn
\n",
"
Washington
\n",
"
Waukesha
\n",
"
Waupaca
\n",
"
Waushara
\n",
"
Winnebago
\n",
"
Wood
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.014671
\n",
"
0.012164
\n",
"
0.063659
\n",
"
0.012197
\n",
"
0.186703
\n",
"
0.013653
\n",
"
0.012156
\n",
"
0.039235
\n",
"
0.081403
\n",
"
0.034584
\n",
"
...
\n",
"
0.017245
\n",
"
0.018460
\n",
"
0.065934
\n",
"
0.013056
\n",
"
0.120046
\n",
"
0.379391
\n",
"
0.034502
\n",
"
0.018779
\n",
"
0.124864
\n",
"
0.068892
\n",
"
\n",
"
\n",
"
1
\n",
"
0.010148
\n",
"
0.000019
\n",
"
-0.026789
\n",
"
-0.001165
\n",
"
0.488611
\n",
"
-0.003380
\n",
"
0.000764
\n",
"
0.125926
\n",
"
0.006220
\n",
"
0.003825
\n",
"
...
\n",
"
0.003234
\n",
"
0.010982
\n",
"
-0.023864
\n",
"
-0.011774
\n",
"
0.010618
\n",
"
-0.225074
\n",
"
0.099319
\n",
"
0.058002
\n",
"
0.452854
\n",
"
-0.009747
\n",
"
\n",
" \n",
"
\n",
"
2 rows × 72 columns
\n",
"
"
],
"text/plain": [
"NAME Adams Ashland Barron Bayfield Brown Buffalo Burnett \\\n",
"0 0.014671 0.012164 0.063659 0.012197 0.186703 0.013653 0.012156 \n",
"1 0.010148 0.000019 -0.026789 -0.001165 0.488611 -0.003380 0.000764 \n",
"\n",
"NAME Calumet Chippewa Clark ... Vernon Vilas Walworth \\\n",
"0 0.039235 0.081403 0.034584 ... 0.017245 0.018460 0.065934 \n",
"1 0.125926 0.006220 0.003825 ... 0.003234 0.010982 -0.023864 \n",
"\n",
"NAME Washburn Washington Waukesha Waupaca Waushara Winnebago Wood \n",
"0 0.013056 0.120046 0.379391 0.034502 0.018779 0.124864 0.068892 \n",
"1 -0.011774 0.010618 -0.225074 0.099319 0.058002 0.452854 -0.009747 \n",
"\n",
"[2 rows x 72 columns]"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"components = pd.DataFrame(pca.components_, columns=dft.columns)\n",
"components"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's visualize those components on a map. We'll use the `coolwarm` colormap (https://matplotlib.org/stable/gallery/color/colormap_reference.html). We want blues to be negative, reds to be positive, and gray to be about zero. This means we'll need to pass in a vmin and vmax that are equally distance from zero. Let's first find a value for this."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.6830682912657138"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cap = components.values.max()\n",
"cap"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's add the rows of the components table as columns to the GeoDataFrame (`wi`)."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(ncols=2, figsize=(10,5))\n",
"wi.plot(column=wi[\"pc0\"], cmap=\"coolwarm\", figsize=(8,8), ax=axes[0], vmin=-cap, vmax=cap)\n",
"wi.plot(column=wi[\"pc1\"], cmap=\"coolwarm\", figsize=(8,8), ax=axes[1], vmin=-cap, vmax=cap)\n",
"axes[0].set_title(\"Component 0\")\n",
"axes[1].set_title(\"Component 1\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although all 72 counties could in theory vary independently, we can take different combinations of just those two maps (added to a base map of the average per county) and explain 95% of all variance.\n",
"\n",
"Let's see what that base map looks like (similar to component 0, it turns out):"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"NAME\n",
"Adams 4.484290\n",
"Ashland 3.334858\n",
"Barron 15.192472\n",
"Bayfield 3.024574\n",
"Brown 85.731733\n",
" ... \n",
"Waukesha 115.279972\n",
"Waupaca 13.570511\n",
"Waushara 5.961534\n",
"Winnebago 48.406193\n",
"Wood 18.998011\n",
"Length: 72, dtype: float64"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean = pd.Series(pca.mean_, index=dft.columns)\n",
"mean"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.0, 'mean')"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAAEeCAYAAADPbQ0xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO29e7Q0WVXg+dsRmRmZed/f+1VFoZTQ6FIe34Jyads0tFAgbTGOOmC3lDTL6l5it65xWkFnhmkfa+GaNb5aZKYW0IKjlgw2TS0txAJhnF4zhVTJGwQ+HsX3ftz3vfmO2PNHnLw3b958RWbkI+Ke3/riu5knIs45mRmxY5999t5HVBWLxWKxxI8z7Q5YLBZLWrEC1mKxWMaEFbAWi8UyJqyAtVgsljFhBazFYrGMicy0O2CxWI42L3TmdEv9yOddovphVb1/DF2KDStgLRbLVNlSn9/JPCPyea9ufOXEGLoTK1bAWiyW6SIgWYl+XiP+rsSNFbAWi2WqiAhOZggBmwDsJJclNkTkmyLy70XksyKyKyLvEpHTIvIhEdkWkY+IyIo59j4R+X9FZENEPiMiL2mp5w0i8iVzztdF5F+37HuJiFwRkV8QkVsicl1E3jCFj2uJCwHJOpG3JJCMXlqSxH8L/CDwHcA/Bz4E/DJwkvB6+3cich74S+DXgWPA/wD8uYicNHXcAl4NLAJvAH5bRF7Q0sYZYAk4D7wReHtTcFsSiICTkchbErAC1hI3/1FVb6rqVeD/AT6hqp9S1QrwAeD5wL8EHlPVx1Q1UNXHgSeBVwGo6l+q6tc05P8G/hr4xy1t1IFfVdW6qj4G7ADPntxHtMSKscFG3ZKAtcFa4uZmy+tyh/fzwDOAHxORf96yLwt8DEBEXgm8lVALdoAi8LmWY1dVtXWKo2TqtSSQNNtgrYC1TIPLwB+p6k+37xARD/hz4PXAB1W1LiL/BUjnHWgZ3osgAVgBa5kG/yfwSRF5BfARQu31PuASsAl4wG2gYbTZlwOfn1JfLePG2GDTiLXBWiaOql4GHiCc/LpNqNH+e8BR1W3g3wHvA9aBnwAenVJXLRNAAHEl8pYExCbctlgs0+Q5xaK+897nRD7vH3/2U0+p6sUxdCk2rInAYrFMGUGcZGikUbEmAovFYhkTVoO1WCzTRUDcdOp6VsBaLJapIoCTkEmrqFgBa7FYpouQWhts6gTsiRMn9J577pl2NyyWI8dTTz11R1VP9j+yHbEabFK45557ePLJJ6fdDYvlyCEiTw93Honxa41K6gSsxWJJHuLYSS6LxWKJH2uDtVgslnFhbbAWi8UyFsRqsBaLxTI+rA3WYrFYxoHVYC0Wi2VcWBusxWKxjAVrg7UcOer1Gtlsbu+9qqKqOC22MlWlVq2SyWRwM/ZSsgyPtcFajhQi4QUfBAEioXbRnpx9Y32NjfV1Tp85S9EKWMuwWA3WcpSo12pkcznK5RK+75PJZMnn8/h+A9d1946rVWsA5AuFaXXVkgpswm3LEUFVuXrlMoHv47oZCoUiuVzOmANqlEulA8d6+fyehmuxWA5iNdgjSrVaYXdnl5Vjx/aG/o7jICKcO38Bx3W5ff0ap0+fQUUol0sUCkUwxwZBQKFYZGFx0QpYy8ikVYO1AvaIsr21RbVa5ZgcB8D3G4gItWoVL59HVTl77jzVahW/0WBhYREIBWtzsmtpeXmaH8GSEkIvgnQOpq2APWIEQYDjOKwcO47jOARBwNrqKpVyibPnzuPl82ysr1Pa3eXMuXM06nXm5ufZ3dlhZ3ubSqXM0vIKyysr0/4olhRh/WAtqcBxHFSVG9evha+DgGq1CsDq6h0KxSKu61KtVlhfW+XY8RNcvXKZei2c0GqaBSyW2JD0TnJZAXvEqJTL3Llze09gtrK7s0O9XufkqdOcOnOGTCaDBgG5bI5jx0+Qy+XwG/ueBKpq7a+WWEiriSCdn8oCcECI1us1arUqq3fudBSuTWrVKrs728zNzVOr1qg3Gpw6c4ZisUgmk8HL5wmCgK2tTWo96rFYBqUZyRV1612n3CUiHxORL4rIF0Tk50z5MRF5XES+av6umHIRkd8TkUsi8lkReUEcn80K2BSiqqytriKOQ7VaoVqpsLG+ztXLl6nVqj3PdV2X+YVFtre2mF9YwPM8Sru7rK+t7R3TjN7yPG/cH8VyRIhbwAIN4BdU9bnAfcCbROS5wJuBj6rqvcBHzXuAVwL3mu0h4B1xfC5rIkgZoXC9w9bmJo1GnXKpRBAEA5+/uLTE1uYmxbli6FVQq3H71k2CIMDLexSLczawwBIzEruJQFWvA9fN620R+RJwHngAeIk57D3Ax4FfMuXv1dBn8QkRWRaRs6aeoen7qUTk2SLy6ZZtS0R+fhhVW0QeNMd/VUQebCl/oYh8zpzze2IMe93asBymXCpx9cplvvn1r7G1uQmENtUowlVEmJubZ3FpCc/Ls762yvraKidOnuKeb/t2isW5cXXfcpQZ3kRwQkSebNke6li9yD3A84FPAKdbhOYN4LR5fR643HLaFVM2En0FrKp+WVWfp6rPA14IlIAPEFHVFpFjwFuBFwMvAt7aIjDfAfx0y3n3m/JubVgMjXqdUmmXG9evUav2Hv734/iJk2RzORzH4fatm2ysr7O0tMzc/LydzLKMkVCDjboBd1T1Ysv28KGaReaBPwd+XlW3WvcZbVXbz4mTqHr5y4CvqerThCr1e0z5e4DXmNd7qraqPgEsi8hZ4BXA46q6pqrrwOPA/Wbfoqo+YT7we9vq6tTGkUZVqVar3L51k8vfepqb10caxQBw+sxZ5hcW2Nrc4PLT30QDZWFxEd/3Y+ixxdIHkehb3yolSyhc/1hV/7MpvmlkDubvLVN+Fbir5fQLpmwkotpgXwv8qXkdVdXuVX6lQ3mvNg5ghgYPAdx9993RPtGYqZR22d3eZHHlBNlcrv8JHajXauxubyBOBj9Qdra3YhV8c/PzFIqhzdX3fYpzc5w4eepAakKLZVyMIx+sMTO+C/iSqv5Wy65HgQeBt5m/H2wp/1kReYRwlL05qv0VIghYEckBPwy8pX2fqqqIjFXV7tWGGRo8DHDx4sWx9mNQfL+B32iwdusGjXqNSqnEmbuegZvJDnS+BgGVconS9ha726FNNVucp1IZzQzQiXyhsGcCWF45Zs0BlokzBj/Y7wN+EviciHzalP0yoWB9n4i8EXga+HGz7zHgVcAlQjPoG+LoRBQN9pXA36vqTfP+ZnOWbUBV+yr7s3fN8o+b8gsdju/VxsyzcfvWnmAE8Bt1bl39Fgsrx5lbWOoqxGqVCtuba5R2ttEIE1SjsLO1tZdroFat4rgu2exgDwKLZWTGEMmlqv8V6Fbpyzocr8CbYu0E0Wywr2PfPAD7qjYcVrVfb7wJ7mNf1f4w8HIRWTGTWy8HPmz2bYnIfUatf31bXZ3amHnmlw87PNRrNdZuXmf1xjUq5RI7Wxus3rzO7etXWL15nZuXv8mNy99gd2uzs3Adk25erVap12qICF4+T61WO5Rc22IZJ0NOcs08A2mwIjIH/CDwr1uKI6naqromIr8GfNIc96uq2vRe/xngD4EC8CGz9WpjJtnZXCeby5PJZrlzvbt9vLSzRWlnq+v+SbOwuHTAPlwsFmk0GmQyGWsusEyEI52LQFV3geNtZatEVLVV9d3AuzuUPwl8V4fyjm3MKjubG9SqlbHVPw6d8tSZsxTNBFeTne1tVu/cJl8ocOr0GTvZZbEMib1zYiSJw+p824oEtWqVjfU1VJVyqcSNa1etq5ZlrIwjF8GsYENlYyRpAnZhcWkvM1a5VGJjfZ1KpXzgmGq1yvWrVzh99pyd+LKMCYGUjpKsgI0R1fHO+sctwOfm5qjX69y+dZNqpbtpo16vc+3KZU6eOk2hWDQLIdpLxxIfabX127skRpKmwa7euU29Xh/o2CAIuHnjOtmcR61WY3l5meWVFUQktTeHZULYJWMsg5A0ATtIuGE7jUYDgI2NDTY3NykUi5w+fZogCGg0GjaFoWUIkmNTjYoVsDEybgEbe/0jVqeqlHZ3uXXrFrVqFVXlLhOqbLVay8AI1gZr6Y2q7i1pPTZiF1rx9Hd3ZwcI1/u6euUKx44fJ5vN2kkxy8BYDdbSk8SZBxiuz7000yAIqNVq3DDZvc5fuGBNBpa+CIJIOjXYdH6qKTARAZswId66zIzF0hUBHIm+JQCrwcbEuF20xsIQJocoD5JSqUS5XKZgl5ix9CGtXgTp/FRTQINkaZfQPdVQnKzeuTOBVixJx0ZyWXoyCQ02iXbeer2O7/t7EWMWyyHCWNlp92IspPNTTYEkCr9J6LCqyubGxtjbsSQbq8FaepJIATvENTrM59ze3mZpedlqsZbupNQGawVsXExAwMbexISeCb7vU6/XrYC1dCTN4dZWwMbEZLwIYpaIE7ymr129Sj6fp1AsUiwWyeVyqb2pLENgNVhLLyZhIoi7jWEEXPggGU4wVioVKpUK62trOI5D0QjbfKFgs3MdcZJiU42KvapjIoluWsMpxPHcCEEQsLOzQ6Ne587NMEtXoVikUCzi5Qt2FQVLKrACNiYmosGOvYXeiEj8dmAjr+u1KvVala2NdUSEfKFAvjhHoVAka80J6SbFblpWwMbERGywU/ZUUFXynseeVBSoVqqxC/7mcjXlUol1wHVdCsU58sUihUIR15oT0oc1EVh6kUwbbPRzqtVqWyXOaIJ/gFN932dne4ud7XAl3pznkS8UKRTn8PJ5a05IAWlN9mIFbEyMPRdsQuqMzBBCvlatUqvumxO8fMHYb+esOSGJNJO9pBArYGNCxSGTP7j8deslI4C2izRlv3T/v31hbf4qigaA45rJNI1FoMsk/bTGhKpSKZeolEusr97BcV0KhWJoTigWyWRsTtrZR1Kb7MUK2JgIgoBqrfv6Vl4+33NhwX6ICBph+ezQedvBcYwTtxH8YkR9XMt0eF6utVWqlcpUNePA99nd2WZ3ZxuAbDa3J2zzhaI1J8wqKR11WAEbE/00ytEvn6ZgHAxVRdUn6DH35nn5kXvV/tAQx0V7NTph6vUa9c0a25thPoS5+QVWTpywmu0sYZeMsfQj6OMHO+qIXmQMTgRjUBpmPSfD7s42pd0dFpdXWFo5ZjXamUCsBmvpjKqysb5OubTb+8BRr5+EXIChr2wEITsFeayqbK6vsbO1ycrxk8wtLNiJsSljbbCWjtRqVTbWB1gaZUTNLim3v+M4BDNkIuiF7/vcuXWDcnmXhcVl8nblhekg2EADyz7lconS7i5LyytUyoNNXOmoInIcGtYYtMfIvZyBJ0e93uD6tasU5+Y4cfLUntnAarWTIjlrbEXFCtiIqCp3bt+mUa+ztbk58HkyojQby+WXzmt6aEq7u1ytXKZQLOKIEGiAl8+Ty3nk8+GEYK1WC1fPrVbJZLMUi8Up9zr5CDbQ4Eijqqyt3iGbzaGqNOrd3bG6Mqo2lBBt6pCvb8Lw/cZexBjAzvY2Xj7PqdNncByHa1ev7plAcrkc9VqNXC5HoVik0Wjguq7VfKNiAw2OLqrKzvZ2JG21E41GY+Tz47ZvaqB4+XybqSD0kT34HlrdxHqJ0ND3NhQyrWfu1WZs0U0ZNI3VeDPZLK5rLn2BRr23f3Gj0eDyt74FHPSSqNVqrK6uAlAoFKgZYXvy1CkymUzLZ02n8IgPm+zlyFGv19ne2mJ3Z3tk4QgxjMZV0Zhv1Fqt2nO/53mHcw+MiOflqI0QcBEHrpuhEuFz+Y1GXxt6uVze+3v16tVQ4FarzM/P47ou2VwOz/OssO1GSr8XK2A7sL21yZ3bt+OtNI4LaMZ9TJPCuM0YfqPBznYYSba2tu9hUigUWDl2zAraTqTUTWugTyUiyyLyfhH5BxH5koh8r4gcE5HHReSr5u+KOVZE5PdE5JKIfFZEXtBSz4Pm+K+KyIMt5S8Ukc+Zc35PzNXXrY1xoarcunkjfuFKUueTxtHrGfgmpvScKpfL3LxxY+aDMSzxMehj43eBv1LV5wDfA3wJeDPwUVW9F/ioeQ/wSuBesz0EvANCYQm8FXgx8CLgrS0C8x3AT7ecd78p79bGWNja2mR3Z2dMtY8uWCZ/Y46hPQ0nh3I5j5y3v3kHtnzLX/M6n99/n8/vv8/vb3mztS8/kztwvjfV6K3mApCWFpoJt6NuCaCviUBEloAfAH4KQFVrQE1EHgBeYg57D/Bx4JeAB4D3aigNnjDa71lz7OOqumbqfRy4X0Q+Diyq6hOm/L3Aa4APmbo6tRE7qjpG4ZpaE1N0RGlUymNtIpMvQovdXIB6eT/SLpsfxrWqffJveFbv3OHsuXPWTNDKGLwIROTdwKuBW6r6XabsGPBnwD3AN4EfV9V1M2r+XeBVQAn4KVX9+1H7MMhj4JnAbeA/icinROSdIjIHnFbV6+aYG8Bp8/o8cLnl/CumrFf5lQ7l9GjjACLykIg8KSJP3h5yeL++tjZStqv+2JspZALfQxDg5XJ4uRz5vIf67ZOU0xuiZzIZTp8+bYVrO+PRYP+Q/dFwk0gj71EZpJcZ4AXAO1T1+cAubUN1o62O9art1YaqPqyqF1X14smTJ4eqf+zXu72fJkajVqFeLVOvlqmVSzTqtQP7p2kCPXP2rF3yphPNlJpRtj6o6t8C7XHsDxCOhjF/X9NS/l4NeQJojrxHYhABewW4oqqfMO/fTyhwbzY7YP7eMvuvAne1nH/BlPUqv9ChnB5txEqj0WBjfX0cVScc+1RoEldy8rXV1T1buqpSq07XZW0mEAm9CKJucKI5cjXbQwO0FnXkPRJ9Bayq3gAui8izTdHLgC8CjwJNT4AHgQ+a148CrzfeBPcBm+YDfRh4uYismMmtlwMfNvu2ROQ+Ywd5fVtdndqIlWZy5nEy8u2ZkiFlYj9FTN9/qVTi2tWrbKyvcvXrX+bm5W8QREiknlqG02DvNEeuZns4SpOTGHkPOlb5t8Afi0gO+DrwBkLh/D4ReSPwNPDj5tjHCA3FlwiNxW8AUNU1Efk14JPmuF9tTngBP0NoLykQTm59yJS/rUsbsRAEAZubG2ysDZANa2RGu0GjpduOicRKw35M102qWq0i6F5i8ltXn8bLFxFHyOY8igtLR89GOzmvgJsiclZVrw848h6JgQSsqn4auNhh18s6HKvAm7rU827g3R3KnwS+q0P5aqc24qBer3Hj2rWBo7TyXi4UciKgLbeohm7rrcM+VSUIgnjdqsaScXvyJPUTRM5z27/GvVf1aoV6i6kgk/Pw8kcodWLTRDAZmqPit3F45P2zIvIIoSvpZospYWiOrLW9UqlECoH1q6WBl0IRwCvMU6nWDhaOwFQ02JQyjJzM5bI9PbUEgWZehfYfy7zfD7dVhO7X0u2rTzO/tMLS8VNHR5Mdw+cUkT8ldPM8ISJXCP3wu42KO468R+XICthaxBj7QKNmdI35gjkqN9oEkCEkbD8XvoKXo16Ox49ag4Dt9VXcTJaF5WOx1DnzjMFEoKqv67Ir0sh7FI6sgK2UIzq7T3l4bsVrfIwjF0HcNTpuZiqZxqaDXZMrVagqfnWXDIIjoctHuIR1+ENL869JuRdGAtX61Dpu0nkBWjqTL86xuHJi2t2YDEJqk70cSQErIrjiEAQ+gQb0MIcNjRPUyfqVUECLEJTreCKoODgiqBHgTYGu0u5p2fZUFyHI9vi5dO+/w+ld28tby8xr7XBsWkX6ODTYuG2lpe1N5haWyM/Nx1rvLKIQeyrOWeFIClhVDQXrBJCwQVB/T44P4/XoZrJU1Y2vYwOg2exE20sy4zAg3b72Lc7ecy+Z1P8O6U24nc5P1QcRYX5xhUwuByLkxuISE/Mtl9In/FQYhzQcQ52O404189dEOarZtNLKyqkzB5b02F5fZePOzfgasD5Vs0sCfpviwhLZnIfjTnbUMi2siSCFtNrNFlaOU9rejC02PG47X1yx8LNAd3ul7Bl+D3xeOXhM+9n+oWxZ/Zh9C3N5Z5uFu45PuxuWETnSAradTM6LL/lG3FpSip7wQS/hpodedJq1O0A2m4220m8CNFjVgHq1Qs7LT7sr40esDfZIsHLyDBKTzSt2H8b0yNcxEO3LScrS4pXSbv+D0sIY0hXOAlbAtuC47tHxPRyAZFzCGB/m9DGeydcZZbh0hTOPNRG0sXjsBOI4bNy+MVI9cWuwod1ysppXMvQ8CB8Fsjfjru3mhb1cAGqS88x+JJdXmGNuYSnmWqPRTFw0fk8GSe0kVzIeAxMmjgtKg9iNsDHXN4stDkeYOkVoKDQUfJX9jf2/AU64qZko02Bvm6WkKl6hyMnzd0/Fg6Bc2sX3fUq7OwRBwOrNkRNK9UewblpHiUGzZvWuJDn6XzeS8gkkik2ueZz2nGqbOtMQ+KrK+p3b1Gs1QJlfXOLEmXOTaTshAjMqVsC2UauUWY/BHzZ24TRDGtYojGPp8ThqnPyS6N3p1BcNAur12li9ClSVRqNO8xvNefkJCfrkTFpFxQrYFhr1Oqs3r8aifeqMa0hTYxx5P0f8pgf5tTsJmtYyEfaG9KphRvZh7fDLJw4unqyq3L5+mVq5xMnzz8ArDLPseH9830eDADeTQYOA+cXJ2YCtBnsE8Bt1GrV4smYlQbjmslmCwN8TFOHf/Z4Hvh8mmjZlcuA/I5ZMgpggXNZhr/yg66pZAkIEZwbjVFs/tbbMiu1F+gH5bI56tdS1jkbjoP90NudRr0XLOQxhzolKaRcRoVou4bguldIuVeOytbV+h5OFuyPXOwiZTIbjp8/iiOD7/mTNFFaDTTeqyu7WRix1iVlWJl7ivwA1CPAjrOoQB5nMGCZuYvhqDmqb+4J10viNOltrt9lev9PRVJDNjddEkMt5BL5PcX5hbO0cIsWBBlbAAkHgs3H7ZmwCNilP46Q43B9FOglXN5NhYSXeFQ5Kuzs44rCztcHC8go5L08QBNRqNXK5XKxtdcOmK0w51VIpPuEKe4m6LYdJ7G0UsePjeHgtHT+F68Z3y6oqa7dv7o1isvk8Xr6A67q4k3YRsxpseonf3zB+MeKoT76+c3AtPcehmSBFMEK9mbh7z57aXJ0h7NdeYm+F6hRmzsfT4vjF9rQfmCJCvhhf8u0gCELh6u9nJ85mJ6OxdiLqindJwQpYiN0OORYTbLPu1jct/rqt7XVqu9N+8SZoZ0s40779MzkPNxPf7Vou7bK7vQVAoThHvV6bYrCFWC+CNCOOkDFP79DFJgynbIYKRnfbGseFOgb/0WmrZZaBmV9cjnxOEPg4TufRWWswTb5QZPn4yVgFeGSsgE0v5Y1VtM0F56DrTjgcF5zwOhCHA4sjwv7QWwQRh8DNdJjr0ua/A+9Dl6Dm++b6WE0BH4xREloJOzZi/GpFBK84F+mcjdU7OI7DYodJsUq5RLklU5fjuuQ8b+R+Do3YSa5Uk59fpLKz1XW/AAQBSsAgvuOZnEetMczKWy19yrg0apU9Id8oRch3OiBWvA5OveGTKczTKO8MdHyjXsPN7K+l5TjOUH6xEI6qWusahCDw91JvBkGABgHVSoV8ocDNq5cPHDtVzTXl2G+WMSiIMTyNJxK6aW0EA+P73YfbXc9p7D8U/RGvibWb1zh++tzA+YqPnTxNEDQf8kqtVuX2jasU21apPXn2PKWdbQoRNeQ40RTbYNP5qSLiFeMOPYzF8330OmaR2QvkGpwpjmLLO1tcf/prBBESETUfCI7jmgQuod/rgXp3d1g5cSq+jg6LTbidXkJZFt8PpoEfg7/i+KXGVBKcSEofHBPAb9QpbW1E+t2ax25vrnfcv7MV3zp0o6DiRN6SwJE3EVRLO6xd/SZxCjS/XiNTmBtiMb4Wsh7U48mL0A3HcSJpRFHJZbNkZH+BQRXFUcXN9bAnNtMctD/vevw8go+by7QW7B1faQQ9P2Mun++ZLKZ9QjKTnzvUHb9a6j/iUCXr5RFCm2hjyN92/fYNttZXKcwvsHT8VN/cxSISTpZ2yE+c8/LMLy3NwLpfYv1g00jgN1i78o3xaHIjVjmZy228rWQlwN2cQMLmHjgL53o/RBQqQ04+Ncky2M9dN5pidkSB5jfqiDhcu3aNCxcu9D1+c3310MM+k8kyt7iE5+UnsGJBf5KikUblSAtYcVy8uYWeHgTDkoyBcDJ6OQp9TXVxmMtHryIyC8srLBvPglBDDTpGJO5ub1Etlw+Vu5kMC4tLs7GSg5AYm2pU0vnYGBARYeXcM3AiusAMQhISqcxSkunx0efG1WQuxeO02PhFpKNwVVXW7tyiUj6cZjFfKM6GcAVCE4ETeUsCA2mwIvJNYBvwgYaqXhSRY8CfAfcA3wR+XFXXJfzVfhd4FVACfkpV/97U8yDwP5pqf11V32PKXwj8IVAAHgN+TlW1WxsjfeI2Ar9B0IjfxzQJXgCe5xkh2xpS0UEj2zum5bGxF29rEl6LsaS1fG7xpz95knEcyOX2hVq7UIlDxkT8qRXI5YttJwqBX6dRj+9a3FxfJfAP+2MvLC2zdOx4bO2Mis2mFfJPVfVOy/s3Ax9V1beJyJvN+18CXgnca7YXA+8AXmyE5VuBi4Tf6VMi8qgRmO8Afhr4BKGAvR/4UI82YiOW9bc6MvsXTL28G6uePZ8BdlZjrHF0pFEmaHT/jd3C5PMx1LvM2ufyRbKOuxfN14wMDHz/gE+tDOCP6/s+m+trh8pFhMWV4zOkvYak1QY7yqd6AHiPef0e4DUt5e/VkCeAZRE5C7wCeFxV14xQfRy43+xbVNUnNFSn3ttWV6c2YsPNZMmMJYnx7GuwcTNbt2zILPapG7VKiXq1Qr1WpVGrUq9VqVcrOJkMkslSXFph5dRZlk+d7SsgK6XdjqOobM4jM4ORW2o8CaJsSWDQb1qBvxYRBf4PVX0YOK2qzSniG0BzIaHzQGss3hVT1qv8SodyerQRG+I45ApFGrV4h7NHT7xCEj/1LN+mXqHI0onTePlCpPNUle3NzvmN67Uq1XIZrxCtzvGS3kiuQQXs96vqVRE5BTwuIv/QutPYS8d6d/VqQ0QeAh4CuPvu6OsVNUZ00+nE6BNIyRNWs0lvESoo2Vz2gC/sfkqVLXQAACAASURBVD7d8Hes9VunrcXvNi4cx+X4mQtD5Qmo12pUK4c9BwAWl4+RnWZily6k1QY70GNDVa+av7eADwAvAm6a4T3m7y1z+FXgrpbTL5iyXuUXOpTTo432/j2sqhdV9eLJkycH+Uh7lLc3qJV3+x8YlZEF7ASSSCdgIm7c1Ms7aGWXoLKzt/mVHfzyNn55G0f7J+0Zxy/luO5QieBVlaDbgoUiLK4cmwm/11bCpTLTaSLo+02LyJyILDRfAy8HPg88CjxoDnsQ+KB5/Sjwegm5D9g0w/wPAy8XkRURWTH1fNjs2xKR+4wHwuvb6urURmxEzVI0KKMKr2QuOjOLF/1o3+NgZ8f/uRv1GuXdwTJ37Z3TqLN+5xY5zyPfIXmLwEyExR5C5EiHyp4GPmCeiBngT1T1r0Tkk8D7ROSNwNPAj5vjHyN00bpE6Kb1BgBVXRORXwM+aY77VVVtTnP+DPtuWh8yG8DburQRG74f4Jjwx+aK1O3DxYORm+HzVlUhCEwqOP+QN8KoAjaJ4jWNDPSgG4OJwCvMUZgbbImYRqNBo15j/c4tatUqIIe01Gwux7GTp8kX4k5sFA9J0Uij0lfAqurXge/pUL4KvKxDuQJv6lLXu4F3dyh/EviuQduIE4X+Nra+uGDca6SZ6WeQxLE9cMQxvpIclPLmRj4QyKBtJdriz9ou6Fv3xe6hNouPhdFu3HAp68NrVcnef+Bq+Nt3i0Zq1KqRcz4s9HGl8hsN1m7fJOt5VCuV0GvAsLO1wdzC4oHjT529QCY7ntFaHCRFI43K7PlrTJg4V5OluZqsBiM/jzXw0d3RYir69iG32O+I5DPi3KtfqwK9J0H7WWkz2RxBMPhD/NSFe/B6aJqqyu2b16mWSziV8qFgAlWltLPD/OISO1ub5AvFmRau40JE7icMenKBd6rq2ybdh3Q+NgYk8P2OYYQWy7QozC/0FK4ApZ1tqua67RSpBeGKBjtbm+SLc5w6e77jMbNE3JNcIuICbycMfHou8DoRee4EPsoBjrQGu37n1lgiuSQWo9wsDrd7I5qEDAyzzdxC78UNw+W2OzrTdMRxnIFXQZgWY1rR4EXAJWPiREQeIQxc+mLcDfXiSAvYbk//kREZ3U3LSqpYkCkkcznMYD+mm8mQ7zOxtXrzestSMP1pt8XOKkNOcp0QkSdb3j9sgqCgc2DTi4fs3tAcaQF7/Mw5ttbvsLU2W/HzIeOWsLMgeCbB9J9UYe6A3klcHNfl1IV7ek5slUu7h5Z86UfUKLBpMWSgwR1VvRh3X+LkSAtYx3FYPn6K4twCt65ejqQZ9EQkBgvBmAXDUZGvM/A5B0mssrB8jEz2sLdCKxurd3rubyfnebhDBCtMA41/pNEtsGmizLZxZkLk8gVOXbg7NlvVrGUq6sRY+jiDH3s2lgDr/8XsueT1IGoE1uLysUjHT4+x5IP9JHCviDxTRHLAawkDlybKkdZgW8l5eY6fPsed61f6H9wX2Vu7CWkNW1Czt/lO2zTdppOrQmPcETdjkIYzIczamAGh77Q5HDe1SjHXhkgYhdWLcmkXRXsM+XVv0KMoruNSnJ98KsZhaIbKxlqnakNEfpYwgtQF3q2qX4i1kQGwAraF4vwCpy/cw+3rl0eaABNHqNRGS5y8SDBeeTUGDbaRKaDFDgmtey1Y2HK4k8niH0g4rXsThiaWwpTux9aFAsq8FjkUeaVulqx0mndsfbh1aqdLp1WRTIZ6dfAEQdnt2xR31qBexqnsHBIlsnSK3Nm7Op4bNqlsrN6JFOZaWF5JxEiqyTgiuVT1McLI0qlhBWwbXqGAly9S3t0euo5Y3LQSmIglIAzbHJasm6M+wvkd68x4NKqdM0sNX2c0p30/v0j++pe7HzC3hHiHNVMNAvxvfYF6tUzNXRq4PTeTSZB5IOTIhsoeRZaOn0BEKO9uD5dTQBLgB5sg7SbtSK6INsOsDf6tb9F4+nPo1ipBYQlODSZgRYRTZ88PleZweiQnO1ZUkvQrTIycl+fE2fNsrN5may3azG0nHMchm8sdziPQmkOgLX+A+C6SyYblapLLjJjf4AApjf1OIpLzCKolXGO3DzbvUP/837b83oM/bPPFOXIjLgs+DcbgRTATWAHbg+L8IvniHNXSLtVKOUyoMUAQQbvty3FdqpVok1abeJDZT4zsZTP7Q929pDIHM3+FjYdl2VyOeqPeYq6QPcU6zEfjwMhJbtKPiDAnPqHN1jwEKxt0TVndmnpNBHBwajWCE3fh+nV0/cbB+pdOIsfO4XihF4EGAfWv/N3Bh2mEaMO5hExstTKOSa5ZwQrYHuRM5vdmirdGvYabyaJBQKNRZ+3WDWpdMsfvIRLrSqHNOpvzMIeCU80OV4SGr+xrPwePczNWgx0Mwd24NnItPuBkOvi5Nuq4mez+Q7lWQrcPLlYoAyT9hvBBnhTPgXasgLXsOYKL65JzXRaPn+Lmtavh7HXLVg/2hVkum40hHSJEtsn2GXKN5XIe0a47m2bh8drCdXeD6if/ksxd/whZWKHxjc8eNgUNOA+Q87xEeQ60YgWs5RD5fH5Pm1TVzjdCTBd81Lm2/q2Oww82eZ4P/Yh1WZ1u10K1ROPSU91PCwYbAXkdPBGSgVgbrOUwjuNQKBQp90h5eMhGOiQS8Ubvd736fsNEDxmjoUBtyqkbgyAga8wxQpiMZxwLUiaOAX2yk5rzVYHRMyjPJlbAjsjJ06e5dfMGlXJnW6zvx+PXqVFjPvscrqpUWxzXZ2EhvEabI312gPDRSXBwyaDJM2jbs7oczCCk1UQw/bsq4biuy5mz5yh2WGQuk8mM5Hh/gPhtBG3Hz+AFPjNdiqkjI1QziG01Wb6vRwP7i8SAiHD67Fnq9Rql3V2qlSq7uzu4MQrY8Zs3RxcicXdxZuSr40AsmdaG/0Qi0tMenM0ld4Ir9IBLaN/7YAVsjGSzOZaWQ0+DarXK7s42jXodP4bE3uO+/DQIyBWKoZ/sYeevAZCB3YkG7tPMiNiYJirH+ImSkve1G7PzW8eLFbBjwvM8PM9jfn6BO7dvUY2QHKQTUWezo16uqkHkYIh2itl4c4/OzC03A5phP+1UnOn3cXisF4FlSHKex9nzF9jd3WF9bW3ooIOoAjZ9DlNTJDYBO4KJoM/+uflkLA3TCRvJZRkJEWF+foFcNse1q1fi9a3sxhQkbPxNpu0xMYbPI4LruntRh0nFarCWkcl5HqfOnOHm9euRz41sFZ3C9Rp/k2m76Qb/PJVzzyWQ/SxTrl/FBZP4J9gLbBHAyy0kd4LLEP/azrOBFbATxnGGs1OqRhxEpUKDnQ0if/fd6mFwEVsR9+CwuW3Fg9Z6GtXKoXSHSSOtGqz1g50wjgjZYSJuEhGGmoQ+DkNMN3+E3zBKXInfqFPd3RqiQ7OBmnywUbckYAXshMl5HgtLg2enH57kq7DDOIulhfalb/pRWh89b/E0UZXIWxKwJoIpkHGjfe2ZTAaNuMbXNIaLIoLT87PpASWucx+bx2hseRwGQ/adBUQQZC/nrpPzkCBjvAmax8ne+70TBVQyaPvnMu99BU6GJiJtOVdbjgFBRSLH5ld3twn8Rp/vf3ZJikYalWT+GgknSlKOXC6H1muRdTk/5rWtBmrTzVGt9wk26HAfefkilda8uuaYas0H12TnVyWbzZnVV2VPtoX7Di7R07pCRL3aAGc/D2tTeO+3dbBDK7vXoC1/RJQJGH/lPOV623ff+uPlVw7sctwM1ONJaBP4fjIFrEKQ0sFKAn+N5ON5HsW5OUq7uz2Py3se9fLOhHo1HdxMhnp9gHy5EsaXNRfW0dZc4oMuWwtGKHfXlkSiDs7biertEY9kcdwMbrZDQu8EkGY/WGuDnRKnTp9heWWl/4FDkpQJZVXFHVDrmoiSM+JaZZHFc0y2xPkTZxLvRWBtsJbYEBGWV44RBMr21mbX4IPhU+XN/gXouBkcx8GPJZFKXIz4vUV9Coy6kKUICyfOUlw+Plo9UyYRTjJDYAXsFBERjp84weLiIqurd6iUywcEbaVaJesV8aulWMVlJpvFdV3qlTL7EmHQFlrvhDCKSBxnKO0p/KyC3xhsAm8ij4xRtcCoIc0j+Nh680ssnbmAm0lmou19ok/qJYWBx0Mi4orIp0TkL8z7Z4rIJ0Tkkoj8mYjkTLln3l8y++9pqeMtpvzLIvKKlvL7TdklEXlzS3nHNtJGNpfjzNlznDh56tC+er0e6w3kui5+rUK9vBtqT82lbjQYcNMD5/iNOo1alXq1wqDqmxKm1wuCYDD7a8t5Y2fkYXZEATvk7E5+YYmV8/ekQLgaG2xKTQRRDE4/B3yp5f1vAr+tqs8C1oE3mvI3Auum/LfNcYjIc4HXAt8J3A/8gRHaLvB24JXAc4HXmWN7tZFKCsUiS8vLh8rd7BBx5l2uP5HJOj91IpfzYloIcgxMWIMd9rGRX1hOtM31qDCQgBWRC8APAe807wV4KfB+c8h7gNeY1w+Y95j9LzPHPwA8oqpVVf0GcAl4kdkuqerXVbUGPAI80KeNVOK6LseOn2Bufv7gjoh5VotZl0zX9HXTvymDtBrcgOheBDKkUJ/+7xgnBwZGA25JYFAb7O8Avwg0F10/DmyoatPh7wpw3rw+D1wGUNWGiGya488DT7TU2XrO5bbyF/dp4wAi8hDwEMDdd9894EeaXU6cPEU2m2VjfZ1sLketstvzdirmMtCootlCOIzfvgMIc24mdNdvdYwPQg1W3QzlYDpOJLOw/ld3RhNcURenhP6rFRw63nXJdViiKMmk1U2rr4AVkVcDt1T1KRF5yfi7FB1VfRh4GODixYsJebZ1x3EcVo4dR0SolEo9Hd0LXg62bgEg5e29ckHBr3e8bBUg64GTj7PbA+OkeGg7VCrKiN/H/LFTqbC97nHEAw2+D/hhEXkVkAcWgd8FlkUkYzTMC8BVc/xV4C7giohkgCVgtaW8Ses5ncpXe7RxJFheOUawtEyjXqNWrbC9vka9bRlrqXdezbY/0xNylUoZx3EJghlMUjei8JchEu+JOCiDm4G8uYX+ByWI5iRXGuk7VlPVt6jqBVW9h3CS6m9U9V8AHwN+1Bz2IPBB8/pR8x6z/280fKw/CrzWeBk8E7gX+Dvgk8C9xmMgZ9p41JzTrY0jg+M45Lw884vLnH3Gt3HunmcdWD00yCZvLSZBIoUL75EELWcYVSyiUM94yfvN+3HUbbCd+CXgERH5deBTwLtM+buAPxKRS8AaocBEVb8gIu8Dvgg0gDephrM3IvKzwIcBF3i3qn6hTxtHlkw2y/GVFXaufT2MOmokMZRWqVcreF6BajXKOmCKVyi0vg1ppiI4mJJg7xABcBw0CA4eqB1kdmNE+/AQgQNRxGvGy6fSeyCtfrCRBKyqfhz4uHn9dUIPgPZjKsCPdTn/N4Df6FD+GPBYh/KObRx1vKXjqF9n9/o3R6hl+irAMHKiUhkuMYqXz4+88ORAjNkGO3/8dPT6E0BSNNKo2EiuBCIiFI6fpb6zSW17fbhKZuCKjmqDHanHg37eUW2wQQNP6y3j2DA4o3WJl3phmaAlY1cUjTSbRvOAXVXWMosUT11IrIAVx0GiumuN1OUBb+BRBawGHdMPttZ6uAk5sFNa886KIBKGIotIYjNm9eSIexFYZpRMYZ75c98W2mMjcigp9KRRpRZ5yD78XTjomb43h9slo9ahOg7YfcO9AUCt9+fytEGg1T0tN2jskm3RcLvhesXoD6WEMAMDqrFgBWzCyR87TenONYJalMmi6ax40IqqksvnqVai9XsYMtnswHdwQx38SBNvQ+A3oCUHw6Ai06+V0cBHhlw4c5ZJa6BBOh+HR4zFu589RB7T6V/QGgQTEfSu41Lro1XuM/3vpSuqxhMiXSihiSDqlgSsgE0BmXyR7NxipHOmrcEC1GtVMplJDKJm7G4c4buPOlJJCmn1g7UCNiVk56OtVDt1G6yhXq+R88YbshuoDh7YMIHvZZSHW2X9Vow9mR0mLWBF5MdE5AsiEojIxbZ9kdKq9sIK2JTgLR5PpG1OgHqtQnag2fHhBFO9Vhs4wcysK0bVjdsEEXLoJgFVCFQibyPyeeBHgL9tLRwyrWpX7CRXSnBzHm6+SKO03f/gWUOVRqOG47oEfo+Y/A6RWv3IeR7o4AlmJmI6GakNTeSDtB+THvKr6peg4++9l1YV+IaJSG0GO10ywU+IyCPm2C/2ascK2JSgQRBRuI5BkIxyk6iSy2ap9BKwQ1Cv1SJluNoLrR0rw7cgbgZx0ydgh+SEiDzZ8v5hk1lvFKKmVe2JFbApoVatwsJxUKMgHRIqxhHGSBB1XHLS3y657+qpHaWPIHtLYY/qajNUqr/YmQ3bdDe0UcevVXFzQ6xyMcMM+dPfUdWL3XaKyEeAMx12/YqqTiRxlBWwKaFU2mXnwNqB7YKi/b0C8dryvPxoN32tWo2cfLovIrM35Tzhdb+SwDjcrlT1nw1xWtS0qj2xk1wpodFo9D9o3Ix8k2hoM42RWdRHR+lTdm4JNzedROnjopkPdkYWPYyUVrVfZVaDTQkz4YAewzVfq1TIeV7HMNpELNs9AKM8h/InzsXWj5lhCn6tIvLfAP8ROAn8pYh8WlVfMWRa1a5YAZsS8oUCpd3p5oYVccjlvN6SsEve1lYcx8HroMmKCJkDyV87VRy+FvO/tDVUbfgEvr+XPKW1bjE1OM1VfM3+/eOkpdlDxmjzR/CrYUhrN4YV4ZniItn5w6sOp4FJR2ap6geAD3TZFymtai+sgE0JjXq9/0HjJmjgV3ZHrkayOfwx+Xo6ZmsKuXw2g5Y2Dxwzqh9DJufhV7sv5TOsLCmeecZMRODFTWgimHYvxoO1waaEpZVjwy3DEicx3SSNeo1coUhmDDPlTS11vMTfQnZhhWwxXWtxtWJDZS0zjZvJsLxybLqdiFG7qpVLOBPw95RxiNsxVOktn4q/0hkirclerIkgRZTLpVjry2SyGA/YrjQFlKKIOOFy0kbQauD3jszqQ71aJZvLEwS+WehRqFXi/YzjmIwWJ4OTyxv34IPJtEEQxyFTmDM2YkCc0Gps/ImhdTAg4S+QpmW620mQRhoVK2BTRJzamIhQjygcG34NcPekQzGXIygPP/GmgU+9FvbBb9SRyCkZp0TQgOpu1wfToW+1uEi13tsLZG4mHc7iQYFZcIIZBwm5Yi2DsHLiZConQfYYy2h+HJVGq1MGUN9mwg1vjKTVBms12BThui7iOGgM8fyqGqtNdVaZift0AGmRyadvscNWkiIwo2IFbIpo1OsEfoyazswJ2WRMSEWttJ9syeaLZNK42KFBEzRpFRUrYFNEOMmV0it1bEz/ASIa0Mtat3gqhdFbbcxGop/4sQI2RcwvLFLa3aFea3PS1/A/7VWmiqruJVuJ5YKfvuxKBL2+6+LycXLF+Qn2ZjqkVL5aAZsmRITi3DyruzeHrqN5s2tY4Wgdiv2mSchdGPVr04AwvL29HmHh5Nk4ejTzpHUOz3oRpAxNqzFrbMxC4vHOfSgsLOO4VgdKMlbApoyd7c3+Bw2IDLiO1SQQxyU7lsUR430gOVkPv7wV6ZxGprOHQDblngNNhnHRSopJYXbuIEssHD95mvmFaEt4d2JWzKduJouIgwY+tZgj1ULiv1MH8Ws9eELnb9vJpNdzoB0bKmtJBF4+j5c/w9LKMe7cvEG1Whm+shlQExzXxW/MQKawQRniO+v2MHPTHB7bxgxcamPBarApJZvLcfr8heln2BqBXL5AfZQHxCDMwo3dQYN1MlmyheIUOjMdNNDIWxKwGmyKcRyHYydOcet636WDZpJ6rcrSra+FQ25VwqwgQfjeb4AqEvjG8yHcV/n2i/hOxiTaVkTVDNnDOpqpVEQD8ANww7T18RHPjZ8rzKU77LkFG2hgSSyFYpFsLnfYN3YARveFHe18DyV7/VKkc+Y+/deRjq995/fHKmBVR11bN8RNceRWJ9JqIrACNuUMKyTjud6jixoRIeO65HdWyX35if4njEio3caoKWo8Dp218u5e4MdRIEipCtvXBisieRH5OxH5jIh8QUT+gyl/poh8QkQuicifmZUWMasx/pkp/4SI3NNS11tM+ZdF5BUt5febsksi8uaW8o5tWAZDNeD2jetDaa/TYmF3laWn/hLvy0/MjCdDFIaxDXYSovXy7ki5dJOEkl43rUE02CrwUlXdEZEs8F9F5EPAfw/8tqo+IiL/O/BG4B3m77qqPktEXgv8JvDfichzCZe6/U7gHPAREfkO08bbgR8ErgCfFJFHVfWL5txObVgG4NaN65RLo6+RlcvlcNUn8BuI4yBOBhWzspVA4PuoBjiOQ4DQqDfQPpqc42YI/P3BuZvJkkXJffWTI/c3Ck6jRmbx5P5CjM2k1wqI0j7gF9UwSbcKiCJBQHbztrnjAxTBQaHpQxwoqG9ClXTfhaspJQKfK3/5FJufuYRfKhNUqgS1BlqrU/3JH+Gef/MTE/ompkiCBGZU+gpYDceYzazJWbMp8FKg+eu/B/hfCIXfA+Y1wPuB35fwEf0A8IiqVoFviMgl4EXmuEuq+nUAEXkEeEBEvtSjDUsfqpUy5d3hhWvGdRE3Q8aBenm7o6lBvCL1+mELpuO4SKb7pZXJz1Op1UCyYYpFcag26mRrU1gVt1YZyZXN9RsULj01Uhc2//az3Hr8sDmk8pLvHane5KAEKZWwA9lgRcQFngKeRahtfg3YUNXm3XUFOG9enwcuA6hqQ0Q2geOmvPUqaj3nclv5i8053dpo799DwEMAd9999yAfKfVsbWwMfa7nefjlbaiHM+xRh+pB4KOBTzafI5OfQwFHBA0aiJOhXDUmCxH8IACC6aVFHLHd9mXB4+TCgz8ytrpnjZhM1zPHQAJWVX3geSKyTLiW+HPG2quIqOrDwMMAFy9eTOejMAJrt2+xu7M91Lm5nEejvD2QUO2ldAhQqRzO6iXOrAUNjCrYY3gwdPkiNc7cvjNMaINN520bKdBAVTeAjwHfCyyLSFNAXwCazpZXgbsAzP4lYLW1vO2cbuWrPdqw9GCYIW8mkyWf9whquxFERtSQ0Nm7kXQGZum7LSv5qZ/4+dQvFQOAhibqqFsSGMSL4KTRXBGRAuFk1JcIBe2PmsMeBD5oXj9q3mP2/42x4z4KvNZ4GTwTuBf4O+CTwL3GYyBHOBH2qDmnWxuWGMnmcmhtl0ZpOzlXblzMgMDvllRn+/NfYeMTn5lwb6ZDMwdxlC0JDGIiOAu8x9hhHeB9qvoXIvJF4BER+XXgU8C7zPHvAv7ITGKtEQpMVPULIvI+4IuEpr03GdMDIvKzwIcJk2K+W1W/YOr6pS5tWHrQb9kYEcFxHLLZLAQ+jcrOkAPd6Wt/oyIjG/9Gv9G7Cdi5Zz+The+6d+T6Zx3lCEdyqepnged3KP86+14AreUV4Me61PUbwG90KH8MeGzQNiy9Eae34Ms5il8v06iXR2snbvmaRHkdw5cgmQ7JtoHv+J/+LZmF9K9mEEY5p1PC2mQvKSSb6xyPkc3lyGdd/PqMBh5M5R4bddWG0TvtuJ1vw8b26D7MluliBWwKyReKxgwQakaumwknsCo7NKqjaa2pY2QngtE12G72xK/9b++kfOXGyPUngbRGclkBm0I8z8PVOtIo42UEx6/SKA3mehWJhFzk4ySWr6CLtChdeponXvYv8SvVOFqZaYJAI29JwArYFJLNebiOiwB+rYoGRyOmfShGVoVieGz10IJLX7/M1mf+YfQ2ZphhPAjS5EVgSRhhZHIyLkALSBcb7N7+PpOWaeBIR3JZkoebyeI34k0lfYiY7/t4MqlGbXS0B5HGocB28SJoUr52k+XRm5lpjnQuAkvyyOWL1Cp2Qqsfo+YS0BgkrNNHwOZPnxy5jVknKUP+qFgbbEpZPnWWkxeeieP2vnlnimmErY54X8eiwfYxEdz52P83eiMzjGp6J7msBptivOIchfkldjfXpt2VgVAE8nP7BSKEOWeNFNvLpdrBYNc8trITcdg/yRtVIJMBJwOuC064Fc6XWH7xd+PkskjGxcllQJzQ9uoI1WvXUr+6QUoVWCtg006SJkjq8yvces4/GamOY6VbSKOxZx/WFgHdTOwS3ssmWbiTxRM9IMTVEaR5npslCE9uScqiqJrk2xqw84++r0X7FlTCdpvtzPkltLzdVUM//cBpTj/wwp6fq3b9a3jnnjXMV5II0hrJZQWsZWiaQqhVbBzUskJvhnb7mkJHlUVimOSq5YogEbKJqQ/VwxFTzd5Jfo6G3+fmd/ssjS4ysvmjdvUrSNYje+JC6jRZ1SOecNuSXNx+N/8ImIWxDw6yO77pIBA6CIluafsiMZYECTPgKwtUv/k53MIC7nz6fArSqsHaSa6Uk/W8aXdhoqTzNg1xF46nUrhCKGCjbqMgIv+riPyDiHxWRD7QTMlq9kVanLUXVsCmnFyh2DUdXiqJXYOdIa06SR4hUdAwXWHUbUQeB75LVb8b+ArwFoC2xVnvB/5ARFyTrvXtwCuB5wKvM8f25AjdeUcTx3E5fu7uZAjZWGRZzJ8zDttgXAI2pSHPyuQ1WFX965b1/p4gXDEFWhZnVdVvAM3FWV+EWZxVVWvAI+bYniTgrrOMSr44z+m7v5383GznFo3DpzTuJWBisQ3GJPQzS2kNOBg6F8EJEXmyZXtoyA78K+BD5vXeoq2G5mKr3cp7Yie5jgiZnMexMxfYuH2Tys4WwSxqQ7EYUOMVsEJAuNDGKJVE6JM4LV4H0lLskj15V/fzkowJNBiCO6p6sdtOEfkIcKbDrl9R1Q+aY36FcIWVPx6mA/2wAvYI4bgZjp0JH8bV8i67m+uUd7aHzrYV/4RSDMuvxCxgNfDDgIBRaBOwzsJx6o0GTeOjahCaInoI4sz8MjJGj5BpM45QWVX9Z732j3p8mgAACI5JREFUi8hPAa8GXqb7Hei2CCs9yrtiBewRxSvM4RXmUFUatSrVcolqeZfK7s4U0xvGkLxa4haxcdDeI0H9tkQ8fbRctzDb5p2kISL3A78I/BNVLbXsehT4ExH5LeAc+4uzCmZxVkLB+lrgJ/q1YwXsEUdEyHp5sl6e+eVjqCrrN69R2lof5Oyx9y8ysdtgg5FnKvaiycQBcYayNcuoWvQM05zkmjC/D3jA4yZw4wlV/TdDLs7aFStgLQcQERaOnRhQwMbd+ExUcZAYEpX62TyaLe4XDJHlzC9tj9yPmWUKix6qate446iLs/bCCljLIbI5Dzebxa/X+xw5e279GrsVNgZi8CII2k0KqSK9obLWTcvSkRPnnjHtLgzHDMbpxyM60imAmkzaD3ZSWA3W0pGslyc/N09ld2dyjc6km9ZsoP4MutXFhJLehNtWwFq6Ulxc6S1gY78p0nmTxUFQK6N+A3FTeMsO7wc786Tw17LEhZcvsnL6nJGjuvc3/KeouAR00/LC0vYRe7imfdASkWOOUZPqsFDoUJdpe++c4IDG0/o6Iz4SLOyXNSepdK+mtl7Kfrnq3kNjzzcVcLI5BEEcB3Hc8G9LMED4GU3qxr2yMN8sIriOgxTnTJkThi2LE4b1Nuva8zRofm9taR8l/AizolHHTVKG/FGxAtbSFTebZW7p2LS7YUk9yVmGOypWwFoslqmiavyNU4gVsBaLZepYG6zFYrGMCWsisFgslnGgyfFrjYoVsBaLZapMKRfBRLAC1mKxTJ0ghpwPs4gVsBaLZbpMIdnLpOibi0BE7hKRj4nIF0XkCyLyc6b8mIg8LiJfNX9XTLmIyO+ZlRc/KyIvaKnrQXP8V0XkwZbyF4rI58w5vyfGy7pbGxaLJT0o0fMQJEUgD5LspQH8gqo+F7gPeJNZTfHNwEdV9V7go+Y9hKsu3mu2h4B3QCgsgbcCLyZcQOytLQLzHcBPt5x3vynv1obFYrHMPH0FrKpeV9W/N6+3gS8RLvb1APAec9h7gNeY1w8A79WQJ4BlETkLvAJ4XFXXVHWdcNnc+82+RVV9wizb8N62ujq1YbFYUsSQix7OPJFssCJyD/B84BPAaVW9bnbdAE6b11FXZTxvXreX06ON9n49RKgtA+yIyJcjfKwTwJ0Ix08T29f4SUo/Yfb7OlyOS4XgqEdyicg88OfAz6vqVmsyClVVERnrI6VXG6r6MPDwMPWKyJO9VqacJWxf4ycp/YRk9TUqSbGpRmWghNsikiUUrn+sqv/ZFN80w3vM31umvNuqjL3KL3Qo79WGxWJJCUozQ1q0LQkM4kUgwLuAL6nqb7XsehRoegI8CHywpfz1xpvgPmDTDPM/DLxcRFbM5NbLgQ+bfVsicp9p6/VtdXVqw2KxpAU92isafB/wk8DnROTTpuyXgbcB7xORNwJPAz9u9j0GvAq4BJSANwCo6pqI/BrwSXPcr6rqmnn9M8AfAgXgQ2ajRxtxMpRpYUrYvsZPUvoJyeprJJIiMKMiSZmNs1gs6WTpxHP1vh/6k8jn/fV7n//UrNukbSSXxWKZKpriSC4rYC0Wy9RJa8LtxC7bLSJ5Efk7EfmMCeH9D6b8mSLyCRN2+2cikjPlnnl/yey/p6Wut5jyL4vIK1rK7zdll0TkzS3lHdsYoM+uiHxKRP5ilvsqIt80ocufFpEnTdnMhUaLyLKIvF9E/kFEviQi3zuj/Xy2+S6b25aI/Pws9nUqpHiSK7ECFqgCL1XV7wGeRxgVdh/wm8Bvq+qzgHXgjeb4NwLrpvy3zXFIGPb7WuA7CUN0/0BCQegCbycM/X0u8DpzLD3a6MfPEUbCNZnlvv5TVX1ei41rFkOjfxf4K1V9DvA9hN/tzPVTVb9svsvnAS8knPz9wCz2dTocYTetWcWE4jbXlM6aTYGXAu835e0hvM2w2/cDLzNP+QeAR1S1qqrfIPR+eJHZLqnq11W1BjwCPGDO6dZGV0TkAvBDwDvN+171TLWvXZip0GgRWQJ+gNCFEFWtqerGrPWzAy8DvqaqTyegrxNBCZeMibolgcQKWNgbcn+aMADhceBrwIaqNswhrWG3e6G6Zv8mcJzoob3He7TRi98BfhFoPnp71TPtvirw1yLylIRhyDBDodGGZwK3gf8kodnlnSIyN4P9bOe1wJ/2qWdW+joZzKKHUbckkGgBq6q+GXZdINTinjPlLnVERF4N3FLVp6bdlwH5flV9AeFQ9U0i8gOtO42WNPbQ6D5tZIAXAO9Q1ecDu7QNf2ekn3tIaP/+YeD/GqWeYZlEG8NxtNMVzjxmaPgx4HsJh1NN74jWsNu9UF2zfwlYJXpo72qPNrrxfcAPi8g3CYfvLyW0H85iX1HVq+bvLUJb4YuYvdDoK8AVVf2Eef9+QoE7a/1s5ZXA36vqzT71zEJfJ4q1wc4YInJSRJbN6wLwg4STHB8DftQc1h7C25x1/VHgb8wT/VHgtRLO3D+TcILg7wgjzu6VcBY+Rzi0e9Sc062NjqjqW1T1gqreY+r5G1X9F7PYVxGZE5GF5mvCkObPM2Oh0ap6A7gsIs82RS8Dvjhr/WzjdeybB3rVMwt9nRwp9iIYKg/jLGzAdwOfAj5LKAD+Z1P+bYRC5xLhUMwz5Xnz/pLZ/20tdf0Kof32y8ArW8pfBXzF7PuVlvKObQzY75cAfzGrfTXHf8ZsX2jWRWjP/SjwVeAjwDFTLoQeDF8DPgdcbKnrX5l2LwFvaCm/aH6zrwG/z35EYcc2evT1ecCT5hr4L8DKLPbTnDNHOKJYaimbyb5OeptferZ+/w//beQNeHLafe+32VBZi8UyVUTkrwhz3Ubljqre3/+w6WEFrMVisYyJxNpgLRaLZdaxAtZisVjGhBWwFovFMiasgLVYLJYxYQWsxWKxjAkrYC0Wi2VM/P/Cr3AdR8RRVwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"wi[\"mean\"] = mean\n",
"cap = mean.max()\n",
"ax = wi.plot(column=wi[\"mean\"], cmap=\"coolwarm\",\n",
" figsize=(5,5), vmin=-cap, vmax=cap, legend=True)\n",
"ax.set_title(\"mean\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How do we weight the combinations of component rows to approximately construct a given row in our original dataset? We can do this with component weights, as given by the `PCA.transform` method."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
0
\n",
"
1
\n",
"
\n",
"
\n",
"
DATE
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
2020-03-16
\n",
"
-363.67
\n",
"
-2.28
\n",
"
\n",
"
\n",
"
2020-03-17
\n",
"
-360.99
\n",
"
-3.83
\n",
"
\n",
"
\n",
"
2020-03-18
\n",
"
-357.78
\n",
"
-5.66
\n",
"
\n",
"
\n",
"
2020-03-19
\n",
"
-356.65
\n",
"
-5.82
\n",
"
\n",
"
\n",
"
2020-03-20
\n",
"
-355.23
\n",
"
-6.59
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1\n",
"DATE \n",
"2020-03-16 -363.67 -2.28\n",
"2020-03-17 -360.99 -3.83\n",
"2020-03-18 -357.78 -5.66\n",
"2020-03-19 -356.65 -5.82\n",
"2020-03-20 -355.23 -6.59"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pc_coef = pd.DataFrame(pca.transform(dft), index=dft.index).round(2)\n",
"pc_coef.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Even though we have 2 number per day instead of the 72 per-county numbers per day, we can use those two numbers as weights on the component rows to approximately reconstruct those 72 numbers. Let's try it for a specific date: Mar 1st, 2021."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 -231.30\n",
"1 4.62\n",
"Name: 2021-03-01, dtype: float64"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mar1 = \"2021-03-01\"\n",
"pc_coef.loc[mar1]"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"