Skip to content

Instantly share code, notes, and snippets.

@darribas
Last active August 29, 2015 14:01
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save darribas/847138dced15727f9fcf to your computer and use it in GitHub Desktop.
Save darribas/847138dced15727f9fcf to your computer and use it in GitHub Desktop.
Contiguity and block -based weights in PySAL
Display the source blob
Display the rendered blob
Raw
{
"metadata": {
"name": "",
"signature": "sha256:b5ea7e1186ce3c370493393ef3359dde69ce0ad9a98dab2bdf8e3af59345b6ea"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"%matplotlib inline"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Combining contiguity with block structures in creating `W` objects\n",
"\n",
"* [Dani Arribas-Bel](http://darribas.org)\n",
"\n",
"In this notebook, we will work through an example on how to combine to rules to build spatial weights. This is inspired by a note by Yamil Velez to the `openspace` list on May 10th. 2014. His [question](https://groups.google.com/d/msg/openspace-list/zgZf_3Ijj0g/WcyaLqZ3wYwJ) was:\n",
"\n",
"```\n",
"I have read the relevant documentation and I still have not found a way to solve this problem. Basically, my goal is to compute the average value in contiguous cells for each cell in a raster file but only for valid contiguous cells. My raster file is made up of both land and water cells but I only want to work with the land cells. When I compute spatial lags, however, values from water cells are introduced. Is there any way to modify the weights so that water cells are automatically given a zero and not factored into the spatial lag calculation?\n",
"\n",
"Thanks so much!\n",
"```\n",
"\n",
"Here I present a potential answer with an example using PySAL data examples."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import pysal as ps\n",
"import numpy as np"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"db = ps.open(ps.examples.get_path('columbus.dbf'))\n",
"db.header"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 10,
"text": [
"['AREA',\n",
" 'PERIMETER',\n",
" 'COLUMBUS_',\n",
" 'COLUMBUS_I',\n",
" 'POLYID',\n",
" 'NEIG',\n",
" 'HOVAL',\n",
" 'INC',\n",
" 'CRIME',\n",
" 'OPEN',\n",
" 'PLUMB',\n",
" 'DISCBD',\n",
" 'X',\n",
" 'Y',\n",
" 'NSA',\n",
" 'NSB',\n",
" 'EW',\n",
" 'CP',\n",
" 'THOUS',\n",
" 'NEIGNO']"
]
}
],
"prompt_number": 10
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The variable `EW` is a binary that differenciates eastern from westerd polygons."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"ew = db.by_col('EW')\n",
"ew[:5]"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 11,
"text": [
"[1.0, 0.0, 1.0, 0.0, 1.0]"
]
}
],
"prompt_number": 11
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from pysal.contrib.viz import mapping as maps\n",
"maps.plot_choropleth(ps.examples.get_path('columbus.shp'), np.array(ew), 'unique_values')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADtCAYAAAAcNaZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd8FHX+/58zk+xuskk2DVJIA4JJCL0qAlJEKYIiYrmz\nYMF6Kj+97xXEs93Z9exyevZ2egqKhyBIP1APpRMQQVoghBDSt2Vn5vfHbNYEkpBszcI8H4957O7s\nzOfzmWT3te95f96f91tQVVVFR0dHRycoiKEegI6Ojs6ZhC66Ojo6OkFEF10dHR2dIKKLro6Ojk4Q\n0UVXR0dHJ4jooqujo6MTRCJCPQAdnfZQXFxMfm4uQxwOFPBscqPnLb1WgRKjkaUrV3L22WeHZPw6\nOrro6oQV6enp1CsKRrz78NYajUiS5O9h6ei0Gd29oBNWiKJIl9RUrF6er4IuujohRRddnbCje/fu\n1Hl5ri66OqFGF12dsKOgVy/vLV1VPUl0VVXF5XLhcDjQV8XrBBrdp6sTdhT07MmiqCiw2dp9rupy\nMXjQIFBVZEVBVhQURUEUBARBQFFVjAYDUUYjJpMJc1QU0dHRmM1mzDExxMbGEhcfzx/+9Cd69uwZ\ngKvTOd3RRVcn7MjNzcVpMHgluoNsNhRAOHFTVVBVLdLB4dC26mpkwIUW/VABHANKIiMZPHSoLro6\nXqGLrk7Y0b17d2pcLq/OldxbS4juLbKVY+pUFavVWweHzpmO7tPVCTuysrKosduRQzUAl0sXXR2v\n0UVXJ+yIiIggtVMn2u9c8A8iUFdbG6LedcIdXXR1wpJu3bp5HTbmKxJQW10dot51wh1ddHXCkoLC\nwtCKbl2oetcJd3TR1QlLCgoLqTcaQ9K3BFh194KOl+iiqxOW5ObmUm8yhaRv3aer4wu66OqEJbm5\nudTIoYlfkECPXtDxGl10dcKSnJwcqtwLHYKNBNi8WJihowO66OqEKUajkaT4+JCEjUmAVRddHS/R\nRVcnbOmakxOSCAYJsOuiq+MluujqhC35hYVeZxvzBRGwOxwh6FnndEAXXZ2wpaCwEEdka1kSAoOE\nLro63qOLrk7Y0qNHD1xRUUHvVwIcTmfQ+9U5PdBFVydsyc3NpTYEScd10dXxBV10dcKWbt26UWmz\nEWzZFYF6lwtFCUXAmk64o4uuTthiNpuJNZuxB7lfAYiMiMBuD3bPOqcDuujqhDU5WVkhCRuLlCR9\ngYSOV+iiqxPW5BUUhEx09aXAOt6gl+vROSWKolBaWsrBgwc5ePAgxcXFnucHDx6ksrKSxx9/nClT\npgR9bD179+aHzz6DIOdhiBBF3dLV8QpddHVa5IYbbmDFihUcPnwYi8VCZmZmk23AgAFkZmZis9m4\n8cYbWb9+Pffddx+mRtm/HA4Hmzdv5vDhw57jBUFoU/+qqmKz2aiurvZsNTU1OJ1OnE4nDoeD0tJS\n6oxGCJDVqQA/GY2oRiMKoAgCClBdV4fLyzptOmc2uujqtMjw4cP5/PPP+f777+nXr1+rx/7www/c\ncMMNJCYm0rVrVxISEqipqWH37t3k5ubSpUsXNmzYQEVFhUd0BUEgMjKSiIgIIiIimoix0+mkpqYG\ng8FAXFycZ4uJicFoNGIwGDAYDFitVuoiAvcxrgcOCwLvvfkmJpMJk8lEVFQUMTExFBQUBKxfndMX\nQVVDEOioEza88847zJ49myVLllBYWHjK4x0OB0VFRdTW1hIdHU1+fj5msxnQLFdHo5VciqLgcrmo\nr68/yWqMiIggLi6OyFOsOKusrKRz586Mq68PyASFDdgUH8+xiooAtK5zJqJbujqtct111xEREcH5\n55/PunXr6Nq1a6vHG41G+vfv3+x7giA0cT34g/j4eIwGAzX19Vj82rKGAhhCsNRY5/RFj17QOSW/\n/e1vGTJkCBs3bgz1UJolOzubQNmhCpzS2tbRaQ+66Oq0CZfLhcFgCPUwmiU/P5+qALWtAMYOet06\n4YkuujptwuFwYAxRIchTUdirFzYxMB9lBTrsj41OeKKLrk6bcDqdHVZ0c3NzEePiAtK2Ah32unXC\nE110ddqEw+HosBZfbm4uchtjf9uLbunq+BtddHXaREd2L+Tm5lIbwMURxhCVetc5PdFFV6dNdGT3\nQnJyMoIgBCQHg4zuXtDxL7ro6rSJjuxeEASBzMxMygPQtgJ+jy3WObPRRVenTXRk0QXIy8sLSNiY\nAphCUBJI5/RFF12dNtGnTx/eeuutUA+jRTIyMigBNgBbgO3AT8BuYB9QDBwByoEqwAo40US1NRQg\nShddHT+iLwPWaRNvvPEGAwYM4LzzzmPUqFGhHs5J1NbWEiEIRIgiMpqgyoCsqk0elYZHQHVvQqNN\nRHNXNDyqqkr6kSOhuCSd0xRddHXaRFpaGm+//TZXX301GzdupFOnTqEeUhNycnIwSxJ92pluUUUT\nYNm9KWgC3fB4GEhPT/f3cHXOYHT3gk6bufDCC7nssst49NFHQz2Uk8jLy8PmRaFIAa26rwGIAsxA\nHBAPJLlfN2RJ09HxB3pqR512sX37dsaPH8+ePXs61MTahg0bOHvgQC70c7u/APHnnMN1M2aQkJBA\nYmJik8e4uLg2J2XX0QFddHW8YPz48UyaNIk777wz1EPxYLfbiYqKYiL+vX2rBQ5JEqrBgBIZiUsQ\ncKoqDlnGVl9PvctFdFQUlthYLHFxJCYmkpicTOfOnemUmkpycvJJYp2UlERqaqofR6kTTuiiq9Nu\nNm/ezAUXXMA///lPJk6ciCRJoR4SAJGiyHmqSjBjDRS06hKNN2ej52pkJKrBgCxJ1AsCTkWhwmZj\nxcqVnHvuuUEcqU5HQZ9I02k3ffv25ZVXXuGRRx7hxhtvxGQyoSgKBQUFDBgwgIsvvphhw4YFfVzG\nyEgcTmdQRVcEjO6tWerrta0R2+LiOH78eIBHptNR0UVXxyumTZvGtGnT2L9/Pw03S0VFRfz4449c\ndtll3HXXXfzxj38Mqr/TZDJhczqJD1qP3iG6C27qnJnooqvjE9nZ2Z7nOTk5TJw4keuvv55LLrmE\n6Oho7rrrLr/0c/z4ccrKynA4HIiiSFZWFnFxcWzZsoVHH32UNWvWcLy6GhOQ5pceA4egKFgDlKBH\np+Ojh4zp+J2MjAweeugh5s+f71M7f/nLX8jPzycxMZGUlBQGDhzIiBEjGDZsGImJiRgMBs4++2xk\nWeaVV17h7lmzcJy62ZAjyrJu6Z7B6JauTkAYPXo0V111FVVVVVgs7SsZqSgKY8eOpaioiEcffZTB\ngwfTs2dPIhqVWldVFavVSkREhCcL2L59+xCjo6GjW5G66J7R6KKrExCio6MZPnw4S5YsYfr06S0e\npygK1dXVlJWVcezYMfbt28dDDz2E0Whk69atdO7cudnzBEE4adFC586dEYzGji+6LpfuXjiD0UVX\nJ2BccMEFLFu2zCO6K1eu5NJLL0VRFOrr63G5XNTX1yOKIkajEZPJhMViYcyYMbzwwgtER0e3q7/O\nnTujhMFCBVFVqautDfUwdEKELro6AWPkyJG8/vrrntdmsxmHw8GWLVuIjo7GbDYTHR3dxG3gC507\nd6Zelv3SViAR0RL06JyZ6KKrEzD69u3LoUOHKCsrIyEhwRN9EBcXF5CEOSkpKTgcHX8qTQJqa2pC\nPQydEKGLrk7AkCSJ7Oxs+vTpQ0VFBfHx8UyZMiVgORuSkpKwOxz8hPbBlhptEc08RqBZncEO4ZFA\ndy+cweiiqxNQ+vXrhyRJPPnkkyQnJwe0L0mSEIFjoiajCk1TNSqq6knl2PAITfPpejZBaPpIozy7\n7g1AUFUEVSVSVRlE2wRcAqx1gajophMO6KKrE3C6desWcMFtwBgZSS+nk7YGqTWI70lbI6Fu9v0T\nth2ACy1F5KkQgboTRLeyspJpl1xCVUUF9fX1OOvrcTU8ulzUu1y4GjZZJjoqii7p6WRlZVFdXc3+\n/fv5y8MPc8MNN7TxynVChS66OgFj48aNfP3117zwwgtB69NsNmNvh+g25NP1NWXPLrREN20RXQGw\nnxCnu2zZMnb+8APZdXWYaFTFgl9dII33uRwObJWVHCoqIgJN+IuLi328Cp1goIuuTkBQVZV77rmH\nBx98sN2LI3whITERe0VF0PprQELLKnYqFOCA2cxfrr++yf4lixdjqasjqY39RUKTxD6VkkRkZCQO\nhwOr1YrVaqWurs7zvPGmKAqiKCKKouaSaeZ549cWi4UBAwa0cWQ6p0JP7ajThAkTJlBUVERMTAwx\nMTHExsYSExNDbm4uTz/9dJvbWbBgAbNnz2bTpk2thoQpisLKlSvZunUrd999t2f/rbfeyrvvvntS\nwhxJkoiOjkaSJARBwGg0YjAYPMeVlpYSdfw4/dt53b6yTBAoVFVOlSX3F0kiccgQVq9d2+Taumdn\nk37ggNfJenaLIjsVhQhJwhARQaQkESmKRIgiEW5ftOT2PwuqitrQt/tRdT9vqBvXsE8FquvrmX7l\nlbwydy6RkZFejlCnAV10z0BUVWXt2rVERkYSHx/v2YxGI5mZmXz88cfExcVRW1tLbW0tNTU1XH75\n5RQXF5OSktKmPh544AFcLhd/+9vfTnpvx44dPPDAA2zatInDhw8jyzKyLLNixQrOOeccRFEkNzeX\n//u//2Pq1KlNxm2z2Thy5AgOhwOXy0VNTQ01jcKvVq9ezYJ//YshQY4OWCkIdFdVMls5phrYZDaz\ntaiIrKwsz/6ysjKyMzIY7XR6HUnRIJCBiMSoB4qiokjr2ZMvv/qqxVWCOm1Ddy+cgdTU1DBixAiG\nDBlCRUUFlZWVVFRUEBkZidPppKCggISEhCbnXHzxxeTn5yNJErm5uRQUFDB37lxP3oMTyczMZN26\ndU32ff/999x9991s2bKFyZMnM3v2bEaMGAHA9ddfz+TJk6mrqyMhIYGKigomTJjQ7Bc8JyenxWtL\nSEhgwbx57fyL+E6EIFDfiv2iADvMZp557rkmgguwatUqUoxGRKfT6/4bR1T4m0igj83Gni1b6Nur\nF6+98QYXXXSRXqbIS3RL9wxEVVUyMzOZP38+gwcP9uyzWq04HA4SExNbPO/IkSPs3buXP//5zwwf\nPpwRI0ZgMBjIysoiPj7e4we87777eOutt4iJicFsNlNdXY3dbufaa69l9uzZZGRkNNtHRUUFRUVF\nWK1Wxo0b1+5r27BhA+ePHcu5lZXtPtcXvpMkamWZREFggPsrtUwUcakqkruUe25eHluLik4Sq5k3\n3si6N9+kW1BH7B2lwH6zGUtKCg//7W9Mnz69w1QOCRd00T1DmTt3LvPnz+frr7/26vyNGzfywAMP\n4HQ6cTgc7N27l5qaGhRFQZZlVFXl4Ycf5txzz6WkpITOnTszZMiQgH9Bjxw5Qrdu3Rgb5CxexYJA\nJXBQVZng3vcVcA7abX8lEDtwIN/98MNJ53bLyqLLwYMdPvl6AypwFCiOiSHSYmHu668zYcKEU52m\n40YX3TOUBjfCO++8w/Dhw/3attVqJTExkf3797fZB+wvZFnGYDBwgaIE3XemAosAC2AVBFBVzscd\n4gWsNBh4Ze5crFYrNTU1VFVVUVJczEcffcQ4WQ6YeyBQqMBBIKJ3b37csiXUwwkbdNE9g3nxxRf5\n9ttv+fDDD/3a7oUXXojD4WDlypV+bbetJCYm0rOigoRTH+p3dooigqqSrKok0HRia5/BgC0yElGW\nob4eUZYpBQySxJAwSNTTHNVAcVYWe/bvD/VQwgZddM9gKioq6Nq1K7t37/bbirH77ruPuXPnsmnT\nJjIzW5vLDxw9e/YkcscOsk59aMhZLUlkyjJdQz0QL7EBmxMSKNMLbbYZvVzPGUxCQgIjRoxg9erV\nfmnv/fff5/nnn2fx4sUhE1zQygWFQzoZBbDKMv7PtxY8IoEaPXlPu9BF9wwnOTmZqqoqn9v54osv\nuPXWW3nvvfc8ERGhIjs7m3Coy3AM7QtoPtWBHRg74HS5kMPUPRIKdNE9w4mLi/NZdJcuXcpvfvMb\nXn311SaLGUJFTk4OrjAIYzoApEhS2E2gNVAL/BgVxT/mztXDxtqBLrpnOBaLxWfRvfvuu/njH//I\nNddc46dR+UZ6ejqiuePbj7WiSHKYWoilaIL7wssvM/Pmm0M9nLBCX5F2hpOYmMju3bu9Pn/9+vXs\n2LGDQYMGsWPHDiwWC/Hx8URFRYVsxVJ6ejqqn0oABQoXYFUUgpPw0r/sMZmoSUhg4ccfe1YU6rSd\njv3J1Ak4gwYN4v333/f6/JKSEuLj47nxxhtxOByexRKKohAdHU1MTAwWiwWLxUJiYiKJiYlYLBZP\nIp3o6Giio6MZOXIkhYWFfrmm9PT0Dl8rrRiIEgSMYRg8ZDUaee2NN3TB9RI9ZOwMx+FwkJGRwcqV\nK/0meqAVXjx48CDFxcWUlJRw5MgRysrKKCsro6amhtraWmw2Gw6Hg/LyciIjI9mxY4df+i4vL6dL\nly6M68D10r4F4kWRAkU55bEdje1xcfz93Xe5+OKLQz2UsES3dM9wjEYj9913H3/4wx9YuHCh39qN\niYmhoKCAgoKCUx7bsILt6NGjfslglZiYiCzLbU4qHgqsokhuGAouaCWKwqEAaEdFn0jTYcaMGaxZ\nsyZk/UdHR9OlSxeWLFnil/YEQSAhIYFqv7Tmf+yAXVFoPq1Qx0dUVZw+ZEQ709FFV4e4uDisViv1\n9W2pfRAYhg4dyoIFC/zWXlpaGh21yPl+IEEUfS4RFAxUtPHudm9rzWb22mxER0eHdmBhjO5e0EEU\nRcxmM3V1dcTHhybX1a233srkyZM9pWR8JTMzk+0dNAlLmSiSGiauhWKgLjubi6dO5fixY9z2u99R\nWFhITExMqIcWtuiiqwNot/hWqzVkojty5EhEUWTz5s307+97sZ3s7GxOTqLYMbCqql9DxRQ0a1RA\nW5bbsBnRvuAuQHYf1/C84XVLzxteV5hMfL9woV8nWc90dNHVAbQqulZraBfPduvWjcWLF/tNdGWD\nATqY77ECkFW1zdWK20IxsBOIcSdNl919yGjugcZVhAVAdNdME1t7DlQqChdfcokuuH5GF10dQLN0\ni4uLyc3NDdkYLrroIubPn8+f//xnn9tKT09Hio7ucKL7M7/m3W2NtsRxNj6mhyiS11aXReMo0VYi\nRtfHxnJXo2KhgURRFKZPm8bnX3xBgsWCJElIogiCQGREBN26daOwTx/Su3QhLy+PadOmBWVcgUCP\n09UB4Oabb+bjjz/mgQce4J577gnJGEpLS8nKyqK8vNxnn+Hy5cu5Yvp0zu5gKQdXiSKZikIXmq9p\nduI+oZX3Gu/zd400B7AuKorK6upWqzk3R319Pdf85jfs3rWLCVOmkJ2dTf/+/TEajWzdupXDhw+j\nKAp2u52yo0f55eef+XHDBpy1tRTabBigSVViBahzby5RpMRg4Nv16+nVq5cfrzh46JauDgCvvfYa\nc+bMYcCAAVxyySV06xb8il0pKSl06tSJFStWMHnyZJ/aSk9Px+ly+Wlk/kFBW/qbiuZv7cgcA4YP\nG9ZuwZVlmSunT2f9kiWk2mx8vHUralQUVZKEqqrECQKS3a5Z2LJMhKpiAs4CYmn5hyO24YmioAoC\na9asCVvR1UPGdDxkZWVx55138sADD4RsDH379uU///mPz+2kp6djs9noSDECR9AmuMIh2Ko6KorJ\n7cwYpygK1159Nd8tXUovm40UIE9VybdaGVpTw9m1tfSsqSGvvp48l4s8VaU70AWIo+2Wuuh0sm/v\n3vZdUAdCF12dJtxzzz0sWbKErVu3hqT/6667ji+//NLndmJjYxEEgY60bqoYSPFDOFygUYFjgtCu\nasyqqnLrzTezYsEC+litAY1BTpNl/vHqq4SrZ7TjfwJ0gkpsbCyPPPIIv/3tb6mrqwt6/5dddhlV\nVVXs2bPHp3YEQSA5KalDrUqrE0WSwyA+tw4wRkXRo0ePNh2vqir/7+67WfDRR/SxWgPus4wAJEkK\nWRY7X9FFV+ckZs6cSf/+/Zk5c2bQrQlRFMnOzva6NHxj0tLSOkzZHidgUxSSQj2QNlAGXHDhhW0W\ntTmzZ/PhG2/Q12olMrBDA7T44SiTKQg9BQZddHVOQhAE5s6dy86dO3nmmWeC3v+YMWOYP3++z+1k\nZWcTfFu9eQ6ixdF21AQ8jamLjWVSGycyn33mGea+8AL9rNagXZsMflm1GCrCd+Q6ASUqKop58+bx\n6quvcvvtt3Ps2DGqq6tZsWIFY8aM4ejRowHr+3e/+x3//e9/vUqqoigKlZWVHDhwQFvwEYDxeUMp\nkBIGPkgFKHU6GTt27CmPfe+993j4/vvpZ7UGNRqjxGTit9deG8Qe/YseMqbTIjk5OWzYsIFbb72V\n3NxcXC4XXbt25ciRIxw8eNAvaRibIy4uDlEUGTt2LImJiTidTurr66mvr8flLoJYUVHB0aNHkWXZ\nszW8J0kSkZGRqKqKyWiEDpCG0CqKnBUG/txKICsjg06dWq9R/NVXX3HHLbcwwGYLejSGGhFBQc+e\nQe7Vf+iiq9MqFouFjz76qMm+0aNH+6WCcHNMmjSJb775BovFwp49e+jXrx8GgwGDwYDRaMRgMBAZ\nGckbb7zB+PHjefzxxzGZTBiNRkwmEwaDwXPr+fHHH3PXrbeGXHTrAKeikBDSUbSNclFk4ilcC2vX\nruWq6dPpY7P9Gj8bROpFkYSEcPhrNo8uujrtJi4ujupq7+ICvv/+e/bt20dlZSXV1dXU1NR4KklU\nVFSwatUq9uzZw3//+1/mzJnDiy++2Gw7H3/8MYMGDSI7O7vFvtLT01E6wAz3fiBJkpA6eAkhgNqY\nGCZMnNji+9999x0TL7yQAqs1JD8iKlBmt3POOeeEoHf/oIuuTruxWCxeie7PP//M8OHDycrKIioq\nylMfzWw2YzabiYuL4/333ycjI4MBAwZQVlbWYltHjx49ZZ6I9PR0HCHMEdxAuSjSJQwEtx4ot9sZ\nPnx4s+9v2bKFC8eNI6+ujtadD4GjEkhLTSU5ORxLemrooqvTboxGI3a7vd3nvfvuuwwePJh169ad\n8tjc3FycTid79+6la9euJ71fXl5Ofn5+q22kpaVhda9KC+WMcbhU/S0HBvTtS1RU1Env1dbWMnXy\nZHJqa0kJ/tA82IC8vLwQjsB39OgFHa9oawzn888/T1ZWFomJiTz++ONcfvnlbTpPFEXy8/OZN2/e\nSe85HA7sdjvdu3dvtY3o6GgMBkNIIxiOuR9D4ftsL5UGAxc1s/R39erV5OXm8suBA2SEYFyNiYRW\n74DCAV10dQLKs88+y6233sqaNWsoLy9n1qxZbT532LBhrFix4qT9R48eJTo6uk3JWDp36hTSsj0H\ngc6S5NcMYIGiymjkwgsv9Lyuq6vj9ltuYfL48aSXliKhLfIIJRZgm5+qRocK3b2g4zfsdjtPPPGE\np96a3W7n6NGj3HXXXV6lahwyZAjLli07aX9paWmzt8DNkZ6eTvmBA6S1u3f/UCVJdA8Df64dOFZb\ny4IFC9i6dSsmk4nfz5qFoaqKs93pFiMEAYc7K1ioiAQUWaa2tjZsSwbpoqvjFc0tD77xxhtZuXIl\nhYWFREREEBERwWOPPeb1l2PAgAHNLsI4cuRImwsjZmVnU/zdd1717ysKYJXlsPDnGoDeqspHjzyC\nHB1NvSCQWlNDaqNjGkQ3lChoEQymMF4GrIuuTrsZOHAgzz33HBs3bsRkMnm2b775hvPPP5/zzz+f\nmJgYYmJiiI2NZcuWLVgsFiwWi2fhQ1soKCigtraWw4cPk56e7tlfXl5OWVkZffr0ITY2lvT0dD74\n4AMMhpMXoubk5HCyrRwcDgEmQSAqDFaiiUAWgKJAbfMZK6QOkLXNCqQkJ7c7z29HInxHrhMybrrp\nJh588EF+/PFHUlJSsNvtOJ1OsrOzKSoqYsOGDTgcDhwOB06n86QVZZIkIUkSERERTR6be24ymXjj\njTe4//77Pf1PnTqV2NhYjh07RllZGX/729/43//+12yoU0ZGBqLZDCHImHYYSBGEVkvihBOSooRc\ndA1ARWUlx48fJzExMcSj8Q5ddHXaTUREBPfeey9bt27l3Xffbde5iqLgdDqx2+3Y7XZPJMKJm81m\nw+Fw8Prrr/PDD03r+sbFxXHppZd6Xr/88sstpqFMT09HMBpDIrp1okhOGCz9bSuRqoozxD8iRqCz\nKPLZZ58xc+bMkI3DF3TR1fGKG264gdzcXMrKyk65Tr8xoih63BFtoaqqiueff77VYwwGQ4uVjNPT\n03GFQCScgD1MUjm2FRNgF0UI8cSgKAjNupLCBT1kTMcrkpKSmDp1Km+88UZA+xk4cCClpaWtHnMq\n0bWHoCLwPiBOFE8rq8YIIZ9IAxDr69m3b1+oh+E1ejVgHa/58ccfmTZtGnv27EGSAlOgxeFwEBsb\nS0lJCUlJzduNvXr1YtasWdx0000AzJ8/n+eee466ujrq6uo4sG8fOJ1NF3S4P/aNP/zNfRFO3Nfc\n16W5NhRVJQUY0Mq1hRsrRZE0VSUvxJJRDlT16MH2XbtCOg5vOZ1+iHWCzMCBA0lNTWXhwoVMmTIl\nIH0YjUZycnL4/PPPufHGG5s9xmAwYLPZPK8feughunbtyqWXXkp8fDy33347OYqC2f2+0MKjP99b\nLwhYTiN7Zh9aprTW1wAGBzOw/fBhSkpK2LVrF7t372by5MkBSzXqb3TR1fGJO+64g5dffjlgogtw\n9tlns3Tp0jaLbllZGc8++yxjxowB4KEHH8S0b1/QcgbIgFVV6RKk/gKNAvwsCPRU1aALhhOtfFAF\nUAPUSxIORcFZV8dZ3boRbzRSY7ejKErYTKzpoqvjE9OnT+f3v/89u3bt4qyzzgpIHxkZGbz66quk\npaWhqqrnFr/heX19PY899hiPP/44iqJQV1dHr169POenpqVRHkQfYCVgFEVMp0nkwlY0f26gfkQU\noArNbVCJNllXDzgURauHJgjEiSIJsoxZlolBs3YNdjvY7ewC9oZRSXZddHV8wmQyceWVVzJ+/HhS\nU7X1S43xLsRLAAAgAElEQVRFsfHjqfYByLJMVVUV8fHxnn3Hjx9HEATmzZuHIAiIooggCEiS5Hne\nsAJOkiTMZnOTW82MLl0oDsC1t0SZIBB9mrgW7Gjxxmerqs/5I+zAEaAaqEWzWp2KgtNtQceIInFA\noqJ4hDUaEFS11YiJKODnnTt9HF3w0EVXx2euv/56XnzxRa6//noiIyM9E1aCILT4vKV9Bw4c4O23\n3+aee+7xtK+qKvfeey+FhYXExcW1e3yZWVms9Prq2oeKJlLdThPR3SgIpAoCCW202hXguHuroqnV\nqqD5vFMkiaQTrNZI0FbDtYNatCxux6OiiDtFhEtHQhddHZ/p168fQ4YMYfDgwYwfP96nttatW8dH\nH33E7373uyb7//73v7NgwQKuvvrqdrfZpUsXMJnAixzA7aUcLRl4VsB7CjzlQKWqMqqZH5A69/sV\n7udOt9Var6pEArGiSCyQ7LZaY9DiU5cDA/0Q56sCa0SRyy69lAkXXcQFF1zgc5vBQhddH2lcDFGS\npDbnmT3duOqqq/jkk098Fl2n09lsboahQ4eyePFir0Q3LS0NMUiiu18USVaU0yIAfosg0ElVOYBm\ntTokCaeq4lQUVMAsisQKAsmyTIwsY0azWiOgWavVn0uIBSDaaOTvL7xAWlqocsh5x+nw2Qgpffv2\nxWw2YzKZEEURURSJjo5mV5jGEHpL9+7d/VKWvSXRHTZsGNu2bfOqzdTUVJQ2JtnxBRtQqiiEY53a\nGmA38D9gtSiyBC0CowKokCTMokimLNNXURgNTADOUxQGyDJnAelouW6DacVFR0aecuFMR0QXXR9J\nSkpiyZIluFwuVFXF5XLxm9/8hsWLF4d6aEFFEARkP9w2Op3OZhdaDBgwgJKSEq/aTE1Nxely+Tq0\nVlGBHaKIRRRpW6bf0KAApcA2YB2wQpJYBKwBjogi0aJIuqJQD+QD44CzZZlCRSEHSEZbDuzN/Zy/\nvdwmQdBF90ykoKCAHY0y2YuiyNixY5tNvn06061bN37++Wef26mvr29WdPv168fx48dbXO7bGqmp\nqdgC7FrYJwgcU1UGd6AwMTvaooYf0fyf34gii4AtgE2SSAYKZJnz0CzX4YpCoaJQLookCEKHWAjR\nGiqEZYpHXXR9JD8/v4noAowZM4ZVq1bhCrB11ZE466yzKCkp8bo0ewMtWbpms5n09HT+85//tLvN\n+Ph4LbuZTyNrmQOCwE5VZaCq0lHSsLiAlcABUSRSkshUFAYoChegWa+D3W6BVNxhWe7z9gkC1arK\nkDCIvpChzcnsOxLh9zPRwSgoKGDhwoVN9qWkpJCZmcmGDRsYMmRIiEYWXCRJIiMjg/POO++UlSLu\nuecepjZTABFaFl2AwYMHs2jRojYXt2xAEATi4+OpOXbM71m/9gsCO1SVQdChMoptBBJEkaHtsLyr\ngZ2qymDCQxiqbTbuvuMOIiQJm82G3Z0m1OF0IssyMWYz8RYLCYmJJCUn06lzZzqlpDBs2DDGjRsX\nsnGHw9+2Q3Oie6GBMWPGsHz58jNGdAGSk5Pp1q1bqyWy16xZw9NPP+2V6J577rm89957Xo2tc+fO\n7RZdBS0cqiGY34o2WSZLEi7AIcuo7kUDGwUBAXfsMZz8vGFrWGSgKAhoPlJ/r+OrQ4tfHd4OwZWB\nHwSBLqoasPJCKt75glvCUV+PbeNGzGgLJBrC0iR3P67ycurRIi+Oofmxa4DPCwvZ6OWkrD/QRddH\nMjMzqayspKqqCovF4tk/duxYXnzxRf70pz+FcHTBJS4ujunTp7eah+HNN9/k6aefbvH91kR3wIAB\nPPnkk16NLS0tjb1FRU327UX7EjrRbsddgoAsirhUlXr3ElQRMAoC0e4JsiRZxiTLmNBWV9WKIn0V\nBVVVUdCyiynQZFNPeN2w1QEHRZGz/OwH3iSKZAKx7Wi3SBSR0OqkhQuq+weiPRX4jgO1Ifa766Lr\nI6IokpeXx08//dTEqj3vvPO4+uqrOXLkiGd57OlOv3792LhxY6ui26VLF2pbqMEFmui2NDnSv39/\nysvLsdvt7S5MmJGRQeP7kf8JAlWqSrIkYVJVJEXBoKoeQW3YJNDSQDYTmVEFWAWhXV/6xpQDJYrC\ncnc4mwDQKM5bwDvr0C7LDGrH8dVAsaJwnhd9tQcRzaJe647/zQGffOAq7Z+UEtCql4QSXXT9QIOL\nobHoWiwW7rrrLm699Vbmz59/RiyaGDp0KK+99lqrx6Snp7dYWgdaF924uDg6derEokWLWnRPtNav\nA83C/F4QsAEjAJMPYW4SmmXrDVa0qILOQJZ7sQFoQtJcft729GJGS1DTVgS0opOBzhdhBIYDR4FS\nUWS3ohAjSYz08n/Q8KPUXrz9n/kLPXrBD7Tk173//vvZs2cPH3zwQQhGFXx69OhxymxP6enp1NbW\nsn///mYtjtbcCwBDhgxh0aJF7R5bWloaYlQU37kr2p6rqvhaxFvE+9jTzaJIAlqS82Sgk3vrDKQ0\n2lLdW1o7tvZmpzBC0EoaWYAeqspwt2Vd58OPniAItNdmFWg+EX0w0UXXDzQXNgZaAu533nmH//f/\n/h8vv/wy9iAsQw0lXbp04dChQ60ek5iYSEFBAfn5+UiSRI8ePZq8X19fT2RkZIvnDxs2jA0bNrR7\nbMnJyVS7XLiAYaraLkuwJby1tEAre9NRnE6RaNcR7KJGvt5mC9Bu0YXQuxd00fUDBQUF7GwhtdyA\nAQNYtGgRX3/9Nd26dePZZ59t9fY6nLFYLFroTis/LoIgsGnTJmw2GwsXLqS+vr7J+3a7vdWigwMH\nDuTw4cPtHtszTz2FyeXiHD/G0op496VvfH5HQEAT3poQ9O2LzemtpauL7mlAjx492L9/Pw5H8yk9\nBg0axIIFC/jqq6/47rvvyM7O5rrrruOzzz6jpiYUH/XAcODAATp16tTmSS6Xy3VSngWHw9Gq6DZM\nprVn4cnNM2eyfdMmhrozYPkLCe9vVf0dPuULDRIUbEvX1+v31tJVddENfwwGAzk5OadcBtuvXz8+\n+eQTz6KJ1157jS5dujBhwgSefvppVq1a1erMfkdn8+bN9O3bt13nnGh12O12jMaWb/4TEhJITExk\nyZIlbWr/lVde4a1//pOhAVgt5oul25FE9wAgiWLQyhn5igstXE92h+e1h47g09WjF/xEz549KSoq\nalImpiWysrK44447uOOOO6iurubrr79m9erVfPrpp2zdupWuXbsyePBg+vfvT3p6OkOGDCErq+Nn\naN28eTN9+vRp8/H9+/fn6NGjKIrisXgdDkeroguai2HhwoVMnDix1eNWr17N3b/7HQOA2DaPqu34\n4tPtKKJbD+wE+oQgHeWprt+Ktqihgl8rTTgUBZeqYhIEEkQRczsn4hQgspU7qWCgi66fKCwspOiE\n4Pu20LCgYPr06YA2e79t2zbWr1/P5s2b+frrr5kzZw47duzo8GFn33//Pddee22bj8/MzCQmJoan\nnnqKvLw8rFYrRUVFFBQUtHresGHD+OKLL1o9prKykvHnn08+WkRAIPAleiFBFNkgy2QAbf+Z8j97\nRJFoID1Et9wqmqg2ro/mRKs00ZCzN04Q6OzO2dtQaUI8RQmfllAA0yl+1AONLrp+omfPnsybN8/n\ndgwGAwMGDGDAgAGAdiuUnZ3N9u3b22RFr1+/nuLiYrKyssjKyiI5OTkoYm2321m9ejXvvvtuu86b\nNm0ajz76KElJSRgMBiIiIjw/QC0xaNAgXnnllVaPGT1yJBZZJjuAt5K+xOkOkGX2AwckySvx8Ac2\nYK+icHZIev/1B+t//FofLalRpQkTIPj5x0CGdi+s8Te66PqJnj178te//tXv7QqCwLRp0/j0008p\nLCxk9+7d/PLLL5SWliKK4kmVFG6//XbMZjPV1dUcOHAAq9VKRkaGR4RP3DIzM4mK8j0D7Jo1a+jd\nuzeJiYntOu+pp55i2bJlxMTE8N1337U6idbAgAEDKCsrw+VynbSQYsOGDUyfPp29v/zCWAK/wsoX\nSY/Et+gHX/lJFIlXVRJC5ONsuPYLod310Xzp0+yHz7sv6KLrJ/Ly8jh06BC7d+8mNzfXr21fdtll\nXHvttXzxxRccPXqUwsJCFEVh2bJlTUTX5XKxfft2SktLiY3VvJh1dXUcPHiQAwcOeLbVq1d79h08\neNCTNrG1LTU1tdX42Q0bNnD22e23mWJiYli/fj0XXXQRmZmZvP3220yYMKHVc5KTk4mLi2PVqlWM\nHTsW0NwyM2fO5NNPP+XOO+/kiSeeCHiaRb+IbogErwptCfLokPSuEYpZfBn8YmT4gi66fsJkMjFn\nzhxuu+02lixZ4tdb+nPOOYeMjAxuvPFGrrnmGgRB4MCBA5x77rlNjquqqiIyMpIrr7ySq6++mosv\nvhiz2Ux+fj75+fnNtq0oCsePH+fw4cNNtu3bt7N06VLP66NHjxIfH9+sIKelpfHtt98yadIkr67P\nYrGwatUqXnrpJaZPn855553HW2+91aSM+okMHDiQBQsWkJqayu9//3vWrl1Ljx49+Pbbb+nVqxdP\nPPFEwCeqJHwT3Qg0EQg2KrBdFElRFJ9X5YUbCrronlbMmjWL999/n3fffZfrrrvOb+2KosiqVaua\n7Gtu5VZSUhKHDh3i888/59133+X2229n6tSpvPLKKy36sURRJDk5meTk5FYjD2RZpqysjMOHD1NS\nUuIR440bN7Jw4UK+/fZbrrrqKp+u8a677mLKlCnMnDmTjIwM4uPjyczMZMyYMQwaNMjjurBarRw/\nfpz33nuP119/nWnTprF8+XIGDRrkeR8CHx3gq6UbQWgs3aNAjarSEZKOBvvqFSAqxInPddH1IxER\nEbz99ttccMEFWCwWLrnkkoD1ZbFYqKioOCmLWUxMDFdffTVXX301paWlTJo0iZUrV/pcpVeSJFJT\nU1vMmDZx4kRWrVrFFVdc4VM/OTk5LF26FLvdzoYNG1i9ejVLlizh008/9Sw+EUWR/v37c9NNNzF9\n+nQSEhKatGG324MSjuUP0ZWDLLoKsF0QyFHVJl9+J1rsqwBkBmksuntBxy/069ePRYsWMWnSJGw2\nm0/WX2skJyczc+ZMLrvsMj744AOys7MBLdph0aJFbNiwgUsvvZQpU6awfPlyn0X3VEybNo3HHnvM\nb+2ZTCaGDRvGsGHD2p2T2Gq1Bk10fcHAr7l1gyVAB9GS20QA6wGbKGJXVepVlWhBwKGqQRPdUCGF\nuK6aLroBYODAgXzzzTdceOGF2Gw2brjhhoD089hjj/Hss88yaNAgnnnmGZKTk3nwwQexWq2MGjWK\ncePGUV9fT2Zm4L9GV111FXfccQfHjh0jOTlQtQfaRoNFXA9+XfZ7IhH4Zuk25FJzQsB8q1agBC0O\n1iZJWGUZFTjsznKWoihY0BaPCKrKYrRJNkuLLfqPM3U5rC66AaJXr16sWLGCcePG8fPPPzNnzhzM\nZrNf+5Akif/7v//jggsu4JprrkFVVR544AEuvfRSRFHk+eefZ82aNfzyyy9+7bc5oqOjSUtLY968\nedx8880B7681EhISMBoMLHE6EQCTKJKLlrc2EPiyuiwSrYKEP0TXjiawx9AE1qYoyKpKrCiSCKTL\nMha0GNiW4l8TJYlD7uOCRUdZnRcsBDXUC5FPc0pKSrj33ntZu3Ytf//735k6dWqHX1nmLe+88w6z\nZs1i7969xMfHh3o4KIpCcXExEydMwFFURPPxG77xH2Ai3lttKwWBHqpKl3ae12DBHkcTWLuiUN9I\nYC1uC7ahblhb+Rk4KgicGyRZ+A8wieCJ7l5g+MyZzD1Fsv1AcqZa+EEjLS2NDz/8kLfffpv777+f\nCRMmsGvXrlAPKyBcd911ZGRkMGfOnFAPBdAm3LKysqiuqPC6pM6p8CX/AoDBXcWiJRS02/2fgO+A\nlZLE12jl1Q+LIlGiSI4sc7aqMgEYqSj0UhQy0ZKZt/cLngDYgmwUnGlWn+5eCBKjR49m06ZNvPDC\nCwwbNoxbbrmF2bNn+93lEGo++eQThgwZgqIoPPfcc21aYRZojpeXE8h0QQq/+mfbi8ktuk5+rTpc\nDjgAhyRhl2UEwOJePZbpvvU34/8lsgDxaHkPnPhWv0ynZXRLN4hERkZy7733smXLFvbu3evJ13A6\neXgKCgrYsmULCxcuZNiwYRw5ciTUQ0IFjoliQCyq9li6drTogS3AOmCVJFEmyxwElgIb0KoDl6BZ\nqQWyzCi0ZbLnKAoFbjdEDIG7HY8AogWBYP7XTp9Pf9vQRTcEpKen8+GHH/LOO+/wl7/8hVGjRrFu\n3bpQD8tvdO3alT179hAbG0uvXr34xz/+gc3W2k10YFm1di3FkZH8LPr/436i6CpoluoutJCs/4oi\ny0WRRcAyYLcg4JIkkoBcWWYocD6aX/gCYISiYBYELGi10aII/iRTsiBQGqS+Ts/ZjdbRRTeEjBo1\nik2bNjFjxgyuvPJKLr74YrZt2xbqYfmFiIgIVqxYwcMPP8yjjz5KVlYWH3zwQUhKpQwaNIg5Dz5I\nuR99lTY0q1VBy5K1SpJYIggsQqvye8ztb+2iKPRWFEahCetoVWWALJMHdEG7nTfQVHxiRZFKv420\neWqB5cD/JImdgkAxWmpFF5CgKFhbKQ4azqhwUrWSYKOLboiJiIjg+uuvZ9euXYwePZqxY8dy7bXX\nUlVVFeqh+YXbb7+d/fv38/DDDzNr1iz69+/Pd999F/RxZGdn42inG8eFtmT2JzRhXSNJLHNbrcuB\nzWhf4k64rVZV5UI0i3WYe0Krq/v99lisMbJMoOuHRKP9YNTJMpWqyi+SxHpB4GtgG2ALYrrJYLoX\nXEBCUlIQezwZXXQ7CCaTiVmzZvHzzz+za9cuvvnmm1APya/cdtttlJSUMHToUM4//3y6d+/ONddc\nw9y5c/nhhx9wOn+t0KWqKsXFxSxbtsxvMca5ubnNiq6CNoG1F9gIrBUEVrqt1q/RhLVCkjCLIpmy\nTF+31ToeLca2EMhHs1ot+GdmOgpwBdjSFNHKv9uBHsBIWWacOwJiINA7oL2HDjUiIuThjHr0Qgcj\nLi6OlJQUpNPw9i4iIoLXXnuNZ599ls8++4zFixfz0ksvceTIEaqrq+nevTsAe/fuRZIk4uPjqays\npLCwkC+++IK0tDSv+05ISKBeVdmOVvXWKUk4VRWnoiDgTqItCKTIMmZ3hYJo3FEJzVh9h4BIUaRr\nANwlvmYvayuJQJ4gsAE4z11DTkSzzINJUCfSIiOxWIK59ONkdNHtgMiyfFJy7tOJmJgYrrvuuiaZ\n2I4dO8Znn32GIAiMGzeOrl27Alqe3Isuuoi+ffuybt26VnMV19bWsmbNGtatW8emTZs4cOAAx48f\np7KyUkuCIwjUCgLxqopZljGjRQIYoN1JtMslCUuAbsEjCF5y866qSrko8j9BYHgI/O3BnkiTIyKI\ni4sLcq9NOX2/2WGMLMunpaXbGsnJydxyyy0n7TcYDCxZsoSrrrqKcePG8fTTT7Njxw527drFvn37\nOH78OFVVVdTU1FBXV0dycjJnnXUWvXr1Yvz48Zx11ln06NGDzMxMcnJysFVW0qm2Fl+9eoHMl+BL\naff2IgD9FIVVaL7cUxeE8j/BtHQVUdQt3XClurqa/fv3s2/fPg4ePEh1dTVWqxWDwcDs2bN9miFV\nFMXj60xLS/MkCk9LS6NTp05nnCADfPTRR2RkZHD11VfTv39/unXrxujRo8nOzvaUHsrJyWl1MUZR\nURE333wzy+bNI6mRD9kb6lWVQNlLEsEt4xMJDEaLHe4EYVOK3RtcoFu6HZ3q6mpWrVpFUVERRUVF\n7Nixg927d+NwOMjOziYnJ4esrCwsFgvR0dG8/vrrjB8/3pNQ2xsef/xx1q5dS0lJCevWrfMkDS8p\nKaGiooLOnTt7RLixKF9yySWkpJy+X5kPP/yQiRMn8s033xDtRSLq2NhYunfvzlIfBRe0DGaBtHSD\nndzcAvQENgkC56lqUCtKBPNKnaqqW7odEVVVWb16NW+++SZffPEFQ4YMoXfv3owYMYKbb76Zs846\nq8Uqu9XV1Xz11Vc+iW7//v3p379/s+85nU5KS0spKSnxbIcPH+att96itLSUv/zlL17329EZOXIk\nycnJfPjhh9x0001etVFaWuqXdI/17omnQBCqFVpZwHFB4HtBYISiBCW0KVg+XRfaBGqdw6Fbuh2J\n4uJi3nnnHd566y1MJhM33ngjTz31VKu1uk5k4sSJ3HfffQETP4PBQGZm5kk5cvPy8vjPf/4TkD47\nErfccgtPPvmk16JbW1tLXXQ0G61WssBr365LVQlU/QEVrQo0QbZ2BaC3orBGENgC9AtSv95epYK2\nSKUGLT2mFS0EzgnIkoSsqrjcm4zmRnE5HCFfdn/Gi67D4eCLL77gzTffZP369Vx++eV89NFHDBo0\nyKsUjMOHD2fHjh2UlZXRqVPwgm8KCgp48skng9ZfqOjcuTO1td4vHXjxxRe58MIL+fLLL1m6dCkG\nQeCc6up2t+MCjF6PonVkNL/+HvdrtdFGo0flhH0nPm/YlEbHNt7nOV4QQBB+fa2qlKD9IIWiioST\nX4W0Dk1Y7WiTYC5B0IRUUZDRfiiMgoBJEIgSBKJVlQRFwSjLGNFcQEZ+XfW3IS6OkpIST6WVUHDG\niu6mTZt48803+fDDD+nXrx833HAD8+fP97l+ktFoZMyYMXz99ddNyqMHmvz8fH7++WdcLtdpHW72\nzDPPMGvWLK/PT05OZsaMGcyYMQO73U5MTAwu2vdFaBCuQLkXatFE/Yh7MlZwbwjCr88b7T/x9Ynv\nRQCiqmqvVRXR/eg5RlU974vufeXAXlEkM8BhZCrwnSiigMcqVdD+tkZRJEoQMAFJsoxRUTxCanA/\nSqDdEbTRejWpKgcOHODss88OxOW0idP329kCa9as4f7772fv3r1cf/31rF+/3hMT6i8mTpzIV199\nFVTRjY6OJjU1lb1799KjR4+g9RtMfvrpJ/bu3eu38kcmk4lOnTpRfuRIu2bsXWjiFCifpwEt5eO5\nIYibbaATsEJRsBO4CUPQBD5bUYgBj6BGuve3N3a6LUg2G/v27fN7u+3hjFkG/P3333PBBRdw7bXX\nMmPGDPbs2cODDz7od8EFmDBhAkuWLEEO4vp10FwMO3bsCGqfweTee+9l6tSpfq3BlpeXx/F2ntMg\nuoEiCm2iLpREAcmSxM4A9yOgZVPrhJbO8sTkP/7G4HKxO8RFBE570d20aROTJ0/msssuY9q0afz0\n00/MmDEjoLfgGRkZdOnShf/9738B66M5evbsSVFRUVD7DBZOp5M1a9b45Fpojn79+tFej64LEANY\nXSHK3UeoyZZljgWhikQwf17igLffeYe0Tp0YOWwYl1x0Eeedey7TLrmE3997L5s3bw74GE5b90J1\ndTVz5szhk08+Yfbs2fz73//GZApe9GGDi+Gcc84JWp/dunXj22+/DVp/weTWW2+lS5cuDB482K/t\nFhYWIsTGQk1Nm8/xiG6ArFEDmhDJeF+Rwh90AmT3pJr3WS86FknAOJcL67Fj1B47xgE0EdyDlvDo\n3//6Fzt37/Z5bqc1TktL9/PPP6ewsJC6ujq2b9/OXXfdFVTBhV9FNxjU1NRw//33M2fOHEaNGhWU\nPoPJH/7wB+bNm8fHH3/s96KeeXl5yO2866kHpABagCKa2NYHrIe2jyNTENgXwD58rTHnDSJazo1U\ntOxwKWhRGoWAWFHBfX/+c8D7P204ePAgl1xyCX/60594//33eeONN0gKUe7Mc845h71791JSUhLQ\nfj799FPy8vLYt28fmzZt8tskU0fhz3/+M//4xz9YsWIFvXv7P+Fgfn4+tVZru5bd1hN4C1QShJCL\nLkCsqlJ/Bi0772Gz8c/XXuOHH34IWB+nhXtBlmVeeuklHnnkEe68804+/vhjjMZARVG2jYiICMaN\nG8fixYu5/vrrWzxOVVXq6uo4duwY5eXlzT42fj5p0iT+9re/ec7fsmULsbGxPPvss0GNCw4G9913\nHy+//DIrVqxocYWerzTksqgDYtt4jhMtBCuQiMBhNP+ugjt2t9FztdGjQtNY3MbPARAEIlSVzmiW\nXXtC3aIJ8KSeIIR8sUJjjEB3m43fXHEF23bsCEhh1bAX3Y0bN3LzzTdjNpv573//S35+fsjGoigK\nP/30Ez/++CNbt25l165dbNu2jfXr11NdXd3sVlVVRUREBMnJySQlJTV5TE5OpqCggOTkZGpqarjr\nrrt46aWXmvT50EMPoaoqI0eOZOnSpWRkZITo6v3LHXfcwfvvv8+yZcsYOHBgwPoRBIGuXbtSvn17\nm0W3HhCDEL96EC1nb0N4mojmS26IpW3YJwERjeJvT4y5FVUVmyhyAChSFEyiSLyikI9WVbg1ooD6\nEIauhYIuwJYjR5hz3308+dRTfm8/bEW3traWBx54gPfff5/HH3+cGTNm+N3f1xin08muXbvYuXMn\n3377LcePH2ffvn2UlZVRXV1NbW0tNTU1GI1G0tLSyMnJYejQoWRlZREfH09cXFyL26n8zTabjREj\nRvDoo48yfPjwJu8JgsAjjzyCxWJhxIgRLF26tNWcs8FGVVX27t1LUVERkyZNOuX/SFEUxo0bx9at\nW1m7di29egU+2WCf3r1ZvX07OW08vl4UiQywEBlEkbNkmTR/9eNuxwUcVxSKJYlVsozZXccth+bF\nQCSwPtdQ+HRPhQDkW628/sorDBk6lMsuu8yv7Yel6C5cuJA77riDkSNHsm3bNq9vq51OJzt37mTn\nzp3s2bOHffv2cejQIcrKyjz5Wa1WK1arFbvdTnR0NMnJyRw8eJDRo0dz7rnnerKMZWVlkZmZidl8\nKtuh/dx+++306NGj1XCp3//+91gsFs477zwWL14cEP9nWykuLmbFihUsX76cFStW4HA4sNlsp3QT\nVGNgSkMAACAASURBVFZWMmTIEEwmE1u2bCE1NdUzGVpWVobJZCIqKuqkR5fL5fnRq62tpba2lm3b\ntlFWVkZcXBzx8fFYLBYSExNJSEggMTGRxMRE4uPjMZlM9O3Xj2/mzYM2Zh+rF4SArUZrQEBzH/ib\nCKAz0FmWcQKHFIUDgsAuVcXsTsx+FppboWEcZyJGoJfVyg3XXUdeXp5fv0+C2pEcKqegpKSEu+++\nmw0bNjB37lzGjBnDkSNHOHToECUlJRw5coSysjKPD/T48eMe8bTb7TgcDux2O3a73SOkZrOZpKQk\nUlNT6dKlCxkZGaSnp5OSkkJqaiopKSmkpKTQqVMnT2xvdnY2zz33HFOnTg34NSuKgtFo5OjRoyQk\nJJzy+H/961/cfffdLFiwgKFDhwZ8fABHjx5l5cqVLF++nOXLl3P8+HFGjRrFmDFjGD16NPn5+cye\nPRuAxx57rNk2li1bxvjx41FVlVmzZrFnzx62bt3KoUOHyMvLIy0tzSPeDf9Dm82GzWYjMjKSmJgY\nzxYbG8vmH3+kqrQUoyhqCU9UFRdaCJTsXmraYEM2LIf13L67b+E9jwCKgqiqSGhVc2PREsK0P8Fk\n21gjSWTLMlkBav9E7EAZUCqKHFUUzKJIqqKQCaxEq2QcCJYIAmcHMDexrxwCjqamsnnbNr9NyoeN\n6D744IM8++yzCG7Hu9PpxOl0EhkZSXR0NLGxsR6rJj4+noSEBJKSkkhISMBisXi2uLg4OnXqREpK\nCsnJyV4tkujTpw+33XYbt912WwCu9GS6du3KN99846khdioWLlzIjBkz+OSTTxg9enRAxzZv3jyu\nvfZaRo8ezZgxYxgzZgy9e/c+KYn7xo0bmTZtGnv27GniYmhIOHTFFVcAMGXKFHr37u3ZevToQWRk\n+5MxPvTQQ3z4yCOcdYpVgY0np2S02+/WHmVB4ICqEgE40KIMokSRKFkmGS2e1R/Bif8VRTLct/3B\nxgGUAAdFEZuq4lRVLgpQXx1ddAF2RUbSaeBAVq5Z45dFVWHhXli2bBmPP/44CQkJLF68mLi4OI+A\nhiK5S3x8PEePHg1afzk5Ofzyyy9tFt1Jkybx73//m8svv5x//vOfTJkyJWBjM5vNDB06lC+//LLV\n4/r160dERAQ//PAD2dnZfPXVV3z55Zd888039OrVi1WrVjFy5Ei/jatv3768bTbDKTKINViybf4U\nqSpHRJFeikISUKuqVMkyVaLIQbSJqki3EEfLMp3Q4kHb644Q3NZ4KDACOWg5EbaJIsXusQQivjQc\n3Bc96uvZtHkzs+68k5defdXn9jq86K5du5arrrqKJ554gr/+9a/07ds31EMiISGBY8eOBa2/KVOm\n8Mc//pGhQ4e2OQHzqFGjWLhwIZMnT/a0EQh69uzZpnwPgiBwxRVXMHnyZOx2O+effz5Tpkxh7ty5\nAQl169OnD5WuwC2mVdFEKM69NWTjUoAaVaVSlqmSpP/f3rmHRVWu/f+z5gDMcFLUQEzFPB9AIcmz\nhW01D71XqWlbs3Npdtj+1Mze9lY7uNuaae4OmuZl2UHDaqftXjXySOIBRARRUUkDJRUBhWEYhpm1\nfn+sYUJEYWbWDKPN57rWNcPMWs96gJnvetbz3Pf35ldRJEuS8LPZDwaJImHIYlPTA1ZUq7FIElU2\npy1JkvCE+aAVuMQfI/2aGxoNfmo1YmUlO1UqtNdZBBWuc7MswLVZe7V+rpIkMrAlm1TbS9oetaJI\nLw9XsagLAeheUcG6tWuJ692bJ596yrX2vHl6IS0tjZEjR/LFF1/QrVs3unbtSpkD6Zru4oknnqCy\nspKvvvrKI+eTJIlp06ZRVFREYmKiQ8c+99xzREREMG/ePLf1LSQkhPz8fJo0aXLDfYuLi8nIyGDA\ngAFuj6MWRZFAvZ7BlZWKL3rtUKvpbrXSUGt7K1CKPBd8Ra2mwGpFA4Sq1egBf6uVALhqszttuZlL\nwNGgIAYPHCjPiQcHExQSQnBICIGBgQQGBvLNhg0cTU4mUpKu8fSt+byu1+p7v/riVdOK0m4vqVJx\nQRQJqfEaNfarRm3bVDWea2w/V8c3m5FD/aqQp4pElQrRJu6BQO8GmFOVAek6HVu3bXMpvd9rR7pb\nt27lscceY9WqVQwbNsy+cCKKoktFH5UgLCzMo8YygiAwZcoUHn/8cYeOKy8vJzExkYyMDPd0DLlv\nXbp04dixY/V+EMPCwhgyZIjb+lITlUpFp/btKTt61OXKv3XhyEhFDTS1bVitlNsWqdp72IXuenTu\n0IEfNm++7vvt2rVjWkYGrZ0we3eFNqLIJWRvX7tJu00opRqjbqsgICILqgnbHL0kUS6KWICmajVa\n5AuZvyShEUU0oogW+X+TgexhHFRPf4KBLhUV3D9yJBlZWU7HxHud6JaWljJr1iy2bt3Kl19+yb33\n3gvI3qcajYaCgoJGSQDYtm0b//3vfzl69CgHDx70+DSHM2XZ161bx8CBA68p7aM01e5mnjT3aQhx\nvXuT6ibRdYXq0Z23UN/Nbp8+fbhoMsllhDzTJTvNbZud6r424Ab9EKBSq+lZz8XtvFpNjtVKQ1Jw\nwoHysjLuGzqU1PR0p4xxvOl/z7Zt24iJiUGSJLKysuyCW01oaCinTp1qlL699dZb/PDDD/To0YP3\n33+fjRs3evT8joquJEl89NFHTJ061Y29kmOddTqdV/r4xvfpQ6Ub3KJcFZ7GEK8bUZ/oNm/enGZh\nYTT+xJ5jVCJP3dTH7VYrxQ4kVrWzWik/c4ZH/vpXp1KYvWKkazAYmD17Nj/88IO9hHldhIWFkZub\n2yhOWtUeue+++67Hzw3UO61y/vx50tLSSE1NtT+WlpYyaNAgxftSVVXF9u3bSUxM5Pvvv6dLly5e\nWYU4JiaGCj8/qKhQvG1XFkK8aaQrIH+26mPS5Ml8+/77hJhM7u+UQqiQ527rq0DRHDlV+keunlNW\nCQKCIBAoSXSSJPsdkwB0NZnYk5TEgrfe4u//+IfD/WpUdu3aJX85KirIysq6ruCCXJQwPz/fg737\ng6ioKM6dO9co54arR7pFRUX89NNPLFiwgAcffJDWrVvTrVs3/v3vf2OxWHj66adJT0+nU6dOLF68\nWJHzWywWkpKSeOaZZ2jZsiXz58+nR48eZGRksGfPHoYPH67IeZSkR48eXDIa3ZJm6kqbos0bwVto\nyGhtzquvUqhW43xJUM/jD5gbMIJVAUOA+4C/AHcDA4A+kkSsKBIiCBwAdqhUHAVygFNaLRqrlbnz\n5pGamupQvxptpFteXs7//u//8u2337JixQpGj64//Do8PJyzZ896oHfXEhkZSVFRUaOcG+TRSHZ2\nNu3bt6ewsJC4uDji4+N5+OGHWbx4MXfcccc1vgZTp07lvffeY/78+U6d02KxsGvXLhITE/nuu++4\n4447GD9+PAcPHmzUaqoNpUmTJoSGhGAsKqrX2MURNJJEpQvHe9tItyGi27RpU2bNns2af/2L7m64\nc3AH/sCVBu4r8EfkQ22aiSKdkWOw9d268eDYsVclXXXt2tWhfjWK6O7Zs4fHH3+cvn37kpmZSVhY\nWIOOi4yMbLQ53cjIyEYNV7vzzjtZuXIl0dHRdO7cuUERHBMmTGDGjBmcPXu2wYuPVquV5ORkvv76\na7777jvatGnD+PHjOXDggFvqybmb7t26UZqcrKjo+okiRhcrR9xsI12A/zdjBkvffZfSigqvziCr\nJgC4qFBErAY5emHAgAG88cYbLrXl0QtuRUUFs2bN4qGHHuKdd97h888/b7DggjzSLS52tIygMrRs\n2RKj0ahIW5cuXeLtt9+mX79+xMfHN2hOTa/XM378eLp27drgkLlmzZoxePDgeke6VquV3bt388IL\nL9CqVStmzJhB27ZtSUlJITU1lZdffvmmFFyAu/r1w6Cw+1wQYHAhbNFd2V3O0NA5XYCgoCD+MX8+\nJ/R6r6jhVh96wKxgGkIAsH/vXpf9fz32v9+3bx+xsbGcO3eOzMxMHnjgAYfbqPaVbQxatmxJeXm5\nU8eaTCZWr17NsGHD7AtyGzZsYOjQoRQVFfHaa68p3Ns/ePbZZ9myZcs1r4uiyC+//MJLL71E69at\neemll4iMjCQ5OZn09HTmzJnT4LRjbyY2Lg5zUH0RmI7RBLhitXIR+fa1AhxO2b0ZR7oAf/vb30i4\n/36y9Hq3uKApidL9iwDyT592uQ6h26cXqqqqeP3111m9ejXvv/++S96ULVq0UGy06SgRERFUVFRg\nNpvrdZMXRZEtW7awZs0aDh48SEFBAa1atWLUqFFMnz6dQYMGERwsW2YPGjSIcePG8eqrrzY4xdcR\nRo8ezWOPPUZqaip33nkn+/btIzExkQ0bNtCsWTPGjx/Pjh076Ny5s+Ln9gZiYmIoVTjp0h85HClT\nEGTXMmTRrWksbncrq+FUJkgSgihikiROIJuUq2sew7WZVTWzrTS1XtPWOtZZHPnrqFQqPv/ySx4a\nM4a0pCSiKyq8ZtRemwKghUoFCiWhCECoSuXyupJbRffUqVNMmjSJ5s2bk5GRQXh4uEvtNW/enIpG\nmsTXarUEBgZy9OhRevXqdc376enprFy5kt27d5Ofn09AQAD33XcfCxYsYMiQIdf93YcOHcpdd93F\nI488wqZNmxTvd0BAAOPGjePhhx/GbDYTEhLChAkT+Pnnnx1eALgZ6dixI2WVlVhQ7sNeBoSoVAyu\ncVteXT6n2pHMYvNQsP9c471SIEQQUAsCYo1sqpqOZ2IN+0n7c1sabs2yPNWCWTM11p5SaxN8bM+r\nLwrYnlcvogU5OGWnVqv5+ptveOD++zmyezc9vFR4DWo1dyic9SdIEpWVriyjukl0JUnis88+4+WX\nX2bu3Lm88MILilR1aEzRBTlkLTs7m169epGXl8fHH3/Mli1bOH36NGazmbvvvpsXX3yRoUOH0r59\n+wb/zh988AGxsbFkZ2fTvXt3xfttNBrp2LEj7777rlva92Y0Gg3t2rShLDeX+t2IG8YVILjW/1ZA\n/jJpkEfCNyIX6CRJ6BQagdeui3aNWCN/J2vWUqvetwr41WBAkiSHvqNarZbvNm5k5LBhHDtwgG4m\nk1dNmVSjtMCpqqq4fPmyS20oLrolJSVMmTKFY8eOsWPHDkXLrbRo0QKTh4OzRVHk9OnTHDp0CIvF\nwmuvvcbMmTO5fPkycXFxPPTQQwwdOpRevXo5nKZbTefOnXnqqaeYOHEihw8fVrT/e/fuZe/eveTk\n5LilqsXNQK/YWI4qKLrlQHMXRlBKZ6TdKNypIX05g1ztw9F0cX9/f/67ZQsJgweTk5VF58pKrxJe\npatvXEQ2eY+MjHSpHUVFd9euXUyePJkHH3yQtWvX1lv7y1GaNm2K2WzGaDSi17vu2W82mzl69CiZ\nmZkcO3aMU6dOUVBQQElJCaWlpZSXl2MwGPDz8+O2226jVatWDBw4kOHDh9O/f39Ff7833niDtm3b\nkpiYyPjx4xVps7CwkClTprBgwYI/reACxPftS/qmTQ0ux1MfFrUavYu3rd4iTgIQ5udHRkaGUx4d\nOp2OpO3bGdivHydPnKBjVZVX/W5KehJfCgnhnXfeYezYsS61o4joVlVVMW/ePD799FNWr17NiBEj\nlGj2GlQqFYGBgeTm5jpUs8hkMpGTk0N2djbZ2dls3ryZ06dPU1paSlBQEBEREbRu3ZqoqChiY2Pt\nEQbVj0EKr37XRdOmTZk9ezYLFy5URHTPnDnDsGHDmDBhAo8++qgCPbx56dmzJxX+/vWKrgjsFwSs\nKhWCKMoLX1xtPagCyq1Wl8r0eJv3gp/BwMG0NLv3sqMEBwezY/du+sbH82teHu3d6GPsCO4wglfC\n+9ll0T158iQTJ04kPDycjIwMbrutoS6jDUMURX799VdSU1PJzMzEYrHw6KOPEhwcTFVVlX2rLt9j\nNpuxWCz2zWw2U1VVRYcOHejevTvdu3cnODiY+fPnM3XqVLf7ujpCnz59WLFihUttGAwGVq9ezcKF\nC3n11Vd58cUXFerdzUd1FMmCN97gbFkZ5wCNINg3tSShFUUCkEuNX0b2XW1rtf6xoIVsI1jzUSWK\nGIGGR5hfizeJbrAosjc52aU2wsLCSE5JYUCfPpwoKPCKEa9gW8j0NpwWXUmSWLNmDa+88grz589n\n2rRpLi+WrVy5km+++YaLFy9y+fJlysrKKCsrQ6PR0LJlS9q2bcukSZMICwtDr9cTEBBg3wIDAwkK\nCrIbL+v1epKSkli+fDnp6elXGWwfPXqUiIgIrxJcgOjoaIqKipzyDL548SLLli3j448/JiEhgY0b\nNxIfH++mnno3JSUlrP7kE95bsgSrwUC4wcAIZAGtlCQqJQkTcthXpSBgUqkostUCi5Gka+d+a9kJ\nmtVqSqxWnDUY9baRbghwODPT5XbCw8M5kJ7OXxISOHbiBF1MpsaNamjEkkc3winRLS4u5tlnn+XE\niRPs3LlTsRXx+fPnM2TIEB544AHatGlD27ZtadOmDaGhoU71cenSpXz99dfXVDRwdKXWU4SHh6PV\naklPT6d3794NOubcuXMsXryYzz77jAkTJrB///5bIqnBGSwWC1OfeYb169cTLghEVVTQhD8Erjq2\n9arJIklyOI4zyGrFVRcOb/r0BQJXysooLi52KEO0LsLCwvhl717uHzmSrAMH6F5R0WgGL970N66J\nwxeiHTt20LNnT1q3bs2BAwcUDUESRZEXX3yRadOmMXr0aKKjo50SXIDExEQGDx5cp7WhVqvlp59+\nwmDwLs8kQRDo1q1bnRlkNZEkiYMHDzJlyhSio6MRBIEjR46wfPnyP63ggvz5+XTtWvqYTHSrqKAp\n7vniBQFVTkaqgPeNdAWgeUCAYhVG9Ho9W5KSuGfMGFL1+pvOh9fdNFh0zWYzc+bM4ZFHHmHVqlUs\nXbpU8egErVbLgQMHFGmrrKzsuiYvCxcupKCgwK3pt84SHx9PSkpKne8VFhaydOlSevbsybhx44iM\njCQnJ4clS5a4HMZyK+Dn50ebVq2ocvN5ggBzA/0KalN9lDeJLkCAyUR6erpi7Wm1WtZ+8QWL3n+f\ndL2e84q1fPPTINHNycmhf//+ZGdnk5GRcUPPW1dYuXIlc+bMISsry+W2jEbjdcOkWrduzYwZMxSP\niVWC2NhYTp8+bf+5uLiYb775hgcffJCOHTty6NAhli1bRm5uLvPmzXNLJd2bmejoaLePrAKRjVSc\nkV1vFd1As9nlxbS6ePLJJ/nPpk3k29LefdQjupIksWrVKgYOHMhTTz3Fpk2b3PolHzFiBNHR0Xz7\n7bcut1VQUHDD+amuXbty/Phxl8+jNDExMVy8eJFZs2YRFxdHVFQUn3zyCaNGjSIvL4+1a9eSkJDQ\n6MU5vZXeffti1Lh3FlGDPDdc4sSx1R4N3kYIKDrSrcndd99NqdmMMlHSjYdSDh43/HSOHTuW//zn\nPzzzzDOUlJQwd+5cjEYj9957LyNHjlSoC1cTFhZGaQOqjhYUFDB37lzUajV+fn7XbJs2bWL69OnX\nPf7kyZOUlZU5VfDRnXTr1o2ysjKCg4P54IMPiI+PR6vVNna3bhpiYmIw6/Xg5sq1gYLABUkiAK5J\nr7WHmtXxXhXKfXmVJBgoOH9escSjmmg0Gu7s1Yvi/fuJULRlz2EBLldWYlYgweaGohsSEsIzzzyD\nXq/HYDCg0+k4fvw4GzZscJvoBgUFNUh0Dx8+zP79+3nhhReuitGtzlibPn06Xbp0qfPYH3/8kSee\neIIffvjBqwQX5EWIqKgoxo0b96fzSVCC6OhoLnsgOF+UJE4Dv1HDYKbmVsNgpuZzCTnO19tQAWF6\nPVlZWfTp00fx9u8fM4ZPMjKIcNEsxtNcBE4j39VYKysV0Ysbiu6nn356zWsff/wxBw8edPnE1yMk\nJKRBnrkGg4HOnTszZcoUh8/xj3/8gzlz5jBkyBBnuuh2wsPDOX/+vE90nSAqKgqTxUIV8hSAu1AB\nnYCO19uhpplNjedlQIogYLCV/DEij6LMtq3KtlkASaVCFASsyIbcfRR2zKpNYFUVhw4dcovo3nff\nfSx6803wctE1A2eB3wVBrgcnSbRSqegsipwJDSUiwvWxusOTXwaDwW23u8XFxSQlJTVoFF1aWmr3\npHWU9957j/Hjx5OTk2NPqIiLi3PKWF1pLBYLhw8fJi4urrG7clNy6dIlgoKCKC8upkn9uzuNBrA4\nWbLHKkkkI8cNV8toU7UaLfKFQieKaCQJjSiiRRaCXCfjys9otZj8/JAEAQQBSRCQsE1x1HguAldM\nJrIUSJKoi+joaPwDAzltNBIlil4xr21EFtgioFKtplIUsUgSQYJAS6C7JBEKCLY7E6skKZJQ5ZDo\npqens2jRItatW+fyiWuTk5PDoEGDGDRoEEuXLq13/7y8PKcMOgAGDx7MunXrOHHiBAaDgU8++YSz\nZ896hegeOXKEVq1a0bSpUp5Yfx4+XrGC2bNmEWE2u72GlxaocsIgOxioOaTYqVbTyWol8gbtGICT\nznQSyNdqmf/WW7Ro0QK1Wo1Go0Gj0Vz1vObPHTted+zuEoIgsGffPkaPGEHa2bOElZcTKUno3HK2\nq7EAl5CnCsoEAZMgYBZFJGRT8maSRKjVShByZIrqOhdSETwvuqtXr6ZDhw7079+/wcfk5uYyYsQI\nRFFEFEUkSZK9PWs9Ly0tZdq0afzrX/9q0Mq8wWCgoKDAqZRZgISEBBISEjh58iRz5sxh7ty5lJeX\nN7ob1969ex36+/qQ2bp1K7NnzCCuogL32xPZRFeBdqrneevdx0nv3SqrlUcffdTlTDMliIqK4vCR\nIyQnJ7Nm9Wr+8913hGg06EwmdGYzwch1yGrPj2u5NuKjelGyHHm+1YA8chVtUzIiUGm1UgqcQPY4\nbqpWE261EipJBNteExyYXxclCY0CkTGC5MB/s7KyksmTJ1NYWMjy5cvR6XRotVqaNWt23StAXl4e\n7dq144svvkCv16NWq+1b9RVWrVYTGhrqkHNYaWkp999/P+3atatz7rmhFBYWsmjRIlJSUsjIyKBL\nly5MnDiRmTNnOt2mK0yePJm7776bp59+ulHOfzOSn59Pzx496FxaSnMPnTMHKFGr6eviPOtutZp2\nVis3umerBH4Ggmos4tQUodqCJIB92uOKKFJuNKLTeWJM6RiVlZUkJyeTmZlJ+oEDZGRkUHjpEqLV\nitU2SLNYLFisVprpdGgEgXKrFaPZTKXFQqBOhySKSCYTTVUqAqxW/PjDTF6DLOKBOOc1XJt9QUGk\npKe7fDfgkOiCXDn25ZdfZuPGjXaHL39/fzZs2HBdg5WIiAjWrVtHQkKCS52tjdFoJCoqitTUVNq2\nbetyeyaTifHjx3P77bfz0UcfKdBDx+nQoQMbN270LaI5wP+MHMnJn36ig5sXmmpyBshXqRjkYiRC\nslpNW6uVNjfYRwIKwV6BV6rj8XqvHRMEyo1GxbNHPUlJSQnHjx+ntLSUqKgoWrRoQZMmTVCpVEx9\n9lmSV63iDg/0IyUoiAMZGS6n2jt8X65Wq1myZAm5ubnk5eXx+++/s2TJEkaNGsXy5cvrvA1q3749\nO3fudKmjdaHX60lISGDXrl2KtHf8+HH279/Pm2++qUh7jvLrr79y5cqVP0XtMiUpLioixIOCC/Kt\naZUC5XYaYrQtALcBkbatlW273ba1BtrYtra2LQpoB/hrtY1a4koJmjZtSr9+/Rg+fDidO3cmLCzM\nPqXYNCzMY/aNWpXK6YrgNVEkrWnMmDHs2bOHFStW8Pe///2a9/v378/XX3/tckG3umjfvj15eXku\ntyOKIs8//zwLFiygWbNmCvTMcZYvX84TTzzhyzZzkACdzuMWfgEoI7rYik26C61afdOL7o1o0rQp\nopszEKvRCoLL9dFAIdEFuerq22+/zf79+695b+HChVRWVrrFUDs8PJwLFy643M7atWuxWCw8+eST\nCvTKcYxGI2vWrOG5555rlPPfzAQEBHhcdA2Anxfag9ZGc4uLrk6nAw9lbGoliStXrrjcjqJDqpiY\nGI4cOXLtSVQqdu3axfr160lMTFTylERERHD+vGseRiUlJbzyyit88MEHjTbKXLt2Lf369aNdu3aN\ncv6bGX9/f4+Lbr4gcJsCoivVKJPuDjSCcMuK7uXLl1n4z38S6qHfT2W1cvHiRdfbUaAvdvz9/bFY\nLHXO67Zp08aeCeZs+EttKisryczMpKTEGeuRP8jPz0ev1zNq1CjGjh3LsmXLOHToEFYPzRP+/vvv\nzJ07l9dff90j57vVKCkp8ahRdgVQLEnoRJFS/ljgchZ3jpcFQfDY59iTSJLEo5MmEXj5ssf8HALL\ny9mogBmXoqIbGBhIaGgoUVFRTJkyhe+///6qlN6ZM2dy7tw5l0Wymvfee4+kpCQWL17sUjsxMTGc\nPn2atLQ0xowZw9GjR5k4cSLNmzdn9OjRLFq0iH379lFVpbxTqyRJTJ06lSlTpviy0Jzk/O+/48nC\nSwLQRBDIU6vZJwhsBbYCp528S3Kn6FpEsdFjz93Bhx9+yP5du+jgwbTicGDb9u0u3zkoOkDQ6/Wc\nOnWKY8eOsWXLFj788EMmT55M7969GTFiBCNGjEClUik60h0+fDi9evVSpL02bdowadIkJk2aBMCF\nCxf45Zdf2L17N1OnTiU3N5c+ffowePBgBg8eTJ8+fVyOf/ziiy/Izc1VfNrlz0JhYSFn8vJoyx/O\nXtR61KCssAUAA2uU+RGBI8BFQcDRySF3V5GoEkXFXcMam4yMDF6dPZveFRWKxN82FH/kcvVbt251\nKXvV4ThdRzEYDOzcuZPNmzezefNmCgoKGDlyJPfccw8DBgwgJibGaS+HhQsXUlRUxKJFixTudd2U\nlJSwZ88edu/eze7du8nKyiI2NtYuwgMGDHDIDyI7O5t77rmHn3/+mZ49e7qx57cuO3bs4N575LI2\nrAAABoRJREFU7wVAJQiyo5dtUwkCoijSRq2mi5tHRL8Bx8F+q1s7eUFCnoaw1tgQBMokiWjkEDB3\n8LOfH79fuHBNncCblbKyMnp07Uqzc+fc9je7EWcAXVwcP23b5vTf1O2iWxNJkjh8+DCHDx8mJSWF\nlJQUzpw5Q+/evRkwYAD9+/enb9++DU5ZXLp0KWfOnGHZsmVu7nndGAwG9u3bZxfhtLQ0unbtahfh\ngQMHXjf8rKysjPj4eObMmcPjjz/u2Y7/iUhKSuLphx4iRoFV5xtRBOwDWtmyxmomLIAsvGpALUmo\nRRE18tyeGjn21s9N/dqsUlFhMt0SnsySJPHQ2LFk/N//0bWR3MpE4KSfH8awMDb9+KNTU4IeFd26\nuHz5Mvv27bOL8IEDB7j99tvtIty/f386depUZ/Xejz76iMzMTFasWNEIPb8Wk8lEamqqXYT37t1L\n27Zt7SI8ePBgWrZsiSRJPPzww4SEhLBq1arG7vYtzYoVK3hn5ky6GY1uPY8J2AaMcutZHEMEtqpU\nWG6RhbRPVq3i1enTiTcaPTqtUBcFwEm9nh+3bKmz+O2NaHTRrY3FYiErK4uUlBT27NlDSkoKBoPB\nLsD9+/cnPj4enU7H6tWr+eWXX1izZk1jd7tOLBYLhw4dsotwcnIyzZo1o0OHDly4cIGUlJSbOj3z\nZmDWzJn8d8mS6/veKsQVIBnZQcxbUlsswE4/Pyq83MO2IWRnZ9PvrruIMxrxlmprhcDxoCB27N5N\nbGxsg4/zOtGti3PnzrF37167CB85coQePXqg1+tp2bIlX331VWN3sUGIokh2djZ79+5l5MiR161W\n7EM51qxZw4IXX6SrAumb1yML2Ze1vUpFJy+qCmEBdmi1mBQoMdPYPPLXv5K+fj2uuR4oz+/A6dBQ\n9qWmNtgI56YQ3dpUVFSQlpbGnj176NSpE2PGjGnsLvnwUnJycrgrNpZBFRVuiRIwA0nAQCDUDe27\nghVIUqup8kD5Incz9dln2f7ZZ3T2wgvIb4KA2KkTh7KyGjR3flOKrg8fjtCtUydCTp7EHXWs84FT\ngkCCF36NRGCLSnVLJEcUFhbSoV07YsvLvWZ6oRoJyNTreWz6dN5csKDe/b1l+smHD7fx/N/+Rr5e\n7xZjmYtAMy81KBLAXizgZqdFixbMe+MNTrnp/+gKAtDZaGTZ0qWkpaXVv79vpOvjVsdqtRIfF4cl\nK4u2Cn3cLUA28pxeD2SLRW8kSavlQmEhoaHeNvnhOFVVVQzq35/CI0foajJ5RZ21mpwFjO3bk338\n+A0rTHjnJdqHDwVRq9Ws37CBMzodrhvzySSrVBhUKnrjvsQGJQjy91fEpMUb0Gq1bN+1i8jYWLJ1\nOq8b8bYCys+f5716ajz6RNfHn4JOnTqx9ssvydTpUCKOwSxJ9BRFmuPeNF5X0anVt4zogmw1kLR9\nO5ExMZzwc1dKiXMIwB3l5SxeuPCG+/lE18efhgceeIB/vvMOh/V6XF0Dt0qS27LI6kNCjkwwIydl\nlANlwGWgGDl+9AJyAL+pslIRv2lvIiAggB+3bMHasiW/edl8eghQWY8hjicd8Xz4aHSef/55jmZl\nsfHzz+lpNDo16rAgRwZoajy3OvGIWg0aDZJKBWq1XMW2erPtYxVFLJKExWrFIoqYLRasVisajQY/\nrRY/rZYAf3/8/PwICAggwN8fnU5HoE6HTqfjjqCgW7L8U5MmTdi+axd39uqFvwftHetDAJrVEy3i\nW0jz8afDYrEwbMgQclJTCVSpkGxluyXbJoJ9s0oSoiTZxc8qilRZLPbYV61N+Pz9/GTh8/fH398f\nXUAAAQEB6PR6dDod+sBA9Ho9gUFB6AMDCQwKQqfTyfs4+Ojv7+8r6WQjLS2NvyQkEG4wEIjsBOaP\n7GWhpXFu5XMEgRM3SJLxia6PPyWlpaWsX78elUqFv22k6G8TzOrnN3pNq9USEBBQpyeID89y7Ngx\nli1dytm8PC6cP8/FwkJKLl/GaDIBcnFOP40GrUqFRqVCC/hZrQgmExqLBX9Ajzw1oMStf45Gw4kb\neG/7RNeHDx+3LGazGYPBQHl5uf2xuLiYCxcucOHCBc7m5ZH/228cO3aM3DNnaBIQgF4QwHYxrcuf\nufbjVa9JEqWVlRhsgl8XPtH14cOHDw/imxjy4cOHDw/iE10fPnz48CA+0fXhw4cPD+ITXR8+fPjw\nID7R9eHDhw8P4hNdHz58+PAg/x+zAzIk/7uIegAAAABJRU5ErkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x4b258d0>"
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We're gonna use this to build block weights:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"block = ps.weights.regime_weights(ew)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 13
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we get the standard contiguity using the shapefile:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"cont = ps.queen_from_shapefile(ps.examples.get_path('columbus.shp'))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 14
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the combination will give the desired structure, one that uses contiguous values, but only \"valid\" ones (here we're assuming \"valid\" means the surrounding observations are also in the same `EW` group). This is basically the intersection of `block` and `cont`, where $j$ is neighbor of $i$ if $j$ is contiguous *and* in the same `EW` group as $i$. This is encoded in the `w_intersection` method of `PySAL`:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from IPython.display import HTML\n",
"HTML('<iframe src=http://pysal.readthedocs.org/en/v1.7/library/weights/Wsets.html#pysal.weights.Wsets.w_intersection width=700 height=350></iframe>')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"<iframe src=http://pysal.readthedocs.org/en/v1.7/library/weights/Wsets.html#pysal.weights.Wsets.w_intersection width=700 height=350></iframe>"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 15,
"text": [
"<IPython.core.display.HTML at 0x4b25b90>"
]
}
],
"prompt_number": 15
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"w = ps.weights.Wsets.w_intersection(cont, block)\n",
"w"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 16,
"text": [
"<pysal.weights.weights.W at 0x4b25ad0>"
]
}
],
"prompt_number": 16
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment