\n",
+ " \n",
+ "BONUS EXERCISE: \n",
+ "\n",
+ "`gtda.time_series.TakensEmbedding` can conduct this transformation. Try avoid using this function and build your own embedding function.\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9780b90c",
+ "metadata": {},
+ "source": [
+ " NOTE:\n",
+ "\n",
+ "From step 2, we present an example originally provided in the appendix of [Khandelwal's and Chandra's paper](https://arxiv.org/abs/2302.09553). This example can be used to verify your code.\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c945e850",
+ "metadata": {},
+ "source": [
+ "### Step 2: Building Laplacian\n",
+ "\n",
+ "**Instruction:**\n",
+ "\n",
+ "In this step, we are going to use `GUDHI`, a Python package specialized in TDA and and higher dimensional geometry understanding. For detailed instructions on the functions we will use, please refer to [this website](https://gudhi.inria.fr/python/latest/rips_complex_user.html#).\n",
+ "\n",
+ "A simplicial complex is constructed on the point cloud using `gudhi.RipsComplex` class, followed by its `create_simplex_tree` method. The resolution threshold `epsilon` is set via the `max_edge_length` parameter. This process identifies the connectivity of the complex within the resolution threshold and produces a simplex tree. The simplex tree serves as a general representation of simplicial complexes. Using its `get_filtration` method, simplicies are retrieved as a collection of lists, where elements are grouped based on their connections. Each dimension up to a specified maximum is represented by its respective collection of lists.\n",
+ "\n",
+ "**Example:**\n",
+ "\n",
+ "Here is an example of a simplex tree $\\mathcal{K}$ with a maximum dimension of 2. In its zeroth dimension, each point is a connected component; In the first dimension, 6 line segments connect 6 pairs of points $[1, 2], [1, 3], [2, 3], [3, 4], [3, 5], [4, 5]$; In the second dimension, a filled trangle is formed among points $[1 ,2 ,3]$.\n",
+ "\n",
+ "$$\n",
+ "\\mathcal{K} = [[[1], [2], [3], [4], [5]],[[1, 2], [1, 3], [2, 3], [3, 4], [3, 5], [4, 5]], [[1, 2, 3]]]\n",
+ "$$\n",
+ "\n",
+ "**Action:**\n",
+ "\n",
+ "Build a simplicial complex by applying functions from `GUDHI` on the point cloud obtained in step 1, and extract its simplex tree. It is recommended to store the simplex tree in a format similar to the example provided, i.e., in the format $[S_0, S_1, S_2, \\dots]$, where $S_i$ represents the set of $i$-simplices.\n",
+ "\n",
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "026d6b61",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import gudhi as gd"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "4db171cb-31db-4aba-8bad-12b4b171e81a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "point_clouds_numpy = np.array(point_clouds)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "6d4bb081-4044-42a1-9272-ba8f94b7ea79",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(64, 5, 4)"
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "point_clouds_numpy.shape"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "aed04d25-c844-433f-bcfa-b7e0b2f9e788",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 21 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 9 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 10 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 9 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 8 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 9 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 8 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 9 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 8 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 7 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 8 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 12 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 18 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 1 - 8 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 10 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 15 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 18 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 16 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 25 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 18 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 16 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 16 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 13 simplices - 5 vertices.\n",
+ "Rips complex is of dimension 2 - 11 simplices - 5 vertices.\n"
+ ]
+ }
+ ],
+ "source": [
+ "simplex_trees = []\n",
+ "for point_cloud in point_clouds_numpy:\n",
+ " rips_complex = gd.RipsComplex(points=point_cloud, max_edge_length=0.1)\n",
+ " simplex_tree = rips_complex.create_simplex_tree(max_dimension=2)\n",
+ " simplex_trees.append(simplex_tree)\n",
+ "\n",
+ " result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \\\n",
+ " repr(simplex_tree.num_simplices()) + ' simplices - ' + \\\n",
+ " repr(simplex_tree.num_vertices()) + ' vertices.'\n",
+ " print(result_str)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "id": "8d4ed0ad-e410-4d1d-a51c-ce3ec1551189",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "64"
+ ]
+ },
+ "execution_count": 17,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "len(simplex_trees)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "id": "442fa468-5d84-4bf9-9d94-bca195da83da",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[0] -> 0.00\n",
+ "[1] -> 0.00\n",
+ "[2] -> 0.00\n",
+ "[3] -> 0.00\n",
+ "[4] -> 0.00\n",
+ "[3, 4] -> 0.01\n",
+ "[2, 3] -> 0.01\n",
+ "[1, 2] -> 0.02\n",
+ "[2, 4] -> 0.02\n",
+ "[2, 3, 4] -> 0.02\n",
+ "[1, 3] -> 0.02\n",
+ "[1, 2, 3] -> 0.02\n",
+ "[0, 1] -> 0.02\n",
+ "[1, 4] -> 0.02\n",
+ "[1, 2, 4] -> 0.02\n",
+ "[1, 3, 4] -> 0.02\n",
+ "[0, 4] -> 0.02\n",
+ "[0, 1, 4] -> 0.02\n",
+ "[0, 3] -> 0.03\n",
+ "[0, 1, 3] -> 0.03\n",
+ "[0, 3, 4] -> 0.03\n",
+ "[0, 2] -> 0.03\n",
+ "[0, 1, 2] -> 0.03\n",
+ "[0, 2, 3] -> 0.03\n",
+ "[0, 2, 4] -> 0.03\n"
+ ]
+ }
+ ],
+ "source": [
+ "fmt = '%s -> %.2f'\n",
+ "for filtered_value in simplex_trees[0].get_filtration():\n",
+ " print(fmt % tuple(filtered_value))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "fa7694cf-0c55-4078-af76-808fa23d1578",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def getSimplexTrees(point_clouds, epsilon):\n",
+ " simplex_trees = []\n",
+ " for point_cloud in point_clouds:\n",
+ " rips_complex = gd.RipsComplex(points=point_cloud, max_edge_length=epsilon)\n",
+ " simplex_tree = rips_complex.create_simplex_tree(max_dimension=2)\n",
+ " simplex_trees.append(simplex_tree)\n",
+ "\n",
+ " return simplex_trees"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "573efefa",
+ "metadata": {},
+ "source": [
+ "**Instruction:**\n",
+ "\n",
+ "Let $S_k$ be the set of $k$-simplicies in the complex $\\mathcal{K}$ with individual simplicies denoted by $s_k ∈ S_k$ written as $[j_0, j_1, \\dots, j_k]$ where $j_i$ is the $i$-th vertex of $s_k$. Note that the vertices are ordered in ascending fashion in the initial point cloud, and this order is kept throughout. The restricted boundary operator $\\partial_k$ is defined on the $k$-simplicies as:\n",
+ "\n",
+ "$$\n",
+ "\\begin{align*}\n",
+ "\\partial_k s_k &= \\sum_{t=0}^{k} (-1)^t [v_0, \\dots, v_{t-1}, v_{t+1}, \\dots, v_k]\\\\\n",
+ "&= \\sum_{t=0}^{k} (-1)^t s_{k-1} (t)\n",
+ "\\end{align*}\n",
+ "$$\n",
+ "\n",
+ "where $s_{k−1}(t)$ is defined as the lower simplex defined from $s_k$ by leaving out the vertex $v_t$. \n",
+ "\n",
+ "**Example:**\n",
+ "\n",
+ "In the simplex tree $\\mathcal{K}$ we have 1 2-simplex and 6 1-simplicies. By leaving out vertice $v_0=1$, $v_1=2$, $v_2=3$, we obtain the lower simplex $s_1=[2, 3]$, $s_2=[1, 3]$, $s_3=[1, 2]$, respectively. Therefore, the boundary operator on the 2-simplicies $\\partial_2$ should be a 6-by-1 matrix:\n",
+ "\n",
+ "$$\n",
+ "\\partial_2 =\n",
+ "\\begin{bmatrix}\n",
+ "1 \\\\\n",
+ "-1 \\\\\n",
+ "1 \\\\\n",
+ "0 \\\\\n",
+ "0 \\\\\n",
+ "0\n",
+ "\\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "In the same way, the boundary operator on the 1-simplicies $\\partial_1$ is:\n",
+ "\n",
+ "$$\n",
+ "\\partial_1 =\n",
+ "\\begin{bmatrix}\n",
+ "1 & 1 & 0 & 0 & 0 & 0 \\\\\n",
+ "-1 & 0 & 1 & 0 & 0 & 0 \\\\\n",
+ "0 & -1 & -1 & 1 & 1 & 0 \\\\\n",
+ "0 & 0 & 0 & -1 & 0 & 1 \\\\\n",
+ "0 & 0 & 0 & 0 & -1 & -1\n",
+ "\\end{bmatrix}\n",
+ "$$\n",
+ "\n",
+ "**Action:**\n",
+ "\n",
+ "Define a function that generates the boundary operator for a specified dimension, using a given simplex tree as input.\n",
+ "\n",
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "20b9d9c3-830a-4f6b-9bcb-6d35e2dd7dde",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def boundary_operator_matrix(simplex_tree, k):\n",
+ " # Get all k-simplices and (k-1)-simplices\n",
+ " k_simplices = [s for s, _ in simplex_tree.get_simplices() if len(s) == k + 1]\n",
+ " k_minus_1_simplices = [s for s in simplex_tree.get_simplices() if len(s[0]) == k]\n",
+ " \n",
+ " # Create a dictionary to map (k-1)-simplices to their index\n",
+ " k_minus_1_dict = {tuple(s[0]): i for i, s in enumerate(k_minus_1_simplices)}\n",
+ " \n",
+ " # Initialize the boundary matrix\n",
+ " boundary_matrix = np.zeros((len(k_minus_1_simplices), len(k_simplices)))\n",
+ " \n",
+ " # Fill the boundary matrix\n",
+ " for j, simplex in enumerate(k_simplices):\n",
+ " for t in range(k + 1):\n",
+ " # Create the (k-1)-simplex by removing the t-th vertex\n",
+ " face = tuple(v for i, v in enumerate(simplex) if i != t)\n",
+ " if face in k_minus_1_dict:\n",
+ " i = k_minus_1_dict[face]\n",
+ "\n",
+ " if k % 2 == 0:\n",
+ " boundary_matrix[i, j] = (-1) ** (t)\n",
+ " else:\n",
+ " boundary_matrix[i, j] = (-1) ** (t + 1)\n",
+ " \n",
+ " return boundary_matrix"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "003595dd-2644-439c-a6df-18b1541871b2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rips_complex_empty = gd.RipsComplex(points=[], max_edge_length=0.5)\n",
+ "simplex_tree_empty = rips_complex_empty.create_simplex_tree(max_dimension=2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "7c4ad409-5a43-4a3d-bd38-6cf32f693a50",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "True"
+ ]
+ },
+ "execution_count": 22,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "simplex_tree_empty.insert([1, 2, 3])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "id": "a5f05d73-3054-4c52-8154-f7977e027e30",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[([1], 0.0),\n",
+ " ([2], 0.0),\n",
+ " ([1, 2], 0.0),\n",
+ " ([3], 0.0),\n",
+ " ([1, 3], 0.0),\n",
+ " ([2, 3], 0.0),\n",
+ " ([1, 2, 3], 0.0)]"
+ ]
+ },
+ "execution_count": 23,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "list(simplex_tree_empty.get_filtration())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "id": "20d7d78a-6f26-491e-a0ea-ae290d36894a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 1., 1., 0.],\n",
+ " [-1., 0., 1.],\n",
+ " [ 0., -1., -1.]])"
+ ]
+ },
+ "execution_count": 24,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "boundary_matrix_empty1 = boundary_operator_matrix(simplex_tree_empty, 1)\n",
+ "boundary_matrix_empty1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "id": "b061587a-fb8a-4db6-8304-574e16252b45",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 1.],\n",
+ " [-1.],\n",
+ " [ 1.]])"
+ ]
+ },
+ "execution_count": 25,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "boundary_matrix_empty2 = boundary_operator_matrix(simplex_tree_empty, 2)\n",
+ "boundary_matrix_empty2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "id": "d0276f5f-69ea-42e8-a8d6-2764c121d153",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[3., 0., 0.],\n",
+ " [0., 3., 0.],\n",
+ " [0., 0., 3.]])"
+ ]
+ },
+ "execution_count": 26,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "laplacian_empty = np.dot(boundary_matrix_empty1.conj().T, boundary_matrix_empty1) + np.dot(boundary_matrix_empty2, boundary_matrix_empty2.conj().T)\n",
+ "laplacian_empty"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "id": "1ab25f8f-f366-4c03-a77c-0f94728dd368",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([3., 3., 3.])"
+ ]
+ },
+ "execution_count": 27,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "eigenvalues_empty = np.linalg.eig(laplacian_empty)\n",
+ "eigenvalues_empty.eigenvalues"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
+ "id": "2ebf4957-a814-47ae-89c3-2a02c088d65a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "simplex_tree_empty.compute_persistence()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "id": "321ebd85-2c15-43d0-bf65-02c6cc0ebc0c",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "[1, 0]"
+ ]
+ },
+ "execution_count": 29,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "simplex_tree_empty.betti_numbers()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 30,
+ "id": "1152db2d-1848-4927-97ba-87dc625b1514",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([], shape=(0, 5), dtype=float64)"
+ ]
+ },
+ "execution_count": 30,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "boundary_matrix = boundary_operator_matrix(simplex_trees[0], 0)\n",
+ "boundary_matrix"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "id": "edb97d46-030e-4dce-bb27-853e4a524da6",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 1., 1., 1., 1., 0., 0., 0., 0., 0., 0.],\n",
+ " [-1., 0., 0., 0., 1., 1., 1., 0., 0., 0.],\n",
+ " [ 0., -1., 0., 0., -1., 0., 0., 1., 1., 0.],\n",
+ " [ 0., 0., -1., 0., 0., -1., 0., -1., 0., 1.],\n",
+ " [ 0., 0., 0., -1., 0., 0., -1., 0., -1., -1.]])"
+ ]
+ },
+ "execution_count": 31,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "boundary_matrix2 = boundary_operator_matrix(simplex_trees[0], 1)\n",
+ "boundary_matrix2"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "id": "246e41ce",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[ 4. -1. -1. -1. -1.]\n",
+ " [-1. 4. -1. -1. -1.]\n",
+ " [-1. -1. 4. -1. -1.]\n",
+ " [-1. -1. -1. 4. -1.]\n",
+ " [-1. -1. -1. -1. 4.]]\n",
+ "\n",
+ "[5. 0. 5. 5. 5.]\n"
+ ]
+ }
+ ],
+ "source": [
+ "laplacian = np.dot(boundary_matrix.conj().T, boundary_matrix) + np.dot(boundary_matrix2, boundary_matrix2.conj().T)\n",
+ "eigenvalues = np.linalg.eig(laplacian)\n",
+ "\n",
+ "print(laplacian)\n",
+ "print()\n",
+ "print(eigenvalues.eigenvalues)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "id": "6c36d802-8ca3-4150-9bc7-abb7c5c31226",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_laplacian(simplex_tree, k):\n",
+ " boundary_matrix = boundary_operator_matrix(simplex_tree, k)\n",
+ " boundary_matrix2 = boundary_operator_matrix(simplex_tree, k + 1)\n",
+ " laplacian = np.dot(boundary_matrix.conj().T, boundary_matrix) + np.dot(boundary_matrix2, boundary_matrix2.conj().T)\n",
+ " eigenvalues = np.linalg.eig(laplacian)\n",
+ "\n",
+ " return laplacian, eigenvalues.eigenvalues"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "81764f0d",
+ "metadata": {},
+ "source": [
+ "**Instruction:**\n",
+ "\n",
+ "The combinatorial laplacian $\\Delta_k$ is defined as:\n",
+ "\n",
+ "$$\n",
+ "\\Delta_k = \\left( \\partial_k \\right)^\\dagger \\partial_k \n",
+ "+ \\partial_{k+1} \\left( \\partial_{k+1} \\right)^\\dagger\n",
+ "$$\n",
+ "\n",
+ "The QPE algorithm will be used to estimate the number of zero eigenvalues of the Laplacian matrix. Since its exponential matrix serves as the unitary matrix in this algorithm, it must have dimensions $2^q \\times 2^q$, where $q$ represents the number of target qubits. It is recommanded to pad the combinatorial laplacian $\\Delta_k$ with an identity matrix with $\\tilde{\\lambda}_{max}/2$ in place of ones, where $\\tilde{\\lambda}_{max}$ is the estimate of the maximum eigenvalue of $\\Delta_k$ using the Gershgorin circle theorem ([details](https://mathworld.wolfram.com/GershgorinCircleTheorem.html)), such that:\n",
+ "\n",
+ "$$\n",
+ "\\tilde{\\Delta}_k =\n",
+ "\\begin{bmatrix}\n",
+ "\\Delta_k & 0 \\\\\n",
+ "0 & \\frac{\\widetilde{\\lambda}_{\\text{max}}}{2} \\cdot I_{2q - |S_k|}\n",
+ "\\end{bmatrix}_{2q \\times 2q}\n",
+ "$$\n",
+ "\n",
+ "where $\\Delta_k$ is the padded combinatorial laplacian and $q = \\lceil \\log_2 |S_k| \\rceil$ is the number of qubits this operator will act on. In QPE, as $2\\pi\\theta$ increases beyond $2\\pi$, the eigenvalues will start repeating due to their periodic form. Thus, $\\theta$ is restricted to $[0, 1)$. As $\\lambda \\to 2\\pi\\theta$ this means that $λ \\in [0, 2\\pi)$. Thus, we need to restrict the eigenvalues of the combinatorial laplacian to this range. This can be achieved by rescaling $\\tilde{\\Delta}_k$ by $\\delta/\\tilde{\\lambda}_{max}$ where $\\delta$ is slightly less than $2\\pi$. Thus, the rescaled matrix $H$ and the unitary marix $U$ for QPE are:\n",
+ "\n",
+ "$$\n",
+ "\\begin{align*}\n",
+ "H &= \\frac{\\delta}{\\tilde{\\lambda}_k} \\tilde{\\Delta}_k\\\\\n",
+ "U &= e^{iH}\n",
+ "\\end{align*}\n",
+ "$$\n",
+ "\n",
+ "**Example:**\n",
+ "\n",
+ "In our example, the combinational laplacian is in the form of a $6\\times6$ matrix:\n",
+ "\n",
+ "$$\n",
+ "\\begin{align*}\n",
+ "\\Delta_1 &= (\\partial_1)^\\dagger \\partial_1 + \\partial_2 (\\partial_2)^\\dagger\\\\\n",
+ "&=\n",
+ "\\begin{bmatrix}\n",
+ "3 & 0 & 0 & 0 & 0 & 0 \\\\\n",
+ "0 & 3 & 0 & -1 & -1 & 0 \\\\\n",
+ "0 & 0 & 3 & -1 & -1 & 0 \\\\\n",
+ "0 & -1 & -1 & 2 & 1 & -1 \\\\\n",
+ "0 & -1 & -1 & 1 & 2 & 1 \\\\\n",
+ "0 & 0 & 0 & -1 & 1 & 2\n",
+ "\\end{bmatrix}\n",
+ "\\end{align*}\n",
+ "$$\n",
+ "\n",
+ "It is padded with $\\tilde{\\lambda}_{max}=6$ and $\\delta=6$ to its nearest power of 2, which is 8 ($q=3$):\n",
+ "\n",
+ "$$\n",
+ "\\begin{align}\n",
+ "H_1 = \\begin{bmatrix}\n",
+ "3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
+ "0 & 3 & 0 & -1 & -1 & 0 & 0 & 0 \\\\\n",
+ "0 & 0 & 3 & -1 & -1 & 0 & 0 & 0 \\\\\n",
+ "0 & -1 & -1 & 2 & 1 & -1 & 0 & 0 \\\\\n",
+ "0 & -1 & -1 & 1 & 2 & 1 & 0 & 0 \\\\\n",
+ "0 & 0 & 0 & -1 & 1 & 2 & 0 & 0 \\\\\n",
+ "0 & 0 & 0 & 0 & 0 & 0 & 3 & 0 \\\\\n",
+ "0 & 0 & 0 & 0 & 0 & 0 & 0 & 3\n",
+ "\\end{bmatrix}\n",
+ "\\end{align}\n",
+ "$$\n",
+ "\n",
+ "**Action:**\n",
+ "\n",
+ "Define a function to build the Laplacian, where Define a function that automatically determines whether the input Laplacian matrix requires padding. If padding is needed, the function will pad and rescale the matrix accordingly. Then, build the unitary based on the padded matrix, in the form of a circuit.\n",
+ "\n",
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "id": "db7a8edd-6f84-4e30-93cf-28217160c951",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_diagonal_value(laplacian):\n",
+ " diagonal = np.diag(laplacian)\n",
+ " row_sums = np.sum(np.abs(laplacian), axis=1) - np.abs(diagonal)\n",
+ " \n",
+ " # Find the maximum upper bound\n",
+ " max_eigenvalue = np.max(diagonal + row_sums)\n",
+ " \n",
+ " return max_eigenvalue"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 35,
+ "id": "8f818f78",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# write your code here\n",
+ "def pad_laplacian(laplacian):\n",
+ " if np.log2(laplacian.shape[0]).is_integer():\n",
+ " return\n",
+ " else:\n",
+ " diagonal_value = get_diagonal_value(laplacian)\n",
+ "\n",
+ " # Calculate new shape\n",
+ " new_shape = int(2 ** np.ceil(np.log2(laplacian.shape[0])))\n",
+ " extra_pad = new_shape - laplacian.shape[0]\n",
+ " new_laplacian = np.pad(laplacian, ((0, extra_pad), (0, extra_pad)), mode='constant', constant_values=0)\n",
+ "\n",
+ " # Padding\n",
+ " rows, cols = np.diag_indices_from(new_laplacian)\n",
+ " new_laplacian[rows[new_shape - extra_pad:], cols[new_shape - extra_pad:]] = diagonal_value / 2\n",
+ "\n",
+ " # Rescaling matrix at the end after padding\n",
+ " new_laplacian = new_laplacian * (6 / diagonal_value)\n",
+ "\n",
+ " return new_laplacian"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "id": "57bd1ba5-78e2-484c-892c-5fa8ceb272f4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[6., 0., 0., 0.],\n",
+ " [0., 6., 0., 0.],\n",
+ " [0., 0., 6., 0.],\n",
+ " [0., 0., 0., 3.]])"
+ ]
+ },
+ "execution_count": 36,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "H_empty = pad_laplacian(laplacian_empty)\n",
+ "H_empty"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "id": "7d39eeca-3d0a-4ba5-bcaa-b65eee05d645",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 3. , -0.75, -0.75, -0.75, -0.75, 0. , 0. , 0. ],\n",
+ " [-0.75, 3. , -0.75, -0.75, -0.75, 0. , 0. , 0. ],\n",
+ " [-0.75, -0.75, 3. , -0.75, -0.75, 0. , 0. , 0. ],\n",
+ " [-0.75, -0.75, -0.75, 3. , -0.75, 0. , 0. , 0. ],\n",
+ " [-0.75, -0.75, -0.75, -0.75, 3. , 0. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. , 3. , 0. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. , 0. , 3. , 0. ],\n",
+ " [ 0. , 0. , 0. , 0. , 0. , 0. , 0. , 3. ]])"
+ ]
+ },
+ "execution_count": 37,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "H = pad_laplacian(laplacian)\n",
+ "H"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "id": "1c256dec-1568-4c39-954e-75632f657577",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[-0.45644749-0.45724905j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0.36411187+0.11431226j, -0.45644749-0.45724905j,\n",
+ " 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " -0.45644749-0.45724905j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, -0.45644749-0.45724905j,\n",
+ " 0.36411187+0.11431226j, 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " 0.36411187+0.11431226j, 0.36411187+0.11431226j,\n",
+ " -0.45644749-0.45724905j, 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , -0.9899925 +0.14112001j,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ,\n",
+ " -0.9899925 +0.14112001j, 0. +0.j ],\n",
+ " [ 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , -0.9899925 +0.14112001j]])"
+ ]
+ },
+ "execution_count": 38,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "from scipy.linalg import expm\n",
+ "\n",
+ "# Calculate e^(i*matrix)\n",
+ "U = expm(1j * H)\n",
+ "\n",
+ "U"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "id": "78b6cb89-42c2-4bec-9f2f-cb38fd09c25c",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 0.96017029-0.2794155j , 0. +0.j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0. +0.j , 0.96017029-0.2794155j ,\n",
+ " 0. +0.j , 0. +0.j ],\n",
+ " [ 0. +0.j , 0. +0.j ,\n",
+ " 0.96017029-0.2794155j , 0. +0.j ],\n",
+ " [ 0. +0.j , 0. +0.j ,\n",
+ " 0. +0.j , -0.9899925 +0.14112001j]])"
+ ]
+ },
+ "execution_count": 39,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "U_empty = expm(1j * H_empty)\n",
+ "U_empty"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 40,
+ "id": "76b24487-b9fa-4483-be70-f2db68aabaed",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def getU(laplacian):\n",
+ " padded_laplacian = pad_laplacian(laplacian)\n",
+ " U = expm(1j * padded_laplacian)\n",
+ " return U"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0acfea73",
+ "metadata": {},
+ "source": [
+ "### Step 3: Applying QPE\n",
+ "\n",
+ "**Instruction:**\n",
+ "\n",
+ "The Betti number is the number of zero eigenvalues in the Laplacian ([ref](https://link.springer.com/article/10.1007/PL00009218)). \n",
+ "\n",
+ "$$\n",
+ "\\beta_k = \\dim (\\ker(\\Delta_k))\n",
+ "$$\n",
+ "\n",
+ "The betti curve is then a series of Betti numbers on different resolution threshold `epsilon`.\n",
+ "\n",
+ "To estimate the number of zero eigenvalues (nullity) in the padded Laplacian matrix (padding didn't add more zero eigenvalues), QPE algorithm is employed. The fundamental concept is that, if the target qubits start out in the maximally mixed state (shown below), which can be thought of as a random choice of an eigenstate, the proportion of all-zero states among all measured states is equal to the proportion of zero eigenvalues among all eigenvalues. Assume the all-zero state show up for $\\{i|\\tilde{\\theta}_i=0\\}$ times in $\\alpha$ shots, the probability of getting all-zero state $p(0)$ is given by:\n",
+ "\n",
+ "$$\n",
+ "\\begin{align*}\n",
+ "p(0) &= \\frac{\\left| \\{i \\mid \\tilde{\\theta}_i = 0\\} \\right|}{\\alpha} = \\frac{\\tilde{\\beta}_k}{2^q} \\\\\n",
+ "\\implies \\tilde{\\beta}_k &= 2^q \\cdot p(0)\n",
+ "\\end{align*}\n",
+ "$$\n",
+ "\n",
+ "Where $\\tilde{\\beta}_k$ is the estimation of $k$-th Betti number. This estimation is then rounded to the nearest integer to obtain the final result."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "72511671",
+ "metadata": {},
+ "source": [
+ "For your reference, the tutorial of QPE in several major quantum computing libraries are listed below:\n",
+ "\n",
+ "* [Qiskit](https://github.com/qiskit-community/qiskit-textbook/blob/main/content/ch-algorithms/quantum-phase-estimation.ipynb)\n",
+ "* [Pennylane](https://pennylane.ai/qml/demos/tutorial_qpe)\n",
+ "* [CUDA-Q](https://nvidia.github.io/cuda-quantum/latest/specification/cudaq/examples.html#quantum-phase-estimation:~:text=Quantum%20Phase%20Estimation-,%C2%B6,-C%2B%2B)\n",
+ "* [Cirq](https://quantumai.google/cirq/experiments/textbook_algorithms#phase_estimation)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "481e000c",
+ "metadata": {},
+ "source": [
+ "**Example:**\n",
+ "\n",
+ "In our example, the probability of measuring all-zero states is approximately $p(0)=0.137 = \\tilde{\\beta}_k / 2^3 \\implies \\tilde{\\beta}_k = 1.096$, which is then rounded to $1$.\n",
+ "\n",
+ "**Action:**\n",
+ "\n",
+ "Utilize your preferred quantum computing library to apply QPE for estimating the number of zero eigenvalues in the Laplacian matrix. Note that the exponential of the Laplacian matrix is used as the unitary operator in QPE.\n",
+ "\n",
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 41,
+ "id": "848ef40c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
+ "from qiskit.quantum_info import Statevector\n",
+ "from qiskit.visualization import plot_histogram\n",
+ "from qiskit.circuit.library import QFT\n",
+ "from qiskit.quantum_info import Operator\n",
+ "from qiskit.primitives import Sampler"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 42,
+ "id": "138e2f9f-2080-427d-b8c6-e12b252152bc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "precision_qubits = 3\n",
+ "target_qubits = int(np.log2(U_empty.shape[0]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 43,
+ "id": "ea94b17f-46bd-4523-a808-3120deb2a00d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# def laplacian_qpe(laplacian, num_qubits):\n",
+ "# # Exponentiate the Laplacian matrix\n",
+ "# exp_matrix = Operator(laplacian)\n",
+ " \n",
+ "# # Create quantum registers\n",
+ "# q_eval = QuantumRegister(num_qubits, 'eval')\n",
+ "# q_vec = QuantumRegister(np.log2(laplacian.shape[0]), 'vec')\n",
+ "\n",
+ "# print(list(q_eval))\n",
+ " \n",
+ "# # Create quantum circuit\n",
+ "# circuit = QuantumCircuit(q_eval, q_vec)\n",
+ " \n",
+ "# # Initialize vector state (assuming uniform superposition)\n",
+ "# circuit.h(q_vec)\n",
+ " \n",
+ "# # Apply QPE\n",
+ "# circuit.h(q_eval)\n",
+ "# for i in range(num_qubits):\n",
+ "# circuit.unitary(exp_matrix, [q_vec[j] for j in range(int(np.log2(laplacian.shape[0])))])\n",
+ "# circuit.append(QFT(num_qubits, do_swaps=False).inverse(), q_eval)\n",
+ " \n",
+ "# # Simulate the circuit\n",
+ "# statevector = Statevector.from_instruction(circuit)\n",
+ " \n",
+ "# # Extract probabilities\n",
+ "# probabilities = statevector.probabilities()\n",
+ " \n",
+ "# # Count non-zero eigenvalues (excluding the zero state)\n",
+ "# non_zero_count = sum(prob for prob in probabilities[1:] if prob > 1e-6)\n",
+ " \n",
+ "# return non_zero_count"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 44,
+ "id": "a34faa96-b901-4999-a6d0-761fd8dcb650",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# def laplacian_qpe(U, num_qubits):\n",
+ "# # Exponentiate the Laplacian matrix\n",
+ "# exp_matrix = Operator(U)\n",
+ " \n",
+ "# # Create quantum registers\n",
+ "\n",
+ "# print(U.shape)\n",
+ " \n",
+ "# q_eval = QuantumRegister(num_qubits, 'eval')\n",
+ "# q_vec = QuantumRegister(int(np.log2(U.shape[0])), 'vec')\n",
+ "# print(\"vec_size\", int(np.log2(U.shape[0])))\n",
+ "# q_aux = QuantumRegister(3, 'aux')\n",
+ " \n",
+ "# # Create quantum circuit\n",
+ "# circuit = QuantumCircuit(q_eval, q_vec, q_aux)\n",
+ " \n",
+ "# # Initialize vector state (assuming uniform superposition)\n",
+ "# circuit.h(q_vec)\n",
+ " \n",
+ "# # Apply QPE\n",
+ "# circuit.h(q_eval)\n",
+ "# for i in range(num_qubits):\n",
+ "# circuit.unitary(exp_matrix, [q_vec[j] for j in reversed(range(int(np.log2(U.shape[0]))))])\n",
+ "# circuit.append(QFT(num_qubits, do_swaps=False).inverse(), q_eval)\n",
+ " \n",
+ "# # Simulate the circuit\n",
+ "# statevector = Statevector.from_instruction(circuit)\n",
+ " \n",
+ "# # Extract probabilities\n",
+ "# probabilities = statevector.probabilities()\n",
+ " \n",
+ "# # Count non-zero eigenvalues (excluding the zero state)\n",
+ "# non_zero_count = sum(prob for prob in probabilities[1:] if prob > 1e-6)\n",
+ " \n",
+ "# return non_zero_count"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 45,
+ "id": "3abfbc52-338e-4542-b2d4-236fe9ffe6ca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Nikhil Magic\n",
+ "\n",
+ "import numpy as np\n",
+ "from math import log2\n",
+ "from scipy.linalg import expm\n",
+ "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
+ "from qiskit.circuit.library import QFT, UnitaryGate\n",
+ "from qiskit.quantum_info import Statevector\n",
+ "\n",
+ "\n",
+ "def project_to_unitary(mat):\n",
+ " \"\"\"\n",
+ " (Optional) SVD-based projection to force exact unitarity.\n",
+ " If 'mat' is close to unitary, it will become exactly unitary \n",
+ " by ignoring singular values.\n",
+ " \"\"\"\n",
+ " U, s, Vdag = np.linalg.svd(mat, full_matrices=False)\n",
+ " return U @ Vdag\n",
+ "\n",
+ "\n",
+ "def laplacian_nullity_qpe(laplacian, num_eval_qubits, shots=100, scale=1.0):\n",
+ " \"\"\"\n",
+ " Estimate the fraction of amplitude in the null space of 'laplacian' by running\n",
+ " QPE on U = exp(i * scale * laplacian). \n",
+ " We return (zero_fraction, circuit) where:\n",
+ " - zero_fraction is the probability that the eval qubits measure |0...0>\n",
+ " - circuit is the full QPE circuit (with measure instructions)\n",
+ "\n",
+ " Args:\n",
+ " laplacian (np.array): Hermitian NxN matrix with N = 2^n\n",
+ " num_eval_qubits (int): Number of qubits in the phase (evaluation) register\n",
+ " shots (int): Unused for local statevector, but kept for compatibility\n",
+ " scale (float): A scaling factor for the exponent in exp(i * scale * L)\n",
+ "\n",
+ " Returns:\n",
+ " zero_fraction (float): Probability that the eval qubits measure |0...0>\n",
+ " qc (QuantumCircuit): The full QPE circuit (with measure)\n",
+ " \"\"\"\n",
+ " dim = laplacian.shape[0]\n",
+ " num_sys_qubits = int(log2(dim))\n",
+ " if 2**num_sys_qubits != dim:\n",
+ " raise ValueError(\"Matrix dimension must be 2^n for standard QPE.\")\n",
+ " \n",
+ " # 1) Build a 'unitary' from the Laplacian via exponentiation\n",
+ " U_matrix = expm(1j * scale * laplacian)\n",
+ " U_matrix = project_to_unitary(U_matrix) # ensure exactly unitary\n",
+ "\n",
+ " # 2) Set up QPE circuit\n",
+ " q_aux = QuantumRegister(num_eval_qubits, 'aux')\n",
+ " q_eval = QuantumRegister(num_eval_qubits, 'eval')\n",
+ " q_sys = QuantumRegister(num_sys_qubits, 'sys')\n",
+ " c_eval = ClassicalRegister(num_eval_qubits, 'meas')\n",
+ " qc = QuantumCircuit(q_aux,q_eval, q_sys, c_eval, name=\"QPE_laplacian\")\n",
+ " \n",
+ " # Put system qubits in a uniform superposition\n",
+ " qc.h(q_sys)\n",
+ " qc.cx(q_sys, q_aux)\n",
+ " # Hadamard on the evaluation qubits\n",
+ " qc.h(q_eval)\n",
+ " \n",
+ " # 3) Controlled powers of U (U^(2^i)), building them as UnitaryGates\n",
+ " # bypassing Qiskit's unitarity check if needed\n",
+ " for i in range(num_eval_qubits):\n",
+ " power = 2**(num_eval_qubits -1 - i)\n",
+ " powered_matrix = np.linalg.matrix_power(U_matrix, power)\n",
+ " powered_matrix = project_to_unitary(powered_matrix) # ensure unitarity\n",
+ " exponent_gate = UnitaryGate(powered_matrix, label=f\"U^{power}\", check_input=False)\n",
+ " ctrl_exponent = exponent_gate.control(1)\n",
+ " qc.append(ctrl_exponent, [q_eval[i]] + list(q_sys))\n",
+ " \n",
+ " # 4) Inverse QFT on the eval qubits\n",
+ " iqft = QFT(num_eval_qubits, inverse=True, do_swaps=True)\n",
+ " qc.append(iqft, q_eval)\n",
+ " \n",
+ " # 5) Measure the eval qubits\n",
+ " qc.measure(q_eval, c_eval)\n",
+ " \n",
+ " qc_no_measure = qc.remove_final_measurements(inplace=False)\n",
+ " # Now we can safely do from_instruction(qc_no_measure)\n",
+ " final_state = Statevector.from_instruction(qc_no_measure)\n",
+ " \n",
+ " probs = final_state.probabilities()\n",
+ " \n",
+ " # 7) Sum probability for all basis states with eval qubits = |0...0>\n",
+ " zero_fraction = 0.0\n",
+ " total_qubits = num_eval_qubits + num_sys_qubits\n",
+ " for basis_index in range(2**total_qubits):\n",
+ " # The lower 'num_eval_qubits' bits correspond to the eval qubits\n",
+ " eval_bits = basis_index % (2**num_eval_qubits)\n",
+ " if eval_bits == 0:\n",
+ " zero_fraction += probs[basis_index]\n",
+ "\n",
+ " return 1 - zero_fraction"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 46,
+ "id": "08697c8c-3700-470b-ac0f-8344ae11bbb6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Terry's version of QPE\n",
+ "\n",
+ "from qiskit.circuit.library import QFT, UnitaryGate\n",
+ "\n",
+ "def laplacian_qpe(U, num_qubits):\n",
+ " # Exponentiate the Laplacian matrix\n",
+ " exp_matrix = Operator(U)\n",
+ " \n",
+ " # Create quantum registers\n",
+ " \n",
+ " q_eval = QuantumRegister(num_qubits, 'eval')\n",
+ " q_vec = QuantumRegister(int(np.log2(U.shape[0])), 'vec')\n",
+ " q_aux = QuantumRegister(3, 'aux')\n",
+ " c_eval = ClassicalRegister(num_qubits, 'meas')\n",
+ " \n",
+ " # Create quantum circuit\n",
+ " circuit = QuantumCircuit(q_eval, q_vec, q_aux, c_eval)\n",
+ " \n",
+ " # Initialize vector state (assuming uniform superposition)\n",
+ " circuit.h(q_vec)\n",
+ " \n",
+ " # Apply QPE\n",
+ " circuit.h(q_eval)\n",
+ "\n",
+ " for i in range(int(np.log2(U.shape[0]))):\n",
+ " circuit.cx(q_vec[i], q_aux[i])\n",
+ " \n",
+ " # for i in range(num_qubits):\n",
+ " # # circuit.unitary(exp_matrix, [q_vec[j] for j in reversed(range(int(np.log2(U.shape[0]))))])\n",
+ " # circuit.unitary(exp_matrix, [q_vec[j] for j in range(int(np.log2(U.shape[0])))])\n",
+ "\n",
+ " for i in range(num_qubits):\n",
+ " power = 2**(num_qubits -1 - i)\n",
+ " powered_matrix = np.linalg.matrix_power(U, power)\n",
+ " # powered_matrix = project_to_unitary(powered_matrix) # ensure unitarity\n",
+ " exponent_gate = UnitaryGate(powered_matrix, label=f\"U^{power}\", check_input=False)\n",
+ " ctrl_exponent = exponent_gate.control(1)\n",
+ " circuit.append(ctrl_exponent, [q_eval[i]] + list(q_vec))\n",
+ " \n",
+ " circuit.append(QFT(num_qubits, do_swaps=False).inverse(), q_eval)\n",
+ "\n",
+ " circuit.measure(q_eval, c_eval)\n",
+ " \n",
+ " # Simulate the circuit\n",
+ " # statevector = Statevector.from_instruction(circuit)\n",
+ " \n",
+ " # Extract probabilities\n",
+ " # probabilities = statevector.probabilities()\n",
+ " \n",
+ " # Count non-zero eigenvalues (excluding the zero state)\n",
+ " # non_zero_count = sum(prob for prob in probabilities[1:] if prob > 1e-6)\n",
+ "\n",
+ " qc_no_measure = circuit.remove_final_measurements(inplace=False)\n",
+ " # Now we can safely do from_instruction(qc_no_measure)\n",
+ " final_state = Statevector.from_instruction(qc_no_measure)\n",
+ " \n",
+ " probs = final_state.probabilities()\n",
+ " # 7) Sum probability for all basis states with eval qubits = |0...0>\n",
+ " zero_fraction = 0.0\n",
+ " total_qubits = num_qubits + int(np.log2(U.shape[0]))\n",
+ " for basis_index in range(2**total_qubits):\n",
+ " # The lower 'num_eval_qubits' bits correspond to the eval qubits\n",
+ " eval_bits = basis_index % (2**num_qubits)\n",
+ " if eval_bits == 0:\n",
+ " zero_fraction += probs[basis_index]\n",
+ " non_zero_count = (1 - zero_fraction)\n",
+ " \n",
+ " # return non_zero_count, circuit\n",
+ " return non_zero_count"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f8724ca3",
+ "metadata": {},
+ "source": [
+ " NOTE:\n",
+ "\n",
+ "The unitary operator can be constructed by converting the exponential matrix of the Laplacian into a quantum circuit. For instance, in Qiskit, this can be implemented using `circuit.unitary(exp_matrix)`. Alternative methods for constructing the unitary operator will be optionally explored in the **BONUS** section at the end of this notebook.\n",
+ "
\n",
+ "\n",
+ " \n",
+ "BONUS EXERCISE: \n",
+ "\n",
+ "1. Why we should measure all-zero states?\n",
+ "2. What is the difference between a maximally mixed state and the $(H\\ket{0})^{\\otimes n}$ state? Two possible aspects are:\n",
+ "\n",
+ "* Plotting of their density matrix.\n",
+ "* Results from QPE.\n",
+ "\n",
+ "3. What parameters affect the accuracy of the estimation before rounding? For Laplacian matrices of varying sizes, how does the accuracy depend on these parameters? Within what range of values do these parameters guarantee (or are highly likely to produce) a correct final result?\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1cbe5cf2",
+ "metadata": {},
+ "source": [
+ "### Step 4: Detecting market crashes"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7f27e81b",
+ "metadata": {},
+ "source": [
+ "**Instruction:**\n",
+ "\n",
+ "At this point, we have betti curves for each window across our dataset, and we wish to use this to detect market crashes. One such way is to take the difference between these curves—or the pairwise distance—for successive windows and look for spikes. This can be done with the $L^p$ norm of the betti curve for each window, defined as follows:\n",
+ "\n",
+ "$$||x||_p = (\\sum_{n}^{i=1} |x_i|^p)^{1/p}$$\n",
+ "\n",
+ "Combining these pairwise distances into a vector, we get a single output curve we can analyze. Experiment with different values of $p$, but a good starting point is the $L^2$ Norm. Using this, it is possible to detect regions where a market crash is occuring. Comparing detected crashes with the price data indicates how accurate the crash detection methodology is ([ref](https://github.com/giotto-ai/stock-market-crashes/blob/master/Stock%20Market%20Crash%20Detection.ipynb)). \n",
+ "\n",
+ "**Action:**\n",
+ "\n",
+ "Use the $L^p$ norm to create pairwise distance curves for successive windows, and then use the results to define when a crash is occuring. Compare this with your data to see how well it performs. You may find the following classical solver is useful.\n",
+ "\n",
+ "**Answer:**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 47,
+ "id": "5ef68f98",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0 : 1\n",
+ "1 : 1\n",
+ "2 : 1\n",
+ "3 : 1\n",
+ "4 : 1\n",
+ "5 : 1\n",
+ "6 : 1\n",
+ "7 : 1\n",
+ "8 : 1\n",
+ "9 : 1\n",
+ "10 : 1\n",
+ "11 : 1\n",
+ "12 : 1\n",
+ "13 : 1\n",
+ "14 : 1\n",
+ "15 : 1\n",
+ "16 : 1\n",
+ "17 : 1\n",
+ "18 : 1\n",
+ "19 : 1\n",
+ "20 : 1\n",
+ "21 : 1\n",
+ "22 : 1\n",
+ "23 : 1\n",
+ "24 : 1\n",
+ "25 : 1\n",
+ "26 : 1\n",
+ "27 : 1\n",
+ "28 : 1\n",
+ "29 : 1\n",
+ "30 : 1\n",
+ "31 : 1\n",
+ "32 : 1\n",
+ "33 : 2\n",
+ "34 : 1\n",
+ "35 : 1\n",
+ "36 : 1\n",
+ "37 : 1\n",
+ "38 : 2\n",
+ "39 : 1\n",
+ "40 : 2\n",
+ "41 : 1\n",
+ "42 : 2\n",
+ "43 : 3\n",
+ "44 : 2\n",
+ "45 : 2\n",
+ "46 : 1\n",
+ "47 : 1\n",
+ "48 : 2\n",
+ "49 : 2\n",
+ "50 : 2\n",
+ "51 : 1\n",
+ "52 : 1\n",
+ "53 : 1\n",
+ "54 : 1\n",
+ "55 : 1\n",
+ "56 : 1\n",
+ "57 : 1\n",
+ "58 : 1\n",
+ "59 : 1\n",
+ "60 : 1\n",
+ "61 : 1\n",
+ "62 : 1\n",
+ "63 : 1\n"
+ ]
+ }
+ ],
+ "source": [
+ "from ripser import ripser\n",
+ "\n",
+ "def classical_betti_solver(point_cloud, epsilon, dim):\n",
+ " '''Return the Betti number on a given point cloud.\n",
+ " Args:\n",
+ " point_cloud: the point cloud after applying the sliding window.\n",
+ " epsilon: resolution threshold.\n",
+ " dim: the dimension on which the Betti number is calculated\n",
+ " '''\n",
+ " result = ripser(point_cloud, maxdim=dim)\n",
+ " diagrams = result[\"dgms\"]\n",
+ " return len(\n",
+ " [interval for interval in diagrams[dim] if interval[0] < epsilon < interval[1]]\n",
+ " )\n",
+ "\n",
+ "# write your code here\n",
+ "count = 0\n",
+ "for pc in point_clouds:\n",
+ " print(count, ':', classical_betti_solver(pc, 0.1, 0))\n",
+ " count = count + 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 48,
+ "id": "264ca9ed-0da9-4616-86bb-5a73ea985673",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "time_series = np.log(pd.read_csv(\"sp500.csv\", header=None).to_numpy().squeeze())\n",
+ "point_clouds_final = getPointClouds(time_series, 5, 4, 5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 49,
+ "id": "5a6bd57c-b653-41e6-82dc-4c10f6c7bd9a",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "KeyboardInterrupt",
+ "evalue": "",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[49], line 11\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[38;5;66;03m# print(laplacian_final)\u001b[39;00m\n\u001b[1;32m 10\u001b[0m U_final \u001b[38;5;241m=\u001b[39m getU(laplacian_final)\n\u001b[0;32m---> 11\u001b[0m non_zero_count_final \u001b[38;5;241m=\u001b[39m \u001b[43mlaplacian_qpe\u001b[49m\u001b[43m(\u001b[49m\u001b[43mU_final\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m3\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 12\u001b[0m betti_values\u001b[38;5;241m.\u001b[39mappend(\u001b[38;5;241m4\u001b[39m\u001b[38;5;241m*\u001b[39m(\u001b[38;5;241m8\u001b[39m \u001b[38;5;241m-\u001b[39m non_zero_count_final \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m8\u001b[39m))\n\u001b[1;32m 15\u001b[0m classical\u001b[38;5;241m.\u001b[39mappend(np\u001b[38;5;241m.\u001b[39msum(eigenvalues \u001b[38;5;241m<\u001b[39m \u001b[38;5;241m1e-6\u001b[39m))\n",
+ "Cell \u001b[0;32mIn[46], line 37\u001b[0m, in \u001b[0;36mlaplacian_qpe\u001b[0;34m(U, num_qubits)\u001b[0m\n\u001b[1;32m 35\u001b[0m \u001b[38;5;66;03m# powered_matrix = project_to_unitary(powered_matrix) # ensure unitarity\u001b[39;00m\n\u001b[1;32m 36\u001b[0m exponent_gate \u001b[38;5;241m=\u001b[39m UnitaryGate(powered_matrix, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mU^\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mpower\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m, check_input\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[0;32m---> 37\u001b[0m ctrl_exponent \u001b[38;5;241m=\u001b[39m \u001b[43mexponent_gate\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcontrol\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 38\u001b[0m circuit\u001b[38;5;241m.\u001b[39mappend(ctrl_exponent, [q_eval[i]] \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mlist\u001b[39m(q_vec))\n\u001b[1;32m 40\u001b[0m circuit\u001b[38;5;241m.\u001b[39mappend(QFT(num_qubits, do_swaps\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\u001b[38;5;241m.\u001b[39minverse(), q_eval)\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/unitary.py:193\u001b[0m, in \u001b[0;36mUnitaryGate.control\u001b[0;34m(self, num_ctrl_qubits, label, ctrl_state, annotated)\u001b[0m\n\u001b[1;32m 185\u001b[0m cmat \u001b[38;5;241m=\u001b[39m _compute_control_matrix(mat, num_ctrl_qubits, ctrl_state\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[1;32m 186\u001b[0m iso \u001b[38;5;241m=\u001b[39m Isometry(cmat, \u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m)\n\u001b[1;32m 187\u001b[0m gate \u001b[38;5;241m=\u001b[39m ControlledGate(\n\u001b[1;32m 188\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mc-unitary\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 189\u001b[0m num_qubits\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnum_qubits \u001b[38;5;241m+\u001b[39m num_ctrl_qubits,\n\u001b[1;32m 190\u001b[0m params\u001b[38;5;241m=\u001b[39m[mat],\n\u001b[1;32m 191\u001b[0m label\u001b[38;5;241m=\u001b[39mlabel,\n\u001b[1;32m 192\u001b[0m num_ctrl_qubits\u001b[38;5;241m=\u001b[39mnum_ctrl_qubits,\n\u001b[0;32m--> 193\u001b[0m definition\u001b[38;5;241m=\u001b[39m\u001b[43miso\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdefinition\u001b[49m,\n\u001b[1;32m 194\u001b[0m ctrl_state\u001b[38;5;241m=\u001b[39mctrl_state,\n\u001b[1;32m 195\u001b[0m base_gate\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcopy(),\n\u001b[1;32m 196\u001b[0m )\n\u001b[1;32m 197\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 198\u001b[0m gate \u001b[38;5;241m=\u001b[39m AnnotatedOperation(\n\u001b[1;32m 199\u001b[0m \u001b[38;5;28mself\u001b[39m, ControlModifier(num_ctrl_qubits\u001b[38;5;241m=\u001b[39mnum_ctrl_qubits, ctrl_state\u001b[38;5;241m=\u001b[39mctrl_state)\n\u001b[1;32m 200\u001b[0m )\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/instruction.py:329\u001b[0m, in \u001b[0;36mInstruction.definition\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 327\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Return definition in terms of other basic gates.\"\"\"\u001b[39;00m\n\u001b[1;32m 328\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_definition \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 329\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_define\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 330\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_definition\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:124\u001b[0m, in \u001b[0;36mIsometry._define\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21m_define\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 121\u001b[0m \u001b[38;5;66;03m# TODO The inverse().inverse() is because there is code to uncompute (_gates_to_uncompute)\u001b[39;00m\n\u001b[1;32m 122\u001b[0m \u001b[38;5;66;03m# an isometry, but not for generating its decomposition. It would be cheaper to do the\u001b[39;00m\n\u001b[1;32m 123\u001b[0m \u001b[38;5;66;03m# later here instead.\u001b[39;00m\n\u001b[0;32m--> 124\u001b[0m gate \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minv_gate\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 125\u001b[0m gate \u001b[38;5;241m=\u001b[39m gate\u001b[38;5;241m.\u001b[39minverse()\n\u001b[1;32m 126\u001b[0m q \u001b[38;5;241m=\u001b[39m QuantumRegister(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnum_qubits, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mq\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:340\u001b[0m, in \u001b[0;36mIsometry.inv_gate\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 336\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Return the adjoint of the unitary.\"\"\"\u001b[39;00m\n\u001b[1;32m 337\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_inverse \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 338\u001b[0m \u001b[38;5;66;03m# call to generate the circuit that takes the isometry to the first 2^m columns\u001b[39;00m\n\u001b[1;32m 339\u001b[0m \u001b[38;5;66;03m# of the 2^n identity matrix\u001b[39;00m\n\u001b[0;32m--> 340\u001b[0m iso_circuit \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_gates_to_uncompute\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 341\u001b[0m \u001b[38;5;66;03m# invert the circuit to create the circuit implementing the isometry\u001b[39;00m\n\u001b[1;32m 342\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_inverse \u001b[38;5;241m=\u001b[39m iso_circuit\u001b[38;5;241m.\u001b[39mto_instruction()\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:160\u001b[0m, in \u001b[0;36mIsometry._gates_to_uncompute\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 155\u001b[0m \u001b[38;5;66;03m# Decompose the column with index column_index and attache the gate to the circuit object.\u001b[39;00m\n\u001b[1;32m 156\u001b[0m \u001b[38;5;66;03m# Return the isometry that is left to decompose, where the columns up to index column_index\u001b[39;00m\n\u001b[1;32m 157\u001b[0m \u001b[38;5;66;03m# correspond to the firstfew columns of the identity matrix up to diag, and hence we only\u001b[39;00m\n\u001b[1;32m 158\u001b[0m \u001b[38;5;66;03m# have to save a list containing them.\u001b[39;00m\n\u001b[1;32m 159\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m column_index \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mm):\n\u001b[0;32m--> 160\u001b[0m remaining_isometry, diag \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_decompose_column\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 161\u001b[0m \u001b[43m \u001b[49m\u001b[43mcircuit\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mq\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdiag\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mremaining_isometry\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcolumn_index\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 163\u001b[0m \u001b[38;5;66;03m# extract phase of the state that was sent to the basis state ket(column_index)\u001b[39;00m\n\u001b[1;32m 164\u001b[0m diag\u001b[38;5;241m.\u001b[39mappend(remaining_isometry[column_index, \u001b[38;5;241m0\u001b[39m])\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:180\u001b[0m, in \u001b[0;36mIsometry._decompose_column\u001b[0;34m(self, circuit, q, diag, remaining_isometry, column_index)\u001b[0m\n\u001b[1;32m 178\u001b[0m n \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mint\u001b[39m(math\u001b[38;5;241m.\u001b[39mlog2(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39miso_data\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m]))\n\u001b[1;32m 179\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m s \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n):\n\u001b[0;32m--> 180\u001b[0m remaining_isometry, diag \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_disentangle\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 181\u001b[0m \u001b[43m \u001b[49m\u001b[43mcircuit\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mq\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdiag\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mremaining_isometry\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcolumn_index\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43ms\u001b[49m\n\u001b[1;32m 182\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 183\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m remaining_isometry, diag\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:237\u001b[0m, in \u001b[0;36mIsometry._disentangle\u001b[0;34m(self, circuit, q, diag, remaining_isometry, column_index, s)\u001b[0m\n\u001b[1;32m 235\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m isometry_rs\u001b[38;5;241m.\u001b[39mucg_is_identity_up_to_global_phase(single_qubit_gates, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_epsilon):\n\u001b[1;32m 236\u001b[0m control_labels \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(\u001b[38;5;28mrange\u001b[39m(target_label))\n\u001b[0;32m--> 237\u001b[0m diagonal_ucg \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_append_ucg_up_to_diagonal\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 238\u001b[0m \u001b[43m \u001b[49m\u001b[43mcircuit\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mq\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msingle_qubit_gates\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcontrol_labels\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtarget_label\u001b[49m\n\u001b[1;32m 239\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 240\u001b[0m \u001b[38;5;66;03m# merge the diagonal into the UCGate for efficient application of both together\u001b[39;00m\n\u001b[1;32m 241\u001b[0m diagonal_ucg_inverse \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mconj(diagonal_ucg)\u001b[38;5;241m.\u001b[39mastype(\u001b[38;5;28mcomplex\u001b[39m, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/isometry.py:281\u001b[0m, in \u001b[0;36mIsometry._append_ucg_up_to_diagonal\u001b[0;34m(self, circ, q, single_qubit_gates, control_labels, target_label)\u001b[0m\n\u001b[1;32m 279\u001b[0m ucg \u001b[38;5;241m=\u001b[39m UCGate(single_qubit_gates, up_to_diagonal\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m 280\u001b[0m circ\u001b[38;5;241m.\u001b[39mappend(ucg, [target_qubit] \u001b[38;5;241m+\u001b[39m control_qubits)\n\u001b[0;32m--> 281\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mucg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_diagonal\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/uc.py:137\u001b[0m, in \u001b[0;36mUCGate._get_diagonal\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 132\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21m_get_diagonal\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 133\u001b[0m \u001b[38;5;66;03m# Important: for a control list q_controls = [q[0],...,q_[k-1]] the\u001b[39;00m\n\u001b[1;32m 134\u001b[0m \u001b[38;5;66;03m# diagonal gate is provided in the computational basis of the qubits\u001b[39;00m\n\u001b[1;32m 135\u001b[0m \u001b[38;5;66;03m# q[k-1],...,q[0],q_target, decreasingly ordered with respect to the\u001b[39;00m\n\u001b[1;32m 136\u001b[0m \u001b[38;5;66;03m# significance of the qubit in the computational basis\u001b[39;00m\n\u001b[0;32m--> 137\u001b[0m _, diag \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_dec_ucg\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 138\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m diag\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/library/generalized_gates/uc.py:178\u001b[0m, in \u001b[0;36mUCGate._dec_ucg\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 171\u001b[0m squ \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 172\u001b[0m HGate()\n\u001b[1;32m 173\u001b[0m \u001b[38;5;241m.\u001b[39mto_matrix()\n\u001b[1;32m 174\u001b[0m \u001b[38;5;241m.\u001b[39mdot(gate\u001b[38;5;241m.\u001b[39mdot(UCGate\u001b[38;5;241m.\u001b[39m_rz(np\u001b[38;5;241m.\u001b[39mpi \u001b[38;5;241m/\u001b[39m \u001b[38;5;241m2\u001b[39m)))\n\u001b[1;32m 175\u001b[0m \u001b[38;5;241m.\u001b[39mdot(HGate()\u001b[38;5;241m.\u001b[39mto_matrix())\n\u001b[1;32m 176\u001b[0m )\n\u001b[1;32m 177\u001b[0m \u001b[38;5;66;03m# Add single-qubit gate\u001b[39;00m\n\u001b[0;32m--> 178\u001b[0m \u001b[43mcircuit\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43munitary\u001b[49m\u001b[43m(\u001b[49m\u001b[43msqu\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43mq_target\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 179\u001b[0m \u001b[38;5;66;03m# The number of the control qubit is given by the number of zeros at the end\u001b[39;00m\n\u001b[1;32m 180\u001b[0m \u001b[38;5;66;03m# of the binary representation of (i+1)\u001b[39;00m\n\u001b[1;32m 181\u001b[0m binary_rep \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mbinary_repr(i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m)\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/quantumcircuit.py:5948\u001b[0m, in \u001b[0;36mQuantumCircuit.unitary\u001b[0;34m(self, obj, qubits, label)\u001b[0m\n\u001b[1;32m 5945\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(qubits, (\u001b[38;5;28mint\u001b[39m, Qubit)) \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(qubits) \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[1;32m 5946\u001b[0m qubits \u001b[38;5;241m=\u001b[39m [qubits]\n\u001b[0;32m-> 5948\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mappend\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgate\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mqubits\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcopy\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
+ "File \u001b[0;32m~/.qbraid/environments/moodys_yhl12u/pyenv/lib/python3.11/site-packages/qiskit/circuit/quantumcircuit.py:2546\u001b[0m, in \u001b[0;36mQuantumCircuit.append\u001b[0;34m(self, instruction, qargs, cargs, copy)\u001b[0m\n\u001b[1;32m 2543\u001b[0m expanded_qargs \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_qbit_argument_conversion(qarg) \u001b[38;5;28;01mfor\u001b[39;00m qarg \u001b[38;5;129;01min\u001b[39;00m qargs \u001b[38;5;129;01mor\u001b[39;00m []]\n\u001b[1;32m 2544\u001b[0m expanded_cargs \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_cbit_argument_conversion(carg) \u001b[38;5;28;01mfor\u001b[39;00m carg \u001b[38;5;129;01min\u001b[39;00m cargs \u001b[38;5;129;01mor\u001b[39;00m []]\n\u001b[0;32m-> 2546\u001b[0m instructions \u001b[38;5;241m=\u001b[39m \u001b[43mInstructionSet\u001b[49m\u001b[43m(\u001b[49m\u001b[43mresource_requester\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcircuit_scope\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mresolve_classical_resource\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 2547\u001b[0m \u001b[38;5;66;03m# For Operations that are non-Instructions, we use the Instruction's default method\u001b[39;00m\n\u001b[1;32m 2548\u001b[0m broadcast_iter \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 2549\u001b[0m operation\u001b[38;5;241m.\u001b[39mbroadcast_arguments(expanded_qargs, expanded_cargs)\n\u001b[1;32m 2550\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(operation, Instruction)\n\u001b[1;32m 2551\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m Instruction\u001b[38;5;241m.\u001b[39mbroadcast_arguments(operation, expanded_qargs, expanded_cargs)\n\u001b[1;32m 2552\u001b[0m )\n",
+ "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+ ]
+ }
+ ],
+ "source": [
+ "betti_values = []\n",
+ "expected = []\n",
+ "classical = []\n",
+ "for epsilon in np.arange(0.05, 0.5, 0.01):\n",
+ " simplex_trees_final = getSimplexTrees(point_clouds_final, epsilon)\n",
+ " laplacian_final, eigenvalues = get_laplacian(simplex_trees_final[43], 0)\n",
+ "\n",
+ " # print(laplacian_final)\n",
+ " \n",
+ " U_final = getU(laplacian_final)\n",
+ " non_zero_count_final = laplacian_qpe(U_final, 3)\n",
+ " betti_values.append(4*(8 - non_zero_count_final * 8))\n",
+ "\n",
+ " \n",
+ " classical.append(np.sum(eigenvalues < 1e-6))\n",
+ "\n",
+ " simplex_trees_final[43].compute_persistence()\n",
+ " expected.append(simplex_trees_final[43].betti_numbers()[0])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9d0a92e8-763f-4342-9b1c-6fc71bfeb7a4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "betti_values, expected, classical"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "31e47a34-6536-4f09-8138-b10067d9c6d7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#scaled betti values\n",
+ "plt.plot(betti_values)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4edd647c-9103-49d5-b6e3-68efa57590db",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(expected)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6b3bae78-ddb2-4f06-a7e5-57ba4b1b210c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(classical)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d1c2e4b5",
+ "metadata": {},
+ "source": [
+ "## *BONUS:* Explore future directions of quantum TDA\n",
+ "\n",
+ "The following is a non-exhaustive list of possible next steps for the quantum TDA pipeline. It is recommended to at least explore 1 option or sub-option.\n",
+ "\n",
+ "- **Find more of applications of TDA in finance**:\n",
+ "\n",
+ " There are several directions where to extend the analysis. Most work on time series analysis has used persistent homology, and more specifically, the $L^{P}$ norm of persistence landscapes, which can be used to detect early warning signals of imminent market crashes. This is precisely studied in the seminal work by Gidea and Katz, see [ref](https://arxiv.org/abs/1703.04385)\n",
+ " - Analyze financial correlation network and their degree of association with Betti curves or other topological features. From the time-series of multiple stock prices, we can build time-dependent correlation networks, which exhibit topological structures and might show some association to the evolution of betti curves or other topological data measures. Generally speaking, the cross correlations in a stock market will be in the form of a high-dimension topological space, with more complicated features. One can also think about other time varying financial graphs (e.g. cryptocurrencies). The following articles can help uncover more applications: \n",
+ " \n",
+ " - [Integral Betti signature confirms the hyperbolic geometry of brain, climate,and financial networks](https://www.arxiv.org/pdf/2406.15505)\n",
+ " - [Using Topological Data Analysis (TDA) and Persistent Homology to Analyze the Stock Markets in Singapore and Taiwan](https://www.frontiersin.org/journals/physics/articles/10.3389/fphy.2021.572216/full)\n",
+ " - Build a ML classifier or regressor on top of vectorized features such as Betti Curves (given their potential to identify trends, patterns or potential turning points in the market) to help with investment or risk management strategies. Show that Betti curves have some predictive skill, as key topological descriptors. See [ref1](https://arxiv.org/abs/2411.13881) and [ref2](https://www.sciencedirect.com/science/article/pii/S2405918823000235) for further information on the topic.\n",
+ "- **A hybrid and more NISQ-friendly quantum TDA pipeline**:\n",
+ " \n",
+ " QPE remains primarily theoretical. Its circuits are simply too deep to run on real hardware. Come up an with iterative or hybrid quantum phase estimation protocol or use tools that increase the algorithmic performance when running quantum circuits on real hardware. Benchmark them against textbook-QPE circuits. Here are some proposals to subtitute the QPE part:\n",
+ " - Variational Quantum Deflation (VQD) Algorithm: VQD is a quantum algorithm that uses a variational technique to find the k eigenvalues of the Hamiltonian H of a given system. [ref](https://quantum-journal.org/papers/q-2019-07-01-156/)\n",
+ " - Variational Quantum Eigensolver (VQE): Using VQE to determine the spectra of adjancency or laplacian matrix. Inspired by: [ref](https://arxiv.org/pdf/1912.12366)\n",
+ " \n",
+ " Finally, run some circuits on simulator + real hardware and compare the performance (runtime, noise effects, # resources) of the new proposal to the QPE solution of the above sections.\n",
+ "\n",
+ "- **A proper procedure of encoding the Laplacian matrix to the unitary**:\n",
+ "\n",
+ " In Step 3, encoding the exponential matrix is recommanded. An alternative approach is to conduct the Paulis decomposition on the Laplacian, then followed by Trotterization. Can you implement this approach in your pipeline? What parameters influence the accuracy? Can you optimize your code to minimize the circuit depth?\n",
+ "\n",
+ "- **Extend the quantum TDA to extract persistent Betti numbers**: \n",
+ " \n",
+ " Implement a quantum TDA algorithm for persistent Betti numbers. Esstimating the persistent Betti numbers is a more general task than estimating the Betti number and it is more practical for TDA. It is an open problem to construct a quantum algorithm for the persistent Betti numbers in a way that is preferable for NISQ devices, and the only current implementation of a quantum algorithm for persistent betti number is shown [here](https://quantum-journal.org/papers/q-2022-12-07-873/pdf/)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0de4a8f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# write your code here\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c098599c-7ab0-48b1-9e33-51d22502310d",
+ "metadata": {},
+ "source": [
+ "# This is the end of the challenge. Good luck!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 50,
+ "id": "e0ce4f0e-d7e5-41cc-aedd-2d9a8a4ef319",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAaJZJREFUeJzt3Xd4U9X/B/B3utJBF5QuKLRQKFD2qmUISKUgKvgTBL4oSxyIE0WpCoooIC4UFUSmAxAHONiUDWXvvSmUtsw2HXQl9/dH6U1ucpM2bdOMvl/Pk+fJPffc25NQmk/O+ByFIAgCiIiIiGyYk7UbQERERFQaBixERERk8xiwEBERkc1jwEJEREQ2jwELERER2TwGLERERGTzGLAQERGRzWPAQkRERDbPxdoNqAwajQbXr1+Ht7c3FAqFtZtDREREZSAIArKyshAaGgonJ9N9KA4RsFy/fh1hYWHWbgYRERGVw9WrV1G3bl2TdRwiYPH29gZQ/IJ9fHys3BoiIiIqC5VKhbCwMPFz3BSHCFhKhoF8fHwYsBAREdmZskzn4KRbIiIisnkMWIiIiMjmMWAhIiIim8eAhYiIiGweAxYiIiKyeWYFLGq1GhMnTkRERAQ8PDzQsGFDTJkyBYIgGL1my5YtUCgUBo+0tDRJve+++w7h4eFwd3dHTEwM9u7dW75XRERERA7HrGXNn376KWbPno3FixcjOjoa+/fvx8iRI+Hr64tXX33V5LVnzpyRLDkODAwUn//2228YN24c5syZg5iYGMycORPx8fE4c+aMpB4RERFVT2b1sOzatQv9+vVD3759ER4ejgEDBqBXr15l6g0JDAxEcHCw+NBNwfvll1/iueeew8iRI9GsWTPMmTMHnp6eWLBggfmviIiIiByOWQFLp06dkJiYiLNnzwIAjhw5gh07dqBPnz6lXtu6dWuEhITg4Ycfxs6dO8XygoICHDhwAHFxcdpGOTkhLi4OSUlJsvfKz8+HSqWSPIiIiMhxmTUkNGHCBKhUKjRp0gTOzs5Qq9X45JNPMHToUKPXhISEYM6cOWjfvj3y8/Mxb948dO/eHXv27EHbtm1x69YtqNVqBAUFSa4LCgrC6dOnZe85bdo0TJ482ZymExERkR0zK2BZvnw5fv31VyxZsgTR0dE4fPgwXn/9dYSGhmL48OGy10RFRSEqKko87tSpEy5cuICvvvoKP//8c7kanZCQgHHjxonHJXsREBERkWMyK2AZP348JkyYgMGDBwMAWrRogStXrmDatGlGAxY5HTt2xI4dOwAAAQEBcHZ2Rnp6uqROeno6goODZa9XKpVQKpXmNJ2IiIjsmFlzWHJzcyWTZQHA2dkZGo3GrB96+PBhhISEAADc3NzQrl07JCYmiuc1Gg0SExMRGxtr1n2JiIiocuUVqvHD1gs4npJp1XaY1cPy2GOP4ZNPPkG9evUQHR2NQ4cO4csvv8SoUaPEOgkJCUhJScFPP/0EAJg5cyYiIiIQHR2NvLw8zJs3D5s2bcL69evFa8aNG4fhw4ejffv26NixI2bOnImcnByMHDmykl4mERERlce6E2mYtqZ4Tunl6X2t1g6zApZZs2Zh4sSJeOmll3Djxg2EhobihRdewKRJk8Q6qampSE5OFo8LCgrw5ptvIiUlBZ6enmjZsiU2btyIHj16iHUGDRqEmzdvYtKkSUhLS0Pr1q2xdu1ag4m4REREVLUu3coRn+cVquHu6myVdigEU2lq7YRKpYKvry8yMzMlyemIiIiofM6mZ+HI1QycTFVh4c7LiG1QC7+OjoGTk6LSfoY5n99m9bAQERGR47udnY9eX22TlEXU9qrUYMVc3PyQiIiIJM6mZxuUWXvSLQMWIiIikriVnW9Q5uvhaoWWaDFgISIiIonkO7kGZT7uDFiIiIjIhny27oxBmbMV568ADFiIiIgIgCAIKFIbTwQbVtOjCltjiAELERER4fmfDyBmaiJOpapkz7/UPbKKWyTFZc1ERETV3K3sfGw4Wbyn37ebzwMAAr2V+HFYe2w4mY6XH4q0WsK4EgxYiIiIqrk3fjssPl91NBUAEOLrjlZhfmgV5medRunhkBAREVE1lleoxvZztwzKva28KkgfAxYiIqJq7N2/jsmW11Da1iAMAxYiIqJq6uqdXPx1KEX2nIuzdZcx62PAQkREVE0NX7jX6LnEUzeqsCWlY8BCRERUDeUWFOHizRzx+POBrSTnFbbVwcKAhYiIqDpauPOy5HhAu7oYHltfPP56cJsqbpFpDFiIiIiqoasy+wU9360hAMDH3QU9mwRWdZNMsq0pwERERFQlsvKLDMrq+Hlg3esPwtfDFU5W3jtIHwMWIiKiaubyrRwxQRwAjOgULj6PCva2QotKx4CFiIiomrh0Kwc/br+IJXuSxbJPn2yBJ9rUtWKryoYBCxERUTWg0Qjo8fkWg/JHW4bCzcX2p7TafguJiIiowsb8ekC23MvGMtoaw4CFiIioGlh3It2g7J3eTazQkvKxj7CKiIiIyk2tEQzKPu7fHE8/UF+mtm1iDwsREZGDS1PlGZS5OdtXCGBfrSUiIiKzHb2aYVDWLNSn6htSARwSIiIicnApGfcAAFFB3mgd5oe6/h5oXsfXyq0yDwMWIiIiB5d9P6ttZFANfDqgpZVbUz4cEiIiInJw2XnFAUtdfw8rt6T8GLAQERE5OFVeIQDA205yrshhwEJEROTA1p1Iw/L916zdjAqz31CLiIiIjNJoBMxYdwZztl4Qy/IKNVZsUcWwh4WIiMgB/XHwmiRYAYBHWoRYqTUVZ1bAolarMXHiRERERMDDwwMNGzbElClTIAiGGfRK/PXXX3j44YdRu3Zt+Pj4IDY2FuvWrZPU+fDDD6FQKCSPJk3sJ10wERGRrTl6LcOgzNPNueobUknMGhL69NNPMXv2bCxevBjR0dHYv38/Ro4cCV9fX7z66quy12zbtg0PP/wwpk6dCj8/PyxcuBCPPfYY9uzZgzZt2oj1oqOjsXHjRm3DXDhaRUREVF6XbuUYlPl7ulmhJZXDrKhg165d6NevH/r27QsACA8Px9KlS7F3716j18ycOVNyPHXqVPz999/4999/JQGLi4sLgoODzWkOERERyXh3xTHsPH9bUvb90Lbw9XS1UosqzqwhoU6dOiExMRFnz54FABw5cgQ7duxAnz59ynwPjUaDrKws1KxZU1J+7tw5hIaGokGDBhg6dCiSk5ON3iM/Px8qlUryICIiomJL9hh+hnaPqm2FllQes3pYJkyYAJVKhSZNmsDZ2RlqtRqffPIJhg4dWuZ7fP7558jOzsZTTz0llsXExGDRokWIiopCamoqJk+ejK5du+L48ePw9vY2uMe0adMwefJkc5pORERULaRm3pMt93C13/krgJkBy/Lly/Hrr79iyZIliI6OxuHDh/H6668jNDQUw4cPL/X6JUuWYPLkyfj7778RGBgoluv20LRs2RIxMTGoX78+li9fjmeffdbgPgkJCRg3bpx4rFKpEBYWZs5LISIickhbz9yULVcoFFXcksplVsAyfvx4TJgwAYMHDwYAtGjRAleuXMG0adNKDViWLVuG0aNH4/fff0dcXJzJun5+fmjcuDHOnz8ve16pVEKpVJrTdCIiomrBx0M6T8XFSYGnH6hvpdZUHrMCltzcXDg5Sae9ODs7Q6MxnYhm6dKlGDVqFJYtWyZO2DUlOzsbFy5cwDPPPGNO84iIiKq9l349KDne/34c/Ox4dVAJswKWxx57DJ988gnq1auH6OhoHDp0CF9++SVGjRol1klISEBKSgp++uknAMXDQMOHD8fXX3+NmJgYpKWlAQA8PDzg61u8tfVbb72Fxx57DPXr18f169fxwQcfwNnZGUOGDKms10lEROTwbmfnS45DfN3h62G/K4N0mbVKaNasWRgwYABeeuklNG3aFG+99RZeeOEFTJkyRayTmpoqWeEzd+5cFBUVYezYsQgJCREfr732mljn2rVrGDJkCKKiovDUU0+hVq1a2L17N2rXtu8ZzURERFWp++dbxOfBPu7YMr673c9dKaEQTKWptRMqlQq+vr7IzMyEj4+PtZtDRERU5fKL1Ih6f614vOfdngjycbdii0pnzuc39xIiIiJyAKdTsyTHth6smIsBCxERkQOYt+OS+HzJ6BgrtsQyGLAQERHZubxCNf49cl087hQZYMXWWAYDFiIiIjuXdEG7b1DvaMfcl48BCxERkZ0pKNJg5aEU3MjKAwCMXLRPPNe2vp+VWmVZDFiIiIjszLebz+P13w7jydm7DM61q+9vhRZZHgMWIiIiO7P2eCoA4Oqde1h7PE1yrl39mtZoksUxYCEiIrIzao02hdqLvxwQn/dtGWKN5lQJBixERER2JrdALVvexQFXB5VgwEJERGRnMnILZcuVLo77se64r4yIiMhB3SuU72FxY8BCREREtk7p4mztJlgMAxYiIiIHwR4WIiIisqiT11UYPDcJKRn3TNbLypOfvwJwDgsRERFZ2CPfbMfui3fQefomk/WGLdhr9BwDFiIiIrIYQRBKr3RfWmae0XMcEiIiIiKLUGsERCSsFo9jIkxnqjUV27i7ctItERERWcB/R69Ljr2ULkbrFqk1SL+/4WGgt9LgvAcDFiIiIrKE/EKN5LhIY7wLJfNeodjDIleNAQsRERFZxMVbOZJjtUZjpCaw88JtAICPuwtuZecbnOeQEBEREVU6QRAwZ+sFSVmh2ngPy6tLDwEAVHlFiA71EcujQ33Qvr4/3F0d92PdcV8ZERGRjXtv5XGDMrWRISH9lUTPP9hAfP7vy13w+4uxUCgUldtAG2J8Zg8RERFZ1JI9yQZl+nNY0lV5GP/HUUQF1ZCU924ejPb1/dEu3B9OTo4bqJRgwEJERGQDFo3sgBEL96FILZ3D8uO2i9h29ia2nb0plr3TuwmULs74Y0ynqm6m1XBIiIiIyAoKijRwud8zsurVLnBxKv5I1h8Sup1TYHDti90aGJQ5OgYsREREVnDuRpY4/NM02AcuzsXBS6FeD0uB2nDVkCPPVTGGAQsREVEVK1Rr0PebHeKxk5NC7G3R72HRHyKqrhiwEBERVbHt524alLk4F38k6y9rLtI7/nNMrOUaZsMYsBAREVWxUYv2i88bBRav/jHWw6I/JFTLyzAlf3XAgIWIiKgK5RWqJcftw4s3OyyZw5KmykNOfpF4Xr+HxdPNcbPZmsKAhYiIqArN3HhOcty/dSgAbQ8LAHy3+bz4XH8Srkc1DViYh4WIiMhK1rzWFU1DilPslyxrBoCrd++Jz1Mz8yTXeLpVz49us3pY1Go1Jk6ciIiICHh4eKBhw4aYMmWKQbpgfVu2bEHbtm2hVCoRGRmJRYsWGdT57rvvEB4eDnd3d8TExGDv3r1mvRAiIiJbU6TWGPSQlOwdNCy2vhisAICzTg+L7udqSoY2eNGvV52YFbB8+umnmD17Nr799lucOnUKn376KWbMmIFZs2YZvebSpUvo27cvevTogcOHD+P111/H6NGjsW7dOrHOb7/9hnHjxuGDDz7AwYMH0apVK8THx+PGjRvlf2VERERW9PfhFES+twaPfL0d9wrUGPfbYfy+/6o4B6VlXT9JfVdn7Uey0/08K6V1CFQnZvUr7dq1C/369UPfvn0BAOHh4Vi6dKnJ3pA5c+YgIiICX3zxBQCgadOm2LFjB7766ivEx8cDAL788ks899xzGDlypHjNqlWrsGDBAkyYMKFcL4yIiMiaXlt2GABw7kY2PvrvBP46lIK/DqWI5zs1rCWpL9dzortCqG09P8wY0MoyjbUDZvWwdOrUCYmJiTh79iwA4MiRI9ixYwf69Olj9JqkpCTExcVJyuLj45GUlAQAKCgowIEDByR1nJycEBcXJ9bRl5+fD5VKJXkQERHZig0n0yXHS/deNagT7OMuOdaddFvSr5Kbr11RtPyFWEQGSjdArE7M6mGZMGECVCoVmjRpAmdnZ6jVanzyyScYOnSo0WvS0tIQFBQkKQsKCoJKpcK9e/dw9+5dqNVq2TqnT5+Wvee0adMwefJkc5pORERUZb5OPFtqHf0dlkuWNQOAWlPcs3L5do7O+eq9sNesV798+XL8+uuvWLJkCQ4ePIjFixfj888/x+LFiy3VPlkJCQnIzMwUH1evGkauREREFWHu/JF52y9i8a7LyM4vwvEU0z3/zz9ouHmh7iqhkuRx32+5YFYbHJlZPSzjx4/HhAkTMHjwYABAixYtcOXKFUybNg3Dhw+XvSY4OBjp6dKusfT0dPj4+MDDwwPOzs5wdnaWrRMcHCx7T6VSCaWyemb6IyIiy7tyOwfdPtsCANj7Xk8EerubrP/XwWv4eNUpAICbS+l9AW3r+RuUSXtYgKQLtw2Glqozs3pYcnNz4eQkvcTZ2RkajfGNmWJjY5GYmCgp27BhA2Jji/dCcHNzQ7t27SR1NBoNEhMTxTpERERV5XZ2PuK+3Coef/jPiVKvefP3I+LzRTsvl1rfRWaCrW6Zm4sCQ37cLR63qONb6j0dnVkBy2OPPYZPPvkEq1atwuXLl7FixQp8+eWXeOKJJ8Q6CQkJGDZsmHj84osv4uLFi3j77bdx+vRpfP/991i+fDneeOMNsc64cePw448/YvHixTh16hTGjBmDnJwccdUQERFRVZn0zwnJBoQ1lKYHIzQaAbqjR2fSs0r9GRG1vQzKFAoFSmKWEF8Pyblrd3NLvaejM2tIaNasWZg4cSJeeukl3LhxA6GhoXjhhRcwadIksU5qaiqSk5PF44iICKxatQpvvPEGvv76a9StWxfz5s0TlzQDwKBBg3Dz5k1MmjQJaWlpaN26NdauXWswEZeIiMjSVh1NlRwv33/N5HLiWZvOGz1nTMPa8qt9XuoeiW83nzfYANHHw9Xsn+FozApYvL29MXPmTMycOdNoHbkstt27d8ehQ4dM3vvll1/Gyy+/bE5ziIiISnUw+S7+OHANb8dHwc/TzWTdrWdvmn3/49czZcu7R9XGljPS+73ft6lBwjhdJblY9LPjulbzFUIANz8kIiIH93/f78KSPcn4fP0Zk/XyCtUYvkCbCLVNPT8ApuePaDSC0Ymx7z7SFLOGtBGPe0TVxuiuDdAxoqbR+7nen3ibdPG2ybZWRwxYiIioWrh0K8fk+bf/OCo5Hts9EgBgauueF385YPScp5szHm0ZIh5/+VTrUtvofH9hy8Wb0rYyRT93ayYiIgeWmVsoPq/lZTodxj9HrkuOS5YZH7mWiV3nb6FduD+ULs7ieUEQsF6nd6VtPT8cTM4Qj73cXKBQKLD+jQdxr0ANfy/Tw1GAtofFsJz9C3wHiIjIYV3SyRRbUKRBSsY9FBRJ54cUqTU4oTcP5Ydn2sFNJ0j437w9GL5gL1R52gDo/2bvklzzUb/mkmOP+5scNg7yRqswvzK119hOzDMGtCzT9Y6MAQsRETmsKzoBy77Ld9B5+iY899N+sUyVV4hun21B3292iGXr33gQ8dHBBqnwd1+8g7gvtPlZDun0pgBA8zq+eLZLhHisLEMCOX1y6fe3je9hcqJudcGAhYiIHFZaZp74/HZOAQDpSqADV+4iJeOe5JrGQd4ApJlnS9zIygcASU8LAPz2/AMG1ygUJia/GCGXUK6GO2dvAAxYiIjIgU1bI7+Jbsmw0Pn0bEn5mw83Fp9n5xUZva9uIARoE725VXCuiVzA4unmLFOz+mHAQkRE1U66Kg9/H07BJ6tPScob6CR0ax9uuN9PiZ+TrkiO3V2LP067RwUCKNt+QnL0e3UaBdaAuysDFoCrhIiIqBr6ff9VfCOTodbfS5tR1tPNBS92a4g5Ww13TP55tzZgCavpgVo1ilcgtavvj+UvxCKspofBNWXhordf33NdDXd1rq4YsBARkUPKK1QbPScXrABARIB0jx/fMqTETxzXXbK6x1RiuNLoDwkZWzVUHXFIiIiIHFKTiWvNqt8xoqbBpoMeroYfk7pp82Mb1Cr38I8c/VVCchN/qyv2sBARkcMxJzPsqw9FYlyvKNlzcvNHUjO0E25n/a+NwfmK0O9hYcI4Lb4TRETkcHSz1kbdX6ZszBs6K4P01ZTJTvvgZ5vF5wE1TGfPNZd+jwqHhLQYsBARkcN5bdlh8fnER5uZrGsqX0rnyAA0Dqph9Hxl0w9Q5JY5V1cMWIiIyKE90KD8k2C9lC5Y/0a3SmyNafpDQHKZb6srvhNERORQinQmxX4zpI3JD/1A77IN6ex4p4dB2fIXYs1vXCnYw2IcAxYiInIYgiDgrd+PACj+sO/bIkRyvmujAPH5T6M6Ys1rXct037r+nhjYrq6kLNjHvYKtNeSql4eFc1i0uEqIiIgcxmvLDosTbuv4e4gf+E+1r4sLN3OwYEQHvLbsECICvPBg49pm3dtHLyeLtwX2+NEPUFy5rFnEgIWIiBzCrex8yeqgG6p88fmMAa3E598PbVeu++vvvlyWpHLm0g9QnJ04EFKC7wQRETmEp+ftkRwvHtWxUu+vdNHmZImJqAknCwzXcA6LcexhISIiu5dfpMaZ9Czx+NRHveFRybsc62a0tUTvCiC3SogBSwkGLEREZPcu3MhBSXLbS9MeMZlbpbx0h4Q8KzkYKqEfoLCHRYtDQkREZNcEQcAj32wHALQK87NIsAIASp19hTzcLPN933BIiB/TJfhOEBGRXdt3+a74vK6fh4maFeOmM1yz7kSaRX6G/rJm10rcWNHecUiIiIjszu3sfLT7eKNB+Xt9m1rsZyp1NkK8k1NgkZ/hrDckZIlcL/aKoRsREdmdYQv2ypaHWrCHRXcOi5uFej50e1g8XJ2ZOE4He1iIiMhu5BWq0WryeuQXaQzOTSplk8OK0g1Yvv9fW4v8DN0A5bFWISZqVj8MWIiIyC6oNQKaTFxr9PzIzuEW/fm6vSoNantZ5GforgrSCBb5EXaLQ0JERGQXLt3KNnneUquDSugmjnN3tcyyZt1kdAIDFgkGLEREZPPuFagR9+U2g/Jfno1B10YB2DbecDflyqbb++FhoYBFl8CIRYJDQkREZHPu5hTgyw1n0TrMD4+2CkHTSdqhoKggb4yPj0KjoBqoX8sLXXR2YLYk3fChsrPoykm+k2vxn2FPGLAQEZHN+efIdfy8+wp+3n0FS/cmS8692L0B4poFVXmbNDo9HvobIVqCpZZO2yuz3vHw8HAoFAqDx9ixY2Xrd+/eXbZ+3759xTojRowwON+7d++KvSoiIrJrH/xzQny+/8pdybmqGI6R0zjIGwBQQ+li8fkyAOBtof2K7JVZPSz79u2DWq0Wj48fP46HH34YAwcOlK3/119/oaBAGyHevn0brVq1Mqjfu3dvLFy4UDxWKpXmNIuIiKqR5nV8rfJzayhdcGRSL4vlYDHAOSwSZgUstWvXlhxPnz4dDRs2RLdu3WTr16xZU3K8bNkyeHp6GgQsSqUSwcHB5jSFiIgclLGhkC6RAXi7dxTq+ntWcYu0fD3Z62Et5Q4TCwoK8Msvv2DUqFFl7hqbP38+Bg8eDC8v6fr1LVu2IDAwEFFRURgzZgxu375t8j75+flQqVSSBxEROYarRiabTn2iBVrW9avaxpDNKPek25UrVyIjIwMjRowoU/29e/fi+PHjmD9/vqS8d+/e+L//+z9ERETgwoULePfdd9GnTx8kJSXB2Vl+nHLatGmYPHlyeZtOREQ27HrGPQBA23p+mDGgJdIy89Es1Ac1vdys3DKyJoVQzoXe8fHxcHNzw7///lum+i+88AKSkpJw9OhRk/UuXryIhg0bYuPGjejZs6dsnfz8fOTn54vHKpUKYWFhyMzMhI+PT9lfBBER2ZQ7OQVoO2UDACAmoiZ+eyHWyi2qeuETVgEAWoX54e+xna3cGstSqVTw9fUt0+d3uYaErly5go0bN2L06NFlqp+Tk4Nly5bh2WefLbVugwYNEBAQgPPnzxuto1Qq4ePjI3kQEZH9013CXN3nnHLbQ6lyDQktXLgQgYGBkuXJpvz+++/Iz8/H008/XWrda9eu4fbt2wgJ4aZPRETVRdKF2xjy425J2fBO4dZpDNkks3tYNBoNFi5ciOHDh8PFRRrvDBs2DAkJCQbXzJ8/H/3790etWrUk5dnZ2Rg/fjx2796Ny5cvIzExEf369UNkZCTi4+PNbRoREdkp/WAFAB5pUb1Xj1ZBqhe7YnYPy8aNG5GcnIxRo0YZnEtOToaTkzQGOnPmDHbs2IH169cb1Hd2dsbRo0exePFiZGRkIDQ0FL169cKUKVOYi4WIqBrb827PKknOZouahfjgZKoK/9emjrWbYlPKPenWlpgzaYeIiGxP5+mbkHJ/ddCvo2PQObJq9geyRVl5hThxXYWO4TUluzc7InM+v7mXEBERWVWhWiMGK3+P7YxWYX7WbZCVebu74oEGtUqvWM1UUX5hIiIieVNXnxKfNwnxtmJLyJYxYCEioiqRea8QMVM34r0Vx6DKKxTLF+68LD5XulhnY0OyfQxYiIjI4gRBQKvJ65Guyseve5LR8sP1uHwrB0VqjVinhpKzFMg4/nYQEZFFCYKAiITVBuXdP9+CVx+KFI/3vCuf3ZwIYA8LERFZ2MHkDKPnvtmkzWruxR4WMoEBCxERWdS59Czx+fj4KNk6bev5VVFryF4xYCEiIny76RxeXnIQ+y/fwcpDKRi3/DCSb+dWyr1vZhVvVjuofRjG9oiUrRPbkMt4yTT2vxERVTNHr2Vg4t8nsHBEB2gEATFTE6HWFOcQ/e9oqljvbHoW/nula7l/TlZeId5dcRz7L98BAAR4uwEAekTVxuYzNyV1uTqISsOAhYiomnn8250AgLZTNpisdzxFheX7ruKpDmEAgP2X7+DFXw7i1Z6RGBYbbvS62VsuoKaXK5Lv5OLfI9fF8ppexVuufDOkDW5m5UMA0POLrQCAvEJ1BV4RVQcMWIiIqoE9F28j+U4uHm4WZNZ1b/95FK3r+aFxkDemrj6FW9n5mPT3CUz6+wROT+kNd1dpz0jy7Vx8uva07L0u3coGUJzJ1dvdVXIut4ABC5nGOSxERA7u5HUVBs3djfF/HMVfB1PMvv7s/Umz+qt9Zm+5YFB369kbRu8TWbuG0XPhtTzNbhdVL+xhISJyYDM3nsXMjefE4wPJd2Xrhfi6IzUzT/acu5H5JV8nnkORRoPx8U0AABqNgIl/nzDalv4yuw//OaYTtp69if/F1Dd6HRHAHhYiIoemG6wAwCqdSbW6khJ6Yu978onbcgvVEARB9tx3m7W9LCsPm+698XAzDHza1ffHuIcbw82FH0dkGn9DiIgcVH5R2eaFrH/jQQBAoLc7loyOAVCcJv+hJoEAgLwCNXLKMMdk3PIjJs+7OfMjh8qPvz1ERA5qzpaLpdZZ+3pXNA7S7pDcKTIAO97pgZ3vPAQXJwUAoFCjQZqR4SJzKBSKCt+Dqi/OYSEicjAFRRq4uTjhq41njdZZ9vwD8Pd0Q1Swt8G5uv7FE2BdnIsDDLVGwPebzxvU01Wos4nhgHZ1kdCnCSb/exJ3cwuw/dyt8rwMIgn2sBAROQi1RsCoRfvQ+P01+OPANZN1w2t5yQYrulycij8iCtUC/jqknZ8y8dFmknqCIEjmxnzcvzlq1VDimyFt0PP+sBJRRTFgISJyEA3fXY1Np4uXFb/1u+n5JN7upXewlwwJTfnvpKR8aEw9yXF+kQbbzmkz1+rmZhnQPgxdGwXg4/7NS/15RKZwSIiIyI5dvpWDofP2ICuvUPa8q7MChWrtCp/oUB+80K1hmXZGllsX9OmTLeDu6owLUx9Bw3dXAwBy8ouM5nepoXTBz8/GlP5CiErBgIWIyI51/3yLyfOB3u54Pa4Rxv9xFJ8NaImB7cPKfO/NZwyTwD3SIgQA4OykgNLFCflFGkmW2shA48nhiCqCAQsRkZ3aevamQZlCAeimTAn1c8fA9mHo3TzYIB1+aTQawz4WLzftx4aX0gX5RQW4V6hGWE0PXL1zD1OfaGHWzyAqK85hISKyQ5m5hRi+YK9BeVxT6V5Bgd7uAGB2sAIASlfDRG9OTtqlyR73z2flFeLqnXsAgJpe5v8corJgwEJEZIeW7UuWLa9fU7onT21vZbl/hrur6Y8Iz/uZa99fqU3HX0PJgIUsgwELEZEdOn5dJT4f0SlcfJ6VVyQGEgBQ08ut3D/DtZTMtOduFO++fCpV25aKBEhEpjBgISKyQ/d0Jro20cmn0is6CF8PbiMel2U1kFHy2weZ5OzEbLZkGQxYiIjsTJFag53ni7PH/vNyZzyls/InOtRXnFsCAN4VCFg0RjY8LOFThlwuRJWFv21ERHYiI7cACijwyrJDuFdY3MPSPNQXTk4KbB3fHfcK1Qj2dUdKRq54TUV6WHTzt8gZ1SVCshu07tAUUWVjwEJEZOMEQcD1zDw8OGMz1HpLjUtW7dSv5SWW6Waa9VIarvQpq5Gdw/HxqlNGzzvpbWbYOsyv3D+LqDQMWIiIbNysTefx5QbjGxnq89TJlVKWFPzGjOocgfbhNbHuRBpmb7mAIR2lKfn1h4w83MofHBGVhgELEZGNMxasRIf6yJZ7SHpYyv9n3slJgdZhfmhZxxePNA9BkxDpZon6PSzuMnlbiCoLAxYiIhuUrspDQZEGwb7usuefal/XaFZZT51hIN3MtOXl5KRAi7q+BuX6vTeuXCFEFmTWKqHw8HAoFAqDx9ixY2XrL1q0yKCuu7v0P58gCJg0aRJCQkLg4eGBuLg4nDt3TvZ+RETVwYWb2YiZmoiuMzaj0XtrZOvMGNAKLkbypPi4u+Kl7g0xLLY+6vp7WKydgzqESeatFMqk8ieqLGYFLPv27UNqaqr42LBhAwBg4MCBRq/x8fGRXHPlyhXJ+RkzZuCbb77BnDlzsGfPHnh5eSE+Ph55eXnleDlERPbvt31XK3yPt3s3wUf9mkOhsFyvh6ebC1aO7SweqzUai/0sIrP6CmvXri05nj59Oho2bIhu3boZvUahUCA4OFj2nCAImDlzJt5//33069cPAPDTTz8hKCgIK1euxODBg81pHhGRQ9CfGwIA/4uphxuqPGw8ZbiDsrU93ioUx69nolPDAGs3hRxYuQc3CwoK8Msvv2DcuHEmI/js7GzUr18fGo0Gbdu2xdSpUxEdHQ0AuHTpEtLS0hAXFyfW9/X1RUxMDJKSkowGLPn5+cjPzxePVSqVbD0iInvzzh9H8dt+wx6WJ9rUQaFag5SMPEx6tJkVWmbcN0PaQBAEi/bmEJU70+3KlSuRkZGBESNGGK0TFRWFBQsW4O+//8Yvv/wCjUaDTp064dq1awCAtLQ0AEBQkHR30aCgIPGcnGnTpsHX11d8hIWFGa1LRGQvDly5KxusAICfhys6NQzAmte6IrZhrSpuWekYrJCllTtgmT9/Pvr06YPQ0FCjdWJjYzFs2DC0bt0a3bp1w19//YXatWvjhx9+KO+PBQAkJCQgMzNTfFy9WvHxXiIiaxswZ5fRc76e3AWZqrdyDQlduXIFGzduxF9//WXWda6urmjTpg3Onz8PAOLclvT0dISEhIj10tPT0bp1a6P3USqVUCq5IygRORZTW/f4ejBgoeqtXD0sCxcuRGBgIPr27WvWdWq1GseOHRODk4iICAQHByMxMVGso1KpsGfPHsTGxpanaUREdumX3VdMnle6MCkbVW9mBywajQYLFy7E8OHD4eIi7aAZNmwYEhISxOOPPvoI69evx8WLF3Hw4EE8/fTTuHLlCkaPHg2geMzz9ddfx8cff4x//vkHx44dw7BhwxAaGor+/ftX7JUREdmR91ceF59/NagV1r/xIHpHy6+wJKqOzB4S2rhxI5KTkzFq1CiDc8nJyXBy0sZAd+/exXPPPYe0tDT4+/ujXbt22LVrF5o1085wf/vtt5GTk4Pnn38eGRkZ6NKlC9auXWuQYI6IyFHdzs6XHD/Rpi4AwJmZY4lECkEwNWpqH1QqFXx9fZGZmQkfH/m9NYiIbFXiqXQ8u3g/AOD0lN7injznb2Sh37c78fQD9ZHwSFNrNpHIIsz5/OZeQkREVrT2eCpWHStO49A4qIZkA8HIQG8c+aCX0RT8RNUJAxYiIivJyivEi78cFI8bBXob1GGwQlSM/xOIiKrYvst38MrSQzh5XZql29+LS5eJjGEPCxFRFcrOL8LAOUkAgH+PXJec83FnwEJkDHtYiIiq0KKdl4ye8/d0q8KWENkXBixERFXo8/VnjZ6r4+9RhS0hsi8MWIiILOyX3VeQ8NdRFKo1Jusx6wqRcZzDQkRkQWqNIGaxLS29vtr+02IRWQx7WIiILChdlSc+X7TrMgBAYaQrJcSXQ0JExrCHhYjIglIz8wzK/D3dcCenQDye0i8aGgFoV9+/KptGZFcYsBARWdCRqxkGZf6erhjYvi5+2HoRADCoQz24ubDDm8gU/g8hIrKgs+lZBmX+nm4Y060hACCspgeDFaIyYA8LEZEFXbmda1C2/8pd+Hm64fL0vlZoEZF9YlhPRGRBd3MLSq9ERKViwEJEZCFzt13A6bTiIaE5T7cTy5c9/4C1mkRktxiwEBFZwInrmZi6+rR4HB7gKT738+SeQUTmYsBCRGQBaXrLmevX9BKfe7lx+iCRufi/hshGnbieiVvZBejWuLa1m0Jmyi0owrOL94vHLer4wsPNGUM6hiGvUIO63DOIyGwMWIhsQFZeIZLv5CI61BcAIAgCnpy9C3mFGvw5phMTitmZztM3SY693Yv/1E77v5bWaA6RQ+CQEJENaPHhevT9ZgfmbS9OJHbt7j3kFRZvlJd04ZY1m0blcDe3UHLs7MRtDYkqigFLOQjcoIwqgVoj4LN1p7HqaKpY9vGqUwCAGevOiGWfrz+LrWdvVnn7qHzyCtUGZU+0qWOFlhA5Fg4JmelMWhYGz03CmO4N8fyDDa3dHLJjC3dewnebL0jKejULAgD8e+S6pHz4gr1MMmYn1p1IAwC4uzph+9sP4dKtHHSMqGnlVhHZP/awmGnmxrO4m1soWa6oLzu/CE/9kCR27xPJKelN0aVQAFdu5xiU8xu6fTh/IwuvLTsMAPBxd0VtbyWDFaJKwh4WM2l0hoM2nkxH3P1vxADw47aLKFBrsPvibey9dAd7L93B6K4NrNFMsnE7z8vPS8nJVyM7v8igfMWhFBQUafBSj4bixFyyPcv3XxOf38jKt2JLiBwPe1jMpHRxFp+P/km7bDE7vwifrD6Fz9adwfZz2g+jIrWmSttH9mHovD2y5Vn5RShSy8+RWnUsFUPm7rZks6gCMnMLMXcbe1WJLIUBi5ncXeXfsltGvk2l81sW6UlX5Rk/l5mH/CLjQa4qrwjHUzIt0SyqoN/2J0uOW9ZlTxhRZWLAYiZPIxkqk+8Y7sgKAEeuZliwNWSPPl2rnf/UMVw6vyFNlYenfkgCANTx80B8dBD0jfn1gGUbSOVyPUMaiD7c1PDfjojKj3NYzJRfpF2yWENZ/PZNX3Mac7ZekK2/cOclPNIipEraRvbhr4Mp4vPfXngACoUCZ9KyED9zm6ReSsY9zBzcGutOpEvKb6jYa2eLTl5XAQBefSgSdfw98ESbulZuEZFjYQ+LmVR52gmR2flFmLbmlNFgBQDq1/Iyeo6qt7r+HlAoihOKebo5y9aRK2caIOu5nZ2Pvw+nIPOeNDFcVl4h9l6+AwDoFBmAQR3qwc2Ff16JKhN7WEzIK1Sj5xdbkZJxDx8+1gwjOkdIknwBwA9bTU+yuyeTRIqqn9lbLqBeTU/0bRmCsJoeuHrnHr4a1Fo8L7d7b4s6vrKb5GkYsVhN/MztuJWt7eE6+mEv+Li7YuLK42JZizqcu0JkCfwKYIKzkwIpGfcAAB/+e9Ksa3s2CQQArD+Rxsy41dyP2y7i07WnMXbJQVy9k4urd4p/p0L9tBvgebsbBizurk5wdzXsYSnS8PfJWnSDFQD4cv1ZAMDB5AyxzEvJ74FElsCAxQQXvf0/VHmFaBBQtiGeknkrhWoBf+rMWaDqRa0R8MlqbYK4rxPPic+DvJUmr32peyS8lIYBSy0vt8prIJklvJan5Pj8jWxk5xeJk+7fiGtsjWYRVQtmBSzh4eFQKBQGj7Fjx8rW//HHH9G1a1f4+/vD398fcXFx2Lt3r6TOiBEjDO7Xu3fv8r+iSlQyv6DEd5vPl3kTswca1hKfv/X7Efxx4BpyCwwTgpFj+37zecnxHwe0icVcnKX//WY8qd3J98m2ddGjSaBsz4uHkfkuQHGARJZz+bZ0NaCbixNmbdIGoSF+7lXdJKJqw6yAZd++fUhNTRUfGzZsAAAMHDhQtv6WLVswZMgQbN68GUlJSQgLC0OvXr2QkiLtcejdu7fkvkuXLi3ny7GsH7ZexPmb2bLnXureEFP6RYvHob7SP1xv/X4Ek/8xb1iJ7Neu87cQPmEVvthwVvb80w/UMyhrGuIjPh/VJVx87qYX2OhP+Cyx5lgqoj9YizXHUmXPU8XcySkwKNt0+oZkHpuSE22JLMas/121a9dGcHCw+Pjvv//QsGFDdOvWTbb+r7/+ipdeegmtW7dGkyZNMG/ePGg0GiQmJkrqKZVKyX39/f3L/4oszNh0FA9XZ3SPKp630iDAy6B3BgB+23/Vkk0jG/I/I5lsS9Sr6WlQ5qkz/KM7d0X/VykrrwgamZ6UMb8eRF6hBmN+PWhma0mtETBv+0WjWyaoNQJGLtT2DjesLT80zOlqRJZT7q8DBQUF+OWXXzBq1CjZD2c5ubm5KCwsRM2a0mRZW7ZsQWBgIKKiojBmzBjcvn3b5H3y8/OhUqkkD2tzd3VGWE1P7Hm3J1a/1hUAMH94eyu3iqxB7pu4/iqg/jKbGequCNL9pr7q1eLfp1Y6mVP/OKgdWgKAzadvlK+xBACY9PdxfLzqFIbO24Pb2YZ5bg5fzcCRa9oMwwPahcne526u4b89EVWOcgcsK1euREZGBkaMGFHma9555x2EhoYiLi5OLOvduzd++uknJCYm4tNPP8XWrVvRp08fqNXGlwNPmzYNvr6+4iMsTP6PR1Vyvz+vIMjHXfx2rD/XwNW5bIEd2bcHpkl7EEd3iUCf5sHi8eOtQhHobTjXQff3xUnnS0BkYA0ceD8Of4zpJJbpLqMFgDeWH5Ycc2WaeX7do02r3+7jjXjnj6OS80kXpD0vDzcLlL3PQ03ky4mo4sodsMyfPx99+vRBaGhomepPnz4dy5Ytw4oVK+Durv1jPXjwYDz++ONo0aIF+vfvj//++w/79u3Dli1bjN4rISEBmZmZ4uPqVesPtbjLjF3rL0nV3TiRHFNeoRoFOnsBbX+7B97r21SScNBYbh4vnYBFv0emVg0lXHXmsuifz8iVzmtZezzN/MZXU7rZq0voD99uOKnNNvzKQ5GIDPSWbJvgrXTBjnd6MFEkkQWVK2C5cuUKNm7ciNGjR5ep/ueff47p06dj/fr1aNmypcm6DRo0QEBAAM6fP2+0jlKphI+Pj+RhKd8PbVumensv3TEo89ALWLLzi/jN14Fl3itEm482SMrCanpCoVCgUCeIkcutAhSvGto2vgc2vdnN6J5VJdJV+fgp6TIEQZAESCU4j6XsbpayQWlaZp44HLTipU54s1cUAOkQ3i+jY1DX33BeEhFVnnIFLAsXLkRgYCD69u1bat0ZM2ZgypQpWLt2Ldq3L31Ox7Vr13D79m2EhNjG/juPtAgxWKUhRy0TiOgHLADw4i8HcCMrDxqNgDxmwXUon607Lek9eUJnnoruh1lADeN5VOrV8kSD2jWMnn+5R6T4fNLfJ/DcTwfwxYYz5W0ywbB3St+Lv2g3m6zlpc2doztJuiZz4xBZnNkBi0ajwcKFCzF8+HC4uEi/BQ4bNgwJCQni8aeffoqJEydiwYIFCA8PR1paGtLS0pCdXbw0ODs7G+PHj8fu3btx+fJlJCYmol+/foiMjER8fHwFX1rl6doowOi5vi2LA6uJfZsZnJP7Jr3uRDomrjyOsUsOosMnG3EjK8+gDtmnX3YnS44/0lnmPrZHQ/H5mG4NUV7urtL/shtPpctuD+HjXj2yrWo0guyQjjleXXbI6L3XHk/FYZ0d1/29tENx/p7aIKUGs9sSWZzZAcvGjRuRnJyMUaNGGZxLTk5Gaqo2B8Ts2bNRUFCAAQMGICQkRHx8/vnnAABnZ2ccPXoUjz/+OBo3boxnn30W7dq1w/bt26FUms4CWpU+G9gKDfSWMf73ShdseONBfPe/trg8vS/8Zb5hyfWwAMVBy5rjacjKK5IkEiP7lZ0vTQq45LkYSdK3WjWUuDy9Ly5OfQSBPuVPLmZsOAkontRdknyuuiSQG75wL7p+uhk5+eVLyngqVYWLN3Nkz+UUFOHFX6RDa7qBSbBOrqUa1SRAJLIms/+X9erVy+g8DP2JspcvXzZ5Lw8PD6xbt87cJlS5ml5ueKl7JN76/YhY1rwMG5y5u5VhKEldPT5YHN2h5Lvi86Ex9RDboJZsPacyZko2RmkiYClUC/C+/8GZU6DGlds5Dj0JtEitwfZzxat35my9IM4tMcdHenuEdW0UIN4z4a9jBvV1UzjoruRyLcOwMRFVDP+XldFBnQ+kuc+0K9M1ZZn7ov/NnOzT9DWnxeefPNGizLmJzFVaJlXdgOh6huMNN2bnFyFdVfy6buhMlp216TzyCtXo991OtP94A67eyTV2C5EgCEi6KM359POzMeLz//R2ZvfWG/bpEhkAV2cF2tTzM/dlEFE5MGApo1ce0k52rFWjbMNVZfnQSlM53odKdXTienHyQlMTaitDaUGw7rd+R5zU3W3GZsRMTcSdnALM2iRdSdjzi604cjUDt7IL0HXG5lLvpT/ZdtHIDgC089L0Failq7HCanoiKaEnlj73gDkvgYjKiQOvZRTi64Hx8VG4cjsHbcL8ynzd8Nj6WJx0xej5S7fkx8/JttwrUMPVWWGwYSEgnS8yqINlkxgKMD2EqDtB3NieQ/YqJ78It+9nEX516SHs0Eujn5Jxz6z76da/PF274lG/J6XEa3GNDMoCyvjlhYgqjj0sZhjbIxIzBrQyax7Ch49HY++7PfFWL/lt549ey2RuFhuXea8QsdMTMWBOkuwePj8lXRafD42pb9G2qA1TrogmPtoM7q7OeKBB8dYXJXMxHMXmM9rtB/SDFTml/b+aubF4Y8rwWtL8Kd4yE2g/H9gKLz5Y/tVdRFRxDFgsTKFQINDHHXHNgozWOZuu3QE6v0iNA1fuoMjUJxNVqeX7riIjtxCHr2bgpt4+M0VqDSbrTNwM9fOwaFvkAqYSz3aJAADsvlicxPDPg461Am3Bjktm1S80MaFdEARsPFUcAF2+LZ3vInfdgHZ1KzxhmogqhgFLFWkc6I0QX/nlrLpDCu/+dRxPzk7CN5uMZ/qlqvFN4jmET1iFT1afEstuZ0s3t5u58Zz43ELzbCX051GUeO+RpuLzdvVtd7fzirhgZPmxMXkm8rPc1Zm/op/07cT1TMnxyM7hZv1cIrIMBixVxMlJgcWjOsqeK9JoP4RKvhXP2nROti5Vjez8Iny54axBue5OzGqNgG83awPL7W/3sHi76vjL9+AMaFdXfF7S09Ih3B+frTuN8AmrcN3M+R22RhAEo3NyRnWOkC03NelY9/1Y+3pXyTn9FPsfPBYNIrI+BixVSD+RXEkyunsFhn9YOa2l6h24chdv/3EEt7Pz0fwD+fxAt3O0Q0L68yiCK5AQrqy6N66N9/s2xW/PS1emKHUy4Jbs+pxboMZ3my8AADpN34S7OdLeIXvS5+vtsuXzhrXHe32byu6Enl9ofFh102ntfBj9nbPHdOdcFSJbxIClCin10qqXBDC5Drj81B49OXsXlu+/hnYfbzRa57Vlh3H+RhYAYNcFbcAyY0BL2RVElU2hUGB01waI0UtMp7vcuWRTvpKl1iV0c8XYm9NpWbLlcc2C4OykkJ13opuy/2x6FmasPS320sj1npVoHORdwdYSkSUwYKlCuj0sPz/bUUyznm8kYHHEPBr2yNlJATedhG29Z27H9DWnxT18+rcOxVPtLbucWY6/pzb1v26w5Okmnw337A35D31bV95tBvJ0eljiZ27D91suiJltmwQXByVPP1Cv4g0koirBgKUK6e4D4+fhJgYwJX9Y9YeG5pu5KoLKZ/OZGwifsEr2XHx0EC5MfQQjO4WLZUUaAXO2XhCPC620b89vL8TKltfXW6Zb4mZWvmy5rfs6UTufK66p8dV2+nQD/pIh1j8PXoNaI4g9NgPayQeaw2OLl6c/aiSJHBFVPQYsVUh3v5F6tTxxN7d4TsHrvx1GXqEaTSetldS/cDMbZFnHrmVi5MJ9suceaxWKL59qDQBwkZkjUaJ1XT8LtKx0jYO8sWR0DFa92kVSrrvpoq7S0vrbqiV7tIkX5w1vj44RxXlmaums7tEN0hoEFM8NyzMyh2XKf9pl6C2M7AmW8EhTzB/eHjMGtCx/w4moUjHTbRXb+15P5Bdq4OvhKplj8Nu+q4aVOfHW4o6mZBg993H/5vC6n/XUy0j20wa1vfBMrGWTxZnSKTKg9Er3mbss2Bao8gpx6/5S8pLl2lOfaIFvN53DyzrbZfz3ShdM+PMY+repg2/vr7AzNqS6aNdl8bmzkdwq7q7O6GlGbw4RWZ59fuWyY4He7giradhlf1tmBYebnX4jtkUFRRpcLuM2CC91b4i5z7SDr4e2p2KQkTkqz3VtIBnqo8p1/Jo2J8qkR5sBACIDa2Dm4DaIDNROjvV2d8V3Q9vi4WZB4o7WpvKwEJH9YQ+LjZi3/aJBWRZ3cq40jd9fAwB495Em2HT6Bt7p3QRt6vnjVlZxoNigtheCvN3RKzoII2XyetSqocTsoW0x5teDknJj80WsrXNkLew8f9ug/OR1FZqF+lihReUz4a9jAIqD91Zl3MNLO5m9eEjIWIr+F7tx+TKRPeFXeCvqEK7NSJqrM+G2R1RtAMCqo6li9zaV3/kb2rlAU1efxu6Ld/DE97sAAF/d30/moahALH3+AdlgpUSjoBoGZbF6y4ttUZCPdoO+iX8ft2JLzJd8pzhtfqB32TcZdL/fM1nSw5JfJD+XRW7PICKyXQxYrOibIW1ky3VXLny+3ni+CComCAJ+3XMF+y7fkT0f9+VW2fLJ/54Qn6eXYQWN/tDPAw1qQlEV+fgraNKj2kytB67ctWJLyu6GKk+S5+aTJ1qU+Vp3vdV3xgIWHwYsRHaFAYsVhfh6iDvr6qqt921y8NwkHEy2jw+aqiIIAgrv76tz+GoG3ltxHAPnJBnUO52mMigrsXDnZfH5O72jSv2ZATW0/y573u2JRSPlt1qwBQpoA6m+NrI0t6BII9nawJRxy4/gfz/uEY9bl3E4CNCuhvp9/1Us33fVaJ4jU5sjEpHtYcBiZUEy6dx9PKTf/HZfvINn5u0xqFedRSSsRqP31iA18x7+OpgiluvmsjmbnoXeM+VTuuvT3z9GjrurM7aN74Htb/dAkI+7TU+2HdyxuJfOVjZCvJNTgMbvr0HbKRuQmln6vkb62x7oToAuTcm/y+m0LLz951F0nJooW89eepuIqBgDFisrSaleMrLQr3Uo6smsIsqR2W+oulp9LFV8/uE/J/Dzbm2eDt0Pw28SpfN/3undRPZ+T7WvK1sup14tT9lVXramb4sQ/PtyF/z8rG30Ak1cqZ07o7uPT1k0lpk7ZIq7a9n+rL0VX3qvGhHZDgYsVuZ6v/u6ZCFDsK87PN1ccHDiw1ZslW17SWelzroT6ZJzx1KKl8FuPn0D/x1NlZwb3qk+4qOD8GDj2pLysDL0rtgbhUKBFnV94elmG/M0VukEmTvO3TJR01ANIzlwjClLz9ffYzsj4n6COSKyDwxYrMxNb8M85f1j3X1iSmislALelpRsXmfMtrO3UKjWYOQiw+y1nm4u+OGZ9vh8oDR7aWn3pMq15niayfP6ewcdTM4w6/7GAhbdnku5oVgism0MWKxMPzlcybHc6pNZm85XSZtsVXZ+EVpNXm+yzp8Hr6HXV9sMymvqpHEP8JJOaq5ucaCxvCSWYm6g/f7KYxX6eS5Gstfq5nFhUkYi+2Mb/cXVmH4Pi6nu7K82nsVrcY0s3SSbNXX1qTLVu6ST0XbVq13g5uyERkHarKhOeh9oz3Y1nnvFEWkEwMTWSJXuTq50ZVDben5G6wqCgKV7ZbapqAQ1lNr/W/a6rxJRdcb/tVamvyJId0nz/7WpY1Bf99vx7ex8rD6WivG/H8GXG7T5WtadSMPWszct0FrryjeymR0AzH2mnWx5eC0vSbBSQndeRB0/j4o3zsa111ktVLIcvKrM0pv8bCooP5WaZVD23ytdZGoaZ6xDR/ffnD0sRPaH/2utLNhX+mGpO3TxxVOtcObj3pIgRjcJ1v9+3IOXfj2I3w9cwzeJ56DKK8TPu6/ghZ8PYPiCvUY3f7NXO3WWuh75oJe4ay8AxETIZ5z1dJP/cPzvlS4YHx+F45PjK7eRNmrusPbic/05Ipa2OOmK5HjXhdvINrLtRE6BtvzvsZ1xadojaG5kR2VjBCO7hupuYGls2IiIbBcDFisL1pv8p/sBq1AooHRxRqu62j/YOTp/6M+kS7+NpmbkSZaPljVJV1U6cOUu3lxe3CNkbkBV0hv1Tu8m8PVwxXuPNEXXRgH49+Uu8PV0xf734wyuMZaJNjzAC2N7RJq9AsVe6b7OnCrcoypXJwDRTYWvuzRd19K9yQCKJ523CvMrVyZhY1N0mgRre9rsIUMxEUkxYLEyP73VQB6uhh+gXwxsLT4v2XPo/A3DrvNrd3MlxxduZhvUsaZf91zBk7N34c+DxT1Cry49hL2X5NPpyylZzdMlMgBA8STKn5+NQYv7AZ1uJlqS0u1RWHEoxUTNilNrBOy+eBuFag2u3dXmxfn2f23F52mZebLXFtzvQZTLRVRWxiYVRwTUwJLnYsweYiIi28CAxcpc9Sbdyg1h+Hq6isuc793vlYj70nAlzPEUlWQ59DPz91ZmUyvkdnY+3lsh3Xhv/cl0PPVDEo7fz52y++JtTPnvJK5nyGdCzc4r/rauP+/HmMQ3u1WgxY5Fd6Kxsb11KkvDd1dj8NzdaPTeGpxKLd4aoW09P3FTQgAokplHczMrX8yd0y0qsNw/P9dIkkU3Fyd0ahhg9hATEdkGBixWpj+WbmzOxd3c4t4F3RUw+r7aeNZmu7r1U63renTWDvyUdBmD5+7G/B2X0Gn6Jsn5tcfT0P+7nWK2X2/30tO0D+kYhoa1zcuQWl1U5YTTc+nFvXxRwT5w0VmaVCQzj2bovN3i8/yi8s+/MrbNAlcGEdm36jGAb8P0Pzy8SplT8cX6M4iPDjZ6Xn/eSk5+Uan3rAq7zt82eX7S3ydky+/mFODFXw6IxwqF6V12z3zcGxm5hUwMZkKUzKqpyqI/HHP5dnGA3SDAS7J6R27i79l07RDmmG4Ny92GAe3q4t0VhrlcuDKIyL7xf7CV6Q8JlRZc9GpmPFiRE/3BOrPbZAm/7Tcvt0ZWXiEEQUCbKRsk5SE+7nBxNv5rq3RxZrBiREn+kwILLmvW3QHbzdkJqffnqtTx90CRzu7IJQGLWiMgfMIqhE9YJbmPn6cbysvNxUk27T57WIjsm1n/g8PDw6FQKAweY8eONXrN77//jiZNmsDd3R0tWrTA6tWrJecFQcCkSZMQEhICDw8PxMXF4dy5c0bu5nhcypjBa0jHegCKA5wjVzMs2KLKd09vTsH/YuqVek1Kxj3ZlOy6S7zJPCU9DAUWnMMya5P2/26hRoN0VXHAEuithEan92XBzkuInZaImRvPGtzjhW4NKtyOn0Z1xMs9IjFzUGuxjD0sRPbNrP/B+/btQ2pqqvjYsKH42+/AgQNl6+/atQtDhgzBs88+i0OHDqF///7o378/jh/XTr6cMWMGvvnmG8yZMwd79uyBl5cX4uPjkZcnv4rA0bg6le2foGRuy71CtWRlTb/WoZj6RAuLtK2ypGRoVy9dnPoIHmwUUOo1E/48hu3nDJPf+Vbgm3d1V9KbZ4nEcbkFRfjfj7vFuVZA8fLiklVCtb2VqF9LO7dEIwCpmXmy200MbFf23bONCavpibfioyQrx/SzShORfTHrf3Dt2rURHBwsPv777z80bNgQ3brJr8b4+uuv0bt3b4wfPx5NmzbFlClT0LZtW3z77bcAintXZs6ciffffx/9+vVDy5Yt8dNPP+H69etYuXJlhV+cPShrN3XJh02RWgN3V+01Mwe1LrXHYvPpG1We3VRXydyEZiE+cHJSGJ00q1BAfG2Hr2Zg5kbDnraMXNvLLWMvnO9P8LZE4rghP+7BrgvG5ynV9lairr8nHmlR+pBmiG/lZR7WTSJnqxPSiahsyv2Vo6CgAL/88gtGjRpl9A9BUlIS4uKkybzi4+ORlJQEALh06RLS0tIkdXx9fRETEyPWcXROTgosGFGchbRBbePb3bvdHzq6V6jGxPsTVFvU8TX63ut2hY9ctM+qGyeWLFMueX2xDbRZab115uz4ergioU9Tk/c6ei3TAi2sHpzv/65oLLD5YWnDlJ5uxf/OXSJrl3ovYyvlyqNlXT8AHEokcgTlDlhWrlyJjIwMjBgxwmidtLQ0BAUFScqCgoKQlpYmni8pM1ZHTn5+PlQqleRhzx5qEoTNb3XHqle6Gq1T0sPy655ksSxFJ19JV71hlkdbhkiOv0m03ryg8zeKe1hCfIsnw+rmBPHx0Pa2ZOQWYoDMcMDmt7qLz8NqOv6+P5ZS8r5XZgfL34dTDCbMmuLhVvqfnMrsCfH1cMXRD3th+9s9Ku2eRGQd5Q5Y5s+fjz59+iA0NLQy21Mm06ZNg6+vr/gICwur8jZUtogAL3iY+GbpKjN05KGziZx+MixTK2mq2vHrxb0ibetpN+B7oVsDKBTAlP7Rkrpyq6QiArywcEQHRIf64Ien2xucp7IpiRMrc0jotWWHDcqG6g1RTn5c+2/sYWLjQ0Ab1FYmH3dXkxsuEpF9KNen2pUrV7Bx40aMHj3aZL3g4GCkp6dLytLT0xEcHCyeLykzVkdOQkICMjMzxcfVq5bZjt6W6C9/BqR7szQN8RGf605utAVpmfkAgHo67ZrQuwkOvv8wHmoShNFdIgBIP9hK9G1R3FPUo0kgVr3aFc1CfQzqUNk4O1X+kFDrMD/J8ZrXuuLj/s0lZbrLzJUupgOHGJ0NLYmIdJUrYFm4cCECAwPRt29fk/ViY2ORmJgoKduwYQNiY2MBABEREQgODpbUUalU2LNnj1hHjlKphI+Pj+Th6Nxklj8rdb41PtYyBBP6NMH/tamDpc89UJVNMymvUI1b2cUBi+4Hl0KhgP/9nanff7QZLk/vi+GdwgEAf7yo/bc3tvMumc9JUbmTbnMLigwSxTUN8YFCocCTbbVDe7obL5YWLJUW0BBR9WV2ClSNRoOFCxdi+PDhcHGRXj5s2DDUqVMH06ZNAwC89tpr6NatG7744gv07dsXy5Ytw/79+zF37lwAxR9ar7/+Oj7++GM0atQIERERmDhxIkJDQ9G/f/+KvzoHIjfEoxvDKBQKvFhKdtB0VV6VJ1Vrq5P4rWYZlyQ31tlVd8PJdBM1yRxOisqbw6LRCGg2yXhSwrimgfjz4DUA0gmvpQVLSlfbGcokItti9l+HjRs3Ijk5GaNGjTI4l5ycjNRU7bbxnTp1wpIlSzB37ly0atUKf/zxB1auXInmzbVdxm+//TZeeeUVPP/88+jQoQOys7Oxdu1auLszW6kuuSGhVnrd8foaB0n30tl+zvh+PpaQmnlPshGdk1PZJlN66vQcPdio9FUlVDbikFAlRCxX7uQalOnOlb2RlS8+b6iz+q20DLZyv+dEREA5elh69epldPv2LVu2GJQNHDjQaGI5oLhn4KOPPsJHH31kblOqFVe9IaHBHcLwVq8ok9d0jKgp2Z+lXs2qm9uy4tA1vPHbEfHYnMmULs5OeKtXY/yUdAUzBrS0RPOqJXFIqBLmsKjuFRqU6d62Q3jxXJRaXm6S3sEO4f76l+H9vk3x8apTAICDyXcr3DYickzW3xWPykQ/S+f0J0v/IG+st8ldRXbANce59CxJsAIAy18wPidJzssPNcLLDzWqzGZVeyW/QpUx6fZ2Tr7J881CffDfK10Q6iddhq5QKDDt/1og4S/t5oS+OkvbD8lsx0BEBHDzQ7uh21XeVy/HijH/61gPr/XUfuj/ffh6pbdLzsNfbTMoC6vC3h2SJ85hqYQhoVOpWaXWaV7HFzW9DIeABncIw/dD24rHJUnlAODVhyIr3DYickwMWOyEbh4WrzJmAnVxdsIbDzcWj/84cM3ocF5l0GgEjF6832L3p4pxElPzV/xeuQVFBmVNgr1lahpSKBSS+Veebs449VFvLH3uAbzak71qRCSPAYud0J3DUlryLX26HyQ7zltu4u2+y3ew8ZThqp6ooLJ9kJFlOVfSHJaUjHv4bvMFSdkL3Rrgl9ExZb6H7h5aSlcneLg5I7ZhLZtKeEhEtoV/HeyE7hwWDzfzph7p5sTIlJksWRkEQcCgubsNyp/tEoHvn24rcwVVtZJFWhXtZfu/73calCX0aSrZGbk0broBSxk3ACWi6o2Tbu2E7hwWc3tY/Dy1kxotsVMvAOQVGo4zzH2mHXpFl747L1UNp0rarTnE1wPpKu2k2/Hxpleryanh5oKIAC/czS0wmBxORCSHAYudkAQsZdhATpfuBoNF6soJWFR5hVi2Nxn9WtdBkI+7bM8NgxXbUllDQrpbQnw9uDX6ta5j9j2cnBT45+XO0GgAb3fX0i8gomqPfbF2ws1FO4fFyczdbL11UqNXVg/L0j3JmLr6NGKmJmLHuVu4k1MgOe9SxiRxVHVKelj+PJBS7nucvK6SJCD0MnN4Upe3uyt8PRmsEFHZsIfFTuj2sNSqUbYU9yW8zNjLpawW77osPn96/h5J+nUAKLLQ0BOV3+GrGQAg7u1UHo98s11yLLeLOBGRJfCvjZ3QDVhq1zBv2wLdHY6vZ+ZVSnv073MzS/ohWNYlrlR15LLTVhR70oioqjBgsRO6AYu5G8S5Ojsh4H6vzDeJ5yrclnsFpjPmPt4qFH+M6VThn0OVy1knuJjy38lyrRaKDJTuT9W6lP2siIgqC4eE7ITusubyfKu9la2dY6LWCJIPL3PdzS0weu6/V7qgeR3fct+bLEf333z+jkuIaxqE2Ia1zLpHXqE2WP3wsWaS4UYiIktiD4udcNWZdFuRYAMAVh1LLb2SCaZyuTBYsV36vzf6E6VLk1eoxrW79wAUT+Qe1KFepbWNiKg0DFjshIuT9p/K3FVC+m6oKjaPxdhcCE5nsG3Oer835i5v1t1JOfGtbvAo4xYRRESVgQGLndBNze9TjrwVjYNqlF6pjFR5xfvIKBTAkI5hYrluUEW2x0kvojQ3vtSduxTobd7EbyKiiuInjJ1QKBT4uH9zjI+PQr1a5u98PGNAK/F5RVc2l/SwdIkMwLT/a4kJfZoAKE4iRrZLf+6Tm5lLknPuByyxDcyb90JEVBk4Y86OPP1A/XJfG+qn/UZcoNZg1KJ9aFXXD6/Fmb87bvKdXADaLQJeeLABnmofhppe5uWHoaqlP4fF3IDl3v0dmj05FEREVsAelmpCd5XR5tM3sOn0DXy18Wy57rX2eBoAIPV+LhaFQsFgxQ7oz31Sm7lNQ05+cQ+LJ1cGEZEV8C9PNaH7bVr3Yyonv8ispal3cwpwJj0LAHNw2Bv9HpayZCMuKNKg8ftrAAD/iyleFeTFHhYisgL2sFQTnjp7vuguZ9XNUHsjKw+T/z2B8zeyjd5nxKJ94vM29fwqt5FkUQY9LGUIWBbsvCQ+X7InGQC4OoiIrIIBSzUyoF1dANIPqlydlR8dP0nEwp2X8dQPSUbvkXI3V3zOpGH2RX/SbZFGgwNX7qLnF1uw7exN2Wsu38oxKKvIhodEROXFgKUaUd4fFtLd/C6/yDDN/p2cAvSeuQ3zd1wyOKebMdebAYtd0R8SOpScgeEL9uLCzRwMW7BX9hq5YT9PJXtYiKjqMWCpRpQuxR80ur0qeYUa2bqn07Iw5b+T0JgYNuDQgH3Rz8OyaNdlZOcXmbxG7jx7WIjIGhiwVCNymyaeSVMBAH7ZfUX2mgydrLa6+8gAgLsrAxZ7Ym4m4sW7LuPjVacMyhmoEpE1MGCpRnSXNpf48N+TAID3Vx6XvUZ3o8MRC7XDBn2aB6NxkHclt5AsST9eeaRFsMn6H/xzQracPSxEZA38y1ONlKdHRDcd++6Ld8Tns59uVyltoqqj0FslVGRmHpYSnMNCRNbAHpZqJCJAPqV/kVp+Hgsgne9CjmX9yfRyXcceFiKyBgYs1YiPh/ymiSkZ94xek1OgnXTZObJ4D5mS5dFkXyprM22m5icia2DAUo2UrBLS1+2zLUavyc3X9rCUJB7j5nf2qWFg5ezYzYCFiKyBAUs14ups/nfsnIIi3M0pwIDZu7D93C0ATBhnr57tElEp9+G/PxFZAwOWakQ/cVhZ3CtQ4+vEc9h/5a5YVtffozKbRVXE3dUZ9WrKz2OSU0MnMJnSL1p8zmXNRGQNDFiqERen0v+5vxnSRnKcU1CEmzqZcQEgIsCrUttFtkejEcSkcdvG95BsmOnJ/DtEZAVmBywpKSl4+umnUatWLXh4eKBFixbYv3+/0fojRoyAQqEweERHa7+xffjhhwbnmzRpUr5XREbpp2FpX99fcrz3vZ54vFUo4qODxLLcfDV83KVDAB78wLJbpnrZNp1OR0FR8Yqx65naidjubk545oH6qOPngYHt6sJFJp8PEZGlmfWX5+7du+jcuTNcXV2xZs0anDx5El988QX8/f2NXvP1118jNTVVfFy9ehU1a9bEwIEDJfWio6Ml9Xbs2FG+V0RGOev1sHjrBSI+7sWriH54pj1e6NYAQPGyZv1dfvVTvJP9MPVPN2rRfnyTeA4A8Ngs7f+/QG93KBQK7JzwED4b2MrSTSQikmXW7LlPP/0UYWFhWLhwoVgWEWF6Ip+vry98fX3F45UrV+Lu3bsYOXKktCEuLggONp15kypGf7de3cmTShcnSWK5klwbC3ZewqjO2n/jgBpKC7eSLOnVno3w2rLDRs8v338Vb8VH4W5uodE6RETWYFYPyz///IP27dtj4MCBCAwMRJs2bfDjjz+a9QPnz5+PuLg41K9fX1J+7tw5hIaGokGDBhg6dCiSk5PNui+VTr9nRBKg6K38+PPgNfF5Vp72w+vJtnUs1DqqCv1a18GCEe2NntcIMLnhJRGRtZgVsFy8eBGzZ89Go0aNsG7dOowZMwavvvoqFi9eXKbrr1+/jjVr1mD06NGS8piYGCxatAhr167F7NmzcenSJXTt2hVZWVmy98nPz4dKpZI8qHT6PSzuOpsh6p/TzXB7I0s76bZerbKvMiHbVK+mqUnTAo6lZFZZW4iIysqsgEWj0aBt27aYOnUq2rRpg+effx7PPfcc5syZU6brFy9eDD8/P/Tv319S3qdPHwwcOBAtW7ZEfHw8Vq9ejYyMDCxfvlz2PtOmTROHmnx9fREWFmbOy6i29Cdc3s7WbmyoH7DMfUa7V9DWszcBFC9zfao932t7ZyofjyDAYM4SEZEtMCtgCQkJQbNmzSRlTZs2LdPwjSAIWLBgAZ555hm4ubmZrOvn54fGjRvj/PnzsucTEhKQmZkpPq5evVr2F1GNOet8EA3pGIYd9xPBAYCz3odYm3qGE6nf7h0FV64QsXumVgoJpZwnIrIWsz59OnfujDNnzkjKzp49azAfRc7WrVtx/vx5PPvss6XWzc7OxoULFxASEiJ7XqlUwsfHR/Kg0ukGJZ0aBkj2FtJNwW+M0oXBiiMwlY8nr1CNHedvisfcN4qIbIVZn0BvvPEGdu/ejalTp+L8+fNYsmQJ5s6di7Fjx4p1EhISMGzYMINr58+fj5iYGDRv3tzg3FtvvYWtW7fi8uXL2LVrF5544gk4OztjyJAh5XhJZIzusI+TQiEJQG7nFMhdIuHGgMUhuJgYEsotUGPq6tPi8QePNTNal4ioKpn1CdShQwesWLECS5cuRfPmzTFlyhTMnDkTQ4cOFeukpqYaDBFlZmbizz//NNq7cu3aNQwZMgRRUVF46qmnUKtWLezevRu1a9cux0siY3TnJjg7mZ9Pxc2ZCeMcgW7g+n+lrPrydpff4ZuIqKqZvYvZo48+ikcffdTo+UWLFhmU+fr6Ijc31+g1y5YtM7cZVA66H1QKhQJPtKmDz9adMVo/KsgbZ9K1K7XYw+IYdOeoPNCgFv46mGLF1hARlQ0/gaoRZ70hof5tTH+71u+B4RwWx6A7cVoQBIyPjxKP6/hxY0sisk38BKpGFHpDQkHeprPWPtwsSHJsau4D2Q/dwFUQgLE9IvHnmFgAQEqGdg+huKZBBtcSEVkLA5ZqKsTXo9RN7Mb2aIgp/bSbVGo0lm4VVQXdocGSnLZyy9Vre5tOP0BEVJUYsFQzXw1qhfceaYqmIaUvBVe6OOPpB7RL1tUCU7Y7At2etpJ/Uvn5SexRIyLbYfakW7JvT7SR5tVwc3ZCgdp414mCWU8dmnC/j0Wuh4X544jIlrCHpZozlaa9xOAOYWhZ1xexDWpVQYuoKok9LDIBS7NQJmQkItvBHpZqztXFCSgwneV2+pMtq6g1VNWE+xGL3JDQ4A71qro5RERGsYelmut1fyVQQA3TK4bIsckNCXFPISKyJexhqebe7t0EQT7u3IW5mtKuEmJwQkS2jQFLNRdQQ4k3e0WVXpEcUskcFqULt10gItvGISGiaiz6/sRa/R4WDhESka1hDwtRNbRxXDdcvpWD9uE1ARguXx/Uoa7cZUREVsOAhagaigysgcjAGkbPc4iIiGwNh4SIyICPO7/LEJFtYcBCRAY83NjDQkS2hQELEQEABrbTzltxceKfBiKyLfyrREQAINm921V2M0QiIuvhXyUiAgC46GS2dWWWWyKyMQxYiAiANBW/XKp+IiJr4l8lIgIgTR7nwlT9RGRjGLAQEQDAWWeirRt7WIjIxvCvEhEBkM5hcWHAQkQ2hn+ViAiAdBiIuzcTka1hwEJEAPRWCbGHhYhsDP8qEREAwIkBCxHZMP5VIiIAwPazt8TnHBIiIlvDgIWIAABJF2+Lz9nDQkS2hn+ViAgA0KC2l/icAQsR2Rr+VSIiAMCkR5uJz5k4johsDQMWIgIA+Hq4is/Zw0JEtoZ/lYgIAOCik+mWk26JyNYwYCEiANz8kIhsm9l/lVJSUvD000+jVq1a8PDwQIsWLbB//36j9bds2QKFQmHwSEtLk9T77rvvEB4eDnd3d8TExGDv3r3mvxoiKrcAbzfxuW4SOSIiW+BiTuW7d++ic+fO6NGjB9asWYPatWvj3Llz8Pf3L/XaM2fOwMfHRzwODAwUn//2228YN24c5syZg5iYGMycORPx8fE4c+aMpB4RWU6gtzt+HNYenm7OUCgYsBCRbVEIgiCUtfKECROwc+dObN++vcw/YMuWLejRowfu3r0LPz8/2ToxMTHo0KEDvv32WwCARqNBWFgYXnnlFUyYMKHUn6FSqeDr64vMzExJUERERES2y5zPb7OGhP755x+0b98eAwcORGBgINq0aYMff/yxTNe2bt0aISEhePjhh7Fz506xvKCgAAcOHEBcXJy2UU5OiIuLQ1JSkjnNIyIiIgdlVsBy8eJFzJ49G40aNcK6deswZswYvPrqq1i8eLHRa0JCQjBnzhz8+eef+PPPPxEWFobu3bvj4MGDAIBbt25BrVYjKChIcl1QUJDBPJcS+fn5UKlUkgcRERE5LrPmsGg0GrRv3x5Tp04FALRp0wbHjx/HnDlzMHz4cNlroqKiEBUVJR536tQJFy5cwFdffYWff/65XI2eNm0aJk+eXK5riYiIyP6Y1cMSEhKCZs2aScqaNm2K5ORks35ox44dcf78eQBAQEAAnJ2dkZ6eLqmTnp6O4OBg2esTEhKQmZkpPq5evWrWzyciIiL7YlbA0rlzZ5w5c0ZSdvbsWdSvX9+sH3r48GGEhIQAANzc3NCuXTskJiaK5zUaDRITExEbGyt7vVKphI+Pj+RBREREjsusIaE33ngDnTp1wtSpU/HUU09h7969mDt3LubOnSvWSUhIQEpKCn766ScAwMyZMxEREYHo6Gjk5eVh3rx52LRpE9avXy9eM27cOAwfPhzt27dHx44dMXPmTOTk5GDkyJGV9DKJiIjInpkVsHTo0AErVqxAQkICPvroI0RERGDmzJkYOnSoWCc1NVUyRFRQUIA333wTKSkp8PT0RMuWLbFx40b06NFDrDNo0CDcvHkTkyZNQlpaGlq3bo21a9caTMQlIiKi6smsPCy2inlYiIiI7I/F8rAQERERWQMDFiIiIrJ5DFiIiIjI5jFgISIiIpvHgIWIiIhsnlnLmm1VyUIn7ilERERkP0o+t8uyYNkhApasrCwAQFhYmJVbQkRERObKysqCr6+vyToOkYdFo9Hg+vXr8Pb2hkKhqNR7q1QqhIWF4erVq8zxUk58DysH38eK43tYOfg+Vhzfw2KCICArKwuhoaFwcjI9S8UhelicnJxQt25di/4M7llUcXwPKwffx4rje1g5+D5WHN9DlNqzUoKTbomIiMjmMWAhIiIim8eApRRKpRIffPABlEqltZtit/geVg6+jxXH97By8H2sOL6H5nOISbdERETk2NjDQkRERDaPAQsRERHZPAYsREREZPMYsBAREZHNY8BSiu+++w7h4eFwd3dHTEwM9u7da+0mWcW2bdvw2GOPITQ0FAqFAitXrpScFwQBkyZNQkhICDw8PBAXF4dz585J6ty5cwdDhw6Fj48P/Pz88OyzzyI7O1tS5+jRo+jatSvc3d0RFhaGGTNmWPqlVZlp06ahQ4cO8Pb2RmBgIPr3748zZ85I6uTl5WHs2LGoVasWatSogSeffBLp6emSOsnJyejbty88PT0RGBiI8ePHo6ioSFJny5YtaNu2LZRKJSIjI7Fo0SJLv7wqM3v2bLRs2VJMuBUbG4s1a9aI5/kemm/69OlQKBR4/fXXxTK+j6X78MMPoVAoJI8mTZqI5/keVjKBjFq2bJng5uYmLFiwQDhx4oTw3HPPCX5+fkJ6erq1m1blVq9eLbz33nvCX3/9JQAQVqxYITk/ffp0wdfXV1i5cqVw5MgR4fHHHxciIiKEe/fuiXV69+4ttGrVSti9e7ewfft2ITIyUhgyZIh4PjMzUwgKChKGDh0qHD9+XFi6dKng4eEh/PDDD1X1Mi0qPj5eWLhwoXD8+HHh8OHDwiOPPCLUq1dPyM7OFuu8+OKLQlhYmJCYmCjs379feOCBB4ROnTqJ54uKioTmzZsLcXFxwqFDh4TVq1cLAQEBQkJCgljn4sWLgqenpzBu3Djh5MmTwqxZswRnZ2dh7dq1Vfp6LeWff/4RVq1aJZw9e1Y4c+aM8O677wqurq7C8ePHBUHge2iuvXv3CuHh4ULLli2F1157TSzn+1i6Dz74QIiOjhZSU1PFx82bN8XzfA8rFwMWEzp27CiMHTtWPFar1UJoaKgwbdo0K7bK+vQDFo1GIwQHBwufffaZWJaRkSEolUph6dKlgiAIwsmTJwUAwr59+8Q6a9asERQKhZCSkiIIgiB8//33gr+/v5Cfny/Weeedd4SoqCgLvyLruHHjhgBA2Lp1qyAIxe+Zq6ur8Pvvv4t1Tp06JQAQkpKSBEEoDhydnJyEtLQ0sc7s2bMFHx8f8X17++23hejoaMnPGjRokBAfH2/pl2Q1/v7+wrx58/gemikrK0to1KiRsGHDBqFbt25iwML3sWw++OADoVWrVrLn+B5WPg4JGVFQUIADBw4gLi5OLHNyckJcXBySkpKs2DLbc+nSJaSlpUneK19fX8TExIjvVVJSEvz8/NC+fXuxTlxcHJycnLBnzx6xzoMPPgg3NzexTnx8PM6cOYO7d+9W0aupOpmZmQCAmjVrAgAOHDiAwsJCyfvYpEkT1KtXT/I+tmjRAkFBQWKd+Ph4qFQqnDhxQqyje4+SOo74e6tWq7Fs2TLk5OQgNjaW76GZxo4di759+xq8Vr6PZXfu3DmEhoaiQYMGGDp0KJKTkwHwPbQEBixG3Lp1C2q1WvKLBABBQUFIS0uzUqtsU8n7Yeq9SktLQ2BgoOS8i4sLatasKakjdw/dn+EoNBoNXn/9dXTu3BnNmzcHUPwa3dzc4OfnJ6mr/z6W9h4Zq6NSqXDv3j1LvJwqd+zYMdSoUQNKpRIvvvgiVqxYgWbNmvE9NMOyZctw8OBBTJs2zeAc38eyiYmJwaJFi7B27VrMnj0bly5dQteuXZGVlcX30AIcYrdmInszduxYHD9+HDt27LB2U+xSVFQUDh8+jMzMTPzxxx8YPnw4tm7dau1m2Y2rV6/itddew4YNG+Du7m7t5titPn36iM9btmyJmJgY1K9fH8uXL4eHh4cVW+aY2MNiREBAAJydnQ1mdKenpyM4ONhKrbJNJe+HqfcqODgYN27ckJwvKirCnTt3JHXk7qH7MxzByy+/jP/++w+bN29G3bp1xfLg4GAUFBQgIyNDUl//fSztPTJWx8fHx2H+iLq5uSEyMhLt2rXDtGnT0KpVK3z99dd8D8vowIEDuHHjBtq2bQsXFxe4uLhg69at+Oabb+Di4oKgoCC+j+Xg5+eHxo0b4/z58/xdtAAGLEa4ubmhXbt2SExMFMs0Gg0SExMRGxtrxZbZnoiICAQHB0veK5VKhT179ojvVWxsLDIyMnDgwAGxzqZNm6DRaBATEyPW2bZtGwoLC8U6GzZsQFRUFPz9/avo1ViOIAh4+eWXsWLFCmzatAkRERGS8+3atYOrq6vkfTxz5gySk5Ml7+OxY8ckwd+GDRvg4+ODZs2aiXV071FSx5F/bzUaDfLz8/kellHPnj1x7NgxHD58WHy0b98eQ4cOFZ/zfTRfdnY2Lly4gJCQEP4uWoK1Z/3asmXLlglKpVJYtGiRcPLkSeH5558X/Pz8JDO6q4usrCzh0KFDwqFDhwQAwpdffikcOnRIuHLliiAIxcua/fz8hL///ls4evSo0K9fP9llzW3atBH27Nkj7NixQ2jUqJFkWXNGRoYQFBQkPPPMM8Lx48eFZcuWCZ6eng6zrHnMmDGCr6+vsGXLFskyyNzcXLHOiy++KNSrV0/YtGmTsH//fiE2NlaIjY0Vz5csg+zVq5dw+PBhYe3atULt2rVll0GOHz9eOHXqlPDdd9851DLICRMmCFu3bhUuXbokHD16VJgwYYKgUCiE9evXC4LA97C8dFcJCQLfx7J48803hS1btgiXLl0Sdu7cKcTFxQkBAQHCjRs3BEHge1jZGLCUYtasWUK9evUENzc3oWPHjsLu3but3SSr2Lx5swDA4DF8+HBBEIqXNk+cOFEICgoSlEql0LNnT+HMmTOSe9y+fVsYMmSIUKNGDcHHx0cYOXKkkJWVJalz5MgRoUuXLoJSqRTq1KkjTJ8+vapeosXJvX8AhIULF4p17t27J7z00kuCv7+/4OnpKTzxxBNCamqq5D6XL18W+vTpI3h4eAgBAQHCm2++KRQWFkrqbN68WWjdurXg5uYmNGjQQPIz7N2oUaOE+vXrC25ubkLt2rWFnj17isGKIPA9LC/9gIXvY+kGDRokhISECG5ubkKdOnWEQYMGCefPnxfP8z2sXApBEATr9O0QERERlQ3nsBAREZHNY8BCRERENo8BCxEREdk8BixERERk8xiwEBERkc1jwEJEREQ2jwELERER2TwGLERERGTzGLAQERGRzWPAQkRERDaPAQsRERHZPAYsREREZPP+H8sZyieONW2PAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "\n", - "**2.** How to implement its quantum counterpart?
\n", - "**3.** Conduct a resource analysis on the quantum algorithm. Which step is most costly?
\n", - "**4.** Get used to Quantum Phase Estimation, QPE.
\n", - "**5.** How is the influence of the classical noise (the random fluctuations in financial markets)?
\n", - "**6.** How is the influence of quantum noise?
" - ] - }, - { - "cell_type": "markdown", - "id": "85d339f8-0a2d-4c84-88d2-ae2f5d7e7597", - "metadata": {}, - "source": [ - "## The challenge" - ] - }, - { - "cell_type": "markdown", - "id": "ddf18fd5", - "metadata": {}, - "source": [ - "### Main idea\n", - "\n", - "Topological Data Analysis (TDA) is a robust and innovative technique that has demonstrated impressive results in detecting financial crashes—critical transitions between financial regimes—and providing early warning signals. By analyzing changes in Betti numbers, TDA can reveal shifts in underlying structures that precede these transitions. In this notebook, we employ a quantum TDA algorithm to calculate Betti numbers, enabling effective early detection of financial crashes.\n", - "\n", - "### What is topological data analysis\n", - "\n", - "Topology studies the properties of geometric objects that are preserved under continuous deformations, such as stretching, twisting, crumpling, and bending. Homology serves as a tool to study and quantify the topological properties of spaces. The homology of an object is often summarized using $k$-th Betti numbers, which count the number of $k$-dimensional holes within the object. For example, the zeroth Betti number, $\\beta_0$, counts the number of connected components; the first Betti number, $\\beta_1$, counts the number of loops; and the second Betti number, $\\beta_2$, counts the number of enclosed volumes, and so on.\n", - "\n", - "In finance, input data is typically represented as time series, which are subsequently transformed into a series of discrete points. To model the underlying structure with them, we construct a simplicial complex (specifically, a Vietoris–Rips complex), a collection of simple shapes that connect the points together. Those simple shapes are called simplex, and they are generalization of the notion of a triangle or tetrahedron to arbitrary dimensions. A $k$-simplex is a collection of $k + 1$ vertices with $k(k + 1)/2$ edges in $k$ dimensions. A resolution threshold, $\\epsilon$, is chosen so that any two points within $\\epsilon$ distance of each other are connected by a line segment in the simplicial complex. As $\\epsilon$ increases, more connections are added, and lower-order components may merge into higher-order components. This results in a decrease in the lower-order Betti numbers and an increase in the higher-order Betti numbers. The study of how topological properties change across a sequence of resolution thresholds is called persistent homology.\n", - "\n", - "
\n",
- " \n",
- "**STEPS**\n",
- " \n",
- "* Input: a time series of stock price; Output: a time-evolution of topological properties.\n",
- "\n",
- "* Preparing point cloud\n",
- " * Apply Taken's embedding.\n",
- " * Apply a sliding window for obatining a time-varying point-cloud.\n",
- "* Building Laplacian\n",
- " * Construct the Vietoris-Rips (VR) complex from the point cloud using [`GUDHI`](https://gudhi.inria.fr/).\n",
- " * Build the boudnary operator of this complex.\n",
- " * Build the Laplacian matrix based on the boundary operators, then pad it and rescale it.\n",
- "* Applying quantum phase estimation\n",
- " * Use Quantum Phase Estimation (QPE) to find the number non-zero eigenvalues of the Laplacian matrix. Round up the results and get the Betti numbers.\n",
- " * Vary the resolution threshold and obtain a series of Betti numbers, which are the Betti curves.\n",
- "* Detecting financial market crashes\n",
- " * Find the relation between Betti numbers and financial market crashes.\n",
- " \n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8a89d505",
- "metadata": {},
- "source": [
- " NOTE:\n",
- "\n",
- "In the coding process, it is recommended to start with the following initial values for the variables: \n",
- "\n",
- "* `N = 4` # dimension of vectors\n",
- "* `d = 5` # time delay\n",
- "* `w = 5` # window size\n",
- "* `epsilon = 0.1` # resolution threshold\n",
- "* `q = 3` # number of precision qubits\n",
- "\n",
- "However, you will be tasked with determining the optimal values later in this challenge.\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1656ba25",
- "metadata": {},
- "source": [
- "### Step 0: Loading data\n",
- "\n",
- "\n",
- "To assess the practical applicability of TDA with quantum techniques, we will analyze a small dataset of the S&P 500 index from the period surrounding the 2008 financial crisis. You may find the associated file: *SP500.csv*"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4ec880cd",
- "metadata": {},
- "outputs": [],
- "source": [
- "import pandas as pd\n",
- "import numpy as np\n",
- "\n",
- "time_series = np.log(pd.read_csv(\"SP500.csv\", header=None).to_numpy().squeeze())"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b4597392",
- "metadata": {},
- "source": [
- "### Step 1: Preparing point cloud\n",
- "\n",
- "**Instruction:**\n",
- "\n",
- "Consider a time series $X = \\{x_0, x_1, \\ldots, x_{L-1}\\}$ of numerical values of length $L$ as input. We choose and fix an embedded dimension $N$ and a time-delay $d \\ge 1$. Taken's embedding theorem convert this time series into a series of $N$-dimensional time-delay coordinate vectors $Z=\\{z_0, z_1, \\dots, z_{L-1-(N-1)d}\\}$:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "z_0 =& (x_0, x_d, \\ldots, x_{(N-1)d}),\\\\\n",
- "z_1 =& (x_1, x_{1+d}, \\ldots, x_{1+(N-1)d}),\\\\\n",
- "\\vdots\\\\\n",
- "z_t =& (x_t, x_{t+d}, \\ldots, x_{t+(N-1)d}),\\\\\n",
- "\\vdots&\\\\\n",
- "z_{L-1-(N-1)d} =& (x_{L-1-(N-1)d}, x_{L-1-(N-2)d}, \\ldots, x_{L-1}).\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "To detect qualitative changes along a time series, we apply a sliding window of size $w$ and assess how topological properties change along the sliding window. For a proper size, $w$ needs to fullfill $N \\ll w \\ll L$. The sliding window gives a time-varying point cloud embedded in $\\mathcal{R}^N$:\n",
- "\n",
- "$$\n",
- "Z^t = \\{z_t, z_{t+1}, \\ldots, z_{t+w-1}\\}, \\quad \\text{for } t \\in \\{0, \\ldots, K-1\\}\n",
- "$$\n",
- "\n",
- "**Action:**\n",
- "\n",
- "Following Takens' embedding theorem, transform the `time_series` into a series of $ N $-dimensional vectors. Afterward, apply a sliding window to these vectors and obtain a time-varying point cloud.\n",
- "\n",
- "**Answer:**"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "7b701be6",
- "metadata": {},
- "outputs": [],
- "source": [
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9376bd0b",
- "metadata": {},
- "source": [
- "\n",
- " \n",
- "BONUS EXERCISE: \n",
- "\n",
- "`gtda.time_series.TakensEmbedding` can conduct this transformation. Try avoid using this function and build your own embedding function.\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9780b90c",
- "metadata": {},
- "source": [
- " NOTE:\n",
- "\n",
- "From step 2, we present an example originally provided in the appendix of [Khandelwal's and Chandra's paper](https://arxiv.org/abs/2302.09553). This example can be used to verify your code.\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c945e850",
- "metadata": {},
- "source": [
- "### Step 2: Building Laplacian\n",
- "\n",
- "**Instruction:**\n",
- "\n",
- "In this step, we are going to use `GUDHI`, a Python package specialized in TDA and and higher dimensional geometry understanding. For detailed instructions on the functions we will use, please refer to [this website](https://gudhi.inria.fr/python/latest/rips_complex_user.html#).\n",
- "\n",
- "A simplicial complex is constructed on the point cloud using `gudhi.RipsComplex` class, followed by its `create_simplex_tree` method. The resolution threshold `epsilon` is set via the `max_edge_length` parameter. This process identifies the connectivity of the complex within the resolution threshold and produces a simplex tree. The simplex tree serves as a general representation of simplicial complexes. Using its `get_filtration` method, simplicies are retrieved as a collection of lists, where elements are grouped based on their connections. Each dimension up to a specified maximum is represented by its respective collection of lists.\n",
- "\n",
- "**Example:**\n",
- "\n",
- "Here is an example of a simplex tree $\\mathcal{K}$ with a maximum dimension of 2. In its zeroth dimension, each point is a connected component; In the first dimension, 6 line segments connect 6 pairs of points $[1, 2], [1, 3], [2, 3], [3, 4], [3, 5], [4, 5]$; In the second dimension, a filled trangle is formed among points $[1 ,2 ,3]$.\n",
- "\n",
- "$$\n",
- "\\mathcal{K} = [[[1], [2], [3], [4], [5]],[[1, 2], [1, 3], [2, 3], [3, 4], [3, 5], [4, 5]], [[1, 2, 3]]]\n",
- "$$\n",
- "\n",
- "**Action:**\n",
- "\n",
- "Build a simplicial complex by applying functions from `GUDHI` on the point cloud obtained in step 1, and extract its simplex tree. It is recommended to store the simplex tree in a format similar to the example provided, i.e., in the format $[S_0, S_1, S_2, \\dots]$, where $S_i$ represents the set of $i$-simplices.\n",
- "\n",
- "**Answer:**"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "026d6b61",
- "metadata": {},
- "outputs": [],
- "source": [
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "573efefa",
- "metadata": {},
- "source": [
- "**Instruction:**\n",
- "\n",
- "Let $S_k$ be the set of $k$-simplicies in the complex $\\mathcal{K}$ with individual simplicies denoted by $s_k ∈ S_k$ written as $[j_0, j_1, \\dots, j_k]$ where $j_i$ is the $i$-th vertex of $s_k$. Note that the vertices are ordered in ascending fashion in the initial point cloud, and this order is kept throughout. The restricted boundary operator $\\partial_k$ is defined on the $k$-simplicies as:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "\\partial_k s_k &= \\sum_{t=0}^{k} (-1)^t [v_0, \\dots, v_{t-1}, v_{t+1}, \\dots, v_k]\\\\\n",
- "&= \\sum_{t=0}^{k} (-1)^t s_{k-1} (t)\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "where $s_{k−1}(t)$ is defined as the lower simplex defined from $s_k$ by leaving out the vertex $v_t$. \n",
- "\n",
- "**Example:**\n",
- "\n",
- "In the simplex tree $\\mathcal{K}$ we have 1 2-simplex and 6 1-simplicies. By leaving out vertice $v_0=1$, $v_1=2$, $v_2=3$, we obtain the lower simplex $s_1=[2, 3]$, $s_2=[1, 3]$, $s_3=[1, 2]$, respectively. Therefore, the boundary operator on the 2-simplicies $\\partial_2$ should be a 6-by-1 matrix:\n",
- "\n",
- "$$\n",
- "\\partial_2 =\n",
- "\\begin{bmatrix}\n",
- "1 \\\\\n",
- "-1 \\\\\n",
- "1 \\\\\n",
- "0 \\\\\n",
- "0 \\\\\n",
- "0\n",
- "\\end{bmatrix}\n",
- "$$\n",
- "\n",
- "In the same way, the boundary operator on the 1-simplicies $\\partial_1$ is:\n",
- "\n",
- "$$\n",
- "\\partial_1 =\n",
- "\\begin{bmatrix}\n",
- "1 & 1 & 0 & 0 & 0 & 0 \\\\\n",
- "-1 & 0 & 1 & 0 & 0 & 0 \\\\\n",
- "0 & -1 & -1 & 1 & 1 & 0 \\\\\n",
- "0 & 0 & 0 & -1 & 0 & 1 \\\\\n",
- "0 & 0 & 0 & 0 & -1 & -1\n",
- "\\end{bmatrix}\n",
- "$$\n",
- "\n",
- "**Action:**\n",
- "\n",
- "Define a function that generates the boundary operator for a specified dimension, using a given simplex tree as input.\n",
- "\n",
- "**Answer:**"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "246e41ce",
- "metadata": {},
- "outputs": [],
- "source": [
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "81764f0d",
- "metadata": {},
- "source": [
- "**Instruction:**\n",
- "\n",
- "The combinatorial laplacian $\\Delta_k$ is defined as:\n",
- "\n",
- "$$\n",
- "\\Delta_k = \\left( \\partial_k \\right)^\\dagger \\partial_k \n",
- "+ \\partial_{k+1} \\left( \\partial_{k+1} \\right)^\\dagger\n",
- "$$\n",
- "\n",
- "The QPE algorithm will be used to estimate the number of zero eigenvalues of the Laplacian matrix. Since its exponential matrix serves as the unitary matrix in this algorithm, it must have dimensions $2^q \\times 2^q$, where $q$ represents the number of target qubits. It is recommanded to pad the combinatorial laplacian $\\Delta_k$ with an identity matrix with $\\tilde{\\lambda}_{max}/2$ in place of ones, where $\\tilde{\\lambda}_{max}$ is the estimate of the maximum eigenvalue of $\\Delta_k$ using the Gershgorin circle theorem ([details](https://mathworld.wolfram.com/GershgorinCircleTheorem.html)), such that:\n",
- "\n",
- "$$\n",
- "\\tilde{\\Delta}_k =\n",
- "\\begin{bmatrix}\n",
- "\\Delta_k & 0 \\\\\n",
- "0 & \\frac{\\widetilde{\\lambda}_{\\text{max}}}{2} \\cdot I_{2q - |S_k|}\n",
- "\\end{bmatrix}_{2q \\times 2q}\n",
- "$$\n",
- "\n",
- "where $\\Delta_k$ is the padded combinatorial laplacian and $q = \\lceil \\log_2 |S_k| \\rceil$ is the number of qubits this operator will act on. In QPE, as $2\\pi\\theta$ increases beyond $2\\pi$, the eigenvalues will start repeating due to their periodic form. Thus, $\\theta$ is restricted to $[0, 1)$. As $\\lambda \\to 2\\pi\\theta$ this means that $λ \\in [0, 2\\pi)$. Thus, we need to restrict the eigenvalues of the combinatorial laplacian to this range. This can be achieved by rescaling $\\tilde{\\Delta}_k$ by $\\delta/\\tilde{\\lambda}_{max}$ where $\\delta$ is slightly less than $2\\pi$. Thus, the rescaled matrix $H$ and the unitary marix $U$ for QPE are:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "H &= \\frac{\\delta}{\\tilde{\\lambda}_k} \\tilde{\\Delta}_k\\\\\n",
- "U &= e^{iH}\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "**Example:**\n",
- "\n",
- "In our example, the combinational laplacian is in the form of a $6\\times6$ matrix:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "\\Delta_1 &= (\\partial_1)^\\dagger \\partial_1 + \\partial_2 (\\partial_2)^\\dagger\\\\\n",
- "&=\n",
- "\\begin{bmatrix}\n",
- "3 & 0 & 0 & 0 & 0 & 0 \\\\\n",
- "0 & 3 & 0 & -1 & -1 & 0 \\\\\n",
- "0 & 0 & 3 & -1 & -1 & 0 \\\\\n",
- "0 & -1 & -1 & 2 & 1 & -1 \\\\\n",
- "0 & -1 & -1 & 1 & 2 & 1 \\\\\n",
- "0 & 0 & 0 & -1 & 1 & 2\n",
- "\\end{bmatrix}\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "It is padded with $\\tilde{\\lambda}_{max}=6$ and $\\delta=6$ to its nearest power of 2, which is 8 ($q=3$):\n",
- "\n",
- "$$\n",
- "\\begin{align}\n",
- "H_1 = \\begin{bmatrix}\n",
- "3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
- "0 & 3 & 0 & -1 & -1 & 0 & 0 & 0 \\\\\n",
- "0 & 0 & 3 & -1 & -1 & 0 & 0 & 0 \\\\\n",
- "0 & -1 & -1 & 2 & 1 & -1 & 0 & 0 \\\\\n",
- "0 & -1 & -1 & 1 & 2 & 1 & 0 & 0 \\\\\n",
- "0 & 0 & 0 & -1 & 1 & 2 & 0 & 0 \\\\\n",
- "0 & 0 & 0 & 0 & 0 & 0 & 3 & 0 \\\\\n",
- "0 & 0 & 0 & 0 & 0 & 0 & 0 & 3\n",
- "\\end{bmatrix}\n",
- "\\end{align}\n",
- "$$\n",
- "\n",
- "**Action:**\n",
- "\n",
- "Define a function to build the Laplacian, where Define a function that automatically determines whether the input Laplacian matrix requires padding. If padding is needed, the function will pad and rescale the matrix accordingly. Then, build the unitary based on the padded matrix, in the form of a circuit.\n",
- "\n",
- "**Answer:**"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8f818f78",
- "metadata": {},
- "outputs": [],
- "source": [
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "0acfea73",
- "metadata": {},
- "source": [
- "### Step 3: Applying QPE\n",
- "\n",
- "**Instruction:**\n",
- "\n",
- "The Betti number is the number of zero eigenvalues in the Laplacian ([ref](https://link.springer.com/article/10.1007/PL00009218)). \n",
- "\n",
- "$$\n",
- "\\beta_k = \\dim (\\ker(\\Delta_k))\n",
- "$$\n",
- "\n",
- "The betti curve is then a series of Betti numbers on different resolution threshold `epsilon`.\n",
- "\n",
- "To estimate the number of zero eigenvalues (nullity) in the padded Laplacian matrix (padding didn't add more zero eigenvalues), QPE algorithm is employed. The fundamental concept is that, if the target qubits start out in the maximally mixed state (shown below), which can be thought of as a random choice of an eigenstate, the proportion of all-zero states among all measured states is equal to the proportion of zero eigenvalues among all eigenvalues. Assume the all-zero state show up for $\\{i|\\tilde{\\theta}_i=0\\}$ times in $\\alpha$ shots, the probability of getting all-zero state $p(0)$ is given by:\n",
- "\n",
- "$$\n",
- "\\begin{align*}\n",
- "p(0) &= \\frac{\\left| \\{i \\mid \\tilde{\\theta}_i = 0\\} \\right|}{\\alpha} = \\frac{\\tilde{\\beta}_k}{2^q} \\\\\n",
- "\\implies \\tilde{\\beta}_k &= 2^q \\cdot p(0)\n",
- "\\end{align*}\n",
- "$$\n",
- "\n",
- "Where $\\tilde{\\beta}_k$ is the estimation of $k$-th Betti number. This estimation is then rounded to the nearest integer to obtain the final result."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "fc95fe90",
- "metadata": {},
- "source": [
- " NOTE:\n",
- "\n",
- "The unitary operator can be constructed by converting the exponential matrix of the Laplacian into a quantum circuit. For instance, in Qiskit, this can be implemented using `circuit.unitary(exp_matrix)`. Alternative methods for constructing the unitary operator will be optionally explored in the **BONUS** section at the end of this notebook.\n",
- "
\n",
- "\n",
- " \n",
- "BONUS EXERCISE: \n",
- "\n",
- "1. Why we should measure all-zero states?\n",
- "2. What is the difference between a maximally mixed state and the $(H\\ket{0})^{\\otimes n}$ state? Two possible aspects are:\n",
- "\n",
- "* Plotting of their density matrix.\n",
- "* Results from QPE.\n",
- "\n",
- "3. What parameters affect the accuracy of the estimation before rounding? For Laplacian matrices of varying sizes, how does the accuracy depend on these parameters? Within what range of values do these parameters guarantee (or are highly likely to produce) a correct final result?\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1cbe5cf2",
- "metadata": {},
- "source": [
- "### Step 4: Detecting market crashes"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7f27e81b",
- "metadata": {},
- "source": [
- "**Instruction:**\n",
- "\n",
- "At this point, we have betti curves for each window across our dataset, and we wish to use this to detect market crashes. One such way is to take the difference between these curves—or the pairwise distance—for successive windows and look for spikes. This can be done with the $L^p$ norm of the betti curve for each window, defined as follows:\n",
- "\n",
- "$$||x||_p = (\\sum_{n}^{i=1} |x_i|^p)^{1/p}$$\n",
- "\n",
- "Combining these pairwise distances into a vector, we get a single output curve we can analyze. Experiment with different values of $p$, but a good starting point is the $L^2$ Norm. Using this, it is possible to detect regions where a market crash is occuring. Comparing detected crashes with the price data indicates how accurate the crash detection methodology is ([ref](https://github.com/giotto-ai/stock-market-crashes/blob/master/Stock%20Market%20Crash%20Detection.ipynb)). \n",
- "\n",
- "**Action:**\n",
- "\n",
- "Use the $L^p$ norm to create pairwise distance curves for successive windows, and then use the results to define when a crash is occuring. Compare this with your data to see how well it performs. You may find the following classical solver is useful.\n",
- "\n",
- "**Answer:**"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "id": "5ef68f98",
- "metadata": {},
- "outputs": [],
- "source": [
- "from ripser import ripser\n",
- "\n",
- "def classical_betti_solver(point_cloud, epsilon, dim):\n",
- " '''Return the Betti number on a given point cloud.\n",
- " Args:\n",
- " point_cloud: the point cloud after applying the sliding window.\n",
- " epsilon: resolution threshold.\n",
- " dim: the dimension on which the Betti number is calculated\n",
- " '''\n",
- " result = ripser(point_cloud, maxdim=dim)\n",
- " diagrams = result[\"dgms\"]\n",
- " return len(\n",
- " [interval for interval in diagrams[dim] if interval[0] < epsilon < interval[1]]\n",
- " )\n",
- "\n",
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d1c2e4b5",
- "metadata": {},
- "source": [
- "## *BONUS:* Explore future directions of quantum TDA\n",
- "\n",
- "The following is a non-exhaustive list of possible next steps for the quantum TDA pipeline. It is recommended to at least explore 1 option or sub-option.\n",
- "\n",
- "- **Find more of applications of TDA in finance**:\n",
- "\n",
- " There are several directions where to extend the analysis. Most work on time series analysis has used persistent homology, and more specifically, the $L^{P}$ norm of persistence landscapes, which can be used to detect early warning signals of imminent market crashes. This is precisely studied in the seminal work by Gidea and Katz, see [ref](https://arxiv.org/abs/1703.04385)\n",
- " - Analyze financial correlation network and their degree of association with Betti curves or other topological features. From the time-series of multiple stock prices, we can build time-dependent correlation networks, which exhibit topological structures and might show some association to the evolution of betti curves or other topological data measures. Generally speaking, the cross correlations in a stock market will be in the form of a high-dimension topological space, with more complicated features. One can also think about other time varying financial graphs (e.g. cryptocurrencies). The following articles can help uncover more applications: \n",
- " \n",
- " - [Integral Betti signature confirms the hyperbolic geometry of brain, climate,and financial networks](https://www.arxiv.org/pdf/2406.15505)\n",
- " - [Using Topological Data Analysis (TDA) and Persistent Homology to Analyze the Stock Markets in Singapore and Taiwan](https://www.frontiersin.org/journals/physics/articles/10.3389/fphy.2021.572216/full)\n",
- " - Build a ML classifier or regressor on top of vectorized features such as Betti Curves (given their potential to identify trends, patterns or potential turning points in the market) to help with investment or risk management strategies. Show that Betti curves have some predictive skill, as key topological descriptors. See [ref1](https://arxiv.org/abs/2411.13881) and [ref2](https://www.sciencedirect.com/science/article/pii/S2405918823000235) for further information on the topic.\n",
- "- **A hybrid and more NISQ-friendly quantum TDA pipeline**:\n",
- " \n",
- " QPE remains primarily theoretical. Its circuits are simply too deep to run on real hardware. Come up an with iterative or hybrid quantum phase estimation protocol or use tools that increase the algorithmic performance when running quantum circuits on real hardware. Benchmark them against textbook-QPE circuits. Here are some proposals to subtitute the QPE part:\n",
- " - Variational Quantum Deflation (VQD) Algorithm: VQD is a quantum algorithm that uses a variational technique to find the k eigenvalues of the Hamiltonian H of a given system. [ref](https://quantum-journal.org/papers/q-2019-07-01-156/)\n",
- " - Variational Quantum Eigensolver (VQE): Using VQE to determine the spectra of adjancency or laplacian matrix. Inspired by: [ref](https://arxiv.org/pdf/1912.12366)\n",
- " \n",
- " Finally, run some circuits on simulator + real hardware and compare the performance (runtime, noise effects, # resources) of the new proposal to the QPE solution of the above sections.\n",
- "\n",
- "- **A proper procedure of encoding the Laplacian matrix to the unitary**:\n",
- "\n",
- " In Step 3, encoding the exponential matrix is recommanded. An alternative approach is to conduct the Paulis decomposition on the Laplacian, then followed by Trotterization. Can you implement this approach in your pipeline? What parameters influence the accuracy? Can you optimize your code to minimize the circuit depth?\n",
- "\n",
- "- **Extend the quantum TDA to extract persistent Betti numbers**: \n",
- " \n",
- " Implement a quantum TDA algorithm for persistent Betti numbers. Esstimating the persistent Betti numbers is a more general task than estimating the Betti number and it is more practical for TDA. It is an open problem to construct a quantum algorithm for the persistent Betti numbers in a way that is preferable for NISQ devices, and the only current implementation of a quantum algorithm for persistent betti number is shown [here](https://quantum-journal.org/papers/q-2022-12-07-873/pdf/)."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "0de4a8f2",
- "metadata": {},
- "outputs": [],
- "source": [
- "# write your code here"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c098599c-7ab0-48b1-9e33-51d22502310d",
- "metadata": {},
- "source": [
- "# This is the end of the challenge. Good luck!"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "tda",
- "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.11.11"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}