diff --git a/moodys_challenge.ipynb b/moodys_challenge.ipynb index 7ee0e87..8dce1b1 100644 --- a/moodys_challenge.ipynb +++ b/moodys_challenge.ipynb @@ -161,15 +161,59 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "4ec880cd", + "execution_count": 1, + "id": "9fb55a85", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", - "import numpy as np\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "4ec880cd", + "metadata": {}, + "outputs": [ + { + "ename": "FileNotFoundError", + "evalue": "[Errno 2] No such file or directory: 'sp500_full.csv'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[2], line 3\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m# Port the market data to the notebook\u001b[39;00m\n\u001b[1;32m----> 3\u001b[0m time_series \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mlog(\u001b[43mpd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread_csv\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43msp500_full.csv\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mheader\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m[\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mto_numpy()\u001b[38;5;241m.\u001b[39msqueeze())\n", + "File \u001b[1;32mc:\\Users\\JackOberman\\anaconda3\\envs\\moodys\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1026\u001b[0m, in \u001b[0;36mread_csv\u001b[1;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)\u001b[0m\n\u001b[0;32m 1013\u001b[0m kwds_defaults \u001b[38;5;241m=\u001b[39m _refine_defaults_read(\n\u001b[0;32m 1014\u001b[0m dialect,\n\u001b[0;32m 1015\u001b[0m delimiter,\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 1022\u001b[0m dtype_backend\u001b[38;5;241m=\u001b[39mdtype_backend,\n\u001b[0;32m 1023\u001b[0m )\n\u001b[0;32m 1024\u001b[0m kwds\u001b[38;5;241m.\u001b[39mupdate(kwds_defaults)\n\u001b[1;32m-> 1026\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_read\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilepath_or_buffer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\JackOberman\\anaconda3\\envs\\moodys\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:620\u001b[0m, in \u001b[0;36m_read\u001b[1;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[0;32m 617\u001b[0m _validate_names(kwds\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnames\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m))\n\u001b[0;32m 619\u001b[0m \u001b[38;5;66;03m# Create the parser.\u001b[39;00m\n\u001b[1;32m--> 620\u001b[0m parser \u001b[38;5;241m=\u001b[39m \u001b[43mTextFileReader\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilepath_or_buffer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwds\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 622\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m chunksize \u001b[38;5;129;01mor\u001b[39;00m iterator:\n\u001b[0;32m 623\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m parser\n", + "File \u001b[1;32mc:\\Users\\JackOberman\\anaconda3\\envs\\moodys\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1620\u001b[0m, in \u001b[0;36mTextFileReader.__init__\u001b[1;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[0;32m 1617\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhas_index_names\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m kwds[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhas_index_names\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[0;32m 1619\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles: IOHandles \u001b[38;5;241m|\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m-> 1620\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_engine \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_make_engine\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mengine\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\JackOberman\\anaconda3\\envs\\moodys\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1880\u001b[0m, in \u001b[0;36mTextFileReader._make_engine\u001b[1;34m(self, f, engine)\u001b[0m\n\u001b[0;32m 1878\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m mode:\n\u001b[0;32m 1879\u001b[0m mode \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m-> 1880\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles \u001b[38;5;241m=\u001b[39m \u001b[43mget_handle\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m 1881\u001b[0m \u001b[43m \u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1882\u001b[0m \u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1883\u001b[0m \u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mencoding\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1884\u001b[0m \u001b[43m \u001b[49m\u001b[43mcompression\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcompression\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1885\u001b[0m \u001b[43m \u001b[49m\u001b[43mmemory_map\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mmemory_map\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1886\u001b[0m \u001b[43m \u001b[49m\u001b[43mis_text\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mis_text\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1887\u001b[0m \u001b[43m \u001b[49m\u001b[43merrors\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mencoding_errors\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mstrict\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1888\u001b[0m \u001b[43m \u001b[49m\u001b[43mstorage_options\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mstorage_options\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 1889\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 1890\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m 1891\u001b[0m f \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles\u001b[38;5;241m.\u001b[39mhandle\n", + "File \u001b[1;32mc:\\Users\\JackOberman\\anaconda3\\envs\\moodys\\Lib\\site-packages\\pandas\\io\\common.py:873\u001b[0m, in \u001b[0;36mget_handle\u001b[1;34m(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)\u001b[0m\n\u001b[0;32m 868\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(handle, \u001b[38;5;28mstr\u001b[39m):\n\u001b[0;32m 869\u001b[0m \u001b[38;5;66;03m# Check whether the filename is to be opened in binary mode.\u001b[39;00m\n\u001b[0;32m 870\u001b[0m \u001b[38;5;66;03m# Binary mode does not support 'encoding' and 'newline'.\u001b[39;00m\n\u001b[0;32m 871\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ioargs\u001b[38;5;241m.\u001b[39mencoding \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m ioargs\u001b[38;5;241m.\u001b[39mmode:\n\u001b[0;32m 872\u001b[0m \u001b[38;5;66;03m# Encoding\u001b[39;00m\n\u001b[1;32m--> 873\u001b[0m handle \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[0;32m 874\u001b[0m \u001b[43m \u001b[49m\u001b[43mhandle\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 875\u001b[0m \u001b[43m \u001b[49m\u001b[43mioargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 876\u001b[0m \u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mioargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 877\u001b[0m \u001b[43m \u001b[49m\u001b[43merrors\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43merrors\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m 878\u001b[0m \u001b[43m \u001b[49m\u001b[43mnewline\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[0;32m 879\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 880\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m 881\u001b[0m \u001b[38;5;66;03m# Binary mode\u001b[39;00m\n\u001b[0;32m 882\u001b[0m handle \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mopen\u001b[39m(handle, ioargs\u001b[38;5;241m.\u001b[39mmode)\n", + "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'sp500_full.csv'" + ] + } + ], + "source": [ + "# Port the market data to the notebook\n", "\n", - "time_series = np.log(pd.read_csv(\"SP500.csv\", header=None).to_numpy().squeeze())" + "time_series = np.log(pd.read_csv(\"sp500_full.csv\", header=None)[1].to_numpy().squeeze())" + ] + }, + { + "cell_type": "code", + "execution_count": 95, + "id": "f433c5cb", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the hyperparameters\n", + "\n", + "dimension = 4 \n", + "time_delay = 30\n", + "window = 40\n", + "epsilon = 0.0 # will be reset\n", + "precision_qubits = 4\n", + "stride = 1" ] }, { @@ -210,11 +254,73 @@ { "cell_type": "code", "execution_count": null, + "id": "def47fab", + "metadata": {}, + "outputs": [], + "source": [ + "# from gtda.time_series import TakensEmbedding\n", + "from gtda.time_series import SlidingWindow" + ] + }, + { + "cell_type": "code", + "execution_count": 96, "id": "7b701be6", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "# Given a time series and the corresponding choices for dimension and time delay it produces the emebeded data\n", + "def embed_series(time_series: np.ndarray, dimension: int, time_delay: int):\n", + " l = len(time_series)\n", + " num_points = l - (dimension - 1) * time_delay\n", + " embedded = np.array(\n", + " [[time_series[z + m * time_delay] for m in range(dimension)]\n", + " for z in range(num_points)],\n", + " dtype=time_series.dtype\n", + " )\n", + " return embedded[np.newaxis, :]\n", + "\n", + "takens_embedding = embed_series(time_series, time_delay, dimension)" + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "id": "5f3cc31f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.099999296028088\n" + ] + } + ], + "source": [ + "# Based on the pairwise distance of the embeded data points we first produce a reasonable choice for epsilone\n", + "\n", + "estim_epsilon = 0\n", + "for a, b in zip(takens_embedding[0][:-1], takens_embedding[0][1:]):\n", + " estim_epsilon += np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))\n", + "estim_epsilon /= (takens_embedding.shape[1] - 1)\n", + "epsilon = estim_epsilon + 1.1\n", + "\n", + "# Retrieve choice of epsilon\n", + "print(epsilon)" + ] + }, + { + "cell_type": "code", + "execution_count": 99, + "id": "76cd85eb", + "metadata": {}, + "outputs": [], + "source": [ + "# Divide our embeded data into time-dependent windows\n", + "\n", + "WindowSlider = SlidingWindow(size=window, stride=stride)\n", + "embeddings_windows = WindowSlider.fit_transform(takens_embedding.squeeze())" ] }, { @@ -274,11 +380,28 @@ { "cell_type": "code", "execution_count": null, + "id": "2baf3027", + "metadata": {}, + "outputs": [], + "source": [ + "import gudhi" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "id": "026d6b61", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "# Generate the simplex complexes for each time-dependent window\n", + "\n", + "simplex_trees = []\n", + "for i in range(len(embeddings_windows)):\n", + " rips_complex = gudhi.RipsComplex(points = embeddings_windows[i], max_edge_length = epsilon)\n", + " simplex_tree = rips_complex.create_simplex_tree(max_dimension=dimension)\n", + " simplex_trees.append(simplex_tree)\n", + " break" ] }, { @@ -337,12 +460,69 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "246e41ce", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "# Given a simplex tree and a dimention k, compute kth boundry operator for the given tree\n", + "def get_boundary_operator(simplex_tree, k):\n", + " D_k = []\n", + "\n", + " for simplex_k in simplex_tree[k]:\n", + " S_k_sub_one = simplex_tree[k - 1]\n", + " d_k_s_k = np.zeros(len(S_k_sub_one))\n", + " \n", + " t = 0\n", + " for j, s_k_sub_one in enumerate(S_k_sub_one):\n", + " if set(s_k_sub_one).issubset(set(simplex_k)):\n", + " d_k_s_k[j] += (-1)**t\n", + " t += 1\n", + " D_k.append(d_k_s_k)\n", + "\n", + " return np.stack(D_k, axis=1) if len(D_k) > 1 else D_k[0][:, np.newaxis]\n", + "\n", + "# Given a simplex tree and a window indice it computes all boundry operators from the window indice - 1 to 0\n", + "def get_boundary_operators(simplex_tree, window):\n", + " D = []\n", + " for k in range(window - 1, 0, -1):\n", + " D_k = get_boundary_operator(simplex_tree, k)\n", + " D.append(D_k)\n", + " return D\n", + " \n", + "# Given a simplex tree and a window indice it returns a filtered tree\n", + "def create_filtered_simplex_tree(simplex_tree, window):\n", + " filtered_tree = [[] for i in range(window)]\n", + " for filtered_value in simplex_tree.get_filtration():\n", + " arr, weight = filtered_value\n", + " filtered_tree[len(arr) - 1].append(tuple(arr))\n", + " return filtered_tree\n", + "\n", + "# Generates the filtered tree corresponding to the toy example given in the notebook\n", + "def create_fake_filtered_simplex_tree(simplex_tree, window):\n", + " filtered_tree = [[] for i in range(window)]\n", + " for filtered_value in simplex_tree:\n", + " arr = filtered_value\n", + " filtered_tree[len(arr) - 1].append(tuple(arr))\n", + " return filtered_tree" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "75ef4871", + "metadata": {}, + "outputs": [], + "source": [ + "# Generates the appropriate filtered trees corresponding to each of the time-dependent window's simplex tree\n", + "filtered_trees = []\n", + "for simplex_tree in simplex_trees:\n", + " filtered_tree = create_filtered_simplex_tree(simplex_tree, window)\n", + " filtered_trees.append(filtered_tree)\n", + "\n", + "# Generates the \"fake tree\" - the one given as a toy example\n", + "fake_tree = [[1], [2], [3], [4], [5],[1, 2], [1, 3], [2, 3], [3, 4], [3, 5], [4, 5], [1, 2, 3]]\n", + "fake_tree = create_fake_filtered_simplex_tree(fake_tree, 3)" ] }, { @@ -423,12 +603,52 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "id": "8f818f78", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "# Given a filtered tree and a dimension k it returns the kth laplacian corresponding to the given simplex tree\n", + "def get_laplacian(filtered_tree, k) -> np.matrix:\n", + " # Approach 1\n", + " boundry_k = get_boundary_operator(filtered_tree, k)\n", + " boundry_k_1 = get_boundary_operator(filtered_tree, k + 1)\n", + "\n", + " return (boundry_k.conj().T @ boundry_k) + (boundry_k_1 @ boundry_k_1.conj().T)\n", + "\n", + "# Given a filtered tree and a dimension k it returnes the kth laplacian padded as shown, as well as the number of qubits and max eigenvalue\n", + "def get_padded_laplacian(filtered_tree, k) -> (np.matrix, int, int):\n", + " laplacian = get_laplacian(filtered_tree, k)\n", + " # Get #S_k and \n", + " s_k, _ = laplacian.shape\n", + " q = int(np.ceil(np.log2(s_k)))\n", + "\n", + " # Compute the maximum possible eigenvalue\n", + " max_eigen = float('-inf')\n", + " for i in range(0, s_k):\n", + " R_i = np.sum(np.abs(laplacian[i,:])) - np.abs(laplacian[i,i])\n", + " \n", + " cur_eigen_max = laplacian[i,i] + R_i\n", + " max_eigen = cur_eigen_max if cur_eigen_max > max_eigen else max_eigen\n", + " new_dim = 2**q\n", + "\n", + " padded_laplacian = np.pad(laplacian, ((0, new_dim - s_k),(0, new_dim - s_k)))\n", + " padded_laplacian[s_k:,s_k:] = (max_eigen / 2) * np.eye(new_dim - s_k)\n", + " \n", + " return (padded_laplacian, q, max_eigen)\n", + "\n", + "# Given a filtered tree and the dimension k it returns the H matrix corresponding to the kth laplacian, as well as the number of qubits\n", + "def get_h_mat(filtered_tree, k, precision=0.999) -> (np.matrix, int):\n", + " delta = precision * 2 * np.pi\n", + " padded_laplacian, q, max_eigen = get_padded_laplacian(filtered_tree, k)\n", + "\n", + " return ((delta / max_eigen) * padded_laplacian, q)\n", + "\n", + "# Computes and saves the H matrices and qubit counts for each window retrieved from our embedings\n", + "H_mats = []\n", + "for filtered_tree in filtered_trees:\n", + " H_mat_real, q_real = get_h_mat(filtered_tree, 1)\n", + " H_mats.append((H_mat_real, q_real))" ] }, { @@ -499,12 +719,154 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "id": "848ef40c", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "# Import qiskit related libraries\n", + "\n", + "from qiskit import QuantumCircuit, transpile\n", + "from qiskit_aer import Aer\n", + "from qiskit.circuit.library import PhaseEstimation\n", + "from qiskit.circuit.library import UnitaryGate\n", + "from scipy.linalg import expm\n", + "from qiskit.visualization import plot_histogram" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "d7d66c4b", + "metadata": {}, + "outputs": [], + "source": [ + "# Given the H matrix corresponding to some window and the number of qubits needed it approximates the kth betti number applying QPE\n", + "def get_betti(H_mat, q):\n", + " H_mat = H_mat_real\n", + " q = q_real\n", + "\n", + " U = expm(1j * H_mat)\n", + " unitaryGate = UnitaryGate(U)\n", + "\n", + " qc = QuantumCircuit(precision_qubits + 2*q, precision_qubits)\n", + "\n", + " qc.h(range(precision_qubits, precision_qubits + q))\n", + "\n", + " for i in range(precision_qubits, precision_qubits+q):\n", + " qc.cx(i, i+q)\n", + "\n", + " pe = PhaseEstimation(num_evaluation_qubits=precision_qubits, unitary=unitaryGate)\n", + "\n", + " pe = qc.compose(pe)\n", + "\n", + " pe.measure(range(precision_qubits), range(precision_qubits))\n", + "\n", + " backend = Aer.get_backend('aer_simulator')\n", + "\n", + " qc_transpiled = transpile(pe, backend)\n", + "\n", + " shots = 2048\n", + " job_sim = backend.run(qc_transpiled, shots=shots)\n", + "\n", + " # # Grab the results from the job.\n", + " result_sim = job_sim.result().get_counts()\n", + "\n", + " num_zeroes = result_sim['0'*precision_qubits] if '0'*precision_qubits in result_sim else 0\n", + " prob_zero = num_zeroes / shots\n", + " betti = 2**q * prob_zero\n", + " betti_num = round(betti)\n", + "\n", + " dims = H_mat.shape[0]\n", + " rank = np.linalg.matrix_rank(H_mat)\n", + " gt_betti = dims - rank\n", + "\n", + " return prob_zero, betti_num, gt_betti" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "e6aa6a77", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(32, 32)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 0%| | 0/1 [00:00" ] }, + { + "cell_type": "markdown", + "id": "736db30d", + "metadata": {}, + "source": [ + "### RESPONSES:\n", + "1. The reason for measuring all-zero state is that it allows us to find the dimention of the kernel of the kth laplacian, which is then used to find the kth betty number. More precisely, the probabilty of measuring the all-zero state is roughly equal to the proportion of zero eigenvalues to all eigenvalues of the kth laplacian. Thus multiplying this probability, times the rows of the kth laplacian (number of possible eigenvalues) which in our case will be $2^q$. We can then use the number of zero eigenvalues, which we just calculated as a rough estimate of the dimention of the kernel of the kth laplacian which rounded to the nearest integer would be our estimate for the kth betty number. \n", + "2. When considering the density matrix, a maximally mixed state has a uniform density. The second matrix has density operators that describe the distribution of the eigenvalues of the combinatorial laplacian. \n", + "3. As we found through testing, one parameter that seemed to affect the accuracy of our estimation before rounding was the number of preicsion qubits used in our circuit. Furthemore, the number of shots taken to perform our experiment also stabelized our results (in a law of large numbers kind of fashion - the probability of the zero state should approach the ideal proportion of zero-eigenvalues). Lastly, we found that tuning the delta parameter that is multiplied to the kth laplacian to compute the $H_0$ matrix also behaved as a threshold for the probabilites of the eigenvalues." + ] + }, { "cell_type": "markdown", "id": "1cbe5cf2", @@ -560,7 +933,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 66, "id": "5ef68f98", "metadata": {}, "outputs": [], @@ -575,12 +948,100 @@ " dim: the dimension on which the Betti number is calculated\n", " '''\n", " result = ripser(point_cloud, maxdim=dim)\n", + " # print(result)\n", " diagrams = result[\"dgms\"]\n", + " # print(diagrams)\n", " return len(\n", " [interval for interval in diagrams[dim] if interval[0] < epsilon < interval[1]]\n", - " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 122, + "id": "b7a9a34d", + "metadata": {}, + "outputs": [], + "source": [ + "# We now scale up to computing our epsilon-dependent betti curves for each of the windows\n", "\n", - "# write your code here" + "epsilon_vals = np.linspace(0.01, 1.0, 25)\n", + "betti_curves = []\n", + "\n", + "for window in embeddings_windows:\n", + " window_bettis = []\n", + " for epsilon in epsilon_vals:\n", + " window_bettis.append(classical_betti_solver(window, epsilon, 0))\n", + " betti_curves.append(window_bettis)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f89f1fbe", + "metadata": {}, + "outputs": [], + "source": [ + "# Given a vector it computes their L^p norm\n", + "def lp_norm(x, p=2):\n", + " return np.sum(np.abs(x)**p)**(1/p)\n", + "\n", + "# Given the window-order epsilon-dependent betti curves, it computes their succesive distance using the L^p norm\n", + "def pairwise_distances(betti_curves, p=2):\n", + " distances = []\n", + " for i in range(len(betti_curves) - 1):\n", + " distance = lp_norm(np.subtract(betti_curves[i+1], betti_curves[i]), p)\n", + " distances.append(distance)\n", + " return np.array(distances) " + ] + }, + { + "cell_type": "code", + "execution_count": 193, + "id": "b86935a5", + "metadata": {}, + "outputs": [], + "source": [ + "# We now compute the pairwise succesive distances in the betti curves and filter only those that are within a certain threshold\n", + "distances = pairwise_distances(betti_curves, p=4)\n", + "\n", + "crash_indices = np.where(distances > 0.99)[0] \n", + "crash_vals = distances[crash_indices]" + ] + }, + { + "cell_type": "code", + "execution_count": 198, + "id": "8612aa08", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 198, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnd0lEQVR4nO3dd3gU1d4H8O+mbRJIIZAKAYJAEBBCkRiaeIkERAQrclEQFRXhKmLFV0G9auwXC4IFxEaxACoCioFQlCIl0nsJJYWWbBLSd94/JtmZ2Z3d7Cbb8/08zz575syZ2bNLyP5yqkYQBAFEREREbszH1RUgIiIiqgsDFiIiInJ7DFiIiIjI7TFgISIiIrfHgIWIiIjcHgMWIiIicnsMWIiIiMjtMWAhIiIit+fn6grYg16vx7lz5xASEgKNRuPq6hAREZEVBEFAUVER4uLi4ONjuQ3FKwKWc+fOIT4+3tXVICIiono4ffo0WrVqZbGMVwQsISEhAMQ3HBoa6uLaEBERkTV0Oh3i4+MN3+OWeEXAUtsNFBoayoCFiIjIw1gznIODbomIiMjtMWAhIiIit8eAhYiIiNweAxYiIiJyewxYiIiIyO3ZFLCkp6fj2muvRUhICKKiojBq1CgcOnTI4jULFiyARqNRPAIDAxVlBEHAjBkzEBsbi6CgIKSmpuLIkSO2vxsiIiLySjYFLOvXr8fkyZOxZcsWrFmzBpWVlRgyZAhKSkosXhcaGoqcnBzD49SpU4rzb731Fj744APMnTsXW7duRZMmTZCWloaysjLb3xERERF5HZvWYVm9erXieMGCBYiKisKOHTswcOBAs9dpNBrExMSonhMEAbNmzcILL7yAkSNHAgC++uorREdHY/ny5bj77rttqSIRERF5oQaNYSksLAQAREREWCxXXFyMNm3aID4+HiNHjsS+ffsM506cOIHc3FykpqYa8sLCwpCcnIzNmzer3q+8vBw6nU7xICIiIu9V74BFr9dj6tSp6NevH7p27Wq2XGJiIubPn4+ffvoJ33zzDfR6Pfr27YszZ84AAHJzcwEA0dHRiuuio6MN54ylp6cjLCzM8OA+QkRERN6t3gHL5MmTsXfvXixevNhiuZSUFIwbNw5JSUm4/vrrsXTpUkRGRuKTTz6p70tj+vTpKCwsNDxOnz5d73sRERGR+6vXXkJTpkzBihUrsGHDhjp3VzTm7++PHj164OjRowBgGNuSl5eH2NhYQ7m8vDwkJSWp3kOr1UKr1dan6kREROSBbGphEQQBU6ZMwbJly7B27VokJCTY/ILV1dXYs2ePIThJSEhATEwMMjIyDGV0Oh22bt2KlJQUm+9PREREdlRZCvz5PpDzj0urYVMLy+TJk7Fw4UL89NNPCAkJMYwxCQsLQ1BQEABg3LhxaNmyJdLT0wEAr7zyCq677jq0b98eBQUFePvtt3Hq1Ck8+OCDAMQZRFOnTsWrr76KDh06ICEhAS+++CLi4uIwatQoO75VIiIistmBFcCaGWL6pUKXVcOmgGXOnDkAgEGDBinyv/jiC9x3330AgOzsbPj4SA03ly9fxsSJE5Gbm4tmzZqhV69e+Ouvv9C5c2dDmWeeeQYlJSV46KGHUFBQgP79+2P16tUmC8wRERGRk106JqUrSwH/IJdUQyMIguCSV7YjnU6HsLAwFBYWIjQ01NXVISIi8nz5B4CzO4HcPcDWOUDbAcC4nwEf++3qY8v3d70G3RIREZEXK7kAfHydMq95e7sGK7bi5odERESklH/ANC8ny+nVkGPAQkREREol+aZ5Qc2cXw8ZBixERESkdPmkaZ7WtWNEGbAQERGRUsYrpnk+rh32yoCFiIiIAEEAqqvMn2/Wxnl1UcGAhYiIiIDFY4F3E4Hcvern+09zbn2McFozERFRY1d8Hjj0q5je+I743DQGGLMQOLgSGPg04O/axVwZsBARETV2SydK6X3LxOfQOKBlL/HhBtglRERE1JhVlgLH15nmB7rXyvEMWIiIiBqzX6aq52tDnFqNujBgISIiaqwunwJ2L1Y/5+Pv3LrUgQELERFRY/XN7ebPHV7tvHpYgQELERFRY1RRAlw8Ih2PmmNUQOPU6tSFAQsREVFjtMUoQEn6N9DnIen49s+dW586MGAhIiJqjApOmeb1fUx8DgwDOg51bn3qwHVYiIiIGqPyItO88Hhg0mYgKBzwca82DQYsREREjc3FY9ICcQDQ52EpHd3Z+fWxAgMWIiKixuLiMeCvD4EdX0h5t3wIdLvbdXWyEgMWIiKixkCvBz7saZrf5TbAL8D59bGRe3VQERERkWN8d696vrapc+tRTwxYiIiIGoODK0zzUl9yejXqi11CRERE3k5fbZo3/D3g2gecX5d6YgsLERGRt9OdM83zdf9xK3IMWIiIiLzduZ2meTHXOL8eDcAuISIiIm9XcFp8juoMtOwFhLcB4pJcWiVbMWAhIiLydhXF4nNkIjDyI9fWpZ7YJUREROTtapfhD2/t2no0AAMWIiIib1dWID5rQ1xajYZgwEJEROTNDqwAdn3j6lo0GMewEBEReSO9Hsh4GfhzlpRXWeay6jQUW1iIiIi80T8LlcEKAHQZ5Yqa2IVNAUt6ejquvfZahISEICoqCqNGjcKhQ4csXvPZZ59hwIABaNasGZo1a4bU1FRs27ZNUea+++6DRqNRPIYOHWr7uyEiIiLRWZW1V/yDnV8PO7EpYFm/fj0mT56MLVu2YM2aNaisrMSQIUNQUlJi9prMzEyMGTMG69atw+bNmxEfH48hQ4bg7NmzinJDhw5FTk6O4bFo0aL6vSMiIiICLh41zQuOcH497EQjCIJQ34vPnz+PqKgorF+/HgMHDrTqmurqajRr1gwfffQRxo0bB0BsYSkoKMDy5cvrVQ+dToewsDAUFhYiNDS0XvcgIiLyGr9MBXZ8ocy780u36xKy5fu7QWNYCgsLAQAREdZHbFeuXEFlZaXJNZmZmYiKikJiYiImTZqEixcvmr1HeXk5dDqd4kFEREQ1jIMVAOhwo/PrYUf1Dlj0ej2mTp2Kfv36oWvXrlZf9+yzzyIuLg6pqamGvKFDh+Krr75CRkYG3nzzTaxfvx7Dhg1DdbXK7pIQx9KEhYUZHvHx8fV9G0RERN6l8Kx6vgePXwEa0CU0adIkrFq1Cps2bUKrVq2suuaNN97AW2+9hczMTHTr1s1suePHj+Oqq67CH3/8gcGDB5ucLy8vR3l5ueFYp9MhPj6eXUJEREQ7vgR+ecw0/6VC59elDg7vEpoyZQpWrFiBdevWWR2svPPOO3jjjTfw+++/WwxWAKBdu3Zo0aIFjh5VGTAEQKvVIjQ0VPEgIiIiAIFhymMfP6DPw66pix3ZtHCcIAj4z3/+g2XLliEzMxMJCQlWXffWW2/htddew2+//YbevXvXWf7MmTO4ePEiYmNjbakeERERfT9eefzUEY+eHVTLphaWyZMn45tvvsHChQsREhKC3Nxc5ObmorS01FBm3LhxmD59uuH4zTffxIsvvoj58+ejbdu2hmuKi8WdI4uLi/H0009jy5YtOHnyJDIyMjBy5Ei0b98eaWlpdnqbREREjUDJBeVxaEsgqJlr6mJnNgUsc+bMQWFhIQYNGoTY2FjDY8mSJYYy2dnZyMnJUVxTUVGBO+64Q3HNO++8AwDw9fXF7t27ccstt6Bjx4544IEH0KtXL2zcuBFardZOb5OIiKgR+KCnlA6JAx7bBWg0rquPHTVoHRZ3wXVYiIio0asqB16Nko6nHQRC3XtohdPWYSEiIiI3kbdXeezmwYqtGLAQERF5g82zpfS4n11XDwdhwEJEROTpKsuAvT9Kx+2ud11dHIQBCxERkac7sUFKXz3CdfVwIAYsREREnqaqAtj9HVCUJx4vvFM616qPa+rkYAxYiIiIPM3Gd4ClE4F5Khsaxic7vz5OwICFiIjI0+yvGVRbcAo48IvyXGsGLEREROQOhGopveQeKd3lVufXxUkYsBAREXmaihL1/HaDnFoNZ2LAQkRE5GlKL6vn+wU6tx5OxICFiIjI01ReUc/3DXBuPZyIAQsREZG3YAsLERERuT0/trAQERGRI+XuARbcDBSctlyuTGf+HFtYiIiIyKHm9gdObgRmdbVc7pvbzJ9jwEJEREQOIwjWl9WdM3+Og26JiIjIIfTVwMvh0nGb/pbLWwpu/IPsUiV3xICFiIjIlfYtUx4HNDFftroKKMoR001jTM8zYCEiIiKHqCxVHuurzJctKwBQ08Ii6E3P+wfbq1ZuhwELERGRK108qjy2FLAczxSfA8OAknzT8xx0S0RERHYnCMCfs5R5lgKWHx8Qn8sKgZhuUn5MNyD+OnYJERERkQOseMI0z1zAYjzYtu9jUvqh9cD9qwGNxn51czN+rq4AERFRo7XjC9M844BFlwP8NBmIulqZ3/kW4O/rgNbJgI/3tz8wYCEiInIHY38Evr0dqK5U5m/+CDiWIT5qpb4E+GmBB35zahVdyftDMiIiIndUVQH41LQbPLwR8PEV0/pqZbmS86bX9pvq0Kq5IwYsRERErnD+oNT9E90V8PUX03qjFpbqCtNrvXisijkMWIiIiJytuhL4ZIB07OMjtbYYj2Ex7iJqpBiwEBEROduxtaZ5PjUtLNVGAYtxAHP/746pk5tjwEJERORsC++S0pGdxGfDGBbjFhajLqEmLRxXLzfGgIWIiMiZKsuUx62vE59rx7AUnQPKi6Xzxl1ClvYa8mIMWIiIiJwpM115fE1Na4uPbKWRje9KaeOAxYv3C7KE67AQERG5yiN/AjFdxbQ8YCk4JaV155TXsIWlbunp6bj22msREhKCqKgojBo1CocOHarzuu+//x6dOnVCYGAgrrnmGqxcuVJxXhAEzJgxA7GxsQgKCkJqaiqOHDli2zshIiJyN9VVpi0ktXsHXTtRClYAZcAiX4a/MFt5fe1Yl0bGpoBl/fr1mDx5MrZs2YI1a9agsrISQ4YMQUlJidlr/vrrL4wZMwYPPPAAdu3ahVGjRmHUqFHYu3evocxbb72FDz74AHPnzsXWrVvRpEkTpKWloayszOx9iYiI3NqeH4D/Ngfm9gcqrgBLHwZ2fQv417SQtOypLF87hgUANDVfz8b7BzViGkGo/6dx/vx5REVFYf369Rg4cKBqmdGjR6OkpAQrVqww5F133XVISkrC3LlzIQgC4uLi8OSTT+Kpp54CABQWFiI6OhoLFizA3XffXWc9dDodwsLCUFhYiNDQ0Pq+HSIiIvt5KUxK9xwP7PxSef6JfUBYK+m4OB94p4OY7noHcMc8oKoceDVKzGvVBxg5G4js6Nh6O5Et398NGnRbWFgIAIiIiDBbZvPmzUhNTVXkpaWlYfPmzQCAEydOIDc3V1EmLCwMycnJhjLGysvLodPpFA8iIiK3cVA59MEkWAGAkDjlsbxLCDVtCRWyHowJq7wqWLFVvQMWvV6PqVOnol+/fujatavZcrm5uYiOjlbkRUdHIzc313C+Ns9cGWPp6ekICwszPOLj4+v7NoiIiOxv/Rt1lzHeYVkesNSuxXLpuJTn27jnydQ7YJk8eTL27t2LxYsX27M+Vpk+fToKCwsNj9OnTzu9DkRE5OVsHTGxeTaw9VOgvAjI+cdy2b7/Mc1TBCw1GyBufM+2OnixeoVrU6ZMwYoVK7Bhwwa0atXKYtmYmBjk5eUp8vLy8hATE2M4X5sXGxurKJOUlKR6T61WC61WW5+qExER1e3SceCDHmL6ycNASLTl8v8sBn57Xkz7BdR9//hk0zz5oFtBD5zYCBz61br6NgI2tbAIgoApU6Zg2bJlWLt2LRISEuq8JiUlBRkZGYq8NWvWICUlBQCQkJCAmJgYRRmdToetW7cayhARETlNyQXgoz7S8aqn675m2SNSeusndZf3UWkvkOf5+gNf3iwdxybVfU8vZ1MLy+TJk7Fw4UL89NNPCAkJMYwxCQsLQ1BQEABg3LhxaNmyJdLTxZX8Hn/8cVx//fV49913MXz4cCxevBjbt2/Hp59+CgDQaDSYOnUqXn31VXTo0AEJCQl48cUXERcXh1GjRtnxrRIREVlh5VOAXrZ2ijbEcnm9HoZBsgCQv7/u12je3jRPoxGnMwt6INSo96Ig27R8I2NTC8ucOXNQWFiIQYMGITY21vBYsmSJoUx2djZycnIMx3379sXChQvx6aefonv37vjhhx+wfPlyxUDdZ555Bv/5z3/w0EMP4dprr0VxcTFWr16NwMBAO7xFIiIiG+xbpjze9Y3l8hvetv01WnRQz+8/TXw23gAxMMy0bCNjUwuLNUu2ZGZmmuTdeeeduPPOO81eo9Fo8Morr+CVV16xpTpERER1O/03kPUtMHgGEGx+GQ4AwNE/bL+/uQG27W8Ejq5R5qW9DsT1MH+v2m4hvdHquL5WjIvxco17jhQREXm/eTXrfGl8gJstzLqpLAW+uV06bnUtcOZvy+NH9HrzA2OH/BdIGgP8cL943GEIkDLZcl1rpy6f2Gi5XCPE3ZqJiKhxuHjU8vmfpiiPBzwpPmssfFV+d6/5cwFNgC63Sce32jAY96Lxfnpcop8BCxERea/Sy1K6SaTlsnt/UB771EwzPrcTOL5eXCZfThCAg9K2M2jVR3k+oKk4kPbRLcDEtXV3R8lf0xi7hBiwEBGRF7soWym2uhwoOA1UVSjLVFcBObuVeaO/Va6L8tUtYndRWaGUN+9G5TU3GQ2+9Q8Wn6OuBlr2sq6+atOdAeCWD6273osxYCEiIu91+YSUPrUZmNUVWDxGyisrFBeI+2SAlPfoFuDqm5UBCwCc3Khcn+XM38rzcUnAdbIxKn71WOBUbfn9x7JMd3ZuhBiwEBGR99KdldJXLojP8plAp7cBhUZrnERdLT6rdc8U1+xxJ29pAYD7ajY7lAccGo3t9VVrYdFa3sW4sWDAQkRE3mvNDPX82m6h8weV+Te8IKXLdebvqzunPA6t2Xm5oWNN1IKkgOCG3dNLMGAhIqLGp+gcsOcH4PcXlPktZCvQtr7O/PV/f6489hdXe0eHIeKzbz33uzNuYYnsJN27kWPAQkREjc+ub4EfHzDND24upQOaAP2mql8vD1jC20gzkOL7ABNWAY/tql+9jMewpExRL9cIceE4IiLyTpVl5s9teEs9P+Iq5XFQeN2vM2U74OMrHbfpW/c15hi3sJibNdQIsYWFiIi802vRtpVv0w8Ia6nM81cZP1ItWza/7QDAz45rpBiPYTGeqdSIMXQjIiLvY8XedwYDnwH+9X/q5/xUNuEtPCOl7/jCtnrVhS0sZrGFhYiIvM8e2aq1UZ0tl73hefPnmrQwzfsgSUo3rWP1XFsZj2FhwGLAgIWIiLzP0geldNrrlstaWi8l4Xog8mr71MkabGExiwELERF5t7YD6i5jjrYpMHmL/epSF5MxLAxYajFgISIi71JdJaVvn2f5S79pjHX3nLrHNG/CKtvqZQ22sJjFgIWIiLyHIADLJ4lpHz+gy63K8+1ukNL3LAUm/WndfcNbA0n3KPNCYutfT3M4hsUsfhJEROQ9fnwQ2Fsz4DYsXlofpcc9wIUjwNjvxQXjmrcH2g+27d6BYZaP7cGkhYXTmmsxYCEiIu9QfF4KVgCgOE9Kj5wtpe/6qn73N959OTC8fvexxDhAYQuLAbuEiIjIO3w1Unl8z4/2vb98TZY2/QEfB3yFGgcoHHRrwE+CiIg8X1U5kL9fOn4+x/67HMtXtLVmyf764BgWs/hJEBGR57twGEDN6rYzCyyvrVJf8hYWtSX77YFdQmaxS4iIiDybIABz+4vplr0cE6wAyjEs9m69qcVpzWYxYCEiIs+WvVlKh8U77nV8ZQHLgV8c9BrGC8fZcWNFD8fQjYiIPE/JBeDtq0zz015z3GvKW1iuXHTMa9ROw64VGueY1/FAbGEhIiLP8/Wt6vlhrRz3mvIxLL5a8+UaQj6GxT/YNIBpxNjCQkREnqOyFHizLVBVZnpu6BuOfW15wHLXl455DfmYla63OeY1PBQDFiIi8gz6auA1C3v/JD/i2NeXT2tu3sExryEPWATBMa/hodglREREnuHiUcvnHTU7qJZiWnOg+XINIV+MjgGLAgMWIiJyfxVXgNl9TPPvXQ5c9S/gsSzH10E+nsRR67DICXrHv4YHYZcQERG5nyuXgHWvAS17i2M5XpftjBzVGRg8A4hMBCLaAVfdYP4+9iRv8PAPcvzrXT7p+NfwIAxYiIjI/ez5Afj7c/GxY4HyXL+pQOIw59dJ3uLh56AuIbkrFxz/Gh7E5i6hDRs2YMSIEYiLi4NGo8Hy5cstlr/vvvug0WhMHl26dDGUeemll0zOd+rUyeY3Q0REXmLV01L69BblOWe0bqiJqvleCghx/HgZAAgMc/xreBCbA5aSkhJ0794ds2fPrrswgPfffx85OTmGx+nTpxEREYE777xTUa5Lly6Kcps2bbK1akRE1BjEJbnmdbUhwLMngacOO+f1OOhWweYuoWHDhmHYMOub4sLCwhAWJkWJy5cvx+XLlzFhwgRlRfz8EBNjYboaERE1DiVmVpFtNwgYPBMIb+3U6igENXPdazdyTp8lNG/ePKSmpqJNmzaK/CNHjiAuLg7t2rXD2LFjkZ2dbfYe5eXl0Ol0igcREXmJgpPq+TfPAlr2dGZNyI04ddDtuXPnsGrVKixcuFCRn5ycjAULFiAxMRE5OTl4+eWXMWDAAOzduxchISEm90lPT8fLL7/srGoTEZEzFZ4Rn1v1AUZ+BOjOATHdgCbNXVsvcimnBixffvklwsPDMWrUKEW+vIupW7duSE5ORps2bfDdd9/hgQceMLnP9OnTMW3aNMOxTqdDfLwDd+gkIiLnKLkIfDdOTPsGiFOXIxNdWydXccbAXg/itIBFEATMnz8f9957LwICLG+XHR4ejo4dO+LoUfVVDbVaLbRaB208RURErrPjCynd6BdOY8Ai57QxLOvXr8fRo0dVW0yMFRcX49ixY4iNja2zLBEReYETG4GXwoC1/5Xykh9yXX3I7dgcsBQXFyMrKwtZWVkAgBMnTiArK8swSHb69OkYN26cyXXz5s1DcnIyunbtanLuqaeewvr163Hy5En89ddfuPXWW+Hr64sxY8bYWj0iIvJEX95smtd5lNOr4VbYJaRgc5fQ9u3bccMN0jLItWNJxo8fjwULFiAnJ8dkhk9hYSF+/PFHvP/++6r3PHPmDMaMGYOLFy8iMjIS/fv3x5YtWxAZGWlr9YiIyBtMO9h4v7BjrgFy9wDdRru6Jm5FIwievzKNTqdDWFgYCgsLERoa6urqEBGRrf7XFSg8LabH/SSuudJYlemA3N1A677K3Zu9kC3f39xLiIiIXKu6UgpWJq4FWvZybX1cLTAUaNvf1bVwO94duhERkfv7/UUpHW06zpEIYMBCRETOUloAvNsJWPEEUFYo5W+dI6X9uGQFqWPAQkREjicIwJttgKIcYPt84I3WwMVjQHWVVCbAdGVzolocw0JERI4lCMDL4ab5H/YEBj4jHT950GlVIs/DFhYiInKsM3+bP7fhLSmtber4upDHYsBCRESOlX9ASg+eoV6mVR/n1IU8FgMWIiICNrwNfD8ByN4C7P4OWPYIcOmEfe5dnC8+97gXGPCkepmEAfZ5LfJaHMNCRNTYnN0JrHwK+Pf34gaD7yYCQrV4bt9SqVz+fuDhDfV/nTIdsGKqGAQBQNMo8bnDEODI78qyfoH1fx1qFBiwEBE1Np/VbK/ydjvL5XL+AXZ+DfS8VzzO3gIsuRe4/hmgz0Tz1236HxDcHLh8Etj7o5Qf3EJ8vn1eTauLAHzUW8yrLK3PO6FGhAELEVFjcPJPMYBIHGbbdT9PAVr1BqKuFhd4K8kXW2dWPgX8Xy7gH6Qsf+kE8MdL6ve6eFR8DgwVH3IVJbbVixodjmEhIvJ2uXuABTcBPz0K7F5i+/W1g2bPbFPmb5plWvboH+bvE5lo/lzzq2yuFjUubGEhIvJmmW8AmenS8emt6uVCWwK6s+rnjFtRaq1/A9BXSjN/9Hqx5cWca+40zXtgjRjk9Jpg/joisIWFiMi7yYMVANi3TL3ctP3Ak4fVz1WUiIu/qdn4rpTe853lugQ0Mc2L7wPc8DzgF2D5Wmr0GLAQEXmrqnLryj1aM4snJBoY97OYDggBOqSJ6cpSoKK47vsse9jyeV8GJVR/DFiIiLyV2hgTY5P+EgfU1mp3PTB1D/DEHsDXX8zTVwK6nIbXR6Np+D2o0eIYFiIib1NVIXaxZL5uvsx9vwJBEUB0Z9Nz4a3FZx9f8Vlfrez6UVNdKaWTxgI3vgKsehYovQQcW2tb/YlUsIWFiMhb6KuBb+8CXo0EshZaLhvRTj1YkfOpaWGprgR2L5by04zGxQiCcmzM8PeAJi2AO+YBHYdaX38iCxiwEBF5i1cigCO/ienlkyyX1YZaPg8APjWN8L9NV+b3vl95XFWmbEXxl61am/RvoN0NYhBD1ADsEiIi8mQXjwFfjQLKC9XP+/iLY1BqxXQD+j1u5c7IKjODbvlQDEhmXBIDJECcRfTPIvVbaEOAccuteC0iyxiwEBF5sg97Wj4fEgsMek5cNG7kx0CPsdbf23i/HwDoPEp89vEV9/+pKlOuUtvCwuJwRA3AgIWIyFOpriqrgaJlJKylGKRcPcJ0Ofy6CHrTvABZy0xAEzFgqbwChLcBCk4BI2bZ9hpEVuIYFiIiT1R6GfjmdtP8xJuUx02jxWdbgxVAfQdlH9nXhn/NQnBlOjFYAcRND4kcgAELEZEn2vmVen5EgvI4JKb+r6EWsMgFBIvPvz4p5WlD6v96RBYwYCEi8kQ5/0jp5EekdFmh1PIBNKzFo66Vac8fFJ/z9kh5tS06RHbGgIWIyBNVlkrpKNl6Kp1uBm7/XDoOsGY2kDlm9g+ypHaxOSI7Y8BCRORpqquA45lieuI6oMe90rnYbsrdlRvSRaM26FYuMKz+9yayEQMWIiJPceWSONh24V3izBwAiE0SB8I+tkvcFyg0DvAPlq6xar0VM+TL7au57lHlsbxrisjOOK2ZiMjdCQJQeAZ4vzsgVCvP1c7aiWgn5clbWAIa0MKS/Ajw+/+ZP68x+pu3Za/6vxZRHdjCQkTk7ja8DczqahqsmBMgG3TbkC6h6x4FHlwL9H9CPO51n/K8cZeRvGWHyM7YwkJE5O7WvaaeH9NNPV8xhqUBXUI+PkCrXkBcD6DzSCC6q/K8cQuLfx3ToIkagC0sRETuSJcDXD5lfhxJj3vEAbdq5C0sDZolVMPHRwxafP2V+cYbKPoYnSeyI5sDlg0bNmDEiBGIi4uDRqPB8uXLLZbPzMyERqMxeeTm5irKzZ49G23btkVgYCCSk5Oxbds2W6tGROQdLhwB3usEvN8N+G8L9TIjZwO+ZhrJA8OA/tOAaycC4a0dV8+e9wIte0vH+joG6RI1gM0BS0lJCbp3747Zs2fbdN2hQ4eQk5NjeERFRRnOLVmyBNOmTcPMmTOxc+dOdO/eHWlpacjPz7e1ekREnm/nlw2/R+pMYPg7gEbT8HuZE9AEmJghHeutHGNDVA82j2EZNmwYhg0bZvMLRUVFITw8XPXce++9h4kTJ2LChAkAgLlz5+LXX3/F/Pnz8dxzz9n8WkREHk2jsvharwlAUS5weJXz61OXrneIK+8mDHR1TciLOW3QbVJSEsrLy9G1a1e89NJL6NevHwCgoqICO3bswPTp0w1lfXx8kJqais2bN6veq7y8HOXl5YZjnU7n2MoTETnLT1OAXV+b5ncbLXa5FJ4BhqY7v16W3DFPnHrtyNYcavQcPug2NjYWc+fOxY8//ogff/wR8fHxGDRoEHbu3AkAuHDhAqqrqxEdrdx/Ijo62mScS6309HSEhYUZHvHx8Y5+G0REjnd6m3qwAgBBzcQWjEmbgIQBzq2XNRiskIM5vIUlMTERiYmJhuO+ffvi2LFj+N///oevvzbzH7MO06dPx7Rp0wzHOp2OQQsReb55Q8yfC2rmvHoQuSGXrMPSp08fbNq0CQDQokUL+Pr6Ii8vT1EmLy8PMTHq26JrtVpotVqH15OIyLksbDYYFO60WhC5I5esw5KVlYXY2FgAQEBAAHr16oWMDGmkuV6vR0ZGBlJSUlxRPSIi5/t7nuXzfvwjjRo3m1tYiouLcfToUcPxiRMnkJWVhYiICLRu3RrTp0/H2bNn8dVXXwEAZs2ahYSEBHTp0gVlZWX4/PPPsXbtWvz++++Ge0ybNg3jx49H79690adPH8yaNQslJSWGWUNERF7vV6mbG7d+Ku66vO414MAvrqsTkRuxOWDZvn07brjhBsNx7ViS8ePHY8GCBcjJyUF2drbhfEVFBZ588kmcPXsWwcHB6NatG/744w/FPUaPHo3z589jxowZyM3NRVJSElavXm0yEJeIyCuVXFAedx8tPqtNbyZqpDSCIFjoNPUMOp0OYWFhKCwsRGhoaN0XEBG5k0OrgUU1Qcr/5Up7AZ0/BHz2L6D3/cCQ/7qufkQOYsv3Nzc/JCJypf0/A/uXi+nIq5UbF0YmAs+eMr8EP1Ejwv8FRESuUqYDvrtXOo5MNC3DYIUIAHdrJiJyvlObgR/uB3L3KPODm7umPkQegKE7EZEzlRcBXwwV03t/VJ4LDHN+fYg8BFtYiIicaetc8+eCI5xXDyIPw4CFiMiZ1r5q/lwYtxghMocBCxGRo/09D/j5MaC60nI5biBIZBbHsBAROZK+WlrF1i+w7rJEpIotLEREjlSUI6W3fVKTMNOSEtbK4dUh8lRsYSEiciTdOdO84AjgykXp+KZ3AEEA4vs4r15EHoYBCxGRI53dYZoX3BzocQ/w5/vicc/xgF+Ac+tF5GHYJURE5Ej5+03zgiKAflPFdHgbBitEVmALCxGRI106YZp3eovYLfRSofPrQ+Sh2MJCRORIpZddXQMir8CAhYjIUf78AMjbK6ZHfyPl3/era+pD5MEYsBAROULObmDNi9JxRDspHdTM+fUh8nAMWIiIHMF4OnOzBCkd0MS5dSHyAhx0S+SucnYDJflA+1RX14RsVVECLBotHccmAQHB4vTlqjJxZhAR2YQBC5E7KNMBl08Csd3EY0EA5g0BqkqBB9ZwQTFP87+uyuPAUPH5lg+cXxciL8EuISJ38EY88MkAYPNs8bjglBisAMCJ9a6rF9VP6SXlscbXNfUg8iIMWOpDEFxdA/IG+mog4xVg3zIp77fnxeeMV6S8ta8CR/9wbt2o/irLTPO63+38ehB5GXYJ2SpvP7BgOND/CaDfY66uDXmyLXOAje8q8zrdLD7v/VGZ/83tZhcZy9eV4dMNxwEADw1sh6jQOnYEJsc68Iv47BcETN0NXDwKtOnr2joReQG2sNgqM11s7pVPVzRWXgR8cRPw10fOqxd5nt//Tz3/0nHTvG6jTfNq5BeV4/NNJ/D5phPILyq3U+WoXs4fApY+KKYDw4CmUQxWiOyELSy2EvRS+tAqIHGYdPzXh0B1BXBiI3DqT/HRd4rz60ju73imen5FMVBebJq/ewlQVQ4MeFIamEvuZ9fXUro413X1IPJCDFhs5Sdrbl90t9RMX14E/P6CafnqKsCXHzMZ+Wqken55EaCvVD+3fzlwfB3wXDbydWU4mFuESyUVOJYvBTi/7cvF0fxiRDTxR6eYUHYPOVPpZfGPFiJyCH6T2srfzBdAcb56flEOEB7vuPqQ59HlWD5XZaFbp6wQyPkH3+4NwvsZR0xOf7j2qCH9+OAOeOLGjg2pKdli59fK47gerqkHkZdiwGKrgKbq+ZdPquef3cGAhZQyXpbSrfsC2X9Jx0XngC9quhnD4oHY7sDBFcrrl9yLsfdvQ682zQwtLB+uEwOV//yrPa6KbGpoYSEn0p1VHicOd009iLwUAxZbVZZK6YAQ8XnNTODPWerlt84FuoxydK3Ik/yzSEpPWAloNOLsszkpynKFp4HbPjMNWIrzEBUaaOju2Xu20BCwpHWJQdeWYY6sPZmTs1t8HviM+EdKN05lJrInzhKyVblOSlcUAWtmmA9WAOWGZ0Ry4a3FYAUwv7eMWr584Dc5V8kFYPf3QGmBMr9MJ7WUtbse6DkO8AtwevWIvBlbWCypLAU+6gMUZgPD3gKSH1Yu8gUAf75fxz2uOK5+5Dk2/Q9o1hbocqu4j0zBKeDWT6Xzarv3xiZZFbBEhWjxYP8EQ5oc6OMUcX+nWs9li9OXf31SyotNcnq1iBoDtrBY4uMnBisAsOoZ267tOFR8PrCCK+M2dn99CPzxEvD9fcDlU2KwAgBhraQygSrjTfyDxIcxfZXiMCo0EC/c3Bkv3NyZs4IcrcRocP2618XnM9ukPK2ZcW5E1CAMWCzxMWqAKisEmre37trOo8RnfaVyzAI1Lvpq5XT39W9K6ZBYy9f2n6Y+yDu4hX3qRrYz7uI9f1Ccil476H7Q806vElFjYXPAsmHDBowYMQJxcXHQaDRYvny5xfJLly7FjTfeiMjISISGhiIlJQW//fabosxLL70EjUajeHTq1MnWqtlf7fiCWhvfNQ1izGnbX0ovnwRkLRS3nKfGZeN7yuOsb6W08fo8t8hWRu7+b6DjEPWWl4Bg86+nr7a9jmQ941WIfbXAhrel47CWzq0PUSNic8BSUlKC7t27Y/bs2VaV37BhA2688UasXLkSO3bswA033IARI0Zg165dinJdunRBTk6O4bFp0yZbq+Z4f74vLr2tpv804KZ3pGN5cz8gBi2rnnVc3ci9HF8PvBQGrHtV/XzvB0zzYrpK6esmSWlfo8Gbpep7CmH/T8DrLcVnsr+Si6Z5R35TjmPz5RgiIkexedDtsGHDMGzYsLoL1pg1a5bi+PXXX8dPP/2EX375BT16SAsr+fn5ISYmxtbquICZ8Sj+wUCHIWK6eXvT1hlAXLZ7JPcXahS+usXy+WZtTfPk3T/ysSsao78rygsBvR7wMcr/bpz0bGajRDJDXy0uQRDdBWg3SP38t3dIxy06AhcOq9yI49WIHMXpY1j0ej2KiooQERGhyD9y5Aji4uLQrl07jB07FtnZ2WbvUV5eDp1Op3i4nH8g0KwNMO0g8EhN69CYJa6tE7mG2l/ixrOA1DYzlM8I8pP9pf7wRvE5rqeU989C5bWHf7etjqS08ingt+fFLRNKLpieP7sDOLdTOk76t/p9rlxyTP2IyPkByzvvvIPi4mLcddddhrzk5GQsWLAAq1evxpw5c3DixAkMGDAARUVFqvdIT09HWFiY4REf7wYrydb+RRwaK6WNxxr4+Du3TuQa712tPE6ZAlwta3HpegcQEm16nb/s50XeqhLZEXj6GPCALCiRT6MFgGUPKY85M8022+dL6bevAn4y2rT0xHrlceJN6vfpOMS+9SIiA6cGLAsXLsTLL7+M7777DlFRUYb8YcOG4c4770S3bt2QlpaGlStXoqCgAN99953qfaZPn47CwkLD4/Tp0856C+b5qUw/Nc7z45RTr1dZBlTL9gJ6/B9gyKvKBQflqyXLybuEjFtkmrQAfP3Nny+9rDw+8LP1dW7s1PZu2mW0L9DBlVJ64NNAZCLQ6WYpTxsKTN3DhSKJHMhpAcvixYvx4IMP4rvvvkNqaqrFsuHh4ejYsSOOHj2qel6r1SI0NFTxcJg7v7Su3Km/TPOM19CoKOJfvt6stAB4K0GZ16ytOJ6pWrYDs7kNNH39gMeygCk7zK98W6soB9j2mfjzVFVher52PAvVrTjP8nndOak76MEM4F8109TlAea45eLKxUTkME4JWBYtWoQJEyZg0aJFGD687g3BiouLcezYMcTG1rFOhTN0GWU6S0ONoDKdVG3RryX3AEV54qBJc39pk2fKeEW5srF8nIr8y6yJ1LpoIiIBaGFhrZ8BT0nplU8Bi/9tfiYSWce4dcrYknuldHBzKS0PKuX5ROQQNs8SKi4uVrR8nDhxAllZWYiIiEDr1q0xffp0nD17Fl999RUAsRto/PjxeP/995GcnIzc3FwAQFBQEMLCxE3annrqKYwYMQJt2rTBuXPnMHPmTPj6+mLMmDH2eI8Nd9W/gMOr1c91uVVcrj/tddNzagFL7UZ2Gh9x6uuUv9XHM5Dn2T5PeSyf5j7gSWDLx2K6/9T6v4Zx68yhlYDaTPvARrIBol4vLs7o14DpxD+oTDGvvffBFcDZ7VKePDAJlk0c0HJnbCJHs7mFZfv27ejRo4dhSvK0adPQo0cPzJgxAwCQk5OjmOHz6aefoqqqCpMnT0ZsbKzh8fjjjxvKnDlzBmPGjEFiYiLuuusuNG/eHFu2bEFkZGRD3599jPwYaN5BmffwBuDRrcCdC8QppMERptepBSyA+EvwwM/i9FT5QmLkucqNBoiP/0W56FuTFuLPyYzLQEgDpu+rjZWq5eMvLT7XWBaQ++Y2YFY3oLy4ftfn7gUuHlE/V1EMfHevMk8bIqVD49TzicghbG5hGTRoEAQL4zAWLFigOM7MzKzznosXL7a1Gs7VpDkwYJq4+Fut2O51X+dvYUXSWo3li8XbnflbSve+H2g7QL2c8doptjI3/gUQWxpqg6SKYnFVVm8eBFpdBRxfJ6b/nCWNLbHF6ueUx+1ukO75y+Om5eXrK8lncvlyBiCRo3G3Zmudlm1udvdC8+XkrBn7Uu4Ga8hQw62ZKaVv/p/jXqeOmWb/WfIPPqz9Ti08C0S0w/5zhXhiSRZyCktRrdejskpAk0B/fDimB/p3cJNWTGuVF4mtKaGxQHGulL/hbXF8z4KbgIJs4ME/1BfnkxME4ORGZd645eIKxQCwb6nynHG3T7sbxFatuKR6vBFqLPafK8R/Fu7EqQtXUAWg31UR+N/oHqoblebryvDsD/9gw+ELiG8ejI/H9kTnuEbSvWsFbn5orYFPS+kmVv6SV1vt1lhRTv3qQ+4ld7f4bO3PRn3VEQSXyiYj1Q7qPpxXjEN5xdCVVaOkQkCFHrh8pRK7sgscV09HeT8JeK+TuDiffA8fAPjoWnGBt5LzwPtWtIAaD7Yd+6P43OVW9fLVRrOxmrUBph0Qu/+IzDicV4xjNcEKAPx57BLyi1Sm0gPILyrHusMXUA3g5MUrOJxXz65OL8UWFmuFtQQGzxCb2Vv2tv66Pg8B2z41f/6i+tRtcjMVV8RgwXjDQkDZrdfTwdOJ65gWv1F/jXRQVuDYujhbeTFwpWYV2h/vB45nKs8Xml8dW1WBrLx8KwNz41GuV9kLrKmHtVAReTAGLLYY8GTdZYwNe0u8btfXwFqV6afndolfQta0xpBrlBYAHyQBEVcBD6wxHYey7TMp3ft+h1Zl16kL6GHm3CuV96IcAdhc3Rkpvvvx89Jv8NqKZtCVqqzTAuCjdUewdNdpdI4Nw5g+rd2/e+iIbKVf42BFTV3/rzLfEJ+Nx/mozfgZNQfodnfdr0kEsRto0bZsrDuYj/yiMpPzoz/ZjOZNAnBDp0iMvrY1Vuw+h01HLuJSibLl5bVf9+PT9UfRpWUYJvRLaPTdQwxYHE2jEWeFJN6kHrAAQP4BILqzmK4qB85lAS17qf81T86362ux++DsdnGRsVDZ+kDVVcBq2V/exrt029nSHdnoYeY7eH61uClpiu9+AMAtWI/HdA+bvVd5lYATF0px4kIpsk4X4M/nBtu9vna1ZY5t5asrAT8zXWiCABxeJaYvHTe9zpi5vYOIVLz8y35sPWF+X6mSimqUVJTiy83ZOJhbbLbs+eIKnC+uwP7cYmRfKsWSh1McVWWPwG9EZ4m8GghtCejOmp6TLzr3y1RxY7vrnwVueN5p1SMV698C1r2mzLtyQRmwZKbLTjq+lezOHlFAlmn+q5VjDent+o7o7SPuJBwdGgBdaYVybEsNrZ8GceGBhhYWt6e6O7IFVaXmAxb5JoXGi77VjkeqlTwJRLaYOaKzooWlwmgyaJMAX7MtLGcKpBaZyKYBiGwaYGhhaewYsDiLjw9wz4/Ax9eZnpP/RVe7C+/6txiwuFJ5kWmwAih38tVXAxtli8NN3W1a3s66demmGrD8UD3QkJ5XNQy9Aw4DrVOwtc1mYOO7SKn8EDlQfjFPuaED/jO4g/Gt3JMgmB+Tc92j0qJ8cpVl5hfQ052R0pM2K8+FtwayZXnD3rCpqkSd48Lw31HieLLlu85i6pIsxfklD6ega8swRflnhgJ7zxbi5g83GfL/b3hnjOrR0il19gScJeRMxgvJ1S5Gp7pEP/cccrrT24CfJotBSbqZrp0rF6V07XodtULi4HDtU8VVle9bqcguh7QOSClqVn2tKAE2vgsA2Bz4H4RDffdzjzCnn3r+mMXi5pJqO6FXWdj64vBvUtp4pen+T9hePyJyOAYszmS8SmltACPff4ZcZ96NwK5vgLevMl/mxweA8zVr4Z/YIOWPnO2cMUcaDZAyGWir/AKvqAlYmmh90DyiZidno66NGYGL0SRAgwAfoFmwP3q0Dnd8fe0lf596fuIwwMdXXDTPmHwX5vwDwB8viwOoAfXWs1pRV9e7mkTGOkY3xVUtgg3dGf2uikBUiPpWElEhWtzQsQV8AbRtHoyO0U1VyzVW7BJyJnkLy73LpFkK5jZBrCw1v7w/OY/GF/DxA6prvgDn9AVSpgB/vi8eX3MX0OMe59crKAIoFcdiHHvjFik/5x/gE9Pit7Uqxm0P3uSkytlRfVeDlv+/+jgFgAAU5QK3zgGiuohBUG8z+wgR2UnnuDBkPHWDVWWjQgPxxf3JDq6R52ILizPJg4+gZtJxVc0gqwqjlpbNs51Tr8buyBppdVNjnW4GZl4CkmWzbfRV4lLwhmOVv+6dYcJK9fxmZgbnFec5ri6OtP5NKZ1oQ8BVJZ9OWtPF+s9CMQCqbbFJGmtyGQBx/SQA6HKb9a9HRA7FgMWZ5PuNNEuQZiosnSj+Nfh6rLI8F5VzvHO7gG/vUD/X9Xbg1pqmCkt7xbTsZf96WSPqamDcz8DDRsvLB5rZObiOZf3d1vb5UnrMIqBNTXdYcAspXx6kNW8vPptrufxNNpjd3LL6N74CjFkCjPzI5uoSkWMwYHG2Jw8Dj/8DBIUrxxjs/Nq0bB2rmpIdnNtl/tzw9wBtTR9ygJm+5OYdgGsftH+9rNXueiC2m3VlbZ0W7A7KCsWl9gEgvqap/OZZYjecvIXp4Q1A51HA3Yukhd+qTBfsAgBsnSulfXzVy/gHAYlDgYAmDak9EdkRAxZnC4lW35TtygXTPHNrSJDtqiqAi8esK9t/mrjBZVC4lGduyf2+UzjOyJHOZUnpoTVr3kR2BG7/DIhMlM4FhgJ3fQl0ukk2mN3CLCEi8jgMWNzFXypNz+UePA3V3bwaCXzYE/jzA2DBzcCZ7WJ+cc1f7807AG0HAEPfBFJnAp2GK69v0gK46yvT+xov6+4uEq5Xz8/d49x6NNQvj4nPvlrru95qu75qZwmZa6nsN7VBVSMi52LA4kqtZcssV5ZI6Q5DxOd9y0x3pCXbnZd1hax5ETi5Efi8Zhn6zNfF545pwH0rgOseMX+fyE6meW0H2K+ejhIiGxv1az32w3KlyyfFZ+O1UiwxDGavaWGRT2+WMzfWh4jcEgMWV7p9nnq+fN8Sc/sPkUQQxIGZpzarn599rXr+quekdFFO3a9j3PXTdoBnbFo5VLZ9wOmtrquHLYpylevc3Pw/66+tbWGprBnDYm4si9omh0TkthiwuFJYS/W/0Jsa/TW54Gbg9N/OqZOnEARpS4OzO4AVTwBfDDUtl2dmwTEA2CrbTC/1pbpfs4lsN+NpB4Gx31tVVZeQB1JdbnVdPeSqKoCSi3WXA4BlDwNfjpCOW/a2/nVqA5asb8XB7OYCFn2V9fckIpdjwOJqITGmecb7n5zcCHw9yinV8RgvhwP/bQEUngX+WSTly9eyyT8gLvJmjXArNv/zDwIeyxJneYXGuvdg257jxed4N1mEquSiOI7o7Xbiv1ldjmcqj+UDoOviXxOw5O0Ffp4CvJuoXi57i/X3JCKXY8Diar61SzTX/EV8zZ3qs4gqip1VI/e3b7mUXvUM8Pfn0rHunJRe/5byOnOtKLasUhuRoP7v42663Ao8lCmuqOwOfp0mpY/8Zr6cmkgbl8q3dr2ZwTNsuy8RuRQDFlczLEhWM5MhNE5c++Hp4y6rktv7fryUPrhCea52XZXDvwP7lirP9XlIXLn2qsHK/PC2dq+iy2k0QFwP91lHZP9yKX1sndliqrQhtpW3puVr4lqguYU9o4jI7TBgcTU/o02waltcgiNMy+r1jq+Pu6vdvM6cY2vFsS0L7zQ9F9AEuPtbYNQcZX5ZHfck+zrws+XzxnsHndlm2/3NtbDIW8ZCYtXLEJHbYsDiasZLvtcuFqc2+6SxT3EuLwLebGO5zD8Lgdkq4zaCm0tp+eBZABAaWSDo7BWUbQ20VzzRsNczt42CfB0XX616GSJyW9yt2dWMf3H6WWjOznwdGPSsY+vjzn5/0bpyl2Qr2j68EfANAKJka6j4GMXp1z3a8Lp5EkEv7kDtLFeMZga16mO+rCAAO790TD3k2ysYt2wSkdtjC4urGc8Iks8a6na3aXn5X8clF8QBqMsnA+tel/IPrACO/mHXaroFcwuAAeJS+moi2imDlVoBsnER4fENq5cniL9OSlc7eXfpDUaDn/0tDIrN22ua9/AG0zxLzLWYycfCMGAh8jgMWFwtNE55LB+7cutc4IV85bos8jUlvhwhDkDN+gZY/6a4UdzfnwNLxgLf3C4tnOUt5FNdnz0l7doLAG3MTF82N+j04fXiLJHpZ+xWPbcmD+icvf7Itk+Vxyc2mN92olw2G27iWmBmARDb3bbXM9flJQ9YfNi4TORpGLC4mnHA4i/7gtVoxL8E43pKeRWyJfzz9yuvLTyrXHpdbUNFVzu9DVg2SWwRsjWgqm2NSn1JXJdjyH+Bq/4lTt8NagY8ddT0GnMr0Ta/ChjwpO0zUDyV/H06c4q8/OdVK2tN3P+Tevna7qCgCHHMSX1WEjbXwhLVWUp7wgrFRKTAgMXVgpopjwOCTcvcKpvVUvsFcP6QabmCbOXxhcOmZVxp+3xg3o3iwNj1bwI/PgCc+sv662tn87QbJD637CWuMxLXQzxuGql2FQHKFoXdSxz7Wvpq4OQmsetJ/jN553wpLV8vR662269Ba92YaWFp3h4Y/4vtXUxE5BYYsLiab4Dy2F8lYAlqJv7FCQCVNSu5zlYZuJjzj1QOAL52kyXZAXG8jfHsj4MrgC+GifUGxC+51c8DhWa6aWq7EYzH/ZgzZXv96uqN5AONLY0FsodXIoAFw8WViGu3RmjVRzmgXG0cTXG+tHZOhxvr//ry1Y7l/LRAwkDbu5iIyC0wYHE14750c2MuSi+JzxePqZ8HxFlE7trUbbzUutwnA4Ftn4lfcltmA//rojx/4Bfgs8FSV4bWioCl53igRYd6V9erGQfJjpR/QHyO7qz8WVcbR/PlLVLa3P4/1jC3zQIH2hJ5NI48czXjLw/51Es1a18Frr7Z/HnjKaTlxYC2jns6g6WABQBWPqWef+USsES+dL7GcgvLC/niNaFcGMws+VgOezMe8HqpZsXm5u2VY0vUApbzB6R0v6n1r0PSWGCFyvVce4XIo7GFxdWMA5a6gotOw227f3pL28o7yq6vbStfphO//N5KUOaHtgR8LcTZfloGK+bUrn9SXeG419g6V0r7BgC6mo0Ow+KVQUpt8KKvBl4KEx9yais9W8svAIhQWXafLSxEHs3mgGXDhg0YMWIE4uLioNFosHz58jqvyczMRM+ePaHVatG+fXssWLDApMzs2bPRtm1bBAYGIjk5Gdu22bgct6ey9OUr1+u+mvIBwNkdDquOQxiPKeg1oe5rCk8DZ/42zW8aZZ86NUa1X9jVDhzDIt9wsroSKMoR0yExgCBbcn/Lx8B7nYHMdNN79Hu84fW4dxkw4Cngts+kPAYsRB7N5i6hkpISdO/eHffffz9uu+22OsufOHECw4cPxyOPPIJvv/0WGRkZePDBBxEbG4u0tDQAwJIlSzBt2jTMnTsXycnJmDVrFtLS0nDo0CFERbn2C2r/uUKMn78N54uVf5UG+mmQ0KIJDuSaThEN8gf8fHzg5+ODQZ2iMH3Y1TicV4QnlmShedMATL6hPb7Zko17rmuNzzL2Qb6zStvnflWtxwt+eXjQD5jzxx70ySlD7SLjy6v7Yqv+aqT7z6vzfTyxJAs5haUQBAHllXpUyFrofQD4+wHlVer1jwq1cgdcNYWnpfSMy+Jg2x1fWL7m58fUB14az6rycPm6Msz64zD+PHoBZy6VorruS0xoIM6L8YU4trbSzKzer/wLMNAXeGbx31hWHYhKQbw2NMgPfa9qjrKKamw8fAHGnTVhQX6Y/e+e6N/B/CysWb/uxLVbpqCf7yVZrmCYJTRwzn5UCz74s/bHSNCLrS9q200kjVUc5uvK8OwP/2DdYfVp+kEBGgT7++HBAQn4bMMJlFZWorQSAHqin88efFvTiPnL3vMYkdTK7HtoTDYdOY/HFu1CZXU1qvT6ms/LvCB/H3w2rrfFnwFXk7+nymo9ymS/ywL8fPH+3T1M6r/pyHk8+u0O6Mqs+58XGuiLj8f2svpzyNeV4cEv/8bus7o677nk72z8sjtXca72/7Yl/j5AaudovHxL14b9nvYANrewDBs2DK+++ipuvdW6GShz585FQkIC3n33XVx99dWYMmUK7rjjDvzvf/8zlHnvvfcwceJETJgwAZ07d8bcuXMRHByM+fPnW7izcxzOKzYJVgCgrEpQDVYAoLQSKCrX43JpFZbtOof8onLsyi7A+eIKHMwtxtbjl7D1xCVsPX4Je/Kt+2u3sia29EM1Dl+S/nNNrZyMRdWDzV1W8yZ+x5GcyziUVwxdWTWKypXBCgDoIQYr5urfILUDL2OuEb9RA0PNFNRIM0nOblf/67v0kmmeB8svKsfCbadxqp7BCiD9QquG+WBFPC/+d9cLelQK0rWFpVVYtTcP61SCFdSc35VdYLEOw3Y9jH6++8yePy+E4Swi8Wu1hWX5a4UquzHzi8rNBisAUFoh4GJJJTYeuYhLVyoVX77yIehbT1yu+7UbiV3ZBbh0pRJF5XUHKwBQWqmv82fA1eTvqczod9nFkkrV+u/KLrA6WAEAXVm1TZ9DflG5xWBFfk+1n09rdv2q1AOr9uY1/Pe0B3D4GJbNmzcjNTVVkZeWlobNmzcDACoqKrBjxw5FGR8fH6SmphrKeDMBPphQ8TQA4Jje/NiLipqAJQjlGHP+fQDAbn0ClL+SZeRN4QvvROLhuerlnKF2HEPzmlk7bQdK57Sy4CUoHLjxFcv3OrfLrlVrTGoDFh/Yf7PHxCrLa/6UQvzLb5P+mrpvZm6mXD3s1rcDAOQL4Xa7JxG5hsNnCeXm5iI6OlqRFx0dDZ1Oh9LSUly+fBnV1dWqZQ4ePKh6z/LycpSXS9GkTmc5grXV/nOFeG3lfvx9/JJJS0R93D7nTwh6KVZesk1sJv9mq/i8Tt8Dg8rfRa5gfqBhpSD+U93jl2HIa6mR/urcUH0NBvruMRyP+TMGi2TXdzo0B8CAetV/4ld/o3NsKAYlRiGtS4ztzY7na/4da1f1la8JEhgGlNf8+5VeBpL+Dax6Wnn9f3YCH9as9htex27NHiBfV4Ytxy9i1d5c/Hn0vNNeV28IWGzfrfndNYfx6YZjSG4Xgfv6JqB/h0i8vfoACrYuxGvCB1bfp1SwYhyJRoN8XRkW/HUCP2edQ26BdVOc/zp20SRPhya4puxzVMAP5VuzsfjvbFzVoimmDG6PEd3dZEC6k2w6ch5Ld57FgRwdjuWZ2RrBgnfXHMaHaw+jY1RTPHdTZ7foHpK/p6O5lt/Th2uPYNXec4hoEoDisiqculiKy9Y0Lxl5d81hfL7xOLq1CsPD119l8jnk68rwyop9yDiQh9JK6/6vvbum4Yt83j7nT7SNaIKRPeJwR694r+we8shZQunp6QgLCzM84uPtu3ndy7/sx59H7ROsAEB5laC4l9rP8EkhFmUw/8u8UiW2LJWV3ysoZ9NsPmX7LyRzcgrLkXHwPF78aR++3Zpd9wUmN6hZGC5e1h3Q73EAGmD4u8qyarOkml8F/Pt7IKYbcPe3tr++m/l2azYeW5yFVXtzbWqObqjagMW3ni0sReXV+OPAeTz7424AwOzM46rByjdVyi7KGZXjDelS1LEGTE130Ldbs/Fx5nGcKShT7aayqd4IRnnN61bpgUP5xXh+6Z46rvI+z/64G0t3ncWB3CJU2B6zAgAqqoG9OcWGnwFXk7+nukKPimoB+3OKsenoJWSd0dUrWKlVWFaFjUcvqn4O327NxorduVYHK/ZSXiXgUH4x3vrtcP1+T3sAhwcsMTExyMvLU+Tl5eUhNDQUQUFBaNGiBXx9fVXLxMTEQM306dNRWFhoeJw+fVq1XH3NHNEZ/dpHIMBOn47WT6O4l3891nZTC1iKBGnl0AN6abGsk/poJLcNt/1FzIgN02Jwp0j8d2QXjE02syiXJbqamSLNZEFV6svAM8eBjmlAyhQxb5jKAMzOo8TnjkOARzaK42A83Njk1vjg7iQM6xqD0EBfp71udU33YX27hEK0vki9OhJv3t4NADB5UDvsgXJxvqHlb+CFqvsVeXmylsNy+Ft+kZoNLccmt8ajg9qhVXigXZuB/XyAxKimeP02z/85stWbt3fDbT1a4uqYEATUc33JAF+ga2xTw8+Aq8nfUx0/WQjw1aBzbFP0bx+BpFahaBZU1xXmhQX6YUD75qqfw9jk1ri5WwyC6vOLvgG0fhokRjXFM2kd6/d72gM4vEsoJSUFK1euVOStWbMGKSkpAICAgAD06tULGRkZGDVqFABAr9cjIyMDU6ZMUb2nVquFVuu4KYqd48Lw7YNi/ZbvOoupS7IadL8fJ/XDuoP5hma/0X1a45ut2bgnubWhW6guFSr/VOWyv1Z/0acgAZfRpvoUut7zDpYkXg281KBqG3w27lp0bWnlcvjGKkuBknwxHSIbo6PRSGttpL0mPmrd/xswP63mwLl/pThDVGggbklqiVuSWmLv2ULc/OEmp7xuQ1pYnryxI/4zWBmcPP2v1kB2CCDbFuigIP6i/KF6IO7wFffsKYbUNF1nd1TN1OOo0EA8M/RqPDP0aqs/o75XNVftFqp1T3JrvHpr4wtUavXvEGnovvgw44jN3RBqPwOuZst7+s+/OpjU3xGfQ1RoID76tziP09qf3Sdv7Iivt5xq0MDZHyf1q//vaQ9hcxtCcXExsrKykJWVBUCctpyVlYXsbPGLd/r06Rg3bpyh/COPPILjx4/jmWeewcGDB/Hxxx/ju+++wxNPSPvKTJs2DZ999hm+/PJLHDhwAJMmTUJJSQkmTLBivY5Gogqmf4lXK/75NPgCI/Fk5aOobBpnUhYAouCCWRJvyRbwCm5u3TVRV0vpgyvNlyObNGQMi+nN9MDrccC5naqn/6iWdhg/LxvwWmew5Od9/e5EZB82Byzbt29Hjx490KOHuEPutGnT0KNHD8yYMQMAkJOTYwheACAhIQG//vor1qxZg+7du+Pdd9/F559/bliDBQBGjx6Nd955BzNmzEBSUhKysrKwevVqk4G4rtAxuikim5r2uwf6aXB1jPqqtEH+QIjWB82C/HBrjzhEhWjRo3U4IpsGoFNMUyS3i0ByQgSS20UgMboprGmZrB10K/ePXrmaZ49WoXh8cAdEhdS0PkVerTh/Z7PDCA30RYjWx6S7yweA1s98/eul8CxQWSJ7ESt/3Pxls0Ta1zFl28NFhWjx7z7xaBMRpBKSWqe24dkX4poM5tQGuAE+ekO3pAbiOivDukbjho4tVJtcw4L80KN1uDLz8gmTcnpBagKP0kjB8TFBCqALhDpWclbZ5ygqRIsbOrYwe0lQgAbNm/hjQIfmiAj2V/3/5O8DJLdrwOq5XqZH63BEBPsjROtj1e+fIH8f058BNyN/T4FGv8uaN/FXrX+P1uE2dcuGBvra9DlEhWjRraW5ZRyU90xOMF1jyppOJX8fYFjX6Pr/nvYgGkEw3vzD8+h0OoSFhaGwsBChoZZ/ODzWnh+AHx+QjnuOA9LSLS/lv2IasF22oNyEVUCbvo6ro9w/S4BlD0nHoS2Bafutv37D28C2z4FJfwJNzH9ZkQ2WTQL+WSiOH+o/tWH3OrsD+OxfpvkvFYrPuXuAuf2B4BbAM7INOwUBeDlcec2Q14Df/09Mt7oWePCPhtWNiDyGLd/fHjlLqFHyNfoz6JYP6953KErZwtKgHXBtkX9QGawAwAQbu3YGPg08dYjBij3VtnAJdpj+VmJ+ITcA4uDoh9YDk4222NBogBHvK/PkqxerbcdARAQGLJ5D3lTexbpVhtFrAnD9c9Lxnh/sWydzPk42zWvW1jmvTeZpagMWO0ylzrViWnBcEtBEZdxSz/HAnV9KxwHBUnrgMw2uGhF5JwYsnkLewmLtSqC+fsAN06XjrG/FJnlH0euBRWMcd39qGE1NX73eDi0slVdM86K6WFkPDdCyl3Ts3wR4PgcYvwK4/tmG142IvBIDFk8hb2HxDzZfTk10Vyl9fJ196qMmezNwSKXrJ6qz416TrOdTE7A0tIWl4DSw0WjBv36PA+N+sv4e8tlA/oFiK0vCAOt3LyeiRocBi6doSMDS/W4pXVpgl+qYEARgwU2m+ddNBu76yjGvSbbR2GkMyzyVXbRvfAVoasNS7X6yn2df75/dQEQNxz9nPIW8S8jWgCVINp1T76Cl4CtLTfPuXgh0Gu6Y1yPbGbqEGvgzENoSKMqRjv/1ou33CAgBIq4Sd982HhxORKSCAYunULSwBJkvpyZQtvqhvqE7s9QoKwR2fAlccycQGguUFZiWYbDiXuzVJaQNkdK3fQ50u7MedfEBHsoU6xLopUsREJFdsUvIU8gDFo2N/2zyLxh7BSw7FgBrXgTe6wQcWwdcMVoS3YexsNvR1CxDlbXIcjlLcvcox0HVNbXeksBQ5ZRmIiIL+K3iKeRdQk1s3NZd/qVijzU4AGDrp1L661FAU6NVie0VGJH9nK1ZRr92b6f6mNtfeWy8PhARkYOwhcVTyFtYmkbZdm2MbEdR3Vn71Ed3RnlcrNxt2+opruQ8ZYX2vydb0ojISRiweAp5wGLrBnG+/lKrzPo3G16XCpU1OOS63gE88FvDX4fsS96VuPr5+q3J0yJRedyyd8PqRERkJf555CnkTe/1aYYvOS+l9dXSAMz6KL1k/tzDG4DY7vW/NzmO/N98y2wgcZi49oktqmSzwYa91bAxLERENmALi6doyKBbY/uWNex6S2u5MFhxX8bdN1fq2A/IWGUpUFCzE7s2FOhxr33qRURkBQYsnsJH1qrSkNYRwHS8ia3MjYVoaCBFjqUx+rmxdT0W+caEU7Yr9wAiInIwfsN4Cnk3kLYe61ZEyhbnauh+QoaARSNuZFeLAzDdm3GgWzvN2VrysUsh0ebLERE5AAMWT6HRAMPfAwbPACISbL9+5GzZgZ0ClqtuAG75AEh9WTy+/fOG3ZccyzhgsXVJ/Ipi8bmtjeNeiIjsgH8Se5JrH6j/tWGtpHR1BfDtXeKOuYPqsTvu5ZPic+0WAf0eF8czNGle//qR4xl3CfnZGLDU7tBs7W7hRER2xBaWxkLepXT4d+DIb0Dm6/W714GfxefCmrVYNBoGK57AuIXF1sX9KkrEZwYsROQCbGFpLBR/Tcu6hMqLbZuaeuUSkL9fTLfiGhwexWTQrRUBS1UF8GrNGj69JojPDFiIyAXYwtJYyL9kSmTTWeUzhorygFXPAecPm7/Pt3dI6VbX2q9+5Hj1aWHZ8rGU3vGF+OzPgIWInI8BS2OSNFZ8lu/WWymb+fFuR2DrHOCLoebvUbsOBwAEcNEwj2ISsFQDp7cBH10LHM1Qv+bSMdM8trAQkQswYGlMaruFimWr3laVm5a7chH4uC+w+WPTc/IVc+W7QJP7M+4SOvM38PVtwIXDwDe3qV+jtvQ+AxYicgEGLI1J7R5ElSVSXmWpetn8fcBv0wG9hd2d+cXlWYxbWLbOBSqKLF9TrnKeLWtE5AIMWBoTtWmstQNo/56nfk3pZSltHNzYugkjuZatKxFv/RT4/f9M87nCLRG5AAOWxkRtobBVz4jPv05Tv0a+0eG3d0rpq28Boq42LU9uzGhl284jLRdf9bR6PlvWiMgFOK25MfGvR4tIhaz76ORGKT3664bXh5zLeCl+W/cSqsUuISJyAbawNCYRV6nnV1uY3iqfRUTe5eCK+l3HFhYicgEGLI1JULh6fmG2ej6gbGFJuF58rp0eTR7Gxs0OzfHnGBYicj4GLI2JuUGyH/Qwf03thneANMuEm995phYd7HMfdgkRkQswYGlMfOoxZKmiRFyOf14acGytmGfLUv7kPlKm2Oc+7BIiIhdgwNKY1CtguQJkvgGc3iLlhcXbr07kPP6BQLO21pcPkC0MeNM7snx2CRGR8zFgaUysCVhuN1qPpaIYKMlX5jVvb786kXvS66VF5R7LUp7jXkJE5AL1Clhmz56Ntm3bIjAwEMnJydi2bZvZsoMGDYJGozF5DB8+3FDmvvvuMzk/dKiF/WyofoxXOo2/Tnn85GHgmjuATjdLeRUlgDZUWY6DLj2X8fL8cod/E3dnBgDdGSnfPxi49kGxZS3pHsCXqyEQkfPZHLAsWbIE06ZNw8yZM7Fz5050794daWlpyM/PVy2/dOlS5OTkGB579+6Fr68v7rzzTkW5oUOHKsotWrSofu+IzDNuYQk0CkQCw8Tnu78F+j0upiuvmAY6PmyY81jG/5ZyC+8C1r8ppj+5XsoPiRbXcHliLzBqtmPrR0Rkhs3fPO+99x4mTpyICRMmoHPnzpg7dy6Cg4Mxf/581fIRERGIiYkxPNasWYPg4GCTgEWr1SrKNWvWrH7viMwz/rKSz/bwC1QuLFd7bsvHytlFTaIcVz9yvOuftXx+1zfis3yFYyIiN2BTwFJRUYEdO3YgNTVVuoGPD1JTU7F582ar7jFv3jzcfffdaNJE2Q+emZmJqKgoJCYmYtKkSbh48aItVSNrGHcH+AdJaeOZH1kLpXSZTkp3v9v+9SLnueYO4N/fmT8v6C1veElE5CI2dUZfuHAB1dXViI6OVuRHR0fj4MGDdV6/bds27N27F/PmKQd2Dh06FLfddhsSEhJw7NgxPP/88xg2bBg2b94MX1/TJuzy8nKUl5cbjnU6nUkZUmHcJSRvOfHxV56Tr3BbnCulIxLsXy9yrmaW/g0FIGeX06pCRGQtp46emzdvHq655hr06dNHkX/33dJf7ddccw26deuGq666CpmZmRg8eLDJfdLT0/Hyyy87vL5exzhgKTlv/tzdC4HPaz77o3+IzwEhQI97HVc/cg5Lg2YFwfZdnYmInMCm30wtWrSAr68v8vLyFPl5eXmIiYmxeG1JSQkWL16MBx54oM7XadeuHVq0aIGjR4+qnp8+fToKCwsNj9OnT1v/Jhoz+RiWnuOB45nq5wCgVW/T61NnAr7+pvnkWSxObxfqt14PEZGD2RSwBAQEoFevXsjIyDDk6fV6ZGRkICUlxeK133//PcrLy3HPPffU+TpnzpzBxYsXERsbq3peq9UiNDRU8SAryIOSdtdLs4IA5Z5B5vhp7V8ncj5LAUllGXBsnXTMfaOIyE3Y3PY7bdo0fPbZZ/jyyy9x4MABTJo0CSUlJZgwYQIAYNy4cZg+fbrJdfPmzcOoUaPQvHlzRX5xcTGefvppbNmyBSdPnkRGRgZGjhyJ9u3bIy0trZ5vi1TJv6g0PsoA5MqFuq/3ZcDiFYzHK8lVlgBrXpSOh77h+PoQEVnB5rbf0aNH4/z585gxYwZyc3ORlJSE1atXGwbiZmdnw8donY5Dhw5h06ZN+P33303u5+vri927d+PLL79EQUEB4uLiMGTIEPz3v/+FVssvSLuSzxLS+FpeREyNX4B960OuIW9p6z4G+MfCmkfGa/UQEblIvTqrp0yZgilT1DdSy8zMNMlLTEyEIAiq5YOCgvDbb7/VpxpkK+MWlu6jgYxXzJeP6gzk75eO2cLiHeQ/B237Ww5YiIjcBKcDNCbyv6w1PsA1d1kub9wCwzEs3kE+cFrQA4NnSMfc2JKI3BQDlsZEo5HSPr5AiPqgZoNONymPOUPIO8hbWAQBGPAkcH9Nd22hbMZdotG/PxGRCzFgaaxCW9a9id2AJ4Gb3pGO9dWOrRM5h2KWUE1XrVow2iTSKdUhIrIGA5bG5tZPgSGvAjFd6y7rpxV36a0lMGDxCvKWttqxZWrdffJyREQuxhWiGpvuo5XHvgFAdYX58vzS8nK1LSwqM8C44i0RuRH+Rmrs1L6ojPUcB8T1ANoOdHx9yLkEC11CMdc4ty5ERBawhaWxs2Yg7S0fOr4e5BpCzc7MalPWe453bl2IiCxgC0tj12m4+NwkyrX1INdSa2kz3l+KiMiF2MLS2A1+SZze3KPuPZ7Ii3HKOhG5OQYsjV3TSOBfL7i6FuQqhllCga6tBxFRHdglRNSYxXYTn41bWNhFSERuhi0sRI3R5L+BS8eA1teJx8bT13ve6/w6ERFZwICFqDGK7Cg+zGEXERG5GXYJEZGpwDBX14CISIEBCxGZ8g9ydQ2IiBQYsBCRKEk2td2H05yJyL0wYCEikXz3bq7LQkRuhgELEYl8GLAQkftiwEJEIkXAYsWmmERETsSAhYhE8oCFY1iIyM0wYCEiEbuEiMiNMWAhIhEDFiJyYwxYiEgkD1I4hoWI3AwDFiIS+fhKabawEJGbYcBCRCKNLGDhoFsicjMMWIhIdGytlGaXEBG5GQYsRCQ6uVFK+3IjdyJyLwxYiEjUvIOUZgsLEbkZBixEJBr6hpTmGBYicjMMWIhIFBQupTlLiIjcDAMWIhJxWjMRuTEGLEQk4uaHROTG6hWwzJ49G23btkVgYCCSk5Oxbds2s2UXLFgAjUajeAQGBirKCIKAGTNmIDY2FkFBQUhNTcWRI0fqUzUiqq8mUVLah7OEiMi92BywLFmyBNOmTcPMmTOxc+dOdO/eHWlpacjPzzd7TWhoKHJycgyPU6dOKc6/9dZb+OCDDzB37lxs3boVTZo0QVpaGsrKymx/R0RUPyHRwN2LgHE/ARqNq2tDRKRgc8Dy3nvvYeLEiZgwYQI6d+6MuXPnIjg4GPPnzzd7jUajQUxMjOERHR1tOCcIAmbNmoUXXngBI0eORLdu3fDVV1/h3LlzWL58eb3eFBHVU6ebgHaDXF0LIiITNgUsFRUV2LFjB1JTU6Ub+PggNTUVmzdvNntdcXEx2rRpg/j4eIwcORL79u0znDtx4gRyc3MV9wwLC0NycrLFexIREVHjYVPAcuHCBVRXVytaSAAgOjoaubm5qtckJiZi/vz5+Omnn/DNN99Ar9ejb9++OHPmDAAYrrPlnuXl5dDpdIoHEREReS+HzxJKSUnBuHHjkJSUhOuvvx5Lly5FZGQkPvnkk3rfMz09HWFhYYZHfHy8HWtMRERE7samgKVFixbw9fVFXl6eIj8vLw8xMTFW3cPf3x89evTA0aNHAcBwnS33nD59OgoLCw2P06dP2/I2iIiIyMPYFLAEBASgV69eyMjIMOTp9XpkZGQgJSXFqntUV1djz549iI2NBQAkJCQgJiZGcU+dToetW7eavadWq0VoaKjiQURERN7L5sUWpk2bhvHjx6N3797o06cPZs2ahZKSEkyYMAEAMG7cOLRs2RLp6ekAgFdeeQXXXXcd2rdvj4KCArz99ts4deoUHnzwQQDiDKKpU6fi1VdfRYcOHZCQkIAXX3wRcXFxGDVqlP3eKREREXksmwOW0aNH4/z585gxYwZyc3ORlJSE1atXGwbNZmdnw8dHari5fPkyJk6ciNzcXDRr1gy9evXCX3/9hc6dOxvKPPPMMygpKcFDDz2EgoIC9O/fH6tXrzZZYI6IiIgaJ40gCIKrK9FQOp0OYWFhKCwsZPcQERGRh7Dl+5t7CREREZHbY8BCREREbo8BCxEREbk9BixERETk9hiwEBERkduzeVqzO6qd6MQ9hYiIiDxH7fe2NROWvSJgKSoqAgDuKUREROSBioqKEBYWZrGMV6zDotfrce7cOYSEhECj0dj13jqdDvHx8Th9+jTXeKknfob2wc+x4fgZ2gc/x4bjZygSBAFFRUWIi4tTLDqrxitaWHx8fNCqVSuHvgb3LGo4fob2wc+x4fgZ2gc/x4bjZ4g6W1ZqcdAtERERuT0GLEREROT2GLDUQavVYubMmdBqta6uisfiZ2gf/Bwbjp+hffBzbDh+hrbzikG3RERE5N3YwkJERERujwELERERuT0GLEREROT2GLAQERGR22PAUofZs2ejbdu2CAwMRHJyMrZt2+bqKrnEhg0bMGLECMTFxUGj0WD58uWK84IgYMaMGYiNjUVQUBBSU1Nx5MgRRZlLly5h7NixCA0NRXh4OB544AEUFxcryuzevRsDBgxAYGAg4uPj8dZbbzn6rTlNeno6rr32WoSEhCAqKgqjRo3CoUOHFGXKysowefJkNG/eHE2bNsXtt9+OvLw8RZns7GwMHz4cwcHBiIqKwtNPP42qqipFmczMTPTs2RNarRbt27fHggULHP32nGbOnDno1q2bYcGtlJQUrFq1ynCen6Ht3njjDWg0GkydOtWQx8+xbi+99BI0Go3i0alTJ8N5foZ2JpBZixcvFgICAoT58+cL+/btEyZOnCiEh4cLeXl5rq6a061cuVL4v//7P2Hp0qUCAGHZsmWK82+88YYQFhYmLF++XPjnn3+EW265RUhISBBKS0sNZYYOHSp0795d2LJli7Bx40ahffv2wpgxYwznCwsLhejoaGHs2LHC3r17hUWLFglBQUHCJ5984qy36VBpaWnCF198Iezdu1fIysoSbrrpJqF169ZCcXGxocwjjzwixMfHCxkZGcL27duF6667Tujbt6/hfFVVldC1a1chNTVV2LVrl7By5UqhRYsWwvTp0w1ljh8/LgQHBwvTpk0T9u/fL3z44YeCr6+vsHr1aqe+X0f5+eefhV9//VU4fPiwcOjQIeH5558X/P39hb179wqCwM/QVtu2bRPatm0rdOvWTXj88ccN+fwc6zZz5kyhS5cuQk5OjuFx/vx5w3l+hvbFgMWCPn36CJMnTzYcV1dXC3FxcUJ6eroLa+V6xgGLXq8XYmJihLffftuQV1BQIGi1WmHRokWCIAjC/v37BQDC33//bSizatUqQaPRCGfPnhUEQRA+/vhjoVmzZkJ5ebmhzLPPPiskJiY6+B25Rn5+vgBAWL9+vSAI4mfm7+8vfP/994YyBw4cEAAImzdvFgRBDBx9fHyE3NxcQ5k5c+YIoaGhhs/tmWeeEbp06aJ4rdGjRwtpaWmOfksu06xZM+Hzzz/nZ2ijoqIioUOHDsKaNWuE66+/3hCw8HO0zsyZM4Xu3burnuNnaH/sEjKjoqICO3bsQGpqqiHPx8cHqamp2Lx5swtr5n5OnDiB3NxcxWcVFhaG5ORkw2e1efNmhIeHo3fv3oYyqamp8PHxwdatWw1lBg4ciICAAEOZtLQ0HDp0CJcvX3bSu3GewsJCAEBERAQAYMeOHaisrFR8jp06dULr1q0Vn+M111yD6OhoQ5m0tDTodDrs27fPUEZ+j9oy3vhzW11djcWLF6OkpAQpKSn8DG00efJkDB8+3OS98nO03pEjRxAXF4d27dph7NixyM7OBsDP0BEYsJhx4cIFVFdXK36QACA6Ohq5ubkuqpV7qv08LH1Wubm5iIqKUpz38/NDRESEoozaPeSv4S30ej2mTp2Kfv36oWvXrgDE9xgQEIDw8HBFWePPsa7PyFwZnU6H0tJSR7wdp9uzZw+aNm0KrVaLRx55BMuWLUPnzp35Gdpg8eLF2LlzJ9LT003O8XO0TnJyMhYsWIDVq1djzpw5OHHiBAYMGICioiJ+hg7gFbs1E3mayZMnY+/evdi0aZOrq+KREhMTkZWVhcLCQvzwww8YP3481q9f7+pqeYzTp0/j8ccfx5o1axAYGOjq6nisYcOGGdLdunVDcnIy2rRpg++++w5BQUEurJl3YguLGS1atICvr6/JiO68vDzExMS4qFbuqfbzsPRZxcTEID8/X3G+qqoKly5dUpRRu4f8NbzBlClTsGLFCqxbtw6tWrUy5MfExKCiogIFBQWK8safY12fkbkyoaGhXvNLNCAgAO3bt0evXr2Qnp6O7t274/333+dnaKUdO3YgPz8fPXv2hJ+fH/z8/LB+/Xp88MEH8PPzQ3R0ND/HeggPD0fHjh1x9OhR/iw6AAMWMwICAtCrVy9kZGQY8vR6PTIyMpCSkuLCmrmfhIQExMTEKD4rnU6HrVu3Gj6rlJQUFBQUYMeOHYYya9euhV6vR3JysqHMhg0bUFlZaSizZs0aJCYmolmzZk56N44jCAKmTJmCZcuWYe3atUhISFCc79WrF/z9/RWf46FDh5Cdna34HPfs2aMI/tasWYPQ0FB07tzZUEZ+j9oy3vxzq9frUV5ezs/QSoMHD8aePXuQlZVlePTu3Rtjx441pPk52q64uBjHjh1DbGwsfxYdwdWjft3Z4sWLBa1WKyxYsEDYv3+/8NBDDwnh4eGKEd2NRVFRkbBr1y5h165dAgDhvffeE3bt2iWcOnVKEARxWnN4eLjw008/Cbt37xZGjhypOq25R48ewtatW4VNmzYJHTp0UExrLigoEKKjo4V7771X2Lt3r7B48WIhODjYa6Y1T5o0SQgLCxMyMzMV0yCvXLliKPPII48IrVu3FtauXSts375dSElJEVJSUgzna6dBDhkyRMjKyhJWr14tREZGqk6DfPrpp4UDBw4Is2fP9qppkM8995ywfv164cSJE8Lu3buF5557TtBoNMLvv/8uCAI/w/qSzxISBH6O1njyySeFzMxM4cSJE8Kff/4ppKamCi1atBDy8/MFQeBnaG8MWOrw4YcfCq1btxYCAgKEPn36CFu2bHF1lVxi3bp1AgCTx/jx4wVBEKc2v/jii0J0dLSg1WqFwYMHC4cOHVLc4+LFi8KYMWOEpk2bCqGhocKECROEoqIiRZl//vlH6N+/v6DVaoWWLVsKb7zxhrPeosOpfX4AhC+++MJQprS0VHj00UeFZs2aCcHBwcKtt94q5OTkKO5z8uRJYdiwYUJQUJDQokUL4cknnxQqKysVZdatWyckJSUJAQEBQrt27RSv4enuv/9+oU2bNkJAQIAQGRkpDB482BCsCAI/w/oyDlj4OdZt9OjRQmxsrBAQECC0bNlSGD16tHD06FHDeX6G9qURBEFwTdsOERERkXU4hoWIiIjcHgMWIiIicnsMWIiIiMjtMWAhIiIit8eAhYiIiNweAxYiIiJyewxYiIiIyO0xYCEiIiK3x4CFiIiI3B4DFiIiInJ7DFiIiIjI7TFgISIiIrf3/656oPwgoyPvAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the models predictions vs the time series itself\n", + "import matplotlib.pyplot as plt\n", + "\n", + "plt.plot(crash_indices, crash_vals, '+')\n", + "plt.plot(time_series - 6)" ] }, { @@ -619,12 +1080,173 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "id": "0de4a8f2", "metadata": {}, "outputs": [], "source": [ - "# write your code here" + "import yfinance as yf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf26ce9c", + "metadata": {}, + "outputs": [], + "source": [ + "start_date = \"2019-01-01\"\n", + "end_date = \"2021-12-31\"\n", + "dimension = 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b117bab", + "metadata": {}, + "outputs": [], + "source": [ + "sector_etfs = [\"XLY\", \"XLK\", \"XLC\", \"XLV\", \"XLF\", \"XLP\"]\n", + "sector_data = {}\n", + "\n", + "overall_max = 0\n", + "for sector in sector_etfs:\n", + " data = yf.download(sector, start=start_date, end=end_date)\n", + " data = data['Close'][sector].to_numpy()\n", + " sector_data[sector] = data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "811e1144", + "metadata": {}, + "outputs": [], + "source": [ + "for sector in sector_etfs:\n", + " plt.plot(sector_data[sector], label=sector)\n", + "plt.title(\"Sector ETF Prices\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5cccd5c", + "metadata": {}, + "outputs": [], + "source": [ + "plt.clf()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff9b65ea", + "metadata": {}, + "outputs": [], + "source": [ + "def verma_embed_series(time_series: np.ndarray, dimension: int):\n", + " embedded = [\n", + " time_series[i:(i + dimension)] for i in range(len(time_series) - dimension)\n", + " ]\n", + " return embedded\n", + "\n", + "for sector in sector_etfs:\n", + " sector_data[sector] = verma_embed_series(sector_data[sector], dimension)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2fca0ed", + "metadata": {}, + "outputs": [], + "source": [ + "def get_sector_complex(timeslice):\n", + " complex = []\n", + " for sector in sector_etfs:\n", + " complex.append(sector_data[sector][timeslice])\n", + " return np.stack(complex).T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "834843be", + "metadata": {}, + "outputs": [], + "source": [ + "sector_complex = get_sector_complex(1)\n", + "sector_complex.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32681a99", + "metadata": {}, + "outputs": [], + "source": [ + "betti_curves = []\n", + "for timeslice in range(len(sector_data[\"XLC\"])):\n", + " betti_curve = []\n", + " sector_complex = get_sector_complex(timeslice)\n", + " for epsilon in range(1, 50):\n", + " betti = classical_betti_solver(sector_complex, epsilon / 75, 0)\n", + " betti_curve.append(betti)\n", + " betti_curves.append(betti_curve)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c9f717f", + "metadata": {}, + "outputs": [], + "source": [ + "betti_curves = np.array(betti_curves)\n", + "distances = betti_curves[1:] - betti_curves[:-1]\n", + "distances = np.linalg.norm(distances, axis=1, ord=6)\n", + "\n", + "max_distance_location = np.argmax(distances)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2da4f748", + "metadata": {}, + "outputs": [], + "source": [ + "spy_data = yf.download(\"SPY\", start=start_date, end=end_date)\n", + "spy_data = spy_data['Close'][\"SPY\"].to_numpy()\n", + "spy_data = spy_data[:-dimension]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "532c1f6b", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(max_distance_location, spy_data[max_distance_location], 'ro', label=\"Max Betti Curve Distance\") \n", + "plt.plot(spy_data)\n", + "\n", + "date_index = pd.date_range(start_date, end_date)\n", + "date_index = [date.strftime(\"%Y-%m-%d\") for date in date_index]\n", + "date_index = date_index[:-dimension]\n", + "date_index = [date_index[i] for i in range(0, len(date_index), 100)]\n", + "\n", + "plt.xticks(np.linspace(0, len(spy_data), len(date_index)), date_index)\n", + "plt.xticks(rotation=70)\n", + "plt.xlabel(\"Time\")\n", + "plt.ylabel(\"S&P Price\")\n", + "plt.title(\"Sector Level Betti Analysis\")\n", + "plt.legend()" ] }, { @@ -638,7 +1260,7 @@ ], "metadata": { "kernelspec": { - "display_name": "tda", + "display_name": "moodys", "language": "python", "name": "python3" }, @@ -652,7 +1274,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.11" + "version": "3.11.9" } }, "nbformat": 4,