{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example usage\n", "\n", "This is an example tutorial of how to `tdads` can be used to analyze persistence diagrams. First, let's make sure we have all the necessary dependencies:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "vscode": { "languageId": "shellscript" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: ripser in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (0.6.8)\n", "Requirement already satisfied: Cython in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from ripser) (3.0.10)\n", "Requirement already satisfied: numpy in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from ripser) (1.26.4)\n", "Requirement already satisfied: scipy in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from ripser) (1.13.1)\n", "Requirement already satisfied: scikit-learn in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from ripser) (1.5.0)\n", "Requirement already satisfied: persim in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from ripser) (0.3.5)\n", "Requirement already satisfied: matplotlib in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from persim->ripser) (3.9.0)\n", "Requirement already satisfied: hopcroftkarp in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from persim->ripser) (1.2.5)\n", "Requirement already satisfied: deprecated in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from persim->ripser) (1.2.14)\n", "Requirement already satisfied: joblib in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from persim->ripser) (1.4.2)\n", "Requirement already satisfied: threadpoolctl>=3.1.0 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from scikit-learn->ripser) (3.5.0)\n", "Requirement already satisfied: wrapt<2,>=1.10 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from deprecated->persim->ripser) (1.16.0)\n", "Requirement already satisfied: contourpy>=1.0.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (1.2.1)\n", "Requirement already satisfied: cycler>=0.10 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (0.12.1)\n", "Requirement already satisfied: fonttools>=4.22.0 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (4.53.0)\n", "Requirement already satisfied: kiwisolver>=1.3.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (1.4.5)\n", "Requirement already satisfied: packaging>=20.0 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (24.1)\n", "Requirement already satisfied: pillow>=8 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (10.3.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (3.1.2)\n", "Requirement already satisfied: python-dateutil>=2.7 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib->persim->ripser) (2.9.0.post0)\n", "Requirement already satisfied: six>=1.5 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib->persim->ripser) (1.16.0)\n", "Note: you may need to restart the kernel to use updated packages.\n", "Requirement already satisfied: matplotlib in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (3.9.0)\n", "Requirement already satisfied: contourpy>=1.0.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (1.2.1)\n", "Requirement already satisfied: cycler>=0.10 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (0.12.1)\n", "Requirement already satisfied: fonttools>=4.22.0 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (4.53.0)\n", "Requirement already satisfied: kiwisolver>=1.3.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (1.4.5)\n", "Requirement already satisfied: numpy>=1.23 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (1.26.4)\n", "Requirement already satisfied: packaging>=20.0 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (24.1)\n", "Requirement already satisfied: pillow>=8 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (10.3.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (3.1.2)\n", "Requirement already satisfied: python-dateutil>=2.7 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from matplotlib) (2.9.0.post0)\n", "Requirement already satisfied: six>=1.5 in /Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "%pip install ripser\n", "%pip install matplotlib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll import all the packages we'll need:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from tdads.machine_learning import *\n", "from tdads.inference import *\n", "from numpy.random import uniform\n", "from numpy import array\n", "from math import cos, sin, pi\n", "from ripser import ripser\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this tutorial we'll create helper functions for generating data and persistence diagrams from circle and spheres:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# function to create a circle dataset and\n", "# compute its diagram\n", "def circle_diagram():\n", " # sample 100 points from the unit circle\n", " theta = uniform(low = 0, high = 2*pi, size = 100)\n", " data = array([[cos(theta[i]), sin(theta[i])] for i in range(100)])\n", " # compute persistence diagram\n", " diag = ripser(data, maxdim = 2)\n", " return [data, diag]\n", "\n", "# function to create a sphere dataset and\n", "# compute its diagram\n", "def sphere_diagram():\n", " # sample 100 points from the unit sphere\n", " phi = uniform(low = 0, high = 2*pi, size = 100)\n", " theta = uniform(low = 0, high = pi, size = 100)\n", " data = array([[sin(theta[i])*cos(phi[i]), sin(theta[i])*sin(phi[i]), cos(theta[i])] for i in range(100)])\n", " # compute persistence diagram\n", " diag = ripser(data, maxdim = 2)\n", " return [data, diag]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our dataset will be 5 circle datasets (and diagrams) and 5 sphere datasets (and diagrams), totalling 10:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "result = [circle_diagram(), circle_diagram(), circle_diagram(), circle_diagram(), circle_diagram(),\n", " sphere_diagram(), sphere_diagram(), sphere_diagram(), sphere_diagram(), sphere_diagram()]\n", "data = [r[0] for r in result]\n", "diagrams = [r[1] for r in result]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will use the bootstrap procedure to determine the significant topological features in each diagram:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def diag_fun(X, thresh):\n", " return ripser(X = X, thresh = thresh, maxdim = 2)\n", "boot = diagram_bootstrap(diag_fun = diag_fun, dims = [0,1,2], alpha = 0.01)\n", "boot_diagrams = [boot.compute(X = d, thresh = 2) for d in data]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The subsetted diagrams show that only the first five diagrams have one loop and only the last five diagrams have one void:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Num clusters:2, num loops: 1, num voids: 5\n", "Num clusters:1, num loops: 1, num voids: 3\n", "Num clusters:1, num loops: 1, num voids: 3\n", "Num clusters:1, num loops: 1, num voids: 3\n", "Num clusters:2, num loops: 1, num voids: 2\n", "Num clusters:4, num loops: 3, num voids: 1\n", "Num clusters:22, num loops: 1, num voids: 1\n", "Num clusters:23, num loops: 2, num voids: 1\n", "Num clusters:5, num loops: 2, num voids: 1\n", "Num clusters:14, num loops: 1, num voids: 1\n" ] } ], "source": [ "for i in range(10):\n", " print('Num clusters:' + str(len(boot_diagrams[i]['subsetted_diagram'][0])) + ', num loops: ' + str(len(boot_diagrams[i]['subsetted_diagram'][1])) + ', num voids: ' + str(len(boot_diagrams[i]['subsetted_diagram'][2])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It seems like we have resolved the two groups of persistence diagrams, but a 2D MDS projection of the 10 diagrams visually confirms as much:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "Array must be symmetric", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31m_RemoteTraceback\u001b[0m Traceback (most recent call last)", "\u001b[0;31m_RemoteTraceback\u001b[0m: \n\"\"\"\nTraceback (most recent call last):\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/externals/loky/process_executor.py\", line 463, in _process_worker\n r = call_item()\n ^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/externals/loky/process_executor.py\", line 291, in __call__\n return self.fn(*self.args, **self.kwargs)\n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py\", line 598, in __call__\n return [func(*args, **kwargs)\n ^^^^^^^^^^^^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py\", line 598, in \n return [func(*args, **kwargs)\n ^^^^^^^^^^^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/utils/parallel.py\", line 129, in __call__\n return self.function(*args, **kwargs)\n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/manifold/_mds.py\", line 104, in _smacof_single\n dissimilarities = check_symmetric(dissimilarities, raise_exception=True)\n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n File \"/Users/jibaccount/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/utils/validation.py\", line 1502, in check_symmetric\n raise ValueError(\"Array must be symmetric\")\nValueError: Array must be symmetric\n\"\"\"", "\nThe above exception was the direct cause of the following exception:\n", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[13], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m mds \u001b[39m=\u001b[39m diagram_mds(p \u001b[39m=\u001b[39m \u001b[39mfloat\u001b[39m(\u001b[39m'\u001b[39m\u001b[39minf\u001b[39m\u001b[39m'\u001b[39m), dim \u001b[39m=\u001b[39m \u001b[39m2\u001b[39m) \u001b[39m# for 2-dimensional homology\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m emb \u001b[39m=\u001b[39m mds\u001b[39m.\u001b[39;49mfit_transform(diagrams)\n\u001b[1;32m 3\u001b[0m plt\u001b[39m.\u001b[39mscatter(emb[:,\u001b[39m0\u001b[39m], emb[:,\u001b[39m1\u001b[39m], color \u001b[39m=\u001b[39m [\u001b[39m'\u001b[39m\u001b[39mred\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mred\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mred\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mred\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mred\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m,\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m])\n\u001b[1;32m 4\u001b[0m plt\u001b[39m.\u001b[39mxlabel(\u001b[39m'\u001b[39m\u001b[39mEmbedding dim 1\u001b[39m\u001b[39m'\u001b[39m)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/src/tdads/machine_learning.py:100\u001b[0m, in \u001b[0;36mdiagram_mds.fit_transform\u001b[0;34m(self, X, y)\u001b[0m\n\u001b[1;32m 98\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(X, \u001b[39mtype\u001b[39m([\u001b[39m0\u001b[39m,\u001b[39m1\u001b[39m])):\n\u001b[1;32m 99\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mException\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mWhen precomputed is True, X must be a ndarray distance matrix of persistence diagrams.\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m--> 100\u001b[0m X_new \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mMDS\u001b[39m.\u001b[39;49mfit_transform(X, y)\n\u001b[1;32m 101\u001b[0m \u001b[39mreturn\u001b[39;00m X_new\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/base.py:1473\u001b[0m, in \u001b[0;36m_fit_context..decorator..wrapper\u001b[0;34m(estimator, *args, **kwargs)\u001b[0m\n\u001b[1;32m 1466\u001b[0m estimator\u001b[39m.\u001b[39m_validate_params()\n\u001b[1;32m 1468\u001b[0m \u001b[39mwith\u001b[39;00m config_context(\n\u001b[1;32m 1469\u001b[0m skip_parameter_validation\u001b[39m=\u001b[39m(\n\u001b[1;32m 1470\u001b[0m prefer_skip_nested_validation \u001b[39mor\u001b[39;00m global_skip_validation\n\u001b[1;32m 1471\u001b[0m )\n\u001b[1;32m 1472\u001b[0m ):\n\u001b[0;32m-> 1473\u001b[0m \u001b[39mreturn\u001b[39;00m fit_method(estimator, \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/manifold/_mds.py:638\u001b[0m, in \u001b[0;36mMDS.fit_transform\u001b[0;34m(self, X, y, init)\u001b[0m\n\u001b[1;32m 635\u001b[0m \u001b[39melif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdissimilarity \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39meuclidean\u001b[39m\u001b[39m\"\u001b[39m:\n\u001b[1;32m 636\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdissimilarity_matrix_ \u001b[39m=\u001b[39m euclidean_distances(X)\n\u001b[0;32m--> 638\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39membedding_, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mstress_, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mn_iter_ \u001b[39m=\u001b[39m smacof(\n\u001b[1;32m 639\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mdissimilarity_matrix_,\n\u001b[1;32m 640\u001b[0m metric\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mmetric,\n\u001b[1;32m 641\u001b[0m n_components\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mn_components,\n\u001b[1;32m 642\u001b[0m init\u001b[39m=\u001b[39;49minit,\n\u001b[1;32m 643\u001b[0m n_init\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mn_init,\n\u001b[1;32m 644\u001b[0m n_jobs\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mn_jobs,\n\u001b[1;32m 645\u001b[0m max_iter\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mmax_iter,\n\u001b[1;32m 646\u001b[0m verbose\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mverbose,\n\u001b[1;32m 647\u001b[0m eps\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49meps,\n\u001b[1;32m 648\u001b[0m random_state\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrandom_state,\n\u001b[1;32m 649\u001b[0m return_n_iter\u001b[39m=\u001b[39;49m\u001b[39mTrue\u001b[39;49;00m,\n\u001b[1;32m 650\u001b[0m normalized_stress\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mnormalized_stress,\n\u001b[1;32m 651\u001b[0m )\n\u001b[1;32m 653\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39membedding_\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/utils/_param_validation.py:186\u001b[0m, in \u001b[0;36mvalidate_params..decorator..wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 184\u001b[0m global_skip_validation \u001b[39m=\u001b[39m get_config()[\u001b[39m\"\u001b[39m\u001b[39mskip_parameter_validation\u001b[39m\u001b[39m\"\u001b[39m]\n\u001b[1;32m 185\u001b[0m \u001b[39mif\u001b[39;00m global_skip_validation:\n\u001b[0;32m--> 186\u001b[0m \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 188\u001b[0m func_sig \u001b[39m=\u001b[39m signature(func)\n\u001b[1;32m 190\u001b[0m \u001b[39m# Map *args/**kwargs to the function signature\u001b[39;00m\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/manifold/_mds.py:369\u001b[0m, in \u001b[0;36msmacof\u001b[0;34m(dissimilarities, metric, n_components, init, n_init, n_jobs, max_iter, verbose, eps, random_state, return_n_iter, normalized_stress)\u001b[0m\n\u001b[1;32m 367\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m 368\u001b[0m seeds \u001b[39m=\u001b[39m random_state\u001b[39m.\u001b[39mrandint(np\u001b[39m.\u001b[39miinfo(np\u001b[39m.\u001b[39mint32)\u001b[39m.\u001b[39mmax, size\u001b[39m=\u001b[39mn_init)\n\u001b[0;32m--> 369\u001b[0m results \u001b[39m=\u001b[39m Parallel(n_jobs\u001b[39m=\u001b[39;49mn_jobs, verbose\u001b[39m=\u001b[39;49m\u001b[39mmax\u001b[39;49m(verbose \u001b[39m-\u001b[39;49m \u001b[39m1\u001b[39;49m, \u001b[39m0\u001b[39;49m))(\n\u001b[1;32m 370\u001b[0m delayed(_smacof_single)(\n\u001b[1;32m 371\u001b[0m dissimilarities,\n\u001b[1;32m 372\u001b[0m metric\u001b[39m=\u001b[39;49mmetric,\n\u001b[1;32m 373\u001b[0m n_components\u001b[39m=\u001b[39;49mn_components,\n\u001b[1;32m 374\u001b[0m init\u001b[39m=\u001b[39;49minit,\n\u001b[1;32m 375\u001b[0m max_iter\u001b[39m=\u001b[39;49mmax_iter,\n\u001b[1;32m 376\u001b[0m verbose\u001b[39m=\u001b[39;49mverbose,\n\u001b[1;32m 377\u001b[0m eps\u001b[39m=\u001b[39;49meps,\n\u001b[1;32m 378\u001b[0m random_state\u001b[39m=\u001b[39;49mseed,\n\u001b[1;32m 379\u001b[0m normalized_stress\u001b[39m=\u001b[39;49mnormalized_stress,\n\u001b[1;32m 380\u001b[0m )\n\u001b[1;32m 381\u001b[0m \u001b[39mfor\u001b[39;49;00m seed \u001b[39min\u001b[39;49;00m seeds\n\u001b[1;32m 382\u001b[0m )\n\u001b[1;32m 383\u001b[0m positions, stress, n_iters \u001b[39m=\u001b[39m \u001b[39mzip\u001b[39m(\u001b[39m*\u001b[39mresults)\n\u001b[1;32m 384\u001b[0m best \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39margmin(stress)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/sklearn/utils/parallel.py:67\u001b[0m, in \u001b[0;36mParallel.__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m 62\u001b[0m config \u001b[39m=\u001b[39m get_config()\n\u001b[1;32m 63\u001b[0m iterable_with_config \u001b[39m=\u001b[39m (\n\u001b[1;32m 64\u001b[0m (_with_config(delayed_func, config), args, kwargs)\n\u001b[1;32m 65\u001b[0m \u001b[39mfor\u001b[39;00m delayed_func, args, kwargs \u001b[39min\u001b[39;00m iterable\n\u001b[1;32m 66\u001b[0m )\n\u001b[0;32m---> 67\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39msuper\u001b[39;49m()\u001b[39m.\u001b[39;49m\u001b[39m__call__\u001b[39;49m(iterable_with_config)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:2007\u001b[0m, in \u001b[0;36mParallel.__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m 2001\u001b[0m \u001b[39m# The first item from the output is blank, but it makes the interpreter\u001b[39;00m\n\u001b[1;32m 2002\u001b[0m \u001b[39m# progress until it enters the Try/Except block of the generator and\u001b[39;00m\n\u001b[1;32m 2003\u001b[0m \u001b[39m# reaches the first `yield` statement. This starts the asynchronous\u001b[39;00m\n\u001b[1;32m 2004\u001b[0m \u001b[39m# dispatch of the tasks to the workers.\u001b[39;00m\n\u001b[1;32m 2005\u001b[0m \u001b[39mnext\u001b[39m(output)\n\u001b[0;32m-> 2007\u001b[0m \u001b[39mreturn\u001b[39;00m output \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mreturn_generator \u001b[39melse\u001b[39;00m \u001b[39mlist\u001b[39;49m(output)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:1650\u001b[0m, in \u001b[0;36mParallel._get_outputs\u001b[0;34m(self, iterator, pre_dispatch)\u001b[0m\n\u001b[1;32m 1647\u001b[0m \u001b[39myield\u001b[39;00m\n\u001b[1;32m 1649\u001b[0m \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_backend\u001b[39m.\u001b[39mretrieval_context():\n\u001b[0;32m-> 1650\u001b[0m \u001b[39myield from\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_retrieve()\n\u001b[1;32m 1652\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mGeneratorExit\u001b[39;00m:\n\u001b[1;32m 1653\u001b[0m \u001b[39m# The generator has been garbage collected before being fully\u001b[39;00m\n\u001b[1;32m 1654\u001b[0m \u001b[39m# consumed. This aborts the remaining tasks if possible and warn\u001b[39;00m\n\u001b[1;32m 1655\u001b[0m \u001b[39m# the user if necessary.\u001b[39;00m\n\u001b[1;32m 1656\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_exception \u001b[39m=\u001b[39m \u001b[39mTrue\u001b[39;00m\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:1754\u001b[0m, in \u001b[0;36mParallel._retrieve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 1747\u001b[0m \u001b[39mwhile\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_wait_retrieval():\n\u001b[1;32m 1748\u001b[0m \n\u001b[1;32m 1749\u001b[0m \u001b[39m# If the callback thread of a worker has signaled that its task\u001b[39;00m\n\u001b[1;32m 1750\u001b[0m \u001b[39m# triggered an exception, or if the retrieval loop has raised an\u001b[39;00m\n\u001b[1;32m 1751\u001b[0m \u001b[39m# exception (e.g. `GeneratorExit`), exit the loop and surface the\u001b[39;00m\n\u001b[1;32m 1752\u001b[0m \u001b[39m# worker traceback.\u001b[39;00m\n\u001b[1;32m 1753\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_aborting:\n\u001b[0;32m-> 1754\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_raise_error_fast()\n\u001b[1;32m 1755\u001b[0m \u001b[39mbreak\u001b[39;00m\n\u001b[1;32m 1757\u001b[0m \u001b[39m# If the next job is not ready for retrieval yet, we just wait for\u001b[39;00m\n\u001b[1;32m 1758\u001b[0m \u001b[39m# async callbacks to progress.\u001b[39;00m\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:1789\u001b[0m, in \u001b[0;36mParallel._raise_error_fast\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 1785\u001b[0m \u001b[39m# If this error job exists, immediately raise the error by\u001b[39;00m\n\u001b[1;32m 1786\u001b[0m \u001b[39m# calling get_result. This job might not exists if abort has been\u001b[39;00m\n\u001b[1;32m 1787\u001b[0m \u001b[39m# called directly or if the generator is gc'ed.\u001b[39;00m\n\u001b[1;32m 1788\u001b[0m \u001b[39mif\u001b[39;00m error_job \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m-> 1789\u001b[0m error_job\u001b[39m.\u001b[39;49mget_result(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtimeout)\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:745\u001b[0m, in \u001b[0;36mBatchCompletionCallBack.get_result\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m 739\u001b[0m backend \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mparallel\u001b[39m.\u001b[39m_backend\n\u001b[1;32m 741\u001b[0m \u001b[39mif\u001b[39;00m backend\u001b[39m.\u001b[39msupports_retrieve_callback:\n\u001b[1;32m 742\u001b[0m \u001b[39m# We assume that the result has already been retrieved by the\u001b[39;00m\n\u001b[1;32m 743\u001b[0m \u001b[39m# callback thread, and is stored internally. It's just waiting to\u001b[39;00m\n\u001b[1;32m 744\u001b[0m \u001b[39m# be returned.\u001b[39;00m\n\u001b[0;32m--> 745\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_return_or_raise()\n\u001b[1;32m 747\u001b[0m \u001b[39m# For other backends, the main thread needs to run the retrieval step.\u001b[39;00m\n\u001b[1;32m 748\u001b[0m \u001b[39mtry\u001b[39;00m:\n", "File \u001b[0;32m~/Documents/shael/PHD/Personal_Research/Miscellaneous_Scripts/tdads/venv/lib/python3.11/site-packages/joblib/parallel.py:763\u001b[0m, in \u001b[0;36mBatchCompletionCallBack._return_or_raise\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 761\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m 762\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mstatus \u001b[39m==\u001b[39m TASK_ERROR:\n\u001b[0;32m--> 763\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_result\n\u001b[1;32m 764\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_result\n\u001b[1;32m 765\u001b[0m \u001b[39mfinally\u001b[39;00m:\n", "\u001b[0;31mValueError\u001b[0m: Array must be symmetric" ] } ], "source": [ "mds = diagram_mds(dim = 1) # for 1-dimensional homology\n", "emb = mds.fit_transform(diagrams)\n", "plt.scatter(emb[:,0], emb[:,1], color = ['red','red','red','red','red','blue','blue','blue','blue','blue'])\n", "plt.xlabel('Embedding dim 1')\n", "plt.ylabel('Embedding dim 2')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we can use a permutation test of our two suspected groups to statistically capture their differences in dimensions 0, 1 and 2 (all are near/below 0.05):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pt = perm_test(p = float('inf'), iterations = 50, dims = [0,1,2])\n", "res = pt.test([[d for d in diagrams[0:5]], [d for d in diagrams[5:10]]])\n", "res['p_values']" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 4 }