{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Welcome and please:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic tricks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Given array $a_0$ that contains integers and indicator array $X$. Return:\n", "\n", "a) boolean array $b_0 = \\left[\\left\\{\\begin{array}{ll} \\text{True,} & x \\in A\\\\ \\text{False,} & x \\notin A \\end{array}\\right. \\quad\\text{for $x$ in $X$} \\right]$\n", "\n", "b) list of boolean arrays $b_0, b_1, ..., b_n$ if there are multiple arrays $a_0, a_1, ..., a_n$ given.\n", "\n", "Source: [my question](https://stackoverflow.com/questions/59656759/what-is-a-best-way-to-intersect-multiple-arrays-with-numpy-array)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[False True True]\n", "[0 0 2]\n", "[0 1 2 3 4 5]\n", "[2 5 4 0]\n", "[array([False, True, True]), array([ True, True, True, True, True, True]), array([ True, True, True, False])]\n" ] } ], "source": [ "X = np.array([2,5,0,4,3,1])\n", "A = [np.array([-2,0,2]), np.array([0,1,2,3,4,5]), np.array([2,5,4,6])]\n", "\n", "print(np.isin(A[0], X))\n", "\n", "Xs, B = np.sort(X), []\n", "for a in A:\n", " idx = np.searchsorted(Xs, a) #locations where to append items of a in X\n", " idx[idx==len(Xs)] = 0\n", " print(idx)\n", " B.append(Xs[idx]==a)\n", "print(B)\n", "\n", "#question: what about np.searchsorted(Xs, a, sorter=np.argsort(X))?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Given array of people descriptions. Replace it by labels 0, 1, ... , each of them corresponds to alphabetical order." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1 1 1 0 1 2 2 1 3 2 2 2 1 1 1 2 2 3 3 1]\n", "[1 1 1 0 1 2 2 1 3 2 2 2 1 1 1 2 2 3 3 1]\n" ] } ], "source": [ "X = np.array(['Female', 'Female', 'Female', 'Baby', 'Female', 'Male', 'Male', 'Female', 'Senior', 'Male', 'Male', 'Male', 'Female', 'Female', 'Female', 'Male', 'Male',\n", " 'Senior', 'Senior', 'Female'])\n", "u, rinv = np.unique(X, return_inverse=True)\n", "print(rinv)\n", "\n", "#Alternative way:\n", "from sklearn.preprocessing import LabelEncoder\n", "le = LabelEncoder()\n", "print(le.fit_transform(rinv)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. An array $X_{x_i}$ of length $n$ is called indicator of term $x_i$ if it's $i$-th item is 1 an any other item is 0. Given an array that contains $n$ distinct countries $c_1$, $c_2$, $c_3$, ..., $c_n$ that might repeat. Replace each of them by indicator array of length $n$. Hint: [numpy broadcasting](https://numpy.org/doc/stable/user/basics.broadcasting.html)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1 0 0]\n", " [0 0 1]\n", " [1 0 0]\n", " [1 0 0]\n", " [0 0 1]\n", " [0 0 1]\n", " [1 0 0]\n", " [0 1 0]\n", " [1 0 0]\n", " [1 0 0]\n", " [1 0 0]\n", " [0 0 1]\n", " [1 0 0]\n", " [1 0 0]\n", " [0 0 1]\n", " [0 1 0]\n", " [0 1 0]\n", " [0 0 1]\n", " [0 0 1]\n", " [1 0 0]]\n" ] } ], "source": [ "countries = np.array(['France', 'Spain', 'France', 'France', 'Spain', 'Spain', 'France',\n", " 'Germany', 'France', 'France', 'France', 'Spain', 'France',\n", " 'France', 'Spain', 'Germany', 'Germany', 'Spain', 'Spain',\n", " 'France'])\n", "countries_indicator = (countries[:,None] == np.unique(countries)).astype(int) # creating indicator array\n", "print(countries_indicator)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4. Given two equal sized 1D arrays $X$ and $Y$ with nonnegative items smaller than $m$ and 2D array $Z$ of size $m \\times m$. Return array that contains $Z_{X_i, Y_i}$ for each $i\\le m$." ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0, 6, 18, 23, 24])" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Z = np.arange(25).reshape((5,5))\n", "X, Y = np.array([0, 1, 3, 4, 4]), np.array([0, 1, 3, 3, 4])\n", "Z[X, Y]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. Given two equal sized 1D arrays $X$ and $Y$ with nonnegative items smaller than $m$ and 2D array $Z$ of size $m \\times m$. Return sub-array of $Z$ that contains rows with indexes taken from $X$ and columns that contains indexes taken from $Y$." ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 1, 3, 3, 4],\n", " [ 5, 6, 8, 8, 9],\n", " [15, 16, 18, 18, 19],\n", " [20, 21, 23, 23, 24],\n", " [20, 21, 23, 23, 24]])" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Z = np.arange(25).reshape((5,5))\n", "X, Y = np.array([0, 1, 3, 4, 4]), np.array([0, 1, 3, 3, 4])\n", "Z[np.ix_(X, Y)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Vectorise Python function `product(*arrays)`" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(4, 2, 0),\n", " (4, 2, 1),\n", " (4, 2, 2),\n", " (4, 4, 0),\n", " (4, 4, 1),\n", " (4, 4, 2),\n", " (4, 0, 0),\n", " (4, 0, 1),\n", " (4, 0, 2),\n", " (3, 2, 0),\n", " (3, 2, 1),\n", " (3, 2, 2),\n", " (3, 4, 0),\n", " (3, 4, 1),\n", " (3, 4, 2),\n", " (3, 0, 0),\n", " (3, 0, 1),\n", " (3, 0, 2),\n", " (1, 2, 0),\n", " (1, 2, 1),\n", " (1, 2, 2),\n", " (1, 4, 0),\n", " (1, 4, 1),\n", " (1, 4, 2),\n", " (1, 0, 0),\n", " (1, 0, 1),\n", " (1, 0, 2)]" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def product1(arrays):\n", " return np.stack(np.meshgrid(*arrays), axis=-1).reshape(-1, len(arrays))\n", "def product2(X):\n", " items = [np.array(item) for item in X]\n", " idx = np.where(np.eye(len(X)), Ellipsis, None)\n", " out = [x[tuple(i)] for x,i in zip(items, idx)]\n", " return list(np.broadcast(*out))\n", "\n", "product2([[4,3,1],[2,4,0],[0,1,2]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "7. Convert iterable `combinations(range(n), k)` to `numpy` array with fastest possible performance. [Source](https://stackoverflow.com/a/63694661/3044825)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2],\n", " [0, 1, 3],\n", " [0, 1, 4],\n", " [0, 2, 3],\n", " [0, 2, 4],\n", " [0, 3, 4],\n", " [1, 2, 3],\n", " [1, 2, 4],\n", " [1, 3, 4],\n", " [2, 3, 4]], dtype=int64)" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from itertools import combinations\n", "\n", "def fromiter(iterable, k):\n", " dt = np.dtype([('', np.intp)]*k) #or np.dtype(','.join('i'*5))\n", " indices = np.fromiter(iterable, dt)\n", " indices = indices.view(np.intp).reshape(-1, k)\n", " return indices\n", "\n", "fromiter(combinations(range(5), 3), 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "8. Given boolean array, return indices where subarray `[True, False, True]` appears." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def rolling(a, window):\n", " shape = (a.size - window + 1, window)\n", " strides = (a.itemsize, a.itemsize)\n", " return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)\n", "\n", "compare_with = [True, False, True]\n", "bool_arr = np.random.choice([True, False], size=20)\n", "rolls = rolling(bool_arr, len(compare_with))\n", "idx = np.where(np.all(rolls == compare_with, axis=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "9. Solve Diophantine equation `3x+7y+10z=60`" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 0 0 6]\n", " [10 0 3]\n", " [20 0 0]\n", " [ 1 1 5]\n", " [11 1 2]\n", " [ 2 2 4]\n", " [12 2 1]\n", " [ 3 3 3]\n", " [13 3 0]\n", " [ 4 4 2]\n", " [ 5 5 1]\n", " [ 6 6 0]]\n" ] } ], "source": [ "l = np.array([3,7,10])\n", "s = 60\n", "unknowns = [range(0, s+1, n) for n in l]\n", "triplets = product1(unknowns) #defined as an exerscise before\n", "suitable_triplets = triplets[np.sum(triplets, axis=1) == s]\n", "solutions = suitable_triplets//l\n", "print(solutions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10. Let $L_i(A)$ define the first item of array `A` which doesn't exceed `i`. A sequence $L_{0}(Y), L_{1}(Y), ..., L_{len(X)-1}(Y)$ is called [`left_argmins`](https://stackoverflow.com/questions/63476467/find-argmin-to-the-left-right-of-each-index-in-array) of `X` if `Y` is a sequence of indices that would sort an array `X`. Find `leftt_argmins` of `X`. What about `right_argmins` of `X`?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = np.array([7, 5, 4, 1, 9, 2, 3, 6])\n", "argidx = np.argsort(x)\n", "u, idx = np.unique(argidx, return_index=True)\n", "left_argmins = argidx[np.minimum.accumulate(idx)]\n", "right_argmins = argidx[np.minimum.accumulate(idx[::-1])][::-1]\n", "print(left_argmins)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "11. * Given an array `X` and replacement array `Y` of shape `(n, 2)`. Replace all the items of `X` that have corresponding mappings available. [Source](https://stackoverflow.com/questions/53155749/replace-elements-in-numpy-array-avoiding-loops)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, 2, 200, 5, 400, 400]\n" ] } ], "source": [ "X = np.array([0,1,2,5,4,4])\n", "Y = np.array([[0,1],[1,2],[3,300],[4,400],[2,200]])\n", "\n", "def replacer(X, Y):\n", " #not vectorized yet\n", " print([Y[Y[:,0]==i, 1][0] if i in Y[:,0] else i for i in X])\n", "replacer(X, Y)\n", " \n", "# def replacer(X, Y):\n", "# sidx = Y[:,0].argsort()\n", "# print('sidx', sidx)\n", "# sorted_indx = np.searchsorted(Y[:,0], X, sorter=sidx) #where to add items from X in sorted Y to keep ascending order\n", "# print('sorted_indx', sorted_indx)\n", "# sorted_indx[sorted_indx==len(sidx)] = -1\n", "# print('sorted_indx', sorted_indx)\n", "# idx_out = sidx[sorted_indx]\n", "# print('idx_out', idx_out)\n", "# out = Y[idx_out,1]\n", "# print('out', out)\n", "# out[Y[idx_out, 0]!=X] = 0\n", "# print('out', out)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Grouping" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Extract groups of array by first column." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[array([219, 275, 441, 455]), array([468, 494, 533]), array([559, 593, 605]), array([613, 679, 686, 692])]\n", "[ 4 7 10 14]\n" ] } ], "source": [ "a = np.array([[1, 275], [1, 441], [1, 494], [1, 593], [2, 679], [2, 533], [2, 686], [3, 559], [3, 219], [3, 455], [4, 605], [4, 468], [4, 692], [4, 613]])\n", "\n", "#Way1\n", "a.sort(axis=0)\n", "u, idx = np.unique(a[:,0], return_index = True) #idx = [0, 4, 7, 10]\n", "print(np.split(a[:,1], idx)[1:])\n", "\n", "#Way2\n", "a.sort(axis=0)\n", "u, cnt = np.unique(a[:, 0], return_counts=True) #cnt = [4, 3, 3, 4]\n", "idx = np.cumsum(cnt)[:-1]\n", "print(np.split(a[:,1], idx)[1:])\n", "\n", "#Way3\n", "import numpy_indexed as npi\n", "print(npi.group_by(a[:, 0]).split(a[:, 1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Group array by first column and leave only the rows where second column attains maximum. [Source](https://stackoverflow.com/a/63758412/3044825)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 2]\n", " [1 4]\n", " [2 6]]\n", "[[0 2]\n", " [1 4]\n", " [2 6]]\n" ] } ], "source": [ "import numpy as np\n", "a = np.array([[2,6],[1,4],[0,1],[1,1],[2,3],[0,2]])\n", "a = a[np.argsort(a[:,0])] #sorting lst by first row\n", "u, idx = np.unique(a[:,0], return_index = True) \n", "print(np.c_[u, np.maximum.reduceat(a[:,1], idx)])\n", "\n", "import numpy_indexed as npi\n", "print(np.transpose(npi.group_by(a[:, 0]).max(a[:, 1])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Group array by first column and sum groups by second column (not using reduceat). [Source](https://stackoverflow.com/a/63750626/3044825)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 1.2, 40. ],\n", " [ 2.3, 27. ]])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def group_by_sum(x):\n", " u, idx = np.unique(x[:,0], return_inverse=True)\n", " s = np.bincount(idx, weights = x[:,1])\n", " return np.c_[u, s]\n", "\n", "x = np.array([[1.2, 10], [2.3, 20], [1.2, 30], [2.3, 7]])\n", "#group_by_sum(x)\n", "np.transpose(npi.group_by(x[:, 0]).sum(x[:, 1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "11. Let `X` and `Y` be two polynomials. Each of the polynomials consists of two rows: powers and coeffs, no zero coeffs. Return a new polynomial `X` $\\times$ `Y`" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 1 2 3]\n", " [1 2 2 1]]\n" ] } ], "source": [ "import numpy_indexed as npi\n", "X = np.array([[0,1],[1,1]])\n", "Y = np.array([[0,1,2],[1,1,1]])\n", "cells = product1([X[0], Y[0]]) #finding pairs of powers\n", "degs = np.sum(cells, axis=1) #summing each pair to get new power\n", "coeffs = X[1][cells[:, 0]] * Y[1][cells[:, 1]] #multiplying corresponding coeffs\n", "terms = np.column_stack([coeffs, degs]) #result as terms\n", "print(np.array(npi.group_by(terms[:, 1]).sum(terms[:,0]))) #reduce similar terms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "15. Given a bunch of groups of ids, return:\n", " - list of connected groups (numpy arrays)\n", " - ids that tells which groups are connected\n", " \n", " Source: [my ccarray compilation](..\\ccarray.pdf)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([array([1, 2, 3, 4, 5])], dict_values([array([0, 1, 2, 3], dtype=int64)]))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "import pandas as pd\n", "import igraph as ig\n", "\n", "def connect_arrays(arrays, return_groups=True):\n", " conc = np.concatenate(arrays)\n", " u, indices = np.unique(conc, return_inverse=True)\n", " bars = [len(n) for n in arrays[:-1]]\n", " bars[0] -= 1\n", " mask = np.ones(len(indices) - 1, dtype=bool)\n", " mask[np.cumsum(bars)] = False\n", " nodes = np.arange(len(u))\n", " edges = (np.array([indices[1:], indices[:-1]]).T)[mask]\n", " graph = ig.Graph()\n", " graph.add_vertices(nodes)\n", " graph.add_edges(edges)\n", " graph_tags = graph.clusters().membership\n", " values = pd.DataFrame(graph_tags).groupby([0]).indices.values()\n", " \n", " if return_groups:\n", " graph_tags = np.array(graph_tags)\n", " pid = [n[0] if len(n) else u[-1] + 1 for n in\n", " arrays] # if items is zero, add non-existent item larger than any\n", " group_ix = np.searchsorted(u, pid)\n", " bad_indexes = group_ix == len(u)\n", " group_ix[bad_indexes] = 0 # stuff I forced to do to avoid index error\n", " group_tags = graph_tags[group_ix]\n", " group_tags = group_tags.astype(float)\n", " group_tags[bad_indexes] = np.nan\n", " groups = pd.DataFrame(group_tags).groupby([0]).indices.values()\n", " return [u[k] for k in values], groups\n", " else:\n", " return [u[k] for k in values]\n", "\n", "connect_arrays([[1,2],[2,3],[3,4],[4,5]], return_groups=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Images\n", "\n", "Let us generate a random `image` using these lines:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0]\n", " [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0]\n", " [0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0]\n", " [0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0]\n", " [0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0]\n", " [0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0]\n", " [0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0]\n", " [0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0]\n", " [0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0]\n", " [0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0]\n", " [0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0]\n", " [0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0]\n", " [0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0]\n", " [0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0]\n", " [0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0]\n", " [0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0]\n", " [0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0]\n", " [0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0]\n", " [0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0]\n", " [0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0]\n", " [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]]\n" ] } ], "source": [ "x, y = np.meshgrid(np.linspace(-np.pi/2, np.pi/2, 30), np.linspace(-np.pi/2, np.pi/2, 30))\n", "image = np.sin(x**2+y**2)[:-1,:-1] > 0.9\n", "print(image.astype(int))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "12. Given a binary `image`, return array of its contour positions (pip instal opencv-python instead of conda install)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[11 2]\n", " [10 3]\n", " [ 9 3]\n", " [ 8 4]\n", " [ 7 4]\n", " [ 6 5]\n", " [ 5 6]\n", " [ 4 7]\n", " [ 4 8]\n", " [ 3 9]\n", " [ 3 10]\n", " [ 2 11]\n", " [ 2 12]\n", " [ 2 13]\n", " [ 2 14]\n", " [ 2 15]\n", " [ 2 16]\n", " [ 2 17]\n", " [ 2 18]\n", " [ 3 19]\n", " [ 3 20]\n", " [ 4 21]\n", " [ 4 22]\n", " [ 5 23]\n", " [ 6 24]\n", " [ 7 25]\n", " [ 8 25]\n", " [ 9 26]\n", " [10 26]\n", " [11 27]\n", " [12 27]\n", " [13 27]\n", " [14 27]\n", " [15 27]\n", " [16 27]\n", " [17 27]\n", " [18 27]\n", " [19 26]\n", " [20 26]\n", " [21 25]\n", " [22 25]\n", " [23 24]\n", " [24 23]\n", " [25 22]\n", " [25 21]\n", " [26 20]\n", " [26 19]\n", " [27 18]\n", " [27 17]\n", " [27 16]\n", " [27 15]\n", " [27 14]\n", " [27 13]\n", " [27 12]\n", " [27 11]\n", " [26 10]\n", " [26 9]\n", " [25 8]\n", " [25 7]\n", " [24 6]\n", " [23 5]\n", " [22 4]\n", " [21 4]\n", " [20 3]\n", " [19 3]\n", " [18 2]\n", " [17 2]\n", " [16 2]\n", " [15 2]\n", " [14 2]\n", " [13 2]\n", " [12 2]\n", " [12 5]\n", " [13 4]\n", " [14 4]\n", " [15 4]\n", " [16 4]\n", " [17 5]\n", " [18 5]\n", " [19 5]\n", " [20 6]\n", " [21 7]\n", " [22 8]\n", " [23 9]\n", " [24 10]\n", " [24 11]\n", " [24 12]\n", " [25 13]\n", " [25 14]\n", " [25 15]\n", " [25 16]\n", " [24 17]\n", " [24 18]\n", " [24 19]\n", " [23 20]\n", " [22 21]\n", " [21 22]\n", " [20 23]\n", " [19 24]\n", " [18 24]\n", " [17 24]\n", " [16 25]\n", " [15 25]\n", " [14 25]\n", " [13 25]\n", " [12 24]\n", " [11 24]\n", " [10 24]\n", " [ 9 23]\n", " [ 8 22]\n", " [ 7 21]\n", " [ 6 20]\n", " [ 5 19]\n", " [ 5 18]\n", " [ 5 17]\n", " [ 4 16]\n", " [ 4 15]\n", " [ 4 14]\n", " [ 4 13]\n", " [ 5 12]\n", " [ 5 11]\n", " [ 5 10]\n", " [ 6 9]\n", " [ 7 8]\n", " [ 8 7]\n", " [ 9 6]\n", " [10 5]\n", " [11 5]]\n" ] } ], "source": [ "#converting extracting contours\n", "import cv2\n", "image = image.astype(np.uint8)\n", "contours, hierarchy = cv2.findContours(image, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)\n", "cells = np.vstack(contours).squeeze() #change shape\n", "print(cells)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "13. Given a binary image, return locations of all the units." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 2 11]\n", " [ 2 12]\n", " [ 2 13]\n", " [ 2 14]\n", " [ 2 15]\n", " [ 2 16]\n", " [ 2 17]\n", " [ 2 18]\n", " [ 3 9]\n", " [ 3 10]\n", " [ 3 11]\n", " [ 3 12]\n", " [ 3 13]\n", " [ 3 14]\n", " [ 3 15]\n", " [ 3 16]\n", " [ 3 17]\n", " [ 3 18]\n", " [ 3 19]\n", " [ 3 20]\n", " [ 4 7]\n", " [ 4 8]\n", " [ 4 9]\n", " [ 4 10]\n", " [ 4 11]\n", " [ 4 12]\n", " [ 4 13]\n", " [ 4 14]\n", " [ 4 15]\n", " [ 4 16]\n", " [ 4 17]\n", " [ 4 18]\n", " [ 4 19]\n", " [ 4 20]\n", " [ 4 21]\n", " [ 4 22]\n", " [ 5 6]\n", " [ 5 7]\n", " [ 5 8]\n", " [ 5 9]\n", " [ 5 10]\n", " [ 5 11]\n", " [ 5 12]\n", " [ 5 17]\n", " [ 5 18]\n", " [ 5 19]\n", " [ 5 20]\n", " [ 5 21]\n", " [ 5 22]\n", " [ 5 23]\n", " [ 6 5]\n", " [ 6 6]\n", " [ 6 7]\n", " [ 6 8]\n", " [ 6 9]\n", " [ 6 20]\n", " [ 6 21]\n", " [ 6 22]\n", " [ 6 23]\n", " [ 6 24]\n", " [ 7 4]\n", " [ 7 5]\n", " [ 7 6]\n", " [ 7 7]\n", " [ 7 8]\n", " [ 7 21]\n", " [ 7 22]\n", " [ 7 23]\n", " [ 7 24]\n", " [ 7 25]\n", " [ 8 4]\n", " [ 8 5]\n", " [ 8 6]\n", " [ 8 7]\n", " [ 8 22]\n", " [ 8 23]\n", " [ 8 24]\n", " [ 8 25]\n", " [ 9 3]\n", " [ 9 4]\n", " [ 9 5]\n", " [ 9 6]\n", " [ 9 23]\n", " [ 9 24]\n", " [ 9 25]\n", " [ 9 26]\n", " [10 3]\n", " [10 4]\n", " [10 5]\n", " [10 24]\n", " [10 25]\n", " [10 26]\n", " [11 2]\n", " [11 3]\n", " [11 4]\n", " [11 5]\n", " [11 24]\n", " [11 25]\n", " [11 26]\n", " [11 27]\n", " [12 2]\n", " [12 3]\n", " [12 4]\n", " [12 5]\n", " [12 24]\n", " [12 25]\n", " [12 26]\n", " [12 27]\n", " [13 2]\n", " [13 3]\n", " [13 4]\n", " [13 25]\n", " [13 26]\n", " [13 27]\n", " [14 2]\n", " [14 3]\n", " [14 4]\n", " [14 25]\n", " [14 26]\n", " [14 27]\n", " [15 2]\n", " [15 3]\n", " [15 4]\n", " [15 25]\n", " [15 26]\n", " [15 27]\n", " [16 2]\n", " [16 3]\n", " [16 4]\n", " [16 25]\n", " [16 26]\n", " [16 27]\n", " [17 2]\n", " [17 3]\n", " [17 4]\n", " [17 5]\n", " [17 24]\n", " [17 25]\n", " [17 26]\n", " [17 27]\n", " [18 2]\n", " [18 3]\n", " [18 4]\n", " [18 5]\n", " [18 24]\n", " [18 25]\n", " [18 26]\n", " [18 27]\n", " [19 3]\n", " [19 4]\n", " [19 5]\n", " [19 24]\n", " [19 25]\n", " [19 26]\n", " [20 3]\n", " [20 4]\n", " [20 5]\n", " [20 6]\n", " [20 23]\n", " [20 24]\n", " [20 25]\n", " [20 26]\n", " [21 4]\n", " [21 5]\n", " [21 6]\n", " [21 7]\n", " [21 22]\n", " [21 23]\n", " [21 24]\n", " [21 25]\n", " [22 4]\n", " [22 5]\n", " [22 6]\n", " [22 7]\n", " [22 8]\n", " [22 21]\n", " [22 22]\n", " [22 23]\n", " [22 24]\n", " [22 25]\n", " [23 5]\n", " [23 6]\n", " [23 7]\n", " [23 8]\n", " [23 9]\n", " [23 20]\n", " [23 21]\n", " [23 22]\n", " [23 23]\n", " [23 24]\n", " [24 6]\n", " [24 7]\n", " [24 8]\n", " [24 9]\n", " [24 10]\n", " [24 11]\n", " [24 12]\n", " [24 17]\n", " [24 18]\n", " [24 19]\n", " [24 20]\n", " [24 21]\n", " [24 22]\n", " [24 23]\n", " [25 7]\n", " [25 8]\n", " [25 9]\n", " [25 10]\n", " [25 11]\n", " [25 12]\n", " [25 13]\n", " [25 14]\n", " [25 15]\n", " [25 16]\n", " [25 17]\n", " [25 18]\n", " [25 19]\n", " [25 20]\n", " [25 21]\n", " [25 22]\n", " [26 9]\n", " [26 10]\n", " [26 11]\n", " [26 12]\n", " [26 13]\n", " [26 14]\n", " [26 15]\n", " [26 16]\n", " [26 17]\n", " [26 18]\n", " [26 19]\n", " [26 20]\n", " [27 11]\n", " [27 12]\n", " [27 13]\n", " [27 14]\n", " [27 15]\n", " [27 16]\n", " [27 17]\n", " [27 18]]\n", "[[ 2 11]\n", " [ 2 12]\n", " [ 2 13]\n", " [ 2 14]\n", " [ 2 15]\n", " [ 2 16]\n", " [ 2 17]\n", " [ 2 18]\n", " [ 3 9]\n", " [ 3 10]\n", " [ 3 11]\n", " [ 3 12]\n", " [ 3 13]\n", " [ 3 14]\n", " [ 3 15]\n", " [ 3 16]\n", " [ 3 17]\n", " [ 3 18]\n", " [ 3 19]\n", " [ 3 20]\n", " [ 4 7]\n", " [ 4 8]\n", " [ 4 9]\n", " [ 4 10]\n", " [ 4 11]\n", " [ 4 12]\n", " [ 4 13]\n", " [ 4 14]\n", " [ 4 15]\n", " [ 4 16]\n", " [ 4 17]\n", " [ 4 18]\n", " [ 4 19]\n", " [ 4 20]\n", " [ 4 21]\n", " [ 4 22]\n", " [ 5 6]\n", " [ 5 7]\n", " [ 5 8]\n", " [ 5 9]\n", " [ 5 10]\n", " [ 5 11]\n", " [ 5 12]\n", " [ 5 17]\n", " [ 5 18]\n", " [ 5 19]\n", " [ 5 20]\n", " [ 5 21]\n", " [ 5 22]\n", " [ 5 23]\n", " [ 6 5]\n", " [ 6 6]\n", " [ 6 7]\n", " [ 6 8]\n", " [ 6 9]\n", " [ 6 20]\n", " [ 6 21]\n", " [ 6 22]\n", " [ 6 23]\n", " [ 6 24]\n", " [ 7 4]\n", " [ 7 5]\n", " [ 7 6]\n", " [ 7 7]\n", " [ 7 8]\n", " [ 7 21]\n", " [ 7 22]\n", " [ 7 23]\n", " [ 7 24]\n", " [ 7 25]\n", " [ 8 4]\n", " [ 8 5]\n", " [ 8 6]\n", " [ 8 7]\n", " [ 8 22]\n", " [ 8 23]\n", " [ 8 24]\n", " [ 8 25]\n", " [ 9 3]\n", " [ 9 4]\n", " [ 9 5]\n", " [ 9 6]\n", " [ 9 23]\n", " [ 9 24]\n", " [ 9 25]\n", " [ 9 26]\n", " [10 3]\n", " [10 4]\n", " [10 5]\n", " [10 24]\n", " [10 25]\n", " [10 26]\n", " [11 2]\n", " [11 3]\n", " [11 4]\n", " [11 5]\n", " [11 24]\n", " [11 25]\n", " [11 26]\n", " [11 27]\n", " [12 2]\n", " [12 3]\n", " [12 4]\n", " [12 5]\n", " [12 24]\n", " [12 25]\n", " [12 26]\n", " [12 27]\n", " [13 2]\n", " [13 3]\n", " [13 4]\n", " [13 25]\n", " [13 26]\n", " [13 27]\n", " [14 2]\n", " [14 3]\n", " [14 4]\n", " [14 25]\n", " [14 26]\n", " [14 27]\n", " [15 2]\n", " [15 3]\n", " [15 4]\n", " [15 25]\n", " [15 26]\n", " [15 27]\n", " [16 2]\n", " [16 3]\n", " [16 4]\n", " [16 25]\n", " [16 26]\n", " [16 27]\n", " [17 2]\n", " [17 3]\n", " [17 4]\n", " [17 5]\n", " [17 24]\n", " [17 25]\n", " [17 26]\n", " [17 27]\n", " [18 2]\n", " [18 3]\n", " [18 4]\n", " [18 5]\n", " [18 24]\n", " [18 25]\n", " [18 26]\n", " [18 27]\n", " [19 3]\n", " [19 4]\n", " [19 5]\n", " [19 24]\n", " [19 25]\n", " [19 26]\n", " [20 3]\n", " [20 4]\n", " [20 5]\n", " [20 6]\n", " [20 23]\n", " [20 24]\n", " [20 25]\n", " [20 26]\n", " [21 4]\n", " [21 5]\n", " [21 6]\n", " [21 7]\n", " [21 22]\n", " [21 23]\n", " [21 24]\n", " [21 25]\n", " [22 4]\n", " [22 5]\n", " [22 6]\n", " [22 7]\n", " [22 8]\n", " [22 21]\n", " [22 22]\n", " [22 23]\n", " [22 24]\n", " [22 25]\n", " [23 5]\n", " [23 6]\n", " [23 7]\n", " [23 8]\n", " [23 9]\n", " [23 20]\n", " [23 21]\n", " [23 22]\n", " [23 23]\n", " [23 24]\n", " [24 6]\n", " [24 7]\n", " [24 8]\n", " [24 9]\n", " [24 10]\n", " [24 11]\n", " [24 12]\n", " [24 17]\n", " [24 18]\n", " [24 19]\n", " [24 20]\n", " [24 21]\n", " [24 22]\n", " [24 23]\n", " [25 7]\n", " [25 8]\n", " [25 9]\n", " [25 10]\n", " [25 11]\n", " [25 12]\n", " [25 13]\n", " [25 14]\n", " [25 15]\n", " [25 16]\n", " [25 17]\n", " [25 18]\n", " [25 19]\n", " [25 20]\n", " [25 21]\n", " [25 22]\n", " [26 9]\n", " [26 10]\n", " [26 11]\n", " [26 12]\n", " [26 13]\n", " [26 14]\n", " [26 15]\n", " [26 16]\n", " [26 17]\n", " [26 18]\n", " [26 19]\n", " [26 20]\n", " [27 11]\n", " [27 12]\n", " [27 13]\n", " [27 14]\n", " [27 15]\n", " [27 16]\n", " [27 17]\n", " [27 18]]\n", "[[ 2 11]\n", " [ 2 12]\n", " [ 2 13]\n", " [ 2 14]\n", " [ 2 15]\n", " [ 2 16]\n", " [ 2 17]\n", " [ 2 18]\n", " [ 3 9]\n", " [ 3 10]\n", " [ 3 11]\n", " [ 3 12]\n", " [ 3 13]\n", " [ 3 14]\n", " [ 3 15]\n", " [ 3 16]\n", " [ 3 17]\n", " [ 3 18]\n", " [ 3 19]\n", " [ 3 20]\n", " [ 4 7]\n", " [ 4 8]\n", " [ 4 9]\n", " [ 4 10]\n", " [ 4 11]\n", " [ 4 12]\n", " [ 4 13]\n", " [ 4 14]\n", " [ 4 15]\n", " [ 4 16]\n", " [ 4 17]\n", " [ 4 18]\n", " [ 4 19]\n", " [ 4 20]\n", " [ 4 21]\n", " [ 4 22]\n", " [ 5 6]\n", " [ 5 7]\n", " [ 5 8]\n", " [ 5 9]\n", " [ 5 10]\n", " [ 5 11]\n", " [ 5 12]\n", " [ 5 17]\n", " [ 5 18]\n", " [ 5 19]\n", " [ 5 20]\n", " [ 5 21]\n", " [ 5 22]\n", " [ 5 23]\n", " [ 6 5]\n", " [ 6 6]\n", " [ 6 7]\n", " [ 6 8]\n", " [ 6 9]\n", " [ 6 20]\n", " [ 6 21]\n", " [ 6 22]\n", " [ 6 23]\n", " [ 6 24]\n", " [ 7 4]\n", " [ 7 5]\n", " [ 7 6]\n", " [ 7 7]\n", " [ 7 8]\n", " [ 7 21]\n", " [ 7 22]\n", " [ 7 23]\n", " [ 7 24]\n", " [ 7 25]\n", " [ 8 4]\n", " [ 8 5]\n", " [ 8 6]\n", " [ 8 7]\n", " [ 8 22]\n", " [ 8 23]\n", " [ 8 24]\n", " [ 8 25]\n", " [ 9 3]\n", " [ 9 4]\n", " [ 9 5]\n", " [ 9 6]\n", " [ 9 23]\n", " [ 9 24]\n", " [ 9 25]\n", " [ 9 26]\n", " [10 3]\n", " [10 4]\n", " [10 5]\n", " [10 24]\n", " [10 25]\n", " [10 26]\n", " [11 2]\n", " [11 3]\n", " [11 4]\n", " [11 5]\n", " [11 24]\n", " [11 25]\n", " [11 26]\n", " [11 27]\n", " [12 2]\n", " [12 3]\n", " [12 4]\n", " [12 5]\n", " [12 24]\n", " [12 25]\n", " [12 26]\n", " [12 27]\n", " [13 2]\n", " [13 3]\n", " [13 4]\n", " [13 25]\n", " [13 26]\n", " [13 27]\n", " [14 2]\n", " [14 3]\n", " [14 4]\n", " [14 25]\n", " [14 26]\n", " [14 27]\n", " [15 2]\n", " [15 3]\n", " [15 4]\n", " [15 25]\n", " [15 26]\n", " [15 27]\n", " [16 2]\n", " [16 3]\n", " [16 4]\n", " [16 25]\n", " [16 26]\n", " [16 27]\n", " [17 2]\n", " [17 3]\n", " [17 4]\n", " [17 5]\n", " [17 24]\n", " [17 25]\n", " [17 26]\n", " [17 27]\n", " [18 2]\n", " [18 3]\n", " [18 4]\n", " [18 5]\n", " [18 24]\n", " [18 25]\n", " [18 26]\n", " [18 27]\n", " [19 3]\n", " [19 4]\n", " [19 5]\n", " [19 24]\n", " [19 25]\n", " [19 26]\n", " [20 3]\n", " [20 4]\n", " [20 5]\n", " [20 6]\n", " [20 23]\n", " [20 24]\n", " [20 25]\n", " [20 26]\n", " [21 4]\n", " [21 5]\n", " [21 6]\n", " [21 7]\n", " [21 22]\n", " [21 23]\n", " [21 24]\n", " [21 25]\n", " [22 4]\n", " [22 5]\n", " [22 6]\n", " [22 7]\n", " [22 8]\n", " [22 21]\n", " [22 22]\n", " [22 23]\n", " [22 24]\n", " [22 25]\n", " [23 5]\n", " [23 6]\n", " [23 7]\n", " [23 8]\n", " [23 9]\n", " [23 20]\n", " [23 21]\n", " [23 22]\n", " [23 23]\n", " [23 24]\n", " [24 6]\n", " [24 7]\n", " [24 8]\n", " [24 9]\n", " [24 10]\n", " [24 11]\n", " [24 12]\n", " [24 17]\n", " [24 18]\n", " [24 19]\n", " [24 20]\n", " [24 21]\n", " [24 22]\n", " [24 23]\n", " [25 7]\n", " [25 8]\n", " [25 9]\n", " [25 10]\n", " [25 11]\n", " [25 12]\n", " [25 13]\n", " [25 14]\n", " [25 15]\n", " [25 16]\n", " [25 17]\n", " [25 18]\n", " [25 19]\n", " [25 20]\n", " [25 21]\n", " [25 22]\n", " [26 9]\n", " [26 10]\n", " [26 11]\n", " [26 12]\n", " [26 13]\n", " [26 14]\n", " [26 15]\n", " [26 16]\n", " [26 17]\n", " [26 18]\n", " [26 19]\n", " [26 20]\n", " [27 11]\n", " [27 12]\n", " [27 13]\n", " [27 14]\n", " [27 15]\n", " [27 16]\n", " [27 17]\n", " [27 18]]\n", "[[ 2 11]\n", " [ 2 12]\n", " [ 2 13]\n", " [ 2 14]\n", " [ 2 15]\n", " [ 2 16]\n", " [ 2 17]\n", " [ 2 18]\n", " [ 3 9]\n", " [ 3 10]\n", " [ 3 11]\n", " [ 3 12]\n", " [ 3 13]\n", " [ 3 14]\n", " [ 3 15]\n", " [ 3 16]\n", " [ 3 17]\n", " [ 3 18]\n", " [ 3 19]\n", " [ 3 20]\n", " [ 4 7]\n", " [ 4 8]\n", " [ 4 9]\n", " [ 4 10]\n", " [ 4 11]\n", " [ 4 12]\n", " [ 4 13]\n", " [ 4 14]\n", " [ 4 15]\n", " [ 4 16]\n", " [ 4 17]\n", " [ 4 18]\n", " [ 4 19]\n", " [ 4 20]\n", " [ 4 21]\n", " [ 4 22]\n", " [ 5 6]\n", " [ 5 7]\n", " [ 5 8]\n", " [ 5 9]\n", " [ 5 10]\n", " [ 5 11]\n", " [ 5 12]\n", " [ 5 17]\n", " [ 5 18]\n", " [ 5 19]\n", " [ 5 20]\n", " [ 5 21]\n", " [ 5 22]\n", " [ 5 23]\n", " [ 6 5]\n", " [ 6 6]\n", " [ 6 7]\n", " [ 6 8]\n", " [ 6 9]\n", " [ 6 20]\n", " [ 6 21]\n", " [ 6 22]\n", " [ 6 23]\n", " [ 6 24]\n", " [ 7 4]\n", " [ 7 5]\n", " [ 7 6]\n", " [ 7 7]\n", " [ 7 8]\n", " [ 7 21]\n", " [ 7 22]\n", " [ 7 23]\n", " [ 7 24]\n", " [ 7 25]\n", " [ 8 4]\n", " [ 8 5]\n", " [ 8 6]\n", " [ 8 7]\n", " [ 8 22]\n", " [ 8 23]\n", " [ 8 24]\n", " [ 8 25]\n", " [ 9 3]\n", " [ 9 4]\n", " [ 9 5]\n", " [ 9 6]\n", " [ 9 23]\n", " [ 9 24]\n", " [ 9 25]\n", " [ 9 26]\n", " [10 3]\n", " [10 4]\n", " [10 5]\n", " [10 24]\n", " [10 25]\n", " [10 26]\n", " [11 2]\n", " [11 3]\n", " [11 4]\n", " [11 5]\n", " [11 24]\n", " [11 25]\n", " [11 26]\n", " [11 27]\n", " [12 2]\n", " [12 3]\n", " [12 4]\n", " [12 5]\n", " [12 24]\n", " [12 25]\n", " [12 26]\n", " [12 27]\n", " [13 2]\n", " [13 3]\n", " [13 4]\n", " [13 25]\n", " [13 26]\n", " [13 27]\n", " [14 2]\n", " [14 3]\n", " [14 4]\n", " [14 25]\n", " [14 26]\n", " [14 27]\n", " [15 2]\n", " [15 3]\n", " [15 4]\n", " [15 25]\n", " [15 26]\n", " [15 27]\n", " [16 2]\n", " [16 3]\n", " [16 4]\n", " [16 25]\n", " [16 26]\n", " [16 27]\n", " [17 2]\n", " [17 3]\n", " [17 4]\n", " [17 5]\n", " [17 24]\n", " [17 25]\n", " [17 26]\n", " [17 27]\n", " [18 2]\n", " [18 3]\n", " [18 4]\n", " [18 5]\n", " [18 24]\n", " [18 25]\n", " [18 26]\n", " [18 27]\n", " [19 3]\n", " [19 4]\n", " [19 5]\n", " [19 24]\n", " [19 25]\n", " [19 26]\n", " [20 3]\n", " [20 4]\n", " [20 5]\n", " [20 6]\n", " [20 23]\n", " [20 24]\n", " [20 25]\n", " [20 26]\n", " [21 4]\n", " [21 5]\n", " [21 6]\n", " [21 7]\n", " [21 22]\n", " [21 23]\n", " [21 24]\n", " [21 25]\n", " [22 4]\n", " [22 5]\n", " [22 6]\n", " [22 7]\n", " [22 8]\n", " [22 21]\n", " [22 22]\n", " [22 23]\n", " [22 24]\n", " [22 25]\n", " [23 5]\n", " [23 6]\n", " [23 7]\n", " [23 8]\n", " [23 9]\n", " [23 20]\n", " [23 21]\n", " [23 22]\n", " [23 23]\n", " [23 24]\n", " [24 6]\n", " [24 7]\n", " [24 8]\n", " [24 9]\n", " [24 10]\n", " [24 11]\n", " [24 12]\n", " [24 17]\n", " [24 18]\n", " [24 19]\n", " [24 20]\n", " [24 21]\n", " [24 22]\n", " [24 23]\n", " [25 7]\n", " [25 8]\n", " [25 9]\n", " [25 10]\n", " [25 11]\n", " [25 12]\n", " [25 13]\n", " [25 14]\n", " [25 15]\n", " [25 16]\n", " [25 17]\n", " [25 18]\n", " [25 19]\n", " [25 20]\n", " [25 21]\n", " [25 22]\n", " [26 9]\n", " [26 10]\n", " [26 11]\n", " [26 12]\n", " [26 13]\n", " [26 14]\n", " [26 15]\n", " [26 16]\n", " [26 17]\n", " [26 18]\n", " [26 19]\n", " [26 20]\n", " [27 11]\n", " [27 12]\n", " [27 13]\n", " [27 14]\n", " [27 15]\n", " [27 16]\n", " [27 17]\n", " [27 18]]\n" ] } ], "source": [ "x, y = np.indices(image.shape)\n", "mask = image == 1\n", "location = np.array([x[mask], y[mask]]).T\n", "print(location)\n", "\n", "#Alternative way 1:\n", "\n", "from scipy import sparse\n", "sparsed_image = sparse.csr_matrix(image)\n", "print(np.array(sparsed_image.nonzero()).T)\n", "\n", "#Alternative way 2:\n", "\n", "print(np.array(image.nonzero()).T)\n", "\n", "#Alternative way 3:\n", "\n", "print(np.array(np.where(image)).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "14. An image is made of $29 \\times 29$ array. Convert it to a more efficient memory data structure." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "from scipy import sparse\n", "sparsed_image = sparse.csr_matrix(image) #this is a data structure that stores rows and columns only, 30% of cells only\n", "image2 = sparsed_image.todense() #convert back to array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Save and retrieve it" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "9. Save array and retrieve it" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.random.randint(0,10, size=(100,100))\n", "np.save(file='data.npy', x)\n", "retrieved_array = np.load('data.npy')\n", "\n", "np.savez_compressed('data.npz', array=x) #not compressed: data.npz\n", "retrieved_array = np.load('data.npz')['array']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10. Dowload a huge .npz file with 3 column array straight from github and perform dimensionality reduction on it using `numexpr`." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(19488504,)\n" ] } ], "source": [ "import tensorflow as tf\n", "import numexpr as ne\n", "path = tf.keras.utils.get_file('cubes.npz', 'https://github.com/loijord/lidar_home/raw/master/cubes.npz')\n", "cubes = np.load(path)['array']\n", "cubes = cubes.astype(np.int64)\n", "s0, s1 = cubes[:,0].max()+1, cubes[:,1].max()+1\n", "d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}\n", "c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)\n", "print(c1D.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "11. Now you've got `C1D` array, can you show a list of groups as like as you do with `pd.DataFrame(cubes).groupby([0]).indices`? [Source](https://stackoverflow.com/q/59239886/3044825)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Future asks:\n", "* [About justification](https://github.com/numpy/numpy/issues/15863)\n", "* [Fill missing values according to frequency](https://stackoverflow.com/questions/63548053/)\n", "* [Link lists that share common elements](https://stackoverflow.com/questions/63583725/)\n", "* [Replace groups of cells by their averages](https://stackoverflow.com/questions/63618938)\n", "* [how to replace all the items of array by indexes of specified lists](stackoverflow.com/questions/63623869/) also used in tagging cc in networkx\n", "* [`np.isin` in 2D case](https://stackoverflow.com/questions/51352527) [update: it is needed to solve problem in conjunction with product]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# More future asks:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Kerneling image\n", "* Zoom of image\n", "* Visualiser of numpy array dimensions, [like in this source](https://papers.bastolasushan.com/what-are-tensors/), or [this](https://beerensahu.wordpress.com/2018/03/21/pytorch-tutorial-lesson-1-tensor/)\n", "* Construct a kernel from given values of its radius.\n", "* Interactive (clickable) graphs.\n", "* Create graph from 3D sparse matrix, actually pointcloud.\n", "* How to intersect 3d matrices faster that it does `numpy` indexed (assuming I understand what is view, ravel, etc.)\n", "* Vectorise `[A.index(b) for b in B]`\n", "* Review of nxviz library\n", "* Review of numpy_indexed library\n", "* Study @Divakar's post with `numba` and `np.ravel()` [Source](https://stackoverflow.com/q/59239886/3044825)\n", "* Study [answers](https://stackoverflow.com/q/63623869/3044825) of how to replace all the items of array by indexes of specified lists? They contain `np.apply_along_axis`, `np.searchsorted` and similar tricks." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A wonderful [explanation](https://towardsdatascience.com/numpy-indexing-explained-c376abb2440d) of `numpy` indexing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Illustration of 3D arrays" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZxcZZnvv7V0VXdXr+kt6azdnc7SSToJSUiCCMyAKA6DF0VmFNdR0VFcPo4z44xzZwYZARnFi44zgg6iFxeE60VAUSIG4bIkgaQTQkKS3tf0XntV13Lq/vHWe/rUfqr7pJOG+n0+9QmcOvWeU9Xn/M7zPu/v+T2mWCxGAQUUUEABCwPz+T6BAgoooIA3EwqkW0ABBRSwgCiQbgEFFFDAAqJAugUUUEABC4gC6RZQQAEFLCCsOd4vSBsKKKCAAvKHKdMbhUi3gAIKKGABUSDdAgoooIAFRIF0CyiggAIWEAXSLaCAAgpYQBRIt4ACCihgAVEg3QIKKKCABUSBdAsooIACFhAF0i2ggAIKWEAUSLeAAgooYAFRIN0CCiiggAVEgXQLKKCAAhYQBdItoIACClhAFEi3gAIKKGABkctlrIACMiIWi6EoCjMzM0QiEaxWK2azGYvFgtlsxmw2YzJlNFsqoIA3JUw5GlMWrB0LSEEsFiMajRKJRBL+W76nJVpJwvJVIOMC3iTIeIEXSLcA3UgmW5PJRDAYpKenh0gkQllZGaWlpTgcDmw2m/oZ+err66O+vp7S0tICGRfwRkfGC7mQXiggJ2KxGJFIhMHBQSoqKnA4HPj9frq7u/H7/axatQqLxUIgEGBycpL+/n5CoRAWiwWHw4HD4aC0tJRAIEAsFsNsFksJkUiEcDiccKwCGRfwRkch0i0gIyTZytTBqVOnKCsrY2JigkgkQlNTEzU1NWr0K8lUIhKJ4Pf78fl8+Hw+RkdHMZvNFBUVJZCxw+GguLg44bjyJcnWZDJhsVjUvLEk5wIZF3CBohDpFqAfiqIk5GkBnE4nY2NjuFwuNmzYQFVVlfpeJuKzWq1UVFRQUVEBCBJubGxUI2Wfz4fL5WJ4eJhgMIjZbFZJWBJySUmJOl40GqW/vx+AxsZGNWrOFBkXCLmACxEF0i1AhaIoRCIRotGoum1iYoKenh6Ki4upqamhoaEhgXDzgclkIhaLYbFYKC8vp7y8POH9aDRKIBDA5/PhdrsZGRkhGAwCqGQcDAax2+2YTCY1sk5+QMhjmc1mrFZrgYwLuKBQIN03OeQ0PhwOoyiKum10dJS+vj4qKirYsmULpaWldHZ2kiMdlRWSdDPBYrFQVlZGWVlZwnZFUVQydrlceDweJiYmACgpKUmIjEtLS1VSlXK25HMwmUzqwl+BjAtYaBRI900KqbGNRCKcPHmS9evXE4vFGB4epr+/n5qaGrZv356QazWZTCoxzwW5SDcTzGazSqxy4W358uUoikIwGFRzxuPj4/j9fmKxWAIZOxwOSkpKsFgsap766NGjbN++PeE42jSFjI4LZFyA0SiQ7psMWrJVFAWTycTU1BR9fX0MDQ2xdOlSdu3apUq+tJgraZ4ryBxwaWkpdXV16vZYLEYgEFDzxpOTk/j9fhRFobi4mJKSEsLhMH6/n9LSUpWMZS47WWucLmdcUFQUMFcUSPdNgnQa23A4zMDAAD6fD5PJxJ49e7BaM18SZrP5nKYX9CLXGCaTSSXj2trahM8Fg0E8Hg/j4+MMDg7i9/uJRqPY7faEyLi0tBSr1ZqVjGdmZrDZbGoUXSDjAvSgQLpvcKQj21AoRG9vL5OTk6xatQqHw8GaNWtyjnW+0gtGjWEymSgpKaGoqAi73c7GjRsB8RuFQiG8Xi9+v5/h4WF8Ph/RaBSbzZZCxkVFRcRiMc6ePYvD4aCmpqZQhVeAbhRI9w0KmbvURmiyesztdrNmzRrWrVuH2WxmYGBA15hGkOaFkp7QEqDJZMJut2O326mpqVG3SzL2+Xz4/X7Onj2Lz+cjEolQVFSkFnfIKFmSMRQKPwrIjALpvsEgyXZgYACHw0FlZSU+n4/u7m6CwSBNTU20tbXN6UY3m81pI129Yy02ctGS8ZIlSxLeC4VCdHZ2AjA2NobP5yMcDmO1WhMiY0nGEunIeHp6mpqaGmw2W6Hw402AAum+QZBc0OD3+5mZmaG7uxtFUWhubqa6unpeN/KFkF4wAkacg81mw263U1VVlRAdh8NhVU0xMTFBX19fSkm0fMnFyt7eXqqqqggGg+qspFD48cZFgXQXOdIVNExNTTEyMkJRURFtbW1UVlbmHCd5kSgdkkkzFosxNTWlRtHJi1EOhyNhYW6hFtL0wAjSSvebFRUVUVVVlVJAEolEVDKemppiYGCAmZkZ1bNiaGiIsrIyHA4Hdrtd/Vyh8OONhwLpLkKkK2gAGB8fp6enB4fDQUNDA2VlZboIV6YNLBZLzv3ksScmJuju7qakpIT169erOU5JLCMjI+pilCTjQCAAoHrvzgUXEqnoeVBJWK1WKisrU/4e0WiUQ4cOYbPZcDqdDA0NpS2JTvanSFf4MTY2xrJly1I8KgpkfGGhQLqLCMkaW7nt7Nmz9PX1UVVVxdatWykpKaG/v193KiAfKZjb7ebAgQOUlZWxefNmHA4HiqIQCoXS5j+1ygAZ5U1OTqaVaTkcjpzEbwSMSnHkQ7qZIFMHy5YtS9gejUYT/ClGRkYIBAKqHC658MNkMjEwMEBDQwPRaJRQKJT2OIXCj/OPAukuAkiydTqdeDweli5diqIoDA0NMTg4SF1dHTt27EiYlprN5oSUQzbkytVKYj9z5gxWq5Xt27cnGNHkGluSsc/nw2azsXTpUmKxGDMzM2pkPDQ0hM/nQ1GUBDKWHr1Gk/G5Si8YhUz+FIqiqGTs8XgYHR1VZxDBYJC+vr4EMpYP1GStsfy3UPix8CiQ7gWMZI1tKBRicnKSQCDA8PAwy5Yt4+KLL05YHZcwm80pq+SZkEmVoCgKIyMj9PX1sWTJElpaWvD5fLoJNxnanK7JZKK4uFg10tF+52AwqBLLwMBAQjVZLBbDarXi8XjOCRnng3NJuplgNpsz+lMcPHgQh8ORUBINUFxcnKI1zlaFd/bsWZYtW6bmjAvyNmNRIN0LEOkKGiKRCENDQ4yNjbF27Vr27t2blXAyEamefWUU3d/fT11dHTt37sRmszE5OYnH40k7hlE3oyxgKCkpSUvGAwMDBAKBFDJOjoyTvX21MIoszwfpZoJUOtTX1ydsT/anmJiYIBAIpPxuWjIeGhpi2bJlhMNhQqFQofDDYBRI9wJCuoKGmZkZent7mZ6epqGhgZqaGlavXp1zrHxIV0ag0WiUwcFBBgcHaWhoSImi56s+MKKarLy8nOLiYlatWgXMkrHP58Pr9ao+C+lMb3KRcb44H6S7f/9+PvyxD2OymNiyfgvbt2/nsssu4y1veUva/bP5U2jJeGpqSn2IBQIBuru705ZEQ6HwY74okO4FgOQODSaTSb3wfT4fa9asYcOGDQQCAdxut64x8yFdgP7+fiYnJ2lsbGT37t1p1QXZFtzmIjmbC5KPoY2Mk30WpB2kjPAkGdvtdvx+P6Ojo/Mi44Uk3Y6ODj78sQ/T09MDm4EKeHbqWZ79xbPcc989EASz3cySJUtY37ye7du2c+mll3LFFVdQWlqaMl623+3gwYNUV1eruXatP0XyIl42Mp6amlJnHunkbW9WFEj3PEJGFbIiyWQy4fF46O7uJhQK0dzcnFDXP5+UQTqEw2GVbJcvX54zZZFtwU3vTbRQxRFa05vkCM/pdNLZ2UkgEFDJGEhrB5mLjM81efT29vLBj3yQo8eOwnrgRqA4zY4RUFwKE84JJqYmeP6Xz/Mf9/8HBMBit1BZVamS8Vve8hb+5E/+JCUvLL+P2WympqYmY0l0siSwqKgooSmprMKbnJxUi0i0qTJt4YdW2vZmUVQUSPc8QFvQIA1WLBYL3d3dmEwmtXosGRaLRbciwWKxZCTIUChEX18fY2NjrFy5kqVLl7J06dKci1KZItV8yoDPd0WajPDsdnuCyY/WKN3n8zE2NqaqAjKRsVGFHunGGR8f56Mf+yjP/r9noQl4D+DIMpAVqIm/WjTbIxB1RZlyTvHi9Iu8+KsX+c8H/hP84jMljhK2bNiikvHll1+edvhMJdFSL67tgyf9KUKhENFolMrKyoQqPPl936yFHwXSXSBkKmhwu91qF4T169enSIS0yDfSTSbomZkZenp6mJqaYtWqVezduxez2czJkyd1jXuhWDvOF+nOQWuUrkU2Mg6FQgwNDVFVVZWgl833XLTRtNfr5dOf/jSPPvEoNALvAnLXt2SGlowBvMAfgRlgHQTqAxycPsjBXx/k3gfvBT+YbCYqKytpXdPKtq3buOSSS7jyyivTtmkymUzYbDZsNltKoNDR0UFdXR3hcJjx8XF6e3sT/Cm0kbHWvzm58CMcDuNyuWhoaHhDFH4USPccI11BA8Do6Ci9vb3qRbd169acY+W7OCb3DQQC9PT04HQ6E9zF8h33QvFeWMgy4Gxk3NHRQUlJCR6Ph7Nnz+YsXsj0XaQ65T3veQ/7/7gfzMAORDrBqDs0BDwLDAGrgcuB1OwCRCHmFimYQ9OHOPTkIb7/s++DD8w2M+UV5bSsbmFb+zY1TaHNCSd/t+rq6hRJozSQlwbz/f39Kf4U8jeUqYnJyUnq6upyFn7IoEY2Q70QUSDdc4R0HRpisRgjIyP09/ezZMkStm/fjtVq5ZVXXtE1Zj5PdKnTPX78OB6Ph6amJjZu3Jh2DL0RbCbSnJmZYWJiQncedD64UKIamZOsq6tLKEpJLl5IR8bSY6G4uJhIJMIDDzzAT37xEyLFEdgATAEdwCHABtgRBFkLLAOWIohZDxTgJaALEe3+ObAky/4WoDr+akocR3EruJwuDk8f5vC+w9z/8P3gEediLbayecNmtrXPRsbRaDRtyqqoqChtSXQkElF/u+npaQYHB5mZmVFniYODg+pvWFxcrF6PUmsM8Lvf/Y6jR49y++236/yBFh4F0jUYstIqFApRVFSEyWQiGo0yNDTE0NAQ9fX1qu5V7q83T6sXXq+Xzs5OXC4Xq1evZtOmTVnJSm8Em0zO0p93enqa6urqhKm3NtorKytLuEnON85lGXC24gW/34/X61Xbzj/00EPce/+9hMwhuARYBWiHCwNOYBpBxGeB1xFRqw2xoJaNjI8CrwKlwFXxfeYKM1AVf61BkPkR4ITYFmmL0OHuoOMPHTzwfx4AL5isJsoqymha2cS29m3s3buXq666ioaGhrSHsFqtVFRUpESp4+PjTExMYLPZ1N9O608hrzW/38/09LQuv5HziQLpGgRtQcP4+DjT09M0NzfT39/P2bNnM0qxjIzcPB4PXV1dhEIhVq9ezczMTIpYPh3yTS9oybapqYn169cTDofV75Ip2pORj8lkUsnYZrOdl+h1oYsjtGT8y1/+ks9/6fO4vC7YhVj4She5FgF18ZcWIQQZOxFkPEIiGduBIIK0dwKbMow/V5xBROFW4ApgBYkPCwAFYp4YHqeHY9PHOPbHY/z40R8LMraYcFQ4WLNiDVs3b2Xv3r1ceeWVLF++PO3hFEWhtLSUpUuXJmyPRqNqvt3j8fD1r3+dl156CUVROHbsGG1tbdxyyy3zIuFvfetb/OAHP8BkMrFlyxZ++MMfJhgPzQUF0p0n0hU0SBeuiYkJVq5cyZ49e85puarL5aKrq0v1zV2yZAnRaFRoOnVAb3phZmYGr9fLkSNHaGpqYsOGDWmj10zR3ujoKNPT0xQVFSVMH2Xr9XRes+lwIUTLkL9O949//CMf/+uPMzo2CtsRqYS5XBY2oD7+0mIS2A/4ECRoBY4BryCi4iqgAREZVyGi5HyePWeB5xCLcDuBdWQmczNiAbASkUOWUCDmjeF1ejk+fZzjzx/nJ0/8BD4HWETqYf9T+2lvb1c/ksmVTl438jq79957+frXv05bWxtbt27lxIkTWa+jXBgaGuLb3/42J06coKSkhBtvvJGf//znfOQjH5nzmFAg3TkjXUGDNgIsKipi165d5zS/OT09TVdXF2azmebm5oTV5bkuuqWDdiHOarWyZ8+eOXeesNvtKY5a0hLS6/UmrHIXFRUl5EAdDodhEaoR0Eu6Tz31FJ/74ucYHhqGduBPEFGsUQgChxERqAURGV+MIGeZnpiM/3s0vo8ZEQlLMl6KyPlWI6Rp2q/lBp6Jf34z4jvM9fzNQEX8tSq+LYyaqghHwgwNDaWQrjZvng0ul4va2lo2bdrEpk2b5niSs4hEIgQCAYqKivD7/TQ2Ns57zALp5onkDg0mk0lthxMIBGhqamL16tV0dnaeE8LVetnabDbWrVuXdqU230W3dKQrq+LcbjfNzc2sW7eOl19+ec7Elymnm8lrNp0Yf2ZmRl3tn48L2UKkF/r6+vjQRz/EkSNHxJRfAV4DOhF51loE2S1nbndiGDiOiGaLgHJgd3w8CS25ET8HL4JAtWTcQSoZVwIBBOnWADeQXSucLxRmUxUWhFrjNGnTCHr9l91ud1pp21ywfPlyvvSlL7Fq1SpKSkq4+uqrufrqq+c9boF0dSK5Q4PJZMLlctHd3U00GlWn9TLizXdxLNcNLMk2EAgwODhIW1tb2oqiuSDZkczv99Pd3Y3H46G5uVntqaYoyoLqdNPpP0dHR/F4PGqZqtb4RlvEUFZWds6VFJkwMTHBxz7xMSH/WoMobCgDooCLRMLrRkSqtvgreWEs3R2qAKeBlxEkaUfkhpvJnS7QRppJ0348zOaJB4AJBPEuiZ/3w8zmjCvi59mIiKzz/ZmHgRcQ312OET+PZNLNx/Te7XYbtpA2PT3Nr371K3p6eqiqquK9730vDz74IB/4wAfmNW6BdLNAW9DQ2dlJXV0d5eXlTE9P093djdVqpaWlJeWPnE/lGGTv3BCLxRgbG6O7u5vy8nLKyspoa2ubV64qGZIMJdl6vV6am5tTVA9GqA+M8F6wWq1py1TTeS2AqCjTpijmozXWIvlB6ff7+fRnPs0vH/ulIM3kwgYLgsCSJVsRZsl4CkF2ZxC5U3v85UCQnA2hSJA/41bmnhvWQuZggwjCjQBXIghRYiZ+jlpiPo54mGjJuC7+uVpSydgFvAiMIaLnVs0+MfHSlm1DfqTrcrnSVnPOBb///e9pampSz+fd7343L7zwQoF0zwXSFTSEQiHGx8d5/fXXKSkpYePGjRkjTYvFklLemA2SpLWkK43De3t7qaysZNu2bZSUlPDKK68YRhoS4XCYkZERRkdHaWlpySgxm++U/FyqFDJ5LciKMq/XqyopZOfeEydOJJCx3W7P+xxlquPLX/4yP3jgByjVClyDIBy9SK4ak5CSMSeCaE+gEhNF8WOEEVFjulxsPnAjpvljiEWytWnGsiOi76VJ27VkPAn0x89XS8ZlCHIdjv/3FlLZJyqOmUywkUhEd/rI5XIZll5YtWoVL730En6/n5KSEp5++ml27tw573ELpKuBlH1Fo9GEgobR0VFGR0epqKhgy5YtaV2btMh3Sqv1SVAUheHh4YQCCq1EJd8oOht8Ph9dXV24XC4qKipob28/58S40BVp6SrKvF4vfX19rFy5Eq/Xm6KkSF68yzSrUBSFO++8k7u+dZcobPhTEiPD+aII4ZFwGEFIlyDkZVr9riQ5d3yfMsTCWB2CxKsQ+eNMf9YgYhGrC7FIdjn5s0ImMg4iIuIDiEo4M7Axfj7pEAWTOfVE88npSgMeI7B7925uuOEGLrroIrVjys033zzvcQukS3rTcEl+AwMD1NbWsmLFCrU80WjInOrY2BgDAwMJxuHp9p1vpCvJNhAI0NLSwtKlS3E6nedcL3surB3nAul3kK4djra5ZjolhZaQ//IDf8nk9OTs1HwKEW0akVIcR5Tt+hDyso3M3q2ZJGPaiHMC6EGQcYxZMq5HpDfK4u+/yqyhztwagqRHDBE1v4B4GNTEv0u220dJT7p60wvnQkp46623cuuttxo65puadDOR7cDAgOqeL428+/v780oZ6IUUeB85coTGxsaM7Xck8o10tXlHr9dLV1cXwWCQlpYW1TZycnLS8JTFYoUeJcXZs2fxer0EwgFBiCUIgulGaGJNiNV/O7NyrBXoW/n3IORZU4hc7bb4OHqQKeIMMBsZjyHSCHIyZo1vewlByCuY/0NjEkG2zviYjcAoImrPhihYrOnTCHoettpWUBcy3pSkK0t1JycnVcVBJBJR7Q7TectardaUltfzQSQSob+/n+HhYaxWK21tbRmNQ7SYS0cIn89HZ2cnoVCIlpYW9TvPZcxsx9J7PosRaZ20TIhiA62VYgwhyZpERJujCJnYi4i7zY4g5CXMysWKEVHqHxHVZWsQqQqj5Fkl8dckYppfiUhVVJAYGZ9CKCJARNMlzBZUrEBI0rLBjyD03vgx2pkld4XcCoeouM9cLpdqkJ4PgsHgOZmJGo03FelqCxpCoRA9PT2UlZWpBQ1au8NkWK1WfD5fXsdTFCVlrHA4TF9fH6Ojoyq5d3Z26h4zny6/sViMo0ePEg6H1cg205hzJcNwOExPTw+jo6MpxQxlZWUJN86FQrqG9UhTYqm5UhOCnMoR5CmhIKb6WjI+gqjwkpdIFKGpXYP+6FYP+hCRrALsRaQT5HmXkqjrjSHIU8raJhALeAfi5ynJuBpBxiuZVVQcQzxA2kg1Wo+gi3RtRbYEg3S73U4wGGRkZCShdVA6OJ3OC9pdTOJNQbrpChoikQhutzulpDUT8p3WW61WotGoSrqhUIje3l7Gx8dZuXJlArlnMxxPdx659vV4PGpnhJaWlowGIxJzsWwMh8P09vYyNjbGqlWr2LVrF9FoNMHMuru7W60mkpFLKBRK+zDKBxcCcUP8PPR+Da1hjDYyjiKm4q8jyOosgiD9zMrFSpmVYTWgXx42hYiePYg0SJuOz5oQEbYDEd1KyAheysUmECT7fPy7WePfK1NqIqrj2FEoKytjw4YN4pDxPm5Hjx5VG7NqWwcl978zUrlwLvGGJt10BQ3adjhFRUW6S1qtVuucZGDRaJTe3l6mpqZYvXo1a9euTSGcfDtCZNrX7XbT1dVFNBqlpaUFRVF0CcXzSS/EYjE6OzsZHR1VZwYmk0n9PauqqhIufNnmRTaN9Pl8HD58OKFxpKyfl05k2XAh5etisTSRbr6wIIi1HND2lowwS24yMpbaXa3DWB0iUtVqYv0Ish1D+CO8k/lHztoIfmX8fF5AKCzCCAlYtgdQBH2kWz4rwzSZTFgsFoqLi1m5cqW6XaYH5QN+aGiI73//+zz99NOYTCb+4R/+gU2bNvH2t789RfOrF6dOneIv/uIv1P/v7u7mq1/9Kl/4whfmNJ4WbzjSTdehwWQy4XQ66e7uBlDb4bzwwgu6b+J8tbcg/nA+n0914sp0rHxSBukI0u1209nZiaIotLS0qHlHvWSqJ70gc94+nw+bzZYQqWf7rLbNizSkbmtrS+jIoHUikxIvbYrCyEIQeb6GeTgYUfCWjpCsZFYoSDIeR+SAX4uPYUdEpNJt7E8RBGkkPIhUwzAiat+MKCHO9Tso8XPKhBgwA0pJ4vWaTqNrMpkoLi6muLhYTZndfffdPPnkk+zfv5/LLruM1157DafTOWfSXb9+PR0dHYBY7F6+fDnXX3/9nMZKxhuGdDN1aJicnKSnpwebzUZra+uccz4yXZALsqrL6XTS3NysS/uaT6SrJVLpLhaLxWhpaUmZWulNW2RLL8gFv5GREVU2J9ufJyMXmWlzulr9rNZ+UqYoZGTc19enSrbKysoStNTn0rlNDwyJdEFfFChhR1S7JXvjBhD54t8gUhARhOOY/IxUUuhdFEtGCEGuJxGph00IEg1n+5AG2dILfsTiWwBKlyUuhOXru7Bq1SquueYarrnmGp0nlhtPP/00LS0trF69OvfOOrDoSVd7E3Z0dKjuRGNjY/T29lJWVsamTZtS2q1oP29EesHr9aqmN83NzZhMJqqqqnSNbbFYUlqQZNvX6XSq3SbSka1EPpFu8n7RaFRVVyxfvly1pxwaGkr5zWSfKj251lz7WCyWtEbWUrIlvRc6OjqIRqMUFxcnRMUL6bdgWKSrJ9+ZCyUIMrUi8rfrEdGjD0HGk4h88SngIImLYksQZLycVB2t1ufBGh9XeyvFq8hyIkoq24SBQUTkXgXYYd26dQm75Ou7cC5yuj//+c953/veZ9h4i5Z002lsg8Egg4ODDA4OUl1dzbZt27IaDksi1VPBkikalcbh4XA4wfRmampKdzpCb6TrdDrp7e0lGo3S3t6eM1+rN22hTS9Eo1FVp9zY2JjiBSzJdS7T8/lM6aVkKxqNYrPZaG5uVhdatLaQ2s4VWjLWlvhecOkFI0hXIsbsXW1C5H3LSDS3iZGopDjL7KKYFUHGpQgiHEZEuTJvnAwFfaSraM5LQeSbhxC56U3xf8+Q4puQr+9Csm3ofBEKhXjssce44447DBtz0ZFuOrKNxWIMDAzg9Xrxer0Zq7mSkQ/pJt+kWuNwbR5VO7YR7dIh0Td3+fLlzMzM6Fog06tKkC2Fent7VbJN1+UCZqPiuUSTRpcBy3bqJSUlCRrn5NY4Q0NDCWbpZrOZmZkZNW0x9xPBmPSCkaSr6BjLxKzBeHPSZ6V/wjii0KMcEd1m+nPnG+k6ERK2WPzY1Yn7aNu7Q/6ka3Sk++STT3LRRRflVADlg0VHutFoVG0NI6fAsh1ObW0tK1eu1L3wkq8iARIJMJ3DmES+edp0+8pjWSwW1Td3cnJSjeZyQU9OV1EUBgcHcbvd1NfXZyRbifkQ50KVAWfqXCFLfGUvt1dffZVIJILNZkuIivX68xq6kGaUqbk2oswXZmZd0FoRDmK5bBvziXQHEX4MDQj5W/K4UVIKhPLJ3Z8L0v3Zz35maGoBFiHpSp+C3t7elHY4J0+eTPCFzQW9pBuLxZiamsLn89HX18f69etTavaTMR8Z2NTUFF1dXVit1pRj5SPvyravJNuBgQEaGhpwOBw0Nzen3VeLC6XAYS6QJb7yN2ltbU2QtKXz5/KJFqEAACAASURBVJXkndxKPRaLGRfp6olO9SJm4Fg6ZV5Zf4MIIkWhtaLMxDpRUgp4IpGI7iozo0nX7/ezb98+7r33XsPGhEVIui6Xi46ODtasWUNra2vCNDffyDXX/rFYjPHxcXp6eigpKaG4uJitW7fqXhzTWzYsSVeSbVFRERs2bEhL7HNVOkgoisLQ0BD9/f00NDSoXg9jY2NzHlMvLkTC1kra0vnzypRVsqSttLRU/A4nELKsWvLvOSZxrnK680W6xa9kZIp0Y4gUxVD8/RginZFtvCgpjVTzXUgzyksXxLrA5OSkYeNJLDrSraioUAX5yTCKdKWdY09PD+Xl5aqd44EDB3RLWPSWDcdiMdxut/rHzUS2EnONdKVrWl9fH/X19TmNdTLhfKcXjICehTStP68W0WgUv9+P2+0WpNSNUAQE4zvI8th6hLNWLaklsckwmnSNGksv6SanCdwIu0kZGzQiJG25zitN14h8JGOFirRzBLPZnPGGKSoqmld6IRaLMTIyQm9vL1VVVapxePL+ei6CXBGpTFnIyLa0tJRt27bNe1wtZK54aGiIvr4+amtr2bVr17yKDdIVUvj9fnp6erBYLJSXl6dMxSXOh5+u0ZDfsbS0VJDI9cyWvk4i7BJHEA0gA8y24pEGN3UIIl7C7N2nh9z0Yj453eRxYujyS1Aj3RlENwlP/POydNmKeDhlI92Y+IzFYlFNw2W5vt6cbjgc1t3A8nxi0ZFuNuTrBCb31xqH19TUsGPHjrR/PCPKdWOxGJOTk3R1dVFSUkJbWxulpaUcPHhQ17h6I10ZQU9MTLBs2TLdio5c0CoiAoEAXV1deL1eVq1ahaIoaafiMi9qxA1xoZQBqzldLSml6/6gIHKafQhpluyJNoPQu9YhWthUIlb2K5lfntioSDeM+G56SXcAIQMzI77XKhIj/GiOseLPUZPJpOrdAWZmZiguLqaysjJF+pfw8QtgBqUXi450s910+Ua6ZrOZyclJhoeHqa+vz0lM+aQvkklX28W3pKSEzZs3qwUbsppuLuMmQxut2+12li9fTmtrq66x9cBkMjEzM8PAwABOp1Nt7yMrATNVl42Pj+PxeNTCBknEUi2wkA0k9ep0n3jiCbZs2ZK2Ekm9yfU0glxBonkMCNLtR5DVCKJzwxkESVcj0hOy+0Mt2bs/aGFUpKvHFSxeuouP2fRBE+lNb3KlUOKk3NTUpG5SFIWXX36ZkpISXC4Xw8PDBIPBlFLxkpIS9b69UB7K2bDoSBeyt/LWQ4qyAKCvrw+bzaY7v5mP/4IkR0m2XV1dOByOBLLVfh+9yBTpyp5qPT09LFmyhB07djA9PZ23HWU2hEIh3G43LpeL1tZWXc5s2uqyaDTKkSNH2Lhxo7pANTk5qTaQzFbQsJB46qmn+MInPsHk9DRhoMhkoqqsjOVr17J9504uu+wy3vrWt+qbfmeCHSHLakUQzkngzxHRsIyMTzNrKGNhloDr4//WkOpnYGSkm+2n9yJKd+XEcgXZ5WW5It00rXpkKnHZsmUJ10EkElF12JOTk+zbt4/77rsPr9fLZz/7WTZv3syVV17J2rVrs37FbHA6nXz84x/n+PHjmEwm7r//fvbu3Tvn8bRYlKSbCblIV+sj0NjYSHt7O4ODg7oXlPItePD7/Rw4cICysjLa29sNMVhOfuBoF/2qqqq46KKL1Co8I8zJIdHGUXpYzGeVOJ1aIFNBg9VqTSBio/rDpUNHRwc3f/CD9Pf18SHgBgTnjcRidHk8dB05wqtHjvDo97/PFAhyMyJAl5KxpvhLIopIOUhjm1GEBeTJ+P4hxBS+BpEvrsE4GVumSDeEiNBdiLy0DZHDzVU7kEsWN0PG805+8Fqt1oSHeWtrK+9+97u5+eabufHGGzl+/Dhnz56dF+l+/vOf5x3veAePPPIIoVBIDQyMwBuKdDOlF7TG4StWrFB1vX6/31CJGczKzLq6upiZmWHPnj2GutlrdaKyNXtlZWVKA0vIn3STp93SWezs2bOqjePrr78+73NPh0wFDeFwWNXQjoyM4HK51JtAm6JIt3CXCcnfc2BggI998IMcOXyY/wHcTaIfzPL46zLNthngKiMJLh0hyei2BtG6RyLEbAnvILPRsfS1/Q9mbR8bmF28K8/jfJMjXQVB/GcRqY4tCMJ9HX0PnkykG0HIysYzvK8TLpeLuro63vrWt4pZyDzgdrt59tlneeCBB4DZEnSjsChJV296QWscnq4rxFw9ctNBS4IVFRVs3bqVjo4Ow9uHSNvKl156iYqKihSFhRZzae0jK/0GBgYYHBxMeEhp95sL5vLZoqIiqqur1ch6ampKLYqRKYqxsTH8fn/etpBut5tPfvzjPP3b3/JW4CfkDtgk1CVBo7wX8iFvG7NOY1s026eAnwKfRKgoBhC54mOIp4TCbHPKBgQR15K+IaUk3Riig0Qf4ru2Itr8SOhpw5NuvxjioTGAmE6sgDJv4gM3n2vF7XbrKo3Xg+7uburq6vjoRz/K0aNH2bFjB/fcc09G06x8sShJNxOknCkYDOY0Doe56XqT1RHa6X0uEsyFbAs82txwJBLh4osvznmcfDW90WhUrVJbtmwZe/bsSZHHzbc4Yr6QjmbSc0Hrl5rJFtJmsyWoKEKhEF/72tf49SOPsDEW43skNnPQjQvNe8EWH6cY0T14Y9L704g87BCivY6UtBUx6zQm9cUBxPc7Gd+nkdSGl/mcuzbSlfngMKKwpA6YIGWmlo/Hh9PpNIx0I5EIhw8f5jvf+Q67d+/m85//PHfeeSe33XabIeMvStLNdPMGAgG1s+6aNWuyGodD/r3BtAtpWrLNNL3PBzKKTiY5rerB4XDQ3t5OR0eHLmLXS5CKohAKhThw4ABLly7N6r9woRQ4pEM2W0gZFX/1q1/lwXvvJYCoSH0Xqda0ujGfhTQtci0y6UWust3q+Gu7ZptMG0hJWy+zkjYLIpWwIcv5KeT2jYhrcFEQkrlpBMmv1oyrQFlpqlfG+SiMWLFiBStWrGD37t0A3HDDDdx5552GjA2LlHST4fP56Onpwev1YrVadbfgyRcyMh4ZGaGnp4fq6up5k61EsjmNVs9bWlqqVsXlg1ykq1U8KIrCjh07UnKqychEuheyVMdms7Fv3z7+6QtfIObz8SnEOtBJ4LvArYgUaCWCgNsQ/RuzmWupv6pRke75sog0M5u01uIYoh1PrimAHt8IOSN4FRGFbya1fVAUKpckRqr5kK7T6TSMdJcuXcrKlSs5deoU69ev5+mnn6atrc2QsWGRkq68wZONwzdt2sSBAwdQFMXwrgKxWAyn08nw8DBAgkog1+f0ejXIfLEk2+Li4owSMz1Tr2zyMpl/rqqqYseOHbz22mu6LnCjFBELheeee45bPvpRJsbGuBm4ltSLfgYR4HUjpLKHgF8gZr/liOCwGdiG6Fxeg6ZhwoWUXsinA0Uu6DXh0bOfzMglWzlqEYGqykTSzNd3YT5qhWR85zvf4aabbiIUCtHc3MwPf/hDw8ZelKTr9Xo5efIkkUgkwTgcZqPRfEg3Vy5VFhuUlZVRXV3Nxo3JybL0kASl51zMZjNTU1MMDw9js9loa2vLGHXq9bRNJkhtXri8vDwhSs/He/d8pxf0HP/UqVN8/KabOHX6NO8H/oLUpggSdkRUux7QNnmZQhBxJyIq/ilC2WBHs/b0R8Q0WXZtmAv0LkblgpEeDnprjLLpghVEymIEUaWWTWUYJSUtdD5tHbdt28bLL79s2HhaLErSBVizZk1araiUjektOc1kZK4oCiMjI/T19anFBrFYjNdee033OcroNdeFMz09zdTUFKFQiI0bN+ac4mfK/yZDS7qTk5N0dnZSWlqaVjOsN4JNR7ryt4pEIqr3wrnqX5Zr1jA+Ps7HPvQhXnz+ea4BbiP7vZ4N0lp2p2bbJPAlxKL7nwLOw9B9WEhVS8wQLYWA9Ittjg+QC1Hm360XjI109ZBujMykK83KSxAphVytCaOkkOaF0KrnXGBRkm55eXnGqf1cncYk6WrJNtmHIRwOGyYxA5GH6uzsxGKxsGTJElauXJmTcEE/Qcrea4cOHaKoqIhNmzZljZ71RJDaYydXwVmtVoaGhvD5fCiKolaYGem9kAl+v59PfuITPPX44+wCfkRqmnI+CAH/DjyDIOGvIXhVwgd0K9DthZNeeL0LBp4T2QerFYJLICqbQjaRSLL5SsYywchIN6RjLGnrqGWRICJXE0LMAKoQZui5mCYirqfBwUFV8pfvQppR6oVzjUVJutlgtVrn5DSmtT6sra1N68OQT0UaZCZdl8tFZ2cnZrNZ7Qhx+vTpeZvpaOF2uzlz5gyBQID29vacpuv5pBcURWF8fJzOzk4qKyvZsWOHquyQKQ/ZYj25wiwQCHD69OkEc3AjouJbbrmFXz3+OCWIdMC/AesQBLmL3O6KmRADfgz8HEGy/wvR0isZDoRkdgtCESE/OwZ0RqBrDI6PQ9erIlq2myFWAoEahIQq97M2NyIYd0eH0NdWnfh+UYQUbQKxGtlAgjJBD4GvWLECk8nE6OgoPp8Pv99PUVERwWAwZ1cPl8tlqJfuucSiJN1cpjf5RqPDw8NMTExQV1eX1fow3xX65Khbkq3JZEppB6+3XTpkj3S9Xi9nzpwhGo3S2trKa6+9lpNwc42pRSAQYHR0lOrqarZu3aqmKZJ/c22LdW1/qYMHD1JXV4fX680ZFefze5eWlnIt8GVEkdRrCHfFbyN4oBwRdC0D2hELYk3ph1Lxh/jnzcA/AG8lv4DUhOCeBuAtoDpphYA+BZ7zwY/i1hjKJPAss74K1eQftRppEamHdGWpsAeRt61EKBOSb58ouWVlUUG6y5fPzk+6u7spLi7GbrendPUoLS3F4XBQVFSE2+3G6XQWSPd8QW96QXZQGB0dpaamZt4+s+kgI1K3201nZyexWIy1a9emnQbp7dyrHVcLn89HV1cXwWBwTt4IuUhXRs4zMzM0NDSwfv36vMaXMJlMCRVmkDkqlr4L2qg4GxGbEFr7lcDbNNsDCO+YE4iZ7nOI9AMInliCUEbtBC5G5GtvQyykfQIRuRp5o8wATwC/RqRAzgLKNgSJvYrQsQYQT4paRIgtq8ccZGZ+IyNdabKTDtpqMuk0tpbE2mkt9Gh5o6Q0f4xGo5SUlFBdXZ3S1cPv96vX/O23305PTw9XXnklbW1tXHvttVx//fU5v2ImrFmzhvLyciwWC1ar1fAFtTcc6RYVFWVt3KjtDVZfX6/mUY0mXBDR35kzZ7BaraxduzZron+ubXikp63P52Pt2rUJSo58kEmV4PP5OHPmDJFIhNbWVnw+X1rP4vnodDNFxdJ3QUbFTqdTfU+qSfTkiksQhRBbNdtiCLI7gVAmHAXuBf4VkYqwINITHoSCYS3zFxhEgP8L/AChpHgXgkf/G0T+U5skDjPblXcYEb674u9VIcLnpcy6jRVhbKSbiXS11WS18XPcRPYpQK4WQkExnra6EDIvpJlMJvV6qa+v57HHHuOyyy7j97//PSdPnjRE0rh///6UJplGYVGSbrYbPFOkm9yIUdo59vf3590ROJf21uPx0NXVhdvtprGxUZd+0GKx6M5FyxbiJ06cwOVy0dLSQl1d3byJT3uxJpO5jDQCgcCCScak70IwGORLn/88L770ElagvKSEZU1NrGtrY+vWrWzfvp2JiYm8ZuMmZu0LrtRs/y/gIeAWBBG/jMjnRhAL8NUIAr4IUUCRa1EeBOe8CHwTEcBeGh9DQoHUSLAIQapLmfVXiAF+RJR5FpG8PkhcOoFYnHMhQvSlJFpA5nunh0lME4QRka22mswT//9chJuJdKOo7mk7d+5M6ZySr/TT4XCwc+fO3DueZyxK0oXMkVlyTlf6CQwODrJ06dIU79x8u01kk2t5vV46OzsJh8OsXbsWt9ut+6KxWCwEg8Gc+4VCIaamphgZGWH9+vVs3LjRkGowSbqhUIiuri7VoDyZzPUuuBkBv9/PX3/qUzz56KPsRRjSWIDOQIDOEyd47cQJ/vjII0whotB3GHDMGCKQvCH+ktvGEcHmCeAIIj1xJyJirUBw3CZgd/xfGRV3At9ARMubECmM5Ig5BvruRBMiveBAkJ6EtIB8AkGIaxBk1ktqlwpty6AKMhNmOP7lFISl5DBiCiDdxUBf+3WpzEg2u3GK86uvq+fh/Q+zffv21I/q7I8WDofn1O8vE0wmE1dffTUmk4lPfvKT3HzzzYaNDYuYdDNBqhekU9bQ0FBWPwG9DSS1+ydPe3w+H52dnczMzNDS0qJGhT6fLy/T82xkJj1tx8fHcTgcrFy5MqWJ33wgi0C6urpoamrKaFC+EMURiqLwj//4j/zoe99jraLwXcRUX0JdmIojiEgLGIF0FbkmRNBYT6K9YwhBpicRqdgORPogiOCzYoR6YSVwE6nrSxJpI918IC0gixFkvCfp/RCCgAcQ7mPH4yep7VIhUxW1iIg5ggjLX43v14J4GmmhR+qWTMxxSZklaOG2r97GLbfckvGjeiVjLpcrpbBiPnj++edpbGxkbGyMt73tbWzYsIHLLrss9wd1YtGSbjYPALfbzYsvvkhjY2NW8xaYn72jTOQHAoG0+dR8ouhMC2nJnrZ79uxhcHBwXj65WkSjUfr7+xkYGKCmpibF/jIZ55p0v/vd73LXv/wL5aEQ/0ZiYUImFJPKB3OF3upXECS6If7SLttMIqLiDwL/g9x2kTGYH+lKZMrp2hBPrXVJ26cRT41hxBMjgIiK7Qhy9CJOfjnpE9p6LkFJzNH4ccbgqquu4sEfP5jTSyQWi+lyGTO6Gq2xUSTX6+vruf766zl48GCBdNMhEomokW00GuXSSy/V3So932IK6fng9/vVyDYdqc2nkaUkw6GhIVauXJlAhvkoHbQ+uVpI9UZ/fz+NjY00NTVhsVjm7OcwXzz66KP8/S23MON28zlEtdfCdU2bhRF1CjXMRuK55Poy5WkI6ebzxAAR5e6Iv7RjDCJC9mTv3HTH09O4EuAYNNQ18Iv9v0ibSpgPjCyMkBLG8vJyfD4fTz31FP/8z/9syNgSi5Z0JYloyXb58uXs3r2bQ4cO6a5kyYd0A4EATqcTp9PJ+vXrqa2tzdkjLF/S1S74NTY2snfv3pS8sNlszmvRTevToPWSqKurU3PceqNnoyPdZ555hg+9//24vV5uAj5C5mn4QsAowy85Vq6rUP3FjfJemO8dbUZ08jWRu6Ik148VRLT2AW7/6u1ZUwnzgZEOY6Ojo6rcLBKJ8P73v593vMOI1YJZLFrSjUQi9PT0MDw8zIoVK9KSkx7oId1AIEB3dzcej4eysjKWLVuWIm9Jh3xlYF6vlxdffDGnp+1c5GXaNkLSWUwrtZqP4Y20ofR4PFRUVOBwOHI+9Lq7u/mrm27itddeYweimOlnwC+ZbW7QhkhPtrNwUa9RlbTaYq1c+5nRN1PXdVAjTj6GPgLP9GNpUgm7L97NXV+/i9raWqampnTLM/N5sBtZjdbc3MzRo0cNGSsTFi3pjoyMYDab50y2EtlINxgM0t3drcqy2tra6OnpySt61dNTbWRkhO7ubhRFYe/evTlXYvPtCDE1NUVfXx8OhyNjZwu9Yybv53Q6OXPmDDabjbKyMs6ePYvX61WF7VJLW15ejt1ux+128+7rruO5Z57hKuB/ItZuQNyrg4gV/9OIqrLHEAFTZfwl5VpvIffUfS4wyvArHB9Hz+zbMCfiXHpYvUinOMi0n/bW07T2aahv4OFnHmbLli34fD58Pl9CN4+ioqKEwpfk8t58zW6MXEg711i0pLtq1aqMhKbXbxbSG70Eg0F6enpwOp00NzcnyLLy8V/Itq/sPNHd3c2SJUu46KKLOH78uO5W8HrOweVy4XK5iMViWc1uQL/hjYx0Zbmxoihs2LABh8NBKBRKSGPIKjOPx0Nvby9fv/NOntu3j22IgoDVSWNb4ttWk6iddSKIuBNBxD8C7kJIUysRNQXTiPWe+ZKmkX7ieluHGQajIt1MnYCToSVdqUqYsfC5z36OW2+9Vd0tVzeP5PLesrIy9T7Q40ftcrlYtWqV3m933rFoSTcbZPSab5XZzMwM3d3dTE9P09zcnFY2lY8iIR05Sk9baRYjzdAVRZlTRVo6aP0XKioq2LBhg66OEHoi3VAoxPT0NH6/n9bWVpYsEd6FyZ81mUyUlpZSXFzMfffdx39985s0RqPcjSjPzwdVCBWDVskQRjgHngbuQ7h8DSCUBFI3uxlRiHAF+v1kjEovJDfTzXY8E4iih/n2MDUq0tV78rK8dwAYg6vfdjXfuvtbuu4Pm83GkiVL1OsHEsvB5TV26NAhLBZLQrPRsrKyhCjY6XTS3t6e77c8b1i0pKunKk0v6SqKwuuvv87U1FRWjSrMT5Gg9bTVmsVAfv3aMp1DIBCgs7MzgRBfffXVOaUNkhEOh+nu7mZ8fBybzcbFF1+cMwJ54IEH+Orf/R1FwSBfQZjMGDWVLkJUjP0QEe3+M6L11yQiIj6DiIpvBT6HIN1qhNz0YuAqEqvCJIz0E881Tj/CAx0TojZY6m2XIZLadeRnfGNUpCtzI9kQi+/ngoZGkUrYtm0bY2NjeTnxaaEtB7fb7aoLXyQSURuOjo6O0tXVRTQapbi4mIceeojBwUHa29vzSklkQzQaZefOnSxfvpwnnnhi3uMlY9GSbjZII/NcCIVC9PT04Pf7KS8vz9nIEvJTO0ginZ6eprOzE5vNlrb9Tr5IJkgZoaerIptrrlYiGo3S19fHyMgIq1evZsWKFZw+fTrn77R9/XqGRkaoAG5EEJxRhPsSosprBkGoVzHLEVLfr60PmEHUBMgOEA8jOkCYEJwmO5m/ldnOMvNFtmBxGmEoNo7wOe+OwaeiMBaFvhEYGYEpm5ixE0EYydTHT7QWQcbpouKFinSDiB90Br74xS/yr//6r+pb+XR7yAYtgVqtViorKxNkYTJ9tXXrVjo6OvjZz37GPffcw+rVq/nlL385r2Pfc889bNy4EbfbPa9xMuENSbq5iDEUCtHb28vExASrV6+msrKS+vr6vHuZ5YLL5cLv99Pb28uGDRt0WSzqgaxe01apZYrQ59oRQusv3NjYyJ49e9RSZT0RuW96mhsQJPZb4D8RmvsqBGdsQkS+2pLZXOhBRK5DwIcQZK7HFt3ObBHDtfFt0vBGRsUnEMoJaW37eURPtLb45/JdtEuXFp1BWCWcROSfPwS4EcViZmatFgBRRYaoV+h1wasu8AzZMBWZCXqD4s6tYdZjoQ7jSDdCetLVqBLefvXb+dEDP0opcDCSdLONI9NXN954Iw8//DD33XcfK1eunHOULTE4OMivf/1rvvKVr3D33XfPa6xMWLSkOxdPXUlSY2NjrF69mj179mA2mxkdHU3bsicd9ES6Ho+Hzs5OFEXBbrcbLgaXC1kHDx5Uq9QyLRrmG+lqF/hqa2tTvCr05n7NZjNXA5fH/z+KIBdtyezjCFKpQARwzYi87V4Sc7BOBNm+ivBX+F/MvwJNa3izB/g/iAa40jkxipCwTcaPXwFsRKQwNiHIeBWZHxhaVYKCIPWX4t/rPQi+lO9le9SXxI/rBSztF/HUU0+hKAqHDh1i3759HDhwgI7DHXjcHpSIAo8iiDhXVJwNyU8MjSphacNSHnn2kYw51Gg0aogPgl7fBRDqBSkZmy/hf+ELX+Cuu+7C4/HMa5xsWLSkmw3J3SPC4TB9fX2Mjo6yevXqlFLXfFIG2faVHgyhUIi1a9dSXV3NCy+8ML8vo4EsnOjv7ycWi6nRZzbkQ7qBQICDBw9SVlaWsdtxtuII7YPQYjajjTksiHxqC7PRJghSO4nIv3YA9wO3I8ipCjHTHUOkJx5AdLsxCjGEt+7dCJK8AjFz7gb+Q7NfGOES9kcEcf5fBBFHEEbo0jZyIyIqLmV2hj6IaPETjo+fnEfW4xlD/FilcV212Wxm9+7d7N69G4BDhw6xa9cuxsbG2LdvH8899xwvd7xM/9H+zFFxtlyxNjMXAHrBGrJy++2386lPfSrreRoZ6erptg3CGCmdDDJfPPHEE9TX17Njxw6eeeaZeY+XCYuWdHMtpMl+Zlrfgky+AvmQbrr0gt/vp6urC7/fn2CDKKG3DXsmqZvU8vb09FBfX69W3em5uPV0pHC73Zw8eZJAIMCuXbuy5pxlCfI//dM/8b1vfpMYUFNVRfOWLVy8Zw/XXHMN27dvx2QyoWeiV4NQGFyq2RZkdsr/VYQkrBv4GCLirEMQ3G6EZncuF/EpRM+zQYTKQVbCniA1ei1CGN0kV9/3AfsQKYMfIh4gboTR1woEqT+JiIr3phkX0PUbyf1ykVB9fT033XQTN910k7otOSo+fuw401PTKCFF/JhasxsZFUvSjasSMqUS0p7nOcjpZoMMAPTIQ3Ph+eef57HHHuM3v/kNwWAQt9vNBz7wAR588MF5j63FoiVdyBx1mc1mxsbG0voWpEO+pCv31RZPrF27Nm1ZsN7OvXLf5JJdWUVWXV09p+4W2SJdv9/PmTNnCIVCrFmzhuHh4ZyLfD/4wQ+49e/+jtJQiH9GEOKZ6WlOPfssTz77LN+96y7CiGn6XLNr0kFwC8I+8bsIYhxgNj1xGBERu5n1uW1CkPClZE4/jAPfAV5A6IFvIvEmiKF/wW818PH4S8KPiJ5/FT/PD5K9mlavWiIMaaftufLryVGxxMTEBL/73e949tlneeXoK/R19DHjmxE/hg2YgTrquOu/72LlypWcPHkyQbKVLNuSMFJBkA95G2Fvescdd3DHHXcAojz9G9/4huGEC4ucdJMRiURUx6zi4mJd02/IX5GglZglF0+kGzsf0pX7SnlZtioyveebHJnPzMzQ1dWFy+WitbWV2tpagsEgg4ODGcd54okn+MInPoF3K+tXtwAAIABJREFUeppPA1czOztdjVAQSEwCf83cSVcLk+bfVfHX2zXvuxBEfBJBxA8hzMJLEYtfSxFlxLsQRPsQIhL9C9Jrd+frKlEaP79lwPPkti/IlV6IINIup4ErVq5MeV/vLCoZtbW1CVFxf38/VquVgYEBHn30UZYsWcLf/u3fzp5HFtmWlojzNR7PBL3kbdTxFhJvCNKVjlzDw8MsX76crVu30t/fr/uPoZd0w+GwKjGrqKjQJTGTkbGe9uNmsxmn00l/fz9FRUWGycsk6Uq/Cql20D4sMumEX3nlFT72vvfR39/PR4B3k1sxUIMIlvLrx5EZ2YiwErEQtgf4aHxbGOhCpAqOIbo//AhBgO9EEGK2Yxmh080kAEhGpjLgGCJ98UegrLqaR3/8Yy6//PKU/fRWXuY8j3jZdrqoGDLLtoLBoFpZNjo6ytTUFMFgkPLy8pxRcTboJd1zVQJ8xRVXcMUVVxg+Lixy0lUUhb6+PoaGhhJkTYFAIG+7xmxVNNrc8OrVq3E4HKrnZi7olZjJpoyhUIiNGzcadiHJ1j69vb0MDg5mVDskqxL6+vr48Pvex5HDh7kesdiUj+BNKgDmCxP5l8oWMSsRe3d823vj55ONcEGQnVGVtHqLupL3kzreSbOZv/3KVxIizpTPG0S6iqLkHTGaTCZKSkooKSlRDaAOHz7Mpk2bVDLOFhWXlJRkDFryMTA30kt3IbCoSVc2fUx25Mq3DXum7hHaCFqbG5bqASN0vX6/n87OToLBIBUVFTQ3NxtGuLFYDKfTqT4s9uzZk/FClmmT6elprrv2Wo6+8gpXIFrk5DLhTjsexpAuzH/Kn88Y+eR0s0Hv1aetXAshFuVOAFddfTX3P/BAzvJtIyNdo8ax2WzY7facUXEgEMBsNicQsYyK9eZ0jfTSXSgsatLduHFjWkLT4+6lRXJ6QVEUBgYGGBwcTIigtePnm6dNhjavKhfhTp06lZe4OxPxa/0d7HY7jY2NtLS0ZB1LURTu/uY32ffEEzhiMRSETOokYrFsK6Jia43OczMy0jWKdPWQqZHpBT2QUfxJRA546bJlPH7//Sn66IyfN5B0jcqNZmrzlBwVQ+ZcsTSdyhUVF0j3AkG+CwuSdLXdFLJ52s6HdGVeeGJiImURTo+8S0LmYJO/q9Pp5PTp0xQXF7Nt2zbVUi8TFEXh3//937n7ttuojUS4C1GJFUEUM3QiWs88i8iLmpmtKpN+t9tJ31fsQop0YeFJV884QYQC45DNxh133MGf/dmf4fV6efXVV1WtarYpuZHpBSPGyReZcsUHDx6kvLw8a1Rst9sL6YWFhhEyERBk5/F4ePHFF6mvr88ZZUiS1rM4JveVygrpYZAur5pPGx65rxxDOovFYjE2btyolhwHAoGMRP7Tn/6UL99yC/j9/C1Chyp/UStCyL+W2S67snT2DGI1/VXgN8xWlS1BdHi5OL6/UaRrhP3h+Yh0sx3Pj4hse4Bdu3fzmyefTHmIZ5qSS9etsrIyw+4BIyPd+cJkMmE2m6mrq8saFd91111qMY+iKLS3t3Pdddfpui+TEQwGueyyy5iZmSESiXDDDTck2FMaiUVNurmQK+8qS167uroIh8NccsklunSw+aQvTCYT4+Pj9PT0sGLFiqwytrk4mEWjUTo7O/H5fKxbty7FQT+dTnf//v389Yc/zOTYGDcjKsT0dgCXpbPaQgHpd3sa0Wj2TgRx3YWoMFuP8Fm4mtwLWemOaQTyIV0jjpnJbimKeFgdBFpbW3nuxz9m06ZNaffNNiWXRDwxMYHH4+HgwYOqF600jbfZbLpJeS4Lackwqo2Toihpzzs5Kv7v//5v7rnnHoqLi1m3bh1Hjx7lXe9615yOabfb+cMf/kBZWRnhcJhLL72Ua665hj17klsrzx9vWNLNlgLQFh1IT9uOjg7dhQd6jMxjsRjDw8P09PTgcDiyLmIln7NedHZ24na7Wbt2bYKzmBZa0j1x4gQf+cu/5PTp07wf+EtEbf98If1u7Qg/BStwM6L5bCdiYeh+RHWZHVHIsBqhnX0bIk2RKbo0MqebC6OIVIoP8SDagchltyHKl/NxFEgnBZNWjrbycn7+wx/S1NREbW1t6odzwGq1UlVVRVVVFaWlpUxNTdHS0oLf71dVMENDQ8zMzCR0aCgvL6e0tDRtGsGIhTSjUhT5FFh4PB42b97Mddddx3XXXTfnY5pMJnXRMhwOEw6HDZtFJGNRk24u05twOJzwx5O9vLq6ulKKDvJ5SmfT9cZiMcbGxujq6mLJkiW0trYSCARyXkS/+93vOHz4MO94xzuyuuBLq8Xp6WmamprYtGlT1t9Bku4lu3dz9OhRLgV+gSA+ozAC/CtCG/teRBWWLBjdBMjYI4qoKjuDILd9wL2I9EE1YsGuHeFR8BZmG1Se6/SCFzHV70MUTFQgHkZ/AH6NyLsGgJXx85OmNxvJLKPTXh0uRE58zGzmc3/zN3zlK1/BbDZz6tQpw4jOZDKpXrQNDbN6E22Hhr6+Pvx+P0BCdVl5ebkhhLnQJcAgdLpGLaRFo1F27NhBZ2cnn/nMZ9LqlY3AoibdbEgmxunpac6cOUNxcfG8iw4yRaRSMVBeXq4axkxMTOD1ejOOdfz4cW666WOcOvU6sIRbb/03LBYH9fX17Ny5kcsvv5x3vetdLF++PKFlekNDQ85uxDBLuseOHmUtwrDlfQiSq2e2u8J68s9l+oHbgEMIN7GvMdvvLB0sCPXDGkSEC4IMx5nNE59AlNB6EOQXRUSHJfFz1GeBkop0pBsBjsRftcAHEMY7gyR2qQDxXXvj77+EIGIfIsrfhCg/3oyIipfFx44hFCCvApddfjkv/uQnCXLAuVaTaZGLLDN1aJC50YmJCXp7e/H5fBw7dkwl4Vw62nQwinTzGcfIhTSLxUJHRwdOp5Prr7+e48ePs3lzvn1OcmNRk66eSNflcnHmzBksFkvCAtN8kEzo2uaMW7ZsSSD0TKmIoaEh3ve+D3Pw4CFEndRtCCoMEY32MDJymscfP8Hjj3+XL33py5hMdiorq9m+fR1XXHE5bW1trFiR23NLbcGOKIGtRER0ryE6LxxG+MhGEARSg4jg9iCm1+kuEAXhh/AEYuHsPoQt41xgQpB/PSK6BUFWjwPfir//MMJm0YPQDLcjSG5T/KXnltNGyzFE2uM5RDR9LSLKlvulu6pKEYTaFv//A4Bt717+5m/+hn379vH/DhzgoTNnmPL5MCMeGG5gculSnn7oobT2nkZEl3MZw2w2U15eTnl5OcuWiSz7wYMHWb9+vdrTTrto53A4VCLO1un5fES6RnYClqiqquKKK67gt7/9bYF080E0GuXUqVPYbDbWrVunq+Agn4KHcDiMx+PhzJkzABlNypMX3ZxOJx/5yMf47W+fRsSY/5vE5SUbIqZbD/x5fFuEWGwAp/MM+/efYP/+nyHWvc1UVFTT1raGSy+9hHe+850pqghJuhaE+N6MMIZpItHQexRBxK8iItc7EKRRhVAlNCNysNPAjxHT6q+RGhHOF68h3L/OIlQQXcBfAZ9CkO4fEGT5KILsnYiUQBviIbEl/t/LSSVPE8Im8pn4WDsROVst9C6kKYiH/uWXX87b3va2hOvm+PHjPPLII1RUVPDFL34x4xgLEenqhclkori4mOLi4oQ8s1QMeDweRkZG8Hq9KIqidnqWZGy32w0zuzkfpDs+Pk5RURFVVVUEAgF+//vf8/d///fzHjcd3nCk6/V61QWmZcuW0draqutz+RjTRKNRhoaGGBsbo7W1NesfXaYigsEgn/3s53jwwYeJxdoR2Uy98aGVWaq8Or4tBozgdp/mpZdO8dJL/5tvfOM7QIzS0mrWrl3BJZfs4qqrrqKyshKLycRMJh9cZrsWJHfhfS3+egXhzuVCOIjZEGmAAYSSIdHMMn+cBe5BEH4zs+5fXZp9yhH5Ye36dAQx3d+PMLR5BPFgABGFX4QgVi9iIesQYlHs3WSO4vWqHGKxGCdPniQUCqkt6MvLy2lqauJf/uVfchJqplX6fHCu9bXZWuV4PJ6ERTsQ5D0yMqJGxXM5t3xI1+v15qza04ORkRE+/OEPE41GURSFG2+8kWuvvTb3B+eARU262gs22dM2EAjkpQSQHrzZ/tjSynFycpKysjK2bduW86YxmUx85zvf4aGHHicaXQ5z6oebdmTEpDiAsPd2IZaw3o7f38uxY69z7NgRvve9hwAPZcRkBxjdqEJM+eW0/xEEMX4bQcSHERHntxE51yryr17zxc/+/yLyqu8jsdFBriILK6l+vDC7UHcQER2Px/d9P7k7A+s1qiktLWXbtm2AqDCUU/Px8XH8fr8q5teav2in37FYzJD0wnyjy3ylXrJVTmlpacKi3fDwMC6Xi3A4zMDAgFpan2wJmUslpDdNYaSXbnt7O0eOHJn3OHqwqEkXBBF2dXXh8XhoaWlRF5fOnj2ru1U6ZFckyAaWk5OTNDc3s2zZMoaHh3MS7pe+9CW++93vE4spiBaK78I45ekEoonNSeDPEHQoo5GlJLZm9GDmzwjNU3wlb5Xt8dcH4v8vXb1k94cXECkIc/yM5ILdW5jtiRZB5ITvjY/75/H9khFjbuoFaXgjcQ0iY56LcHOlFxSEFvl1YJNGO2u327Hb7QkG9pmm5qWlpZSXlzMzM0M4HNbdISHt+ZynvHA6SNnVSo0FpXbRbnJykr6+PkKhUEqlXWlpqXo/RSIRXYbp2uMuJix60u3s7KSuro62traEH19vR2CJdKSb7C7W2tqK2WzG6/VmLY746U9/yuc+9494vWFEOuA0YnL+X4hbX7q8Xkb6RuDZ4EdkXF9E9CNIzgmnQzkmzITmWSNWTHpPAa2r13vi22KI1MNxhL3iK8BjiLxyZXycCKL7w/osxzRKp6sX2Uh3CJHGMJeW8j+//GU+85nPZB0r3dRcURR1ah4KhThz5oxa3TgX5cCFJPVKN4520U4iFoupMwOv15syM5Alv7nSDEakZ84HFj3ptre3py1zzceYPHl/RVHo7+9ncHCQFStWpHSeyCQZ27dvH3/1V59lfHwS+ARiqUpK6hXEbXsaEScdRKzJQ6KA660I/UDyjRRBaAZ+gyDq/0JkLfVi/qRrR7+Ri9Z0/J2a7WMINcL3EVP9XHRxPgxvkvfzIBbvhk0mPnnLLdx22210dHTM6YY3m82qnnZ4eJj29nYsFguhUAiPx5Ox3Le8vByHw5FCahdSpKs3Qs20aBeNRtU1GZfLxfj4uOrzq42Ki4uLMZlMeDwewzpsLyQWPelmwlzsHSORCIODg/T19bF06dKMVWTJhP7KK69w002foLe3BzHpfi+ptV5mhLx+JbPLVdLNQBLxUUSWNMJsA6tNiD/TE4is6dw0AybM6E+2pEc+pJsJ9cxWoOm5zefip5sOesfQkm4EkbfuAHbv3s0ffvYzlSSMVB6YTCY1PZGsHJDR4NDQED6fj1gspqYnZMnqfKPU86GvTQeLxUJl5f9v783DoyzP9v/PTDJZyE4gJCSBkIQs7BCCWhW7fN2wtVXbitWioj/7q6J4WK1a+6rUvmCxbiitb0uLb+sX6eKCdV9aSlUIAQQFsi+EhOwkk8kyme35/nHnfvJMMjOZDQJxzuOYgyPD5Jl7njxzPtd9Xed1XglERkYya9YsYmJi1KKdzJefOHGCwcFBdu/ezWeffYbFYuHAgQPMmTPH7zTN8ePHWbVqFS0tLej1em677TbWrl3r9+cYC2c96Xoak+NtekFRFPr7+2loaCAtLW3MWWRSBlZTU8N1193E559/DlyFEDv54oWrdTOQUwEUxMCbCkRU+1dEtOxA1O/fG/r/fASBe3uR630upI1ENMExsfH1a3m6DW90iBz1LiBpyhTe3rZtVA9+MEh3rGNo230lHA4H/f39mEwmOjs76ezspLu724mI4+Li1GjQG5xp9pBaFZG2aJeSMpz1z8nJISYmhr/85S8899xzHDlyhGeffZbzzjvP5/cLDw/nySefZMmSJZhMJoqKirj44ouZM2fO2L/sB8560nUHb9ILsi24urqasLAw0tPTvZKYdXR0cM89P+XgwS8QUevLCLPDYECHiLHeQ6QgYhAK1DxEqaoU+AhR4zcjUhOFCBVtIUIz4MJvgvAzItIF3zrfgpnT9XbU+TGgOTychx97zG3eNhikC74XgbT2hiCKvDNmzMBgMLhMT2iJ2J2E60yJdCW8kYwlJyeTl5dHcXExzz33XEDvl5aWpjaJxMXFUVhYSFNTU4h0fYW7mV8S2rbgBQsWqEYhntDb28uKFVewd+/BoWduwbupYd7CiBjm/Q9E2SoPuBPnOjyItoVKhn29jiCyjpOAQUREXIBIQ8wBslEIO2MiXT2+EWmwIl1PMCP0vtXA4qIi3n7nnYBUBV6tKQiuXIqiEBYWpuZItW5kVqtVTU+4knBJQj7TSNdb851T4aVbX1/PZ599dsp8F2ACkK6vkYK2i0zbFmw2m91GxjabjXvuuYctW/4vDkc2okfqKKIk9DuG3WRzEX1U5+OsNh0Lg4g0wouIVMI0YO3QsVx9vngEoWpzu30MGyweQUjJ9gytw0IfDh4lincxcz7CI9fX2DyK059eONUuYw7EX3IPkJ6RwaubN6utn2NFs2dC5dyTJaPBYCApKcmpeUcr4ZKWo2azWc0tSyL2JT0hcaqnT4xEMM1uQARV11xzDc8888wpGXYpcdaTrrfo7++nqqoKi8XC7NmzR90hXXkkOBwOHn/8cdavfxabLRlRxCoaeWSG7VoOI5Sn/80wEecgtv4XMFolakfkbTcj6GUSsAZhB+Nrji0G0ZawEFHIA9Gb9XvgQ+ysoIvZvMNh3uEIP+cEeiKZjMIi+jkPMTp89JDvYUQRnKjT20i3C9EV9w/ErWg+Qt+Riu9qZ1fm5CcQLcFKdDQvbN6sTm3QFmwMBoPTFl2rJw0UwTiOr5aMrnwXmpub6e/vJy4uTtUUm81mwsPDnRo7xuowO91G6N3d3U6RfSCwWq1cc801XH/99Vx99dVj/0IAOOtJdyxbw76+Purr6zGZTOosMlcYWXj7/e9/z09/uo6BgTDgpwgplzsrlJFkN4CIOisQKtWtCGvvOAQRZyHq+B8hItQw4EbgOwQnVWEH3kEQ7hxEBC0tXaTnqA0H9XRQyYcc4UOO8BgN6DCQCMyjn3MREbFUEgcr0h2LIsyIbHYZ4owPIPYBLzF8tvIRcf4CBBHP9OK4EiaEleNxnY5bfvQjHnvsMZVMEhMTSUpKQqfTodPpsNlsaq60o6OD/v5+wsLCMJvNNDU1uZVynS4EowimKAqRkZGjJjXI9ITJZHJKT2jlW3FxcWr+NRik60uuXHpJBwpFUbjlllsoLCz06JURLJz1pOsOFouFwcFBDh48SG5u7qjmiZGQhbcdO3bwox/9hO5uE/BjRPzn62mKRsRl84HvDj03iCDiPyIaG/QIepH+VQoiWs7Ff2txBeF/9fzQcTcw7Is1EtqBPFJJa0ehkS4q+Q9H+JQveIJ6dOhJRM9shBdrM/5FmxJhuI50HQxbJ8YjztxuRKnwfM1rOhB2P28h2of7EAWwXITXgvS7zWVYJa0gbhilCCvH4mXL+Nff/qbehB0OB4qiOP0L4qYeHx9PfHy8ugW32Wzs378fRVGcpFxaNy4tGZ1KBKNBwJ3niKv0hN1uV9UT7e3t1NbWqlrawcFBtUU+MjLSr3X5ErkHy+zmk08+4c9//jPz589X27rXr1/PihUrxvhN/zDhSNdms1FfX09raysRERHMmzfPKwH13r17WblyNR0dncBNCAlYsApkIBSfTyLitrsQjak2BH1UIPy93kDodpMQti8LEcQ9m7FzxFWIrrfmoeNfgO+0GIaIGWcCFw9FtQoKzXTxCnv5K9HouAiFCESprpjhaDPDy3d0lV44jtjq24Gv4dynp32tnmErSC0OAGUxMSQvXsyuI0fo7u6mX1GYgSDhLoQqISU5mbe2beOCC5zdGuQXXRupSeK12+0oiqKSscUiSpIpKSmkpqaqTm6y02wkGWmJWDu/K1jjbYLRkebtXDGpiBjZYTYwMMDBgwddTq2Qn9/d1IqRa/HFYSwYhbQLLrggaH8Lb3DWk668m9rtdnVsemZmJueddx5lZWVemd5ccME32LdvL6L54AEEYQXr1FQivHJbEYY032eYzMMYbqCV3llWBBFXIlIT7yE62RIRLmNaIo4dOu7vEBHuzUPHCeaftRl4BhEfzmGAMvYOPfs+sA/hNtaFIEe57V+IIOIsRm/7tT93I/SwrUOvP8fF673BJCApPp5333sPi8VCdXU1DQ0NHD58mI8//hhDaSnXX3klq1evJi4ujqamJuLj4z3mKeXz8l+Hw0F9fT1tbW3k5eURHh6Ow+FQrzGpIEhJSVGvy8HBQdWNq7GxUXUki4uLw2az0d/f77NZeLARaIpCamkNBgM5OTnq866mVsgJF9pcsZZkx9tL93TgrCddh8PhtovMW/+F+vpmRM5WQRDMLxDRZjKCEL+CoBJfTlfr0HEqEHYuN+Nd44QBIRXLY9jxVg5ElzN4P0LEhhEMx4ErESY3wfqT9gJ/QJSxpiL8v6KBCizYWYQY1a5FGYKI9w79exJxC8lBlB8XI4hVmtjsQvThpSMy2u7aUXyJQY4dO0ZzczOzZs2isLCQyy67jHvvvVf9f+mDbDKZOHbsGL29vapRS3x8vBrFjcxNyqkgqampLFu2bBRJuUtPyC16cnKyk6FLT08Pzc3N1NTUqOOctFGhv7aI/iBYudiRcDW1wm63qyZAciis3W5Xh2rKFI43ud1gqxdOF8560tXpdFgsFpdj0731X0hMjKOjYzmCGEHEbV8gIs19iCJYD94RcS8il1qCyEQ+jJCABQJt/vVrCHVEE4Jkz0VkQnchzG8MCEuZDATNXYRwQPAWNkTs+jsEyV6Js7hMx4Cb3ywcemhxDFHS2zO0wpMIvUcEQj1wDZ69eN3lf0dCAcyDg9jtdoqLiz3KqFwRgVa1IN3AYmJiiIqKoquri8jISBYtWuRWu+suPSEJWD7ka2XeU9YatGuQRSttVCjJ+FQU7E6nf0NYWJiaI5fQ+vO2t7fT09NDaWmp6lHsTjliMplCpDse0Ov15OTkuLzTeku6SUnRCLJUn0E4gGkHjWuJuITRRJyPoJS9CFL+H0SMFyw4ECY3/xg67gsMZz4vG/rXjpjwJbW6OxHKhXAEEacjNv4XIVIVWigIapR55wtw7YCmx+zDqmcipj78/5rnShGJlpU+HMcdZPnxUyA/LY3sbN8HB8mef+0X2Gq1Ul1dTUtLC/Hx8VgsFj777DO13TYuLo74+HiPuVBXRGy329Wcb2ZmplNULCPcsLAwlYglGbW0tNDb26tGhXINrsyefEUwIl1fcrEjoW31lbuOrKwsJycyqRzR6/WYzWZKSkrQ6XSYzeaATcxXr17Nm2++SUpKCocPHw7oWN7grCddTzAYDAwMuIvLhpGUJCdaeXwVw0S8Zug5LRG/jSA6AyIK/RX+pyZG4g2E/CsWkR8udvM650KYgHQ3q0AkAD5FiK+k220aIvNaicglzxlar7uoxTfSdYVkvG+Q8NQcIUfvDBgM/Pqpp1i9enWAKxNRV3t7OzU1NaSnp5Ofn6+Sp4zIenp66O7u5vjx4wwODhIZGammJuLj4902FvT391NRUYHBYKCoqEglbHcFO4Do6GiioqKYNm2aWrCTeeLOzk7MZjN79+5V/Wn9aW4IRqRrs9mC3gLsyqPYbrfT0NCAoiicPHmSSy65BLPZzD333MMNN9zg7rAecdNNN7FmzRpWrVoV8Pq9wYQgXZkDGglvI12x1eka83WjoSXiIgQZf8Cwi+w+RKrBhH854n3ARoQo6g6EfM3XC1vrbvZ/hp4To36EIOs3Q2tNQJgtenb1B13AHg6Bei/0I24dtcClV1zBxo0bSUpKCjhi6+vro6KigsjISIqKikaZHmkjstTUVGDYG9ZkMql52oGBAbWpQhbr2tra6OzsJC8vb1TxZ2TBTkKmJCQJy4JdREQEycnJJCcn09XVRVFRkVq0GtnYIYnYU2NHsCLdYJHuWGZTs2bN4p577mHHjh18+umn2Gw2r4Ird1i+fDn19fV+/76vmBCk6w5j2TsqiqJeoATsTJCM2JZPZnRq4iTORPw4gojl/N0CRP5Xzt9tQEyFaEBMC1uJ/8PHXcGCUEVsG3r/QUTKYizChWBEur56L8jX2hkenDlnzhwO/e1vJCYmYjKZaGxsVEfdS6KR0edYZGCz2airq6Orq4v8/Hyf8oRab1htY4H0x21paaG8vFz1R2hpaaGvr4/4+Hi1cOQOer3ebcHObDZTXV1NfHy8OrInISFhzMYO2ZGmNcE500g3OnpsnXpvb68qW5Odc2cLJjTpurN31G4hk5KSKCws5B//2BXgu00BtzGgN0S8HkHEBkQxKxrRCfc1gvdncgAfIqaahSGmWmQg0hfe9poFbhHpbXEMhnW/DYhUQmRcHNtffJHLLrtMfc1I68Pe3l416qysrFSLYlp1gsFgQFEUWlpaqK+vJzMzk9zc3KBJt+x2O42Njeh0Os477zyioqJUEpTFspE3CfnwlBvV6XQ0NTXR2NhITk4OKSkpPjV2yPOj9ejt7e2loqJCvRH409hxKmwdPcFoNJ5Sf4RTiQlBup48dUdGutJdLDo6mkWLFhEdHT00XG+snO5YSEQQlx3vUgDuiPj7iK3+JERn2Yahn6cg8q1fYTgi9gWfI/x+OxCqhvma/wvDe9PGwM3QfYl0FYS3ba1Ox90//Sk///nPx4wOR1bHpQdtT0+PerO1WCxYrVaio6PJzs5WI8RA4XA4OHbsGK2trcyePdspHxkeHu7SgEbeJFpaWqiqqlLnqGlvEhEREZhMJsrLy0lgbESlAAAgAElEQVRMTGTZsmUqyY2lnBiZJ46JiWHSpEmkpaWh0+koLS1lxowZPjV2jMTptHUEQbpno3IBJgjpuoNWp9vb20tlZSXg7C4GDJkj97o6hA/QI05nP8JjwR9MRngk5AMPDT3XhYiIDyKUETIi9paImxDa44OIFMJljM6q+ka6wUgvjAUr4tPWAwUFBXz4r3/5HdloPWitVis1NTVq376iKE6NC5JoJOH50s568uRJKisrmTZtmkstr7u1ubtJmEwmOjo6qK2tVX0Ppk2bRmJiIlarVY1e3R0XXHfYaXPFXV1dOBwOp8YOaYvqqbFj5Cy3U1FI84Tu7u6g2zqeLkwI0vV04dlsNr744gsGBgaYPXu2yw4WURQJlHRBpAb68J90QaQVRsrXLhx6SEjVxCFGE3EywnmgGJG2eBOhE/4B7vPCYXjvHxZ4esFTpCslYB8DqWlp7Pr739V++EAg8/cNDQ3MnDmT/Px89brRFsXMZjM9PT0YjUav1Qlms5nKykoURWHhwoVe5SQ9Qd4kpHzs5MmT5ObmkpSUpPo+a9emlbB56m7TFuykJM5sNrN48WLCw8PViFgW7KSmeWRjh0yRSLN0SZLR0dGYTKaAGju8Je9geuled9117Ny5k46ODjIyMli3bh233HJLUI7tChOCdF3BYrFQW1tLf38/eXl5TJ061e3FKL50/lc/hxGOIN1AMMmLY7jSEcscsSTiRxGFsW8jImJP8DbS7QPMvIW4rcxHxM6+xjfucrpSAtYfHs7td97J1VdfTUxMDK2trR6lWGPBaDRSWVlJfHw8S5cuHdVEI6HT6YiOjtaknNyrE6RwXxbMZs+eHTSbQYCBgQHKy8uJiIhwUlLExMSoawPRZtzT06OSYH9/P+Hh4U6pCS0JKopCa2srdXV1ZGVlkZqa6vKcussTy8g8MTFRLdjZ7Xaqq6tRFCXgxg5fcrrBIt2XX3557BcFEROOdO12O/X19bS0tJCVlUVMTIzTbCVXSE9PRzh+eTtJyx0MMOTE5T9iEKY3vmJkjviHCGeDsQgXBA16ytTaEOmJg4CdvYRxFAMKZiwIa56lCIcvScSeLqyRka5WAnbt9dfzm9/8hvDwcJVQRpKdNur0FNlJD4aBgQEKCwv9EtG7Uye0tbWpk0diYmKoqamhvr7eaW3eGLyMhMwJt7W1MXv2bKfOOVdwZ8kobxL19fX09fWh1+vVSHTSpEksXrzY42QMb/LEkohlhJycnMyUKVNUCads93XX2BEbGztKHuatZvhUTI04XZgQpCurso2NjRw/fpz09HR1bHpDQ8OYf0gxNjoC4SawELE1l8YyvhhqyPRCIIglcOIGkUrwNk8bhmv1gixjfYq4GUUg8s219HLd0Gt6qKCWChrZQTsKA1gQ/W5aIs5l+GKTpKuVgBXOmcPhV18lM3PYRt0VoVgsFpWIZWRnMBhGEbGs8M+aNYtp06YFTZVgsViorKzEarWyePFip5HjWl+Huro6ley0239PMrGuri4qKytJSUmhuLjY7y36yFZnadTT3NzMlClTsNlsHDp0CEVRRsnrPEWZrohYNnzIIqGrxg4ZTY9s7Dh27BhWq9WpsUMSuje+C1lZWX6dn/HGhCDd/v5+SktLmTp16igPBqlg8CS4FtAhur32I2RVLyGm7iYgyFhrYuiuoBOM9EIMwSHdaHwj3ZE53TaEW4Jp6DhFiBvRSYQvsEQ8DNnfDGeiTVRSTiX72JMylY6ODgYcDmYNHUW6B28DDLGxbP/f/3WSgHlCREQEU6ZMcTKjl1v8np4empqaMBqNGAwGpk6dqo6niYmJCYh4FUWhsbFRlWq5Sle58nWQY9R7enpGycS0xbra2losFgsLFiwIOCeshdFopKKiguTkZDUQkZDnRt7AqqurR0Wj8fHxLr870ku4sbFxlErDVcFOPicbO2REDDg1dgwODlJaWjpmY0co0h1nREdHO7VVauEt6er1kTgccxierACCbHYhSPhNhOtWF2IrvwBhNrMAUbiKITiRbjwErA8A38ZIhmte24eYq1CPoMZsnGe+eTO1LA5xTj7j89paVVv62muv8c477/DxoUNEGY38+N57+a//+q+AW1BlfrW5uZmwsDDOPfdcVWLV09Oj5valF6wkO28LPpK4kpKSnKRa3sDVGHVpbtPT00N1dTVGo5HIyEgSEhJoa2tz0hL7C5vNRnV1Nb29vcydO5eYmJhRr9FG4RKKoqjyupMnT3Ls2DEsFgtRUVFO0XBNTQ3x8fEuzYU8ddi5Sk9oz1FnZydFRUXqzUoqOAYGBtQZbv/6179oa2sL2s3p3XffZe3atdjtdm699VYeeOCBoBzXHSYE6er1ercaQm/tHcPCDDgcIxUM4cDXhx4SZkS55yPgFYQJTTdilkIvgqQzEPGcP11k8QSnqOdLpKsfeu1+hG8uCN3x1xmdE9bjndJBkPOePXuIjIwkPDychQsXcskll5CdnR20qQoOh4Pjx4/T3NxMTk6OUwTlKuqURCxtHT1t/7U5YXfE5Q/CwsIICwujtbWVuLg4Fi1apI6WMplMtLW1OVkeam8UY+/YUH9/xowZTioNbyALYDExMeoMNa2qo6GhAZPJhMFgoLe3VyVfb+bHecoTm81mqqqqSEhIwG63qwQbHx+vFuykptlkMlFRUcGdd95JeHg43/rWt3jkkUe8/oxa2O127rjjDj744AMyMjIoLi7myiuvPGXj12GCkK43Y3jGgsEQgdVq8uLdohBaV+12uBdBwr9CuIB9OPTcTMTWewkiIs5n7FbbOIIT6Xqb01UQJN859L4RiMg2G9dFRV/6yRTmz59PeXk5g4ODpKSk0NfXR2lpqSrDkl9Yf5QJJ0+epKqqiqlTp3q0c5Rw1Zyg3f43NDSo/rp6vZ6BgQEyMjJUw/JgwG63U1tbS3d3N/n5+U76XEn+06eLeXbaqLOzs5O6ujq1oUObw5YBh9lspqKiAr1e79I7wl9I+9T6+npSUlJYsmQJer1eza/LG4V2N6H1nfC0m9DpdLS1tVFXVzdmhx0I9cbdd9/NBx98wBtvvEFCQgLt7e1+f7a9e/eSm5urutOtXLmSHTt2hEg3EIzlvwCicUJ8Yb0hXVeIRUizDiNmoP0ZoZH9GPgPwoEsErF1z0F0hC1EEPFshid5gYgwA+35AhHpjhWRtiOi9m5EBFuIuEmMpT3wdFwF4aK7Ewjj6NGj5OXljcq/ycjJZDLR1NSE2Wz22q1Lq4sNNAc6cvvf09NDeXm52ijQ29vL/v37Aec8rDeeDiOhdS9bunTpmDcZd1GntHvUOp0pioLVaiUjI4Pp06cHlJrQQsrBTCYT8+bNc4r2XeXX3ZnEa9uc5bkbHBykrKwMg8HgJOXz1NhhNpt56qmnOH78OJGRkURERAypj/xDU1OTU/E2IyODkpISv4/nDSYM6XpyGnOXXpCmIaLQEo3JFGiDhGwFLkTU7KUpuhk4iiDiTxBOZP9guPCWx3ChKp7gkK6nmVd9CEVCHSKiLUBoe+cw9iXhiXS7gH8jWo3nAwcoLi52SS7aDigJs9msbv9PnDgxiohjY2NpbW2lra3N42RnfyA71Xp7e5kzZ84oeZk7o3NJxDI14Soi1kagixcv9noemStonc6mTZuGyWSirKyMuLg4kpKS6Ovro7y8HLPZ7CSv82d8fGdnJ1VVVWq0783vemsSPzg4iM1mY9q0aWpziifo9XoOHjzI2rVrufLKK6mrqwvKjcUVZ5zq0UkThnTdwWAwqIMEJaSrVHt7Ozk5OcydO5eEhEm0tPgb6UrEIMiuH2eFQxQixbAEuG3ouT4EEZciIuLXEDliBUHSlyEUE3LQuC9Dxhlax8gLSuptDyA6165DpDPKEakEbwpvrvrJBhkemp4BrEJcWgd8uoBd6WGlVrelpYUjR44QFhZGbGwsXV1d2O32gJomQHzpmpubOXbs2KhONS1cGZ1rc4zNzc2YTCYnc53Y2FiMRqOqudVW+AOFNk0xsq1dQsqz3DVOuNMSW61WKioqsNlsHqdleAvtuRsYGKCsrIzY2FimT59Of3+/ugOw2WxOygmZdhocHGTjxo3s3LmTP/zhDyxYsCCg9WiRkZHB8ePH1Z8bGxvV9M6pwoQhXU+Rruxbl0UXObzy3HPPVS+4yZNjCLwVOBaREx1Juq4QgyDVYuD2oed6EMY01wPfQUSjOxEyLRsiFbEMERXPR3jkuiObKIbJUUG0HuxC/MlXIMhRImzoON4UyLSRrsKwMXoMzsN3FLzP/bqHw+GgqamJsLAwvvKVr6hfQqnVlRGxjOrkwxsi7u3tpby8nNjYWI+dau6g9U2QW1wpw2ptbeWLL75Ar9djMBhoamrCZDKphBdIlCbntY2VppAm4K62/z09PS61xFarVS1KagdsBgopuWtqaiIvL0+NhLVpp5Em8e+99x4bNmzAYrGQl5fH3Xff7dSNFwwUFxdTVVVFXV0d6enpbN++nW3btgX1PUZiwpCuO8j0QktLC7W1taSkpHDOOeeM2gYmJsYjcpuBYBLilPpL3vGIMTkGxBh47ZCbakRueA/CC/ckgtQKEDN0ZZ44leFGBgWRt/03YERE2otdvK8e3yPdZobmNgytuWDE64Z79f0pQsnOwo6OjlGdWa6aJnwhYpvNRm1tLUajcVQxK1DIqdRms5ni4mJiYmLU7iyty5nNZlPbZL1VJlgsFioqKnA4HH5HoO60xJ2dndTU1KjevPX19XR2djqpOvw1tOnv71ejW08FT23qJCkpiW3btpGWlsa6deswm80cOHCA5OTkoBJveHg4zz//PJdeeil2u53Vq1czd+7coB3f5Xue0qOfRri7I/f399PSItpq3Wl5gaFtozHAVcQw7DQWCCIRZKb9UuUCdw09JI4A7yJmtr2BIOJwRG42dugYryH6w76D+z+3/BJ4S7oORE46H2HE4y7toaO/v98nUpNex7W1tUyfPt3rzqyxiFi2EcvnU1JSmDt3btC0nlpv3pFdcHLulzZP7E6ZoJWIycYErVmPjECDBZleaWpqIj8/X02BaPOwTU1No3LY3nj/KopCQ0MDzc3NFBQUeN3MsH//fu6++26+973vsWvXLvU9vvWtbwX+gV1gxYoVrFix4pQc2xUmDOmOhLRylBfKWBIQcbGVB/iukxARXqANEhEIwhyrBXnu0EPCgTC8eRcxBN2CmDoxFunJSNeb9MKJoX9vYLhhwh10fPrppyQnJzsRibvuMO24nEALTuBMxPLYBoOBrKwsBgYGqKysdPJzkGTiyc/BFWTxKiYmxus0hSdlQk9PD11dXRw7dgyz2azKxGbNmhVUD9m+vj7KyspcNjm4y2G762DT5omlhresrIykpCSv5HwgCo4bNmxg9+7d/OlPfzrlEed4YcKRrlaRkJeXR1xcHAcOHBjz94Sd3RFEMWsBQtrla6dULGLrHaxI11foEemDxQiP3XsYm3Bh2HvBU6RrRMjf2hg2WR8LOhYuXMiUKVNcdodJops0aRKtra10d3e7lJcFArvdTl1dHSdPnnQ7ikfr1OWLsY43x/YF2u11SkoKdXV1anpFURSn+WdRUVFONzJffH+lF0NHRwcFBQVe70TcdbDJpg7p/dvf34/D4SA1NZXExESv7BpLS0u55557uPbaa9m5c2fQdNFnIibMJ7PZbFRVVdHR0aEqEmRxzZvmiBkzZiC2588hJE8WxPZ5GcPNDRl4diGLQUSLgUa6UQTelRaBLx657iNdK0JhcRSh4U1DSM28gRiR7SqPaLVa1fxrR0cHYWFhREdHqx6tkowDUSXINEV6erpb6RqMbazT0tIyiojlVFpZzAq0lVmLrq4uKioqSE1NdUqvjLSblOvzRedsNBopLy8nJSUlKOvWpk6kWiMjI4Np06bR19dHV1cXDQ0N6o1CK/+TRdH169dTWlrKSy+9RGFhYUDrORswYUi3u7ub6OhoJ0UCeK+5y8nJAaYiikMgvAfeREzMfQ/RsaVDGN4sQ0STC4Z+R2ISgqQCJV1/I92Rx/BWPRCGIF5tpKsAVQhdcToiTRGLSF94S+Y6txrpwcFB6uvrmTRpEueffz4REREqEff09KgdTtJBTD682fprR537m6ZwZ6wjozlZIGxvb8dsNvu0PnewWq1UVlZisVg8mqFr7Sbd6ZxHRuwxMTGcPHkSs9k8qskhUDgcDmpra+nq6nLSOMfGxrq8UZhMJt577z1+9atfYTabKSws5Pbbb3c5YGAiYsKQbkpKSkB/NJFX05JlFmKk+hrNc58jFAS7gb8i1A6TEKqBZYgKvpngRLqBtgJH4j05ShNzSbptCMWDA6EXTtO81lvvBQDd0KTlYVitVmpra+np6RmlHDAYDOpocQmtg5jUmrpTJWgVD65GnQcCWSg7ceIEeXl5Khl7spr0log9FeF8gTudc2NjI1VVVURGRqIoiprH1Wp1/b1RyMh52rRpFBUVuY2ctTeKuLg4ysvLyczM5Je//KXa9TdlyhSvGiXOdkwY0h0LY3l0CtIdK7pcMPSQcCA0qu8D7yBagG0IUp6L0NL6bpwtWnhPd3pBQdwsPkTM3j0XcRMZ+SWSqQhvoFMbU0Y2IXjb4SStAEcS8UhVgk4nCD45OZl58+Y5+dwGCtkaPHny5FFFIXcR8chpDloi1naH9ff3q23H/miFPcFisVBVVYXdbuecc85RJWbaG1l7e7tTjl3rmeDp72O329VZc75Ezrt37+a+++7jhz/8Ib/+9a/Vc3n55ZcH/oHPEkwY0h3L9GasMSCiOGBBEJW3eS49QqN6gea5XESU+RBwHJF+WIggsYWIFuGxtrvBIF1vI91exE1CQaQS5iGaM9yt0ZtZvoMIPbGNpKQkenp61BHfwSAWLdFJFYLD4WDGjBnqz9ocpz/FJhhuDe7r63PZGuzN+iRcEbHdbsdms5GRkUFaWlrQikfayDk7O3uUrtXVjUzbNDGWFabMOaenpzN79myvzml/fz+/+MUv+Pzzz9m+fTt5eXlB+axnIyYM6XqCbJBwdVHL7e7Jk1Lj2ot3FX93MAAPI5oUrEAFwjLxE+D3iFE8mYiuMjmhQjtXAUTK4lSTrg1h4/gZwh/YAVzL2DK1scZKlgOfkJGRzl/+sovBwUGqqqr8HpfjDp5GnUtoc5yy2KQt5sjUxKhPoZkj5qk12BdoiVhuyadMmUJiYiK9vb1UV1e7nILh69ZfzlaLjIz06QbnrmlC3ijq6+vp7e3FYrGg0+nIzMwkISFhzB2koihqdHvzzTfz9NNPB2Vq8NmMLw3pjlQwaFuCZ86cyezZsxFb8kBJVzs9woCIHOcBNw49N4CQppUCbwFPInLDsxhu8XUgBtlcgChi+fOFj8L9CB7ZEmwAvokg2j/h3WgidzldMVYyIsLMY489xnnnnUdnZycREREkJCTQ2dmJ1WodU1DvDbwddT4yxzmy6t/Y2OhUVZca09ra2lOy3ZfG4n19fU5b8pGqCXmj0BYTtfIwV0Qsh0LKnPNYs9W8QXh4uErEnZ2dVFZWkp2dTWxsLCaTSbXClNOLtU0T0h943bp1HD16lL/97W/k5uYGvKaJgAlDup7utlp7R0VRVJPnlJQU9Utrt9sJCzNgtwfqvxCGZ51uNMLEZinw46HnjAiS3YMYYnMIQYivIEhuHqLVdxGjFRPuEMFo0u1EqDO6EOQuW4IteK90GElwA4j0RA3XX38t69f/N7W1tcTHx7Nw4ULCwsKcBPVVVVWj3Lm8tUkMdNS5q6q/NNA2Go3U19djMpmIiIhAr9dz/Phxp86wQCCvubEiZ3c5bE9EHB4eTl1dnU+NCN7CarVSVVXF4OCg0zBLbZHSbrerxjqNjY28+uqrvPLKK/T39/OVr3yFX/7yl2RkZLh7i6Di6aefZsuWLeh0OubPn8/WrVsDNuwJNiYM6cLY9o7d3d1UVlYyadIklixZgsFgUOc46XQ6wsIisNt7AlyFAd+9FxJwzg1fiyDlnyDI+G1EeuKvCC1xDIJ8z2HYl3dkdC5J14Eg1T1AJcKt7Ns4/+l9MSaX6QUHImLfQ35+Pn/96z4GBwc5fvz4qCkLUsMp3ZtkZ5PRaOTEiROYTCb1dVoilhGsdjpEsC0ddTodAwMD1NfXk5qaSlFRETqdTvX7lZ1hFovFZYvuWDCbzZSXlxMeHu63sbg7IjYajTQ0NNDT00NERARGo5Hq6mqPEbEvaG9vp7q62uOodhDda9KTuLe3l+7ubjIzM7n99tvp6OjgxRdfxGAwUFxc7PdavEFTUxObNm3i6NGjREdH8/3vf5/t27dz0003ndL39RUTinTdQVEUampqiIiIoLCwkJiYGCeylReTsIEMNNINxhh2WUjTI0h1oeb/HAgCfg+hmNiKSE9MRkSu5zA8QFOPiJr3IUj5e7hOIWiJdKwioh6Rq95OdLTC73+/laKiIk6cOEF2drbLgY2jjuCis0k7N6yxsVE1v5azzpKTkykqKgp6db+yshKbzTYqco6OjiY6OtpJZypbdLWzw2QLrDY9IV8vt/vBtnUEUZiqra1l2rRpLF68GL1e76Rzbm9vp6+vT7Vy9IWIpbGOoigsWbLEK52zoij85z//4YEHHuC2225j8+bNQW0Y8RY2m42BgQEMBgP9/f2n3KbRH0wo0h0Z6crqc2trK1OnTqWwsBC73a7OYBp5UURGRtLXFwzSDVSnOwn3Uyz0CJOZCzXPmRHjgj4CtgObhn7fgSiUfQ1RrHMH3dDDhudxQn2IyNvGj398Ew888AA1NTVYrdaAt7Uje/0tFotqxp2WlobZbGb//v1OFXVvpE2uICfZHj9+3GsDGW2LrtSSStMaOVK8rq4Om82mfuETEhKC4kerhTYvPH/+fCdpnCud80giHtmCLVUJ8hy2trZSW1vrUvXgDiaTiYcffpja2lpef/31cRuNnp6ezr333suMGTOIjo7mkksu4ZJLLhmXtXjChCJdCYfDQUNDA01NTWRlZZGQkEBXVxcWiwW9Xu8U3WoRGRmOsFA0499QSRCkdbrHsEcBVww9JIyIFMWVjB4u6QqeSNeOaAwpZdGihWzb9i5Go5HGxsagjwwfa9S5rKhLaZM2mktISBizGcFkMlFeXk5CQgLFxcUBFfW0pjWpqanqaJuuri4yMjKwWCx88cUXTqYw8uHP+0of3czMTK8VFZ6I2GQyqUSs04nuQbkb9MZHQlEU/v3vf/Pggw/y4x//mN/+9rfjEt1KdHV1sWPHDurq6khMTOR73/seL730EjfccMO4rckVJhTpSn2i3Hadc845an6uvb2d/fv3o9Pp1As/ISHBabtVXDyPN9/8Iw7Hc8AMhk3GFyIMcLyJ5KSJeSCIQfg/BAJpSuPLRGBXLbvHgZ3ExhpYv/4pZs+eTXV1NYmJiaSkpGCxWIiMjAzKl82bUefairqENppz17UWHh6udsIVFBS4nLQQCDo6OtyOttH66Won/coJE9rR5q6g9dENhvualojld6auro7p06ej1+tpaGgYs2HCZDLx85//nIaGBt544w1mzpwZ0JqCgQ8//JBZs2apapCrr76aTz/9NES6pxLV1dWYzWanIpmiKERGRqo2cXa7nZ6eHoxGoyp8l1/Qp57ayCOPiCp2WVkZu3fvZu/e39DR0Y6iSAOccxieZ5bGaDlXFIFPoPA10nUHA6KI5g1Gkq4J+A863QnuvXctt99+u3ozy8zMVEmksbERk8mk5mnlzcyXbX+go85dRXNaaZjUmEZHR5OSkqLm/IKx7R8cHFTzn9rqvhZaUxjtpF9XNokjpVft7e3U19cH3UcXRJGvrKyMyMhIiouLR+XLXTVMbN68GYvFwqFDh7jlllv4zW9+E9Q8uydUVFRw7bXXqj/X1tbyi1/8grvvvhsQplV79uyhv7+f6OhoPvroI5YuXXpa1uYLdK6q/RoEPm/lNMJisWC321WydZdGGInBwUHq6upoaWkhIiJCzd9pt6z19fX8/e9/51//+jeHDtViMnUgotr5iG4zaYBzJ2LQ5M8D+CQvIBQLbwRwDIDzhtaX48Vr/whcjpCjHQAOcv75X2Hr1i20t7ej1+vJy8tzS1TyZiYffX19TrnDhISEUdt+bW41EM8Bd5CNAgaDQbVH1K5xcHBQHWfuqzRMu/bc3Fwnra2/0PrVdnV10d7ejk6nY/LkySQmJvokr/Nm7Y2NjT4V+Xp6enjwwQc5ceIES5cupba2ltraWj799NPT3vBgt9tJT0+npKTEKcp+5JFH+Mtf/kJ4eDiLFy9my5YtAe8M/ITbC3lCke59992nzrsqKioiLi5uzC/xyZMnqa6uJikpiaysLAwGg1qpNhqN6hdURiGShGNiYti7dy87duzg3//eTVVVI2azHKGTgRi3swjR9utrRPV/gRcRU4MDwVcRxj3e2OW9iPBaKCc5OYGXX/5fUlNTXY7L8RbabX9PT4+67U9ISMBgMNDS0kJCQgK5ublB9U/Vdqt5ahSQGl3tGj0pEiTkbLW4uDhycnKCuvaRTQ6JiYmqvE5qYRVFGTWF2FvSGxgY4OjRo8TGxpKbm+vV7ymKwj//+U8eeugh1q5dy8033zyuuVuA999/n3Xr1vHJJ5+M6zo84MtBuhUVFezZs4eSkhIOHDiAxWJh3rx5FBUVUVxczNy5c9UvUFdXF/X19YSFhTF79uwxi0Fy8qtMTUhJkzYa1uv1bNu2jTfffJOjR+tpbm7Dbu/DdX7Y00X7D+BXiCnBgeBSROQ6f4zXdQN/Q6eDRx/9OatWrVLH5WRkZAT1C9bb20tVVRW9vb1ERUVhs9mIiopSz2GgjQhdXV1UVlaSkpLCzJkzfV67doyO9oYrNbp9fX309fVRWFgY1NlqgDptITExkezsbLeEqL0We3p66O3tHUXEWp2z/FySzH0ZnWM0GvnZz35GW1sbL7zwApmZmUH5rE7et98AABrMSURBVIFi9erVLFmyhDVr1oz94vHBl4N0R8JsNnPw4EH27NlDaWkpR44cwWAwYDAYiIyM5IknnqCgoMBvUtFW0o1GI/39/ej1egYHB4mPj2f27NkMDAywY8cO3n//fUpKjrrJDy9ieKAkiBbdnyBahQPBlQjN7xI3/y8Nyg+zaNEC/v73v9Lc3ExkZCS5ublB3ZZpTVhmzpxJWlqaKvFzFW1qi0zeVPulo5bFYiE/Pz+oLmNyRplsD1YUxe+uOldwOBzU1dXR2dnpdpy6N8eQ+VdJxDKXHBkZSXt7O5MnTyYnJ8fr6PaDDz7g4Ycf5p577mHVqlXjHt1KWCwWpk+fzpEjR4I+HTiI+HKS7ki88sorPProo6xYsYKoqCj27dvHsWPHyMjIoLi4mKKiIpYuXUpSUpLPucWBgQGqqqqwWq2kpKQwODiI0WhUt6vaSM59flh2mcUCTyM0sYHgewi51zkjnlcQ0rj/kJycxH//9yPk5ORgt9vJyspi+vTpQc3R9fb2UlFRQUxMDDk5OWMWXmS0KdM7JpPJqcikJTnt0Mbs7Oygjg0HZzIvKChQd0Tarjpvo01X6O7upqKigmnTpjFjxoygEpvUqbe3txMbG6ua1cg1yoLnyPfs7u7mwQcf5OTJk7zwwgvqePkzBTt27GDz5s28//77470UTwiRLog2wcmTJzulEuS8qJKSEkpKSti3bx8mk4nCwkKVhBcuXOi2gGSz2dRx1bm5uaOKEtrtqvyCOhwO4uLiVCKeNGkSpaWlvP766+zatZvy8gYsli5El1kRw74LvuaHf4jQ656veU74L4SHm9i48Zd8+9vfVlUJUVFRKslJApFrjI2N9ZkQgjnqfCTJSSK2Wq3ExsaSnZ1NQkJC0EjLH2NxbVedJGLA6W8tSU7b5FBYWBjUyByEpKusrIzk5GRmzZqlnhfpk6BdozSnOXToEFFRUWzdupX77ruPG2644bRGt93d3dx6660cPnwYnU7HH//4R84777xRr1u5ciWXXnopN99882lbmx8Ika4vsFqtfPHFFyoRf/7554SHh7NkyRKWLFnC0qVLyc7O5rXXXmPGjBlkZmaqGkdvoN0KyvzwyEq/zWbjpZde4sMPP6SsrJ4TJ9qx2XoQ3gnFiMYHaQvp7n3/P6AJuAjhcVsClHPVVd9m06ZnqKurIzo6mtzc3FF5VEkg2khOdoNpbxauiEhrjZiZmUl6enpQo0+73a7acWZkZKjKCe0aXemwvYXWWHz27NkBSaJckZzD4cBisTB16lRmzJhBbGxs0M6PP6kKu93OwYMHWb9+PTU1Nep0hzvvvJOVK1cGZV3e4MYbb+TCCy/k1ltvxWKx0N/fPyr33N/fT2ZmJrW1tUGdjHwKECLdQCAnse7bt4+SkhLeeecdvvjiCwoKCli+fLkaEQcieZKaSKPRSGtrK319fUyaNImpU6eq7bE9PT289tprfPjhh17mh+9C+NumAZ+SlTWTv/1tG2FhYS7H5fiyRq0aQVtMtNvtVFRUEBUV5ZLMA0V7ezs1NTWkp6eTkZEx6nxr8+xSumYwGJzSO64GNsKw6qGtrY38/PygTiWG4SYHu91OWlqa6uXgr0fCSPT09FBWVuZTEVFRFN555x3WrVvH/fffzw9+8AP0ej0mk4n+/v7TljPt6elh4cKF1NbWBvUGPY4IkW6wsHv3bjZs2MATTzxBbGysGg3v3btXlVdJydqSJUt8+vJIF7SkpCRmzZqlEogkOVlg0hJIXV2dh/xwG1BNZGQMmzc/zUUXXaSOy5GFrEAhmxC6u7tpbW1Vi4jJycnqOoMhnjebzVRUVKh6YV+KfNrxPkajUZ0qoT2PAwMDVFRUMHXqVLKysoK6rdaOKnLX5OBOXudqFtxIyMi/u7ubOXPmeN1ccvLkSe6//34GBgbYvHnz0Miq8cHBgwe57bbbmDNnDocOHaKoqIhnn302qAM0TzNCpBsseHLKt9vtlJWVUVJSQmlpKQcOHMBut7NgwQKWLl3K0qVLKSwsHFWJN5vNVFVVYbPZyMvLc3uhabuYpG5TURSn7XR0dDR79+7l9ddf5513PmDOnDx++9vfUlNTQ3x8PNnZ2UHtIBo56jw9PV1VI8ibhWx5lSTnS6Vfa+sYLLcurZl5V1cXra2t2Gw2EhISmDx5slt9rj8YGBigrKxMTeP4ckx5s5Dn0dUIImkdOX36dDIzM726kSqKwltvvcVjjz3Gz372M1auXDnu0eW+ffs499xz+eSTTzjnnHNYu3Yt8fHxPPbYY+O6rgAQIt3xgGyy2L9/P3v37qWkpEQ1WykqKmLhwoXs2bOHefPmceWVV/rV1eQu9yor0ydPnmRwcJCCgoKgjssB51Hns2fPdht9ajut5M0CGNU2PDK6lF4MycnJZGVlBVVRIc3sa2trmTlzJqmpqW5vFlqS83YNiqLQ0NBAc3Mz+fn5QZtMrF1jS0sLVquV+Ph4kpKS1HPpKaXT2dnJfffdh81mY/PmzWeM5KqlpYVzzz2X+vp6AP7zn//w+OOP89Zbb43vwvxHiHTPFMjIcNOmTWzZsoXs7GxMJhMzZ85Uo+ElS5aQkJDgd/RhsViora2ltbVVHbstR9LIaDOQXGswRp1rC0xGo9GpbTgmJoauri4GBwdV/+NgQmssnpeX5/ZcyJ2FVjEhlSeeZGHSyUymiYLdIivHFWVkZDB9+nQnnwlXXXXyJvzGG2+wfv16HnroIa699tpxj25H4sILL2TLli3k5+fz6KOP0tfXxxNPPDHey/IXIdI9k6AoCuvXr+eWW24hNTUVh8NBdXW1mhvev38//f39zJ07VyXiefPmeZXHHJkXDg8PV7fTkjyMRiM2m02N4hISErze8kt7wbS0NDIzM4Oa+7RYLBw7dowTJ06oNwvtdjohISGgho1gGIu7a0KIi4tTZ4f19vb63eTgCTabjcrKSvVm5E7GOLKr7o477qChoQG9Xs/q1av52te+xkUXXRTUtXlCVlaWen2Fh4ezb98+l687ePCgqlzIzs5m69atQdshjANCpHu2wWKxcPDgQZWIDx8+TFRUFIsXL1aJODs7WyU9s9lMTU0Ng4OD5OfnjxkduoriYPSWX0ZDcrS5TqfzaHzjL/r6+igvL2fSpElOuU9tt5q22cSTN4Ir9PT0UF5ezuTJk4Mefdrtdpqamqivr1ej5rHMfnyFtI70pQiqKAqvv/46jz/+OA899BALFixg//791NbW8sgjj/i9Fl+RlZXFvn37gjpm6SxAiHTPdiiKQnd3N6WlpWqhrra2lrS0NKKiomhubmbLli3k5ub6HX3KLb8kYrnlB0F+ubm5Hmdl+fuesrkkPz9/TO2lO28Edy25drudmpoajEZj0MfAg4g+q6qqGBgYoLCwUG288WT2I9cZGRk55rm0Wq2qzKygoMDrSL+trY2f/OQnGAwGnnvuuaC4oPmLEOmO+I8Q6Z69OHToEDfeeCO5ublkZGRw4MABuru7yc/PV01+5Owvf4jy5MmTVFRUEBcXR1RUFCaTCbPZrNohBioJ6+zspKqqKmBjHVfdaiAGOvb29pKamkpOTk7Qu6vk4EZvo09t7tVoNKrj392Z/UjDc1/amxVF4dVXX2Xjxo08+uijXH311eOeu501a5baWv+jH/2I2267bVzXc5oQIt2JiObmZgYGBsjOzlafs9lsHDlyRDX5OXjwIDqdjkWLFqlNHPn5+R6313LUucPhID8/36ltWjuyXJKHNtKU+WFPBKc1/s7Pzw96qmJwcJDy8nKsVitJSUn09/erUbuW4Py9Gcn1AxQUFPhdlHRn9hMdHa0arc+ZM8frFuHW1lZ+8pOfEB0dzbPPPnvGRJYnTpxg+vTptLW1cfHFF/Pcc8+xfPny8V7WqUaIdL+sUBSF3t5e9u/fr6YlKisr1em6RUVFLFu2jNTUVKxWK4cOHcJut6vzybyBq0jT1VgkgMbGRpqamnw6vi+f1ZOxuNzyy3UODAyMapLwtH3XNjkEy7jc1fHr6uqYMmWK2gnpKX0C4vy/8sor/PrXv+YXv/gF3/nOd8Y9unWHRx99lNjYWO69997xXsqpxpeHdN99913Wrl2L3W7n1ltv5YEHHhjvJZ1xkGYue/fuVSPimpoabDYbX//611m5ciVLliwJyBNAOxZJEvHg4CCxsbHMmDGDxMTEoEa4/hiLa5sktF1/I13hwsPDVT8GWegLpnE5DMvYDAYDeXl5Tikbd+mTDz/8EBBdkmlpaWzatCnoo969gd1uZ+nSpaSnp/Pmm286/V9fX58qs+vr6+Piiy/m4Ycf5rLLLjvt6zzN+HKQrt1uJy8vjw8++EC1a3z55ZeZM2fOeC/tjMbjjz/Orl27uOuuu2hubmbv3r189tlnWCwW5s+fr+aH58yZ43P+VtoL9vb2kpOT49TaPDg4OGoskq9kZrfbqaur4+TJkxQUFARsLO7KFc5sNmO325k+fTqpqaleWTb68n7SmjIvL89r0rRarWzatIm3336b6OhojEYjUVFRvPbaa0GfpTYWnnrqKfbt20dPT88o0q2treWqq64CROrrBz/4AQ899NBpXd844ctBurt37+bRRx/lvffeA2DDhg0APPjgg+O5rDMeRqOR+Pj4UVGt2Wzms88+czKBj42NVXPDS5cudesBq3Uac1do8mYskieCk00Cp0IzDMP2iJMnT2bKlClq5592EKfWstEfD+aysjKfo+eWlhbWrl3L5MmTefrpp9VxRD09PX5ZcAaCxsZGbrzxRh566CGeeuqpUaT7JYbbi2FCTQNuampyGieSkZFBSUmJz8dZvXo1b775JikpKRw+fDiYSzwj4U6mFRUVxXnnnad6miqKQmdnJ6WlpezZs4ft27fT0NDAjBkzVJOfoqIiTpw4QXl5OXPmzKGoqMhtoUkOAJ00aZJqtqIdRaOdNDwyGq6ursZisajqjGBCRs9dXV3MmTNHlZklJiaSkZEBiKhNyutqa2tVNzPtOt0Z1CiKoua2fWkRdjgcbN++nU2bNrF+/XquuOIKp+MHe3yQN7j77rvZuHGjmvIIYWxMKNJ1FbX7k5O86aabWLNmDatWrQrGsiYMdDodU6ZM4fLLL+fyyy8Hhv1bS0pKeP/997nrrruwWCxceOGFtLW10dvby4IFC7zO30qCjY+PdyI4ud1vaGigt7dXtb2UJtzBGi3U1dVFRUUF06dPZ+nSpW6vn/DwcJKSkpwIU2tQc+LECcxms9p+LcnYarVSVlZGfHw8xcXFXjdpNDc3s3btWqZOncq///3vM6JTSwYmRUVF7Ny5c7yXc9ZgQpFuRkYGx48fV39ubGxk+vTpPh9n+fLlqvFGCJ6h1+vJyckhJyeHjz/+mDVr1nD77bdTXl5OSUkJW7du5YsvvsBgMLB48WI1P+xLE0d4eDhRUVHU19cTHx/PkiVLcDgcKsEdP35czQ+PLIB5C6vVSlVVFWaz2e/oOSIigilTpqhSLa0kTJL54OAgiYmJREREqGkdT+t0OBxs27aN559/ng0bNrBixYozRpnwySef8MYbb/D222+rn/OGG27gpZdeGu+lndGYUDldaY340UcfkZ6eTnFxMdu2bWPu3Lk+H6u+vp5vfvObX4r0QrDgzvZSURR6enpUE/i9e/dSU1PDtGnTnPLDrhoAvDUW92Yskrt8p2xCyMrKCnrHHQxP+ZV+GNK8XD60Jjqy/TosLIwTJ05w1113kZaWxpNPPhl0U3VPMJvNLF++nMHBQWw2G9/97ndZt26d29fv3LmTX//616Gc7jC+HDnd8PBwnn/+eS699FLsdjurV6/2i3CDgePHj7Nq1SpaWlrQ6/XcdtttrF27dlzWcrrgjqx0Oh0JCQl84xvf4Bvf+AYwrKktKSlhz549bN68mc7OTvLy8tT8sNFopLKyku985zsUFxd7jIx1Oh0xMTHExMQ45YelOY1MS2g9EaKjo6mvrycsLMxj7tlfyPl7HR0dTgY4sbGxxMbGqrsw7TqPHz/Ohg0bOHLkCN3d3fzwhz/klltuOe352sjISP75z38SGxuL1Wrlggsu4PLLL+fcc889reuYiJhQkW4wEWik29zcTHNzM0uWLMFkMlFUVMTrr78ekq95gN1u5+jRo/zrX//ihRdewGQyMWPGDPLz89VouKCgICCNrGyQaGxspLOzE4PB4GSwPpYfrbeQBjtTp071enQOiJTYXXfdxfTp07nmmms4evQopaWlbNy4kaysrIDX5Q/6+/u54IIL+O1vf8s554ycLB2CG3w5It0zCWlpaWrEFRcXR2FhIU1NTSHS9YCwsDDmz5/PW2+9xc9+9jN+8IMfOJnAP/HEE1RUVJCUlKQqJYqLi30afGm1WqmvrycmJobly5cTHh6u5iO7u7tpaGhwORbJl0kXtbW1o5QP3vzen/70J/7nf/6HJ554gosvvhidTscVV1zh1e+fCtjtdoqKiqiuruaOO+4IEW6QEIp0XeC6665j586ddHR0MG3aNNatW8ctt9zi9/Hq6+tZvnw5hw8fHhdZz0SCNIHXzqY7ceIEs2bNcjKBH6k7djgcNDQ00NraOubQSZkfNhqNbsciuZp00d3dTXl5OWlpacyYMcPrG8Hx48e58847yc7OZuPGjWfcNdLd3c1VV13Fc889x7x588Z7OWcLvhzNEWcient7ueiii3jooYe4+uqrff59XwsaX0ZIE/g9e/aoJvBms1k1gY+NjWXnzp3cf//9zJo1y6/mAU9jkWJjY+nq6mJgYMAngxqHw8GLL77I73//e5588km+8Y1vnDHKhJFYt24dMTExXwbPhGAhRLrjAavVyje/+U0uvfRS7rnnHr+OIc3GtQWNZ599NlTQGAODg4Ps2bOHX/7ylxw9epSZM2eiKApLlixRI2J/CVjCarXS1NTEsWPH1Dywt2ORGhoaWLNmDXl5eWzcuDHoPr+Bor29HYPBQGJiIgMDA1xyySXcf//9fPOb3xzvpZ0tCOV0TzcUReGWW26hsLDQb8IFUZWXX0ir1YrVaj1jo6EzCZGRkYSHh3PllVfy7rvvotfr6e7uVgeEvvrqq9TV1ZGenq6ScFFREcnJyV6dX5vNRnV1NQMDAyxbtozo6GinsUhdXV3U19c7jUVqbW2loKCAv/zlL2zdupUnn3ySr3/966f17+mtqqa5uZkbb7wRu92Ow+Hg+9//fohwg4RQpHuK8PHHH3PhhRcyf/58NZpav349K1as8PlYIwsav/rVr4K93C8lZJ5XpiVKS0sxGo0UFBSMMoHXwpfROdqxSP/1X//F7t27MZvNfOtb3+L888/n+uuvD7pUzRNCqprThlB6YSIg0IKGJwu+EASsVusoE3i9Xs/ixYspKCjggw8+YNWqVVx66aVetzbb7Xb+8Ic/8OKLL/LMM89QXFzMoUOH2LdvH2vWrDmtBjUj8e1vf5s1a9Zw8cUXj9saJihCpDtREEhBw5MFXwiuIU3gN23axPPPP8+CBQtoampSPQeWLl1KcXEx06ZNcxnx1tXVceeddzJ//nzWr18f9HHygSCkqjmlCOV0z1aMLGh8+OGH3H///T4fp7Gxkbfeeku14AvBO8icul6v5/PPP2fq1KnqhAdpAv+73/2OtrY2cnNzVSJeuHAhL7/8Mn/+85959tlnufDCC8+oXHxvby/XXHMNzzzzTIhwTzNCke4Zjs8//3xUQePhhx/2+Tjf/e53efDBBzGZTKEe+VMAu91ORUWFqh9+9913WbZsGS+++KLXErLThWCoakIYE6H0wpcZb775Jm+//Ta/+c1vAjYmycrKUmd0hYeHs2/fviCvdmLAnfnPeENRFG688UYmT57MM888M97LmcgIke6XGQ8++CB//vOfnVper776ar8s+LKysti3b98ZM2k2BAFvjfeDqaoJwSNCpBuCQDAi3RDpnnnYtWsXsbGxrFq1KmRHembALemOn1YlhLMSOp2OSy65hKKiIn73u9/5fZzu7m6++93vUlBQQGFhIbt37w7iKr98WL58uTorLYQzGyH1wpcMX/3qV/nqV7/q9+9/8sknTJ8+nba2Ni6++GIKCgpYvny5z8dZu3Ytl112GX//+9+xWCz09/f7vaYQQjibEIp0Q/AJ0ng7JSWFq666ir179/p8jJ6eHnbt2qU6t0VERJzWqQghhDCeCJFuCF6jr69Pnfra19fH+++/71dnXG1tLVOnTuXmm29m8eLF3HrrrfT19QV7uSGEcEYiRLoheI3W1lYuuOACFi5cyLJly7jiiiu47LLLfD6OzWbjwIED/PjHP+azzz4jJiaGxx9/3OfjVFRUsGjRIvURHx8fkkGFcOZDURRPjxBCCDqam5uVmTNnqj/v2rVLWbFiRUDHtNlsyrRp05T6+voAV3fm4Z133lHy8vKUnJwcZcOGDS5fs3LlSiU1NVUJDw9X0tPTlS1btpzmVYYwAm55NVRIC+G0IzU1lczMTCoqKsjPz+ejjz4K2OXqo48+Iicnh5kzZwZplWcG7HY7d9xxBx988AEZGRkUFxdz5ZVXjjpfL7/88jitMARfESLdEMYFzz33HNdffz0Wi4Xs7Gy2bt0a0PG2b9/OddddF6TVnTnYu3cvubm5ZGdnA7By5Up27NgRsmI8ixEi3RDGBYsWLQpaC7HFYuGNN95gw4YNAR3n6aefZsuWLeh0OubPn8/WrVu9tm88VWhqaiIzM1P9OSMjg5KSknFcUQiBIlRIC+GsxzvvvMOSJUuYNm2a38doampi06ZN7Nu3j8OHD2O329m+fXsQV+kfFBcdo2eip0MI3mOsNuAQQjjjodPptgPvKYrid45Cp9OlA3uAhUAP8DqwSVGU94OzSr/XdR7wqKIolw79/CCAoiiBhfUhjBtCkW4IZzV0Ot0k4GLg1UCOoyhKE/BroAFoBozjTbhDKAVm63S6WTqdLgJYCbwxzmsKIQCESDeEsxqKovQripKsKIoxkOPodLok4NvALGA6EKPT6W4IxhoDgaIoNmAN8B5QBvxVUZQj47uqEAJBqJAWQggC/weoUxSlHUCn070KfAXw3f8yyFAU5W3g7fFeRwjBQSjSDSEEgQbgXJ1ON0knKlXfQESWIYQQVIRIN4QQAEVRSoC/AweALxDfDf+9K0MIwQ1C6oUQQgghhNOI/wdlB0snKFmW/gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "# prepare some coordinates\n", "x, y, z = np.indices((8, 8, 8))\n", "\n", "# draw cuboids in the top left and bottom right corners, and a link between them\n", "cube1 = (x < 3) & (y < 3) & (z < 3)\n", "cube2 = (x >= 5) & (y >= 5) & (z >= 5)\n", "link = abs(x - y) + abs(y - z) + abs(z - x) <= 2\n", "\n", "# combine the objects into a single boolean array\n", "voxels = cube1 | cube2 | link\n", "\n", "# set the colors of each object\n", "colors = np.empty(voxels.shape, dtype=object)\n", "colors[link] = 'red'\n", "colors[cube1] = 'blue'\n", "colors[cube2] = 'green'\n", "\n", "# and plot everything\n", "fig = plt.figure()\n", "ax = fig.gca(projection='3d')\n", "ax.voxels(voxels, facecolors=colors, edgecolor='k')\n", "\n", "plt.show()\n", "\n" ] }, { "cell_type": "code", "execution_count": 199, "metadata": {}, "outputs": [], "source": [ "a = np.arange(60).reshape(3,4,5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Measuring performance of combinations" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ " 0%| | 0/11 [00:00" ] }, "metadata": { "needs_background": "dark" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import perfplot\n", "from itertools import chain, combinations\n", "from scipy.special import comb\n", "from scipy.linalg import pascal\n", "\n", "def C(n, k): # huge memory overload...\n", " if k == 0:\n", " return np.array([])\n", " if k == 1:\n", " return np.arange(1, n + 1)\n", " elif k == 2:\n", " return np.less.outer(np.arange(n), np.arange(n))\n", " else:\n", " x = C(n, k - 1)\n", " X = np.repeat(x[None, :, :], len(x), axis=0)\n", " Y = np.repeat(x[:, :, None], len(x), axis=2)\n", " return X & Y\n", "\n", "def C_indices(data):\n", " n, k = data\n", " return np.transpose(np.nonzero(C(n, k)))\n", "\n", "def comb_index(data):\n", " n, k = data\n", " count = comb(n, k, exact=True)\n", " index = np.fromiter(chain.from_iterable(combinations(range(n), k)),\n", " int, count=count * k)\n", " return index.reshape(-1, k)\n", "\n", "def stretch(a, k):\n", " l = a.sum()+len(a)*(-k)\n", " out = np.full(l, -1, dtype=int)\n", " out[0] = a[0]-1\n", " idx = (a-k).cumsum()[:-1]\n", " out[idx] = a[1:]-1-k\n", " return out.cumsum()\n", "\n", "def numpy_combinations(data):\n", " #print(data)\n", " n, k = data\n", " x = np.array([n])\n", " P = pascal(n).astype(int)\n", " C = []\n", " for b in range(k-1,-1,-1):\n", " x = stretch(x, b)\n", " r = P[b][x - b]\n", " C.append(np.repeat(x, r))\n", " return n - 1 - np.array(C).T\n", "\n", "def numpy_combinationsv(data):\n", " #print(data)\n", " n, k = data\n", " x = np.array([n])\n", " P = pascal(n).astype(int)\n", " C = []\n", " for b in range(k-1,-1,-1):\n", " x = stretch(x, b)\n", " r = P[b][x - b]\n", " C.append(np.repeat(x, r))\n", " return np.array(C).T\n", "\n", "def nump(data, i=0):\n", " n, k = data\n", " i = 0\n", " if k == 1:\n", " a = np.arange(i, i+n)\n", " return tuple([a[None, j:] for j in range(n)])\n", " template = nump((n-1, k-1), i+1)\n", " full = np.r_[np.repeat(np.arange(i, i+n-k+1),\n", " [t.shape[1] for t in template])[None, :],\n", " np.c_[template]]\n", " return tuple([full[:, j:] for j in np.r_[0, np.add.accumulate(\n", " [t.shape[1] for t in template[:-1]])]])\n", "\n", "def nump2(data, i=0):\n", " n, k = data\n", " a = np.ones((k, n-k+1), dtype=int)\n", " a[0] = np.arange(n-k+1)\n", " for j in range(1, k):\n", " reps = (n-k+j) - a[j-1]\n", " a = np.repeat(a, reps, axis=1)\n", " ind = np.add.accumulate(reps)\n", " a[j, ind[:-1]] = 1-reps[1:]\n", " a[j, 0] = j\n", " a[j] = np.add.accumulate(a[j])\n", " return a\n", "\n", "\n", "def build_args(k):\n", " return {'setup': lambda x: (k, x),\n", " 'kernels': [numpy_combinations, numpy_combinationsv, nump, nump2],\n", " 'n_range': [x for x in range(1, k)],\n", " 'xlabel': f'N',\n", " 'title': f'test of case C({k}, k)',\n", " 'show_progress': True,\n", " 'equality_check': False}\n", "\n", "outs = [perfplot.bench(**build_args(n)) for n in (12, 15, 17, 20)]\n", "fig = plt.figure(figsize=(20, 20))\n", "for i in range(len(outs)):\n", " ax = fig.add_subplot(2, 2, i + 1)\n", " ax.grid(True, which=\"both\")\n", " outs[i].plot()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0, 1, 2],\n", " [0, 1, 3],\n", " [0, 1, 4],\n", " [0, 2, 3],\n", " [0, 2, 4],\n", " [0, 3, 4],\n", " [1, 2, 3],\n", " [1, 2, 4],\n", " [1, 3, 4],\n", " [2, 3, 4]], dtype=int32)" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def stretch(a, k):\n", " l = a.sum()+len(a)*(-k)\n", " out = np.full(l, -1, dtype=int)\n", " out[0] = a[0]-1\n", " idx = (a-k).cumsum()[:-1]\n", " out[idx] = a[1:]-1-k\n", " return out.cumsum()\n", "\n", "def numpy_combinations(data):\n", " #print(data)\n", " n, k = data\n", " x = np.array([n])\n", " P = pascal(n).astype(int)\n", " C = []\n", " for b in range(k-1,-1,-1):\n", " x = stretch(x, b)\n", " r = P[b][x - b]\n", " C.append(np.repeat(x, r))\n", " return n - 1 - np.array(C).T\n", "\n", "numpy_combinations((5,3))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 4 }