diff --git a/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb b/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb index db085b2..bc1e0d9 100644 --- a/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb +++ b/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb @@ -38,13 +38,50 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "import cereeberus.distance.ilp_solver_iterations as solver_iter\n", + "from cereeberus import Interleave, MapperGraph\n", "import matplotlib.pyplot as plt\n", - "import math" + "import math\n", + "from cereeberus.data import ex_mappergraphs" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "torus = ex_mappergraphs.torus(0, 2, 10, 13)\n", + "line = ex_mappergraphs.line(0, 16)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "ename": "KeyError", + "evalue": "14", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[8], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m myInt \u001b[38;5;241m=\u001b[39m Interleave(torus, line)\n\u001b[0;32m----> 2\u001b[0m \u001b[43mmyInt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:66\u001b[0m, in \u001b[0;36mInterleave.fit\u001b[0;34m(self, pulp_solver, verbose, max_n_for_error)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m min_n \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m checked_results:\n\u001b[1;32m 65\u001b[0m myAssgn \u001b[38;5;241m=\u001b[39m Assignment(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mF, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mG, n\u001b[38;5;241m=\u001b[39mmin_n)\n\u001b[0;32m---> 66\u001b[0m prob_status \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptimize\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 68\u001b[0m checked_results[min_n] \u001b[38;5;241m=\u001b[39m (prob_status, myAssgn)\n\u001b[1;32m 70\u001b[0m prob_status, myAssgn \u001b[38;5;241m=\u001b[39m checked_results[min_n]\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:1816\u001b[0m, in \u001b[0;36mAssignment.optimize\u001b[0;34m(self, pulp_solver)\u001b[0m\n\u001b[1;32m 1804\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21moptimize\u001b[39m(\u001b[38;5;28mself\u001b[39m, pulp_solver \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 1805\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively.\u001b[39;00m\n\u001b[1;32m 1806\u001b[0m \u001b[38;5;124;03m This function requires the `pulp` package to be installed.\u001b[39;00m\n\u001b[1;32m 1807\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1813\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 1814\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1816\u001b[0m map_dict, prob_status \u001b[38;5;241m=\u001b[39m \u001b[43msolve_ilp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1818\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m prob_status \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOptimal\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 1819\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/ilp.py:169\u001b[0m, in \u001b[0;36msolve_ilp\u001b[0;34m(myAssgn, pulp_solver, verbose)\u001b[0m\n\u001b[1;32m 167\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m up_or_down \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mup\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;66;03m#NOTE: the change in block indices\u001b[39;00m\n\u001b[1;32m 168\u001b[0m bou_n \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mB_up(other_map, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m'\u001b[39m)[block]\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[0;32m--> 169\u001b[0m bou_0 \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mB_up\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstarting_map\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m0\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[43mblock\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[1;32m 170\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m starting_map \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mF\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 171\u001b[0m map_V \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mphi(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m0\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mV\u001b[39m\u001b[38;5;124m'\u001b[39m)[block\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mget_array()\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/labeled_blocks.py:539\u001b[0m, in \u001b[0;36mLabeledBlockMatrix.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 529\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m 530\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 531\u001b[0m \u001b[38;5;124;03m Get the i'th block matrix.\u001b[39;00m\n\u001b[1;32m 532\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 537\u001b[0m \u001b[38;5;124;03m The item from the block matrix.\u001b[39;00m\n\u001b[1;32m 538\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 539\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[43mblocks\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n", + "\u001b[0;31mKeyError\u001b[0m: 14" + ] + } + ], + "source": [ + "myInt = Interleave(torus, line)\n", + "myInt.fit()" ] }, { diff --git a/Notebooks_Unused/n_search_fail.ipynb b/Notebooks_Unused/n_search_fail.ipynb new file mode 100644 index 0000000..51d2ac9 --- /dev/null +++ b/Notebooks_Unused/n_search_fail.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "id": "47f4a073", + "metadata": {}, + "outputs": [], + "source": [ + "%reload_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "282a714d", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import Interleave, MapperGraph, ReebGraph\n", + "\n", + "import cereeberus.data.ex_reebgraphs as ex_rg\n", + "import cereeberus.data.ex_mappergraphs as ex_mg" + ] + }, + { + "cell_type": "markdown", + "id": "e191f088", + "metadata": {}, + "source": [ + "# Create torus and line mappers of mismatched heights" + ] + }, + { + "cell_type": "markdown", + "id": "d1a6c618", + "metadata": {}, + "source": [ + "### It works when the differnce is 1" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6c6b8ae9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Interleaving distance: 3\n", + "Interleaving distance (with dist_fit): 3\n" + ] + } + ], + "source": [ + "line = ex_mg.line(a=0, b= 16)\n", + "torus = ex_mg.torus(a=0, b =2, c = 13, d = 17)\n", + "myInt = Interleave(torus, line)\n", + "result = myInt.fit()\n", + "dist_result = myInt.dist_fit() # dist_fit function optimizes with the distance matrix\n", + "print(f\"Interleaving distance: {result}\")\n", + "print(f\"Interleaving distance (with dist_fit): {dist_result}\")" + ] + }, + { + "cell_type": "markdown", + "id": "30459cf0", + "metadata": {}, + "source": [ + "### But fails when the difference is more than 1" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "ee3dea41", + "metadata": {}, + "outputs": [ + { + "ename": "KeyError", + "evalue": "17", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[9], line 4\u001b[0m\n\u001b[1;32m 2\u001b[0m torus \u001b[38;5;241m=\u001b[39m ex_mg\u001b[38;5;241m.\u001b[39mtorus(a\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, b \u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m, c \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m13\u001b[39m, d \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m18\u001b[39m)\n\u001b[1;32m 3\u001b[0m myInt \u001b[38;5;241m=\u001b[39m Interleave(torus, line)\n\u001b[0;32m----> 4\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mmyInt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 5\u001b[0m dist_result \u001b[38;5;241m=\u001b[39m myInt\u001b[38;5;241m.\u001b[39mdist_fit() \u001b[38;5;66;03m# dist_fit function optimizes with the distance matrix\u001b[39;00m\n\u001b[1;32m 6\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInterleaving distance: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mresult\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:66\u001b[0m, in \u001b[0;36mInterleave.fit\u001b[0;34m(self, pulp_solver, verbose, max_n_for_error)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m min_n \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m checked_results:\n\u001b[1;32m 65\u001b[0m myAssgn \u001b[38;5;241m=\u001b[39m Assignment(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mF, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mG, n\u001b[38;5;241m=\u001b[39mmin_n)\n\u001b[0;32m---> 66\u001b[0m prob_status \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptimize\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 68\u001b[0m checked_results[min_n] \u001b[38;5;241m=\u001b[39m (prob_status, myAssgn)\n\u001b[1;32m 70\u001b[0m prob_status, myAssgn \u001b[38;5;241m=\u001b[39m checked_results[min_n]\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:1816\u001b[0m, in \u001b[0;36mAssignment.optimize\u001b[0;34m(self, pulp_solver)\u001b[0m\n\u001b[1;32m 1804\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21moptimize\u001b[39m(\u001b[38;5;28mself\u001b[39m, pulp_solver \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 1805\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively.\u001b[39;00m\n\u001b[1;32m 1806\u001b[0m \u001b[38;5;124;03m This function requires the `pulp` package to be installed.\u001b[39;00m\n\u001b[1;32m 1807\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1813\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 1814\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1816\u001b[0m map_dict, prob_status \u001b[38;5;241m=\u001b[39m \u001b[43msolve_ilp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1818\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m prob_status \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOptimal\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 1819\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/ilp.py:169\u001b[0m, in \u001b[0;36msolve_ilp\u001b[0;34m(myAssgn, pulp_solver, verbose)\u001b[0m\n\u001b[1;32m 167\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m up_or_down \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mup\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;66;03m#NOTE: the change in block indices\u001b[39;00m\n\u001b[1;32m 168\u001b[0m bou_n \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mB_up(other_map, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m'\u001b[39m)[block]\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[0;32m--> 169\u001b[0m bou_0 \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mB_up\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstarting_map\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m0\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[43mblock\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[1;32m 170\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m starting_map \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mF\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 171\u001b[0m map_V \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mphi(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m0\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mV\u001b[39m\u001b[38;5;124m'\u001b[39m)[block\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mget_array()\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/labeled_blocks.py:539\u001b[0m, in \u001b[0;36mLabeledBlockMatrix.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 529\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m 530\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 531\u001b[0m \u001b[38;5;124;03m Get the i'th block matrix.\u001b[39;00m\n\u001b[1;32m 532\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 537\u001b[0m \u001b[38;5;124;03m The item from the block matrix.\u001b[39;00m\n\u001b[1;32m 538\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 539\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[43mblocks\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n", + "\u001b[0;31mKeyError\u001b[0m: 17" + ] + } + ], + "source": [ + "line = ex_mg.line(a =0, b= 16)\n", + "torus = ex_mg.torus(a=0, b =2, c = 13, d = 18)\n", + "myInt = Interleave(torus, line)\n", + "result = myInt.fit()\n", + "dist_result = myInt.dist_fit() # dist_fit function optimizes with the distance matrix\n", + "print(f\"Interleaving distance: {result}\")\n", + "print(f\"Interleaving distance (with dist_fit): {dist_result}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b442c31a", + "metadata": {}, + "source": [ + "I checked the code and I think the problem is in the boundary matrix computation. Which affects the edge-vertex parallelogram checks. Tried fixing, but haven't figured it out yet." + ] + }, + { + "cell_type": "markdown", + "id": "c17be785", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "interleavingenv", + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks_Unused/sandbox_fixing_optimization.ipynb b/Notebooks_Unused/sandbox_fixing_optimization.ipynb new file mode 100644 index 0000000..b7a9323 --- /dev/null +++ b/Notebooks_Unused/sandbox_fixing_optimization.ipynb @@ -0,0 +1,149 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "3d7c8bc8", + "metadata": {}, + "outputs": [], + "source": [ + "%reload_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a632ce73", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import Interleave, MapperGraph, Assignment\n", + "import cereeberus.data.ex_mappergraphs as ex_mg\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import networkx as nx" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "44844a7b", + "metadata": {}, + "outputs": [], + "source": [ + "mg1 = ex_mg.line(0, 20)\n", + "mg2 = ex_mg.torus(0, 2, 17, 18)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "0657e22b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "-\n", + "Trying n = 2...\n", + "n = 2, status = None\n", + "\n", + "-\n", + "Trying n = 3...\n", + "n = 3, status = None\n", + "\n", + "-\n", + "Trying n = 6...\n", + "n = 6, status = 1\n", + "\n", + "-\n", + "Trying n = 4...\n", + "n = 4, status = 1\n", + "\n", + "-\n", + "Trying n = 3...\n", + "n = 3, status = None\n" + ] + }, + { + "data": { + "text/plain": [ + "4" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "myInt = Interleave(mg1, mg2)\n", + "myInt.fit(verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6a5a4786", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import ReebGraph\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88cb17ae", + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "ReebGraph.__init__() got an unexpected keyword argument 'values'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[4], line 8\u001b[0m\n\u001b[1;32m 5\u001b[0m values \u001b[38;5;241m=\u001b[39m points[:, \u001b[38;5;241m2\u001b[39m]\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# Build and compute the Reeb graph\u001b[39;00m\n\u001b[0;32m----> 8\u001b[0m rg \u001b[38;5;241m=\u001b[39m \u001b[43mReebGraph\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvalues\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mvalues\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcoords\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpoints\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 9\u001b[0m rg\u001b[38;5;241m.\u001b[39mcompute()\n", + "\u001b[0;31mTypeError\u001b[0m: ReebGraph.__init__() got an unexpected keyword argument 'values'" + ] + } + ], + "source": [ + "# Example point cloud (you can replace this with your own)\n", + "points = np.random.rand(500, 3) # 500 points in 3D\n", + "\n", + "# Define a scalar function on these points — say, height (z-coordinate)\n", + "values = points[:, 2]\n", + "\n", + "# Build and compute the Reeb graph\n", + "rg = ReebGraph()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "interleavingenv", + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/cereeberus/cereeberus/distance/ilp.py b/cereeberus/cereeberus/distance/ilp.py index f7e805f..575e93a 100644 --- a/cereeberus/cereeberus/distance/ilp.py +++ b/cereeberus/cereeberus/distance/ilp.py @@ -1,7 +1,5 @@ -# from cereeberus import ReebGraph, MapperGraph, Interleave from .labeled_blocks import LabeledBlockMatrix as LBM from .labeled_blocks import LabeledMatrix as LM -# import cereeberus.data.ex_mappergraphs as ex_mg import matplotlib.pyplot as plt import numpy as np @@ -27,9 +25,6 @@ def build_map_matrices(myAssgn, map_results): Returns: dict: the dictionary containing the final map matrices """ - - # get the function values - func_vals = myAssgn.all_func_vals() # create a dictionary to store the final matrices final_LBMs = {} @@ -92,6 +87,324 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): psi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + # create the decision variables (NOTE: these are all the decision variables for all the diagrams) + for thickening in ['0', 'n']: + for obj_type in ['V', 'E']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + # create lp variables for phi + n_rows = myAssgn.phi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.phi(thickening, obj_type)[block].get_array().shape[1] + + phi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('phi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += phi_vars[block][thickening][obj_type][(a,b)] == myAssgn.phi(thickening, obj_type)[block].get_array()[a][b] + + + # create lp variables for psi + n_rows = myAssgn.psi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.psi(thickening, obj_type)[block].get_array().shape[1] + + psi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('psi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += psi_vars[block][thickening][obj_type][(a,b)] == myAssgn.psi(thickening, obj_type)[block].get_array()[a][b] + + # create the other decision variables + z_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + + map_product_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + for obj_type in ['V', 'E']: + for starting_map in ['F', 'G']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + if starting_map == 'F': + n_rows_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.phi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[1] + + # set the z variables + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_1) for c in range(n_cols_1)), cat='Binary') + + else: + n_rows_1 = myAssgn.phi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.psi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_2 = myAssgn.phi('n', obj_type)[block].get_array().shape[1] + + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_2) for c in range(n_cols_1)), cat='Binary') + + + + + # create the constraints + for block in func_vals: + for starting_map in ['F', 'G']: + if block not in myAssgn.all_func_vals(map=starting_map): # If the block is not in the starting map, skip it. This is to account for the graphs with different function values + continue + # set the other map based on starting map + if starting_map == 'F': + other_map = 'G' + else: + other_map = 'F' + + for up_or_down in ['up', 'down']: # deals with 1 (up, down) and 2 (up, down) + + if block == myAssgn.all_func_vals(map=starting_map)[-1]: # skip the last block for this type of diagrams + continue + + #set the matrices + if up_or_down == 'up': #NOTE: the change in block indices + bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block+1].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block+1]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block+1].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block+1]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + else: + bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + + # set the dimensions + shape_m_mix = map_V.shape[0] + shape_n_mix = map_V.shape[1] + shape_o_mix = bou_n.shape[1] + shape_p_mix = map_E.shape[1] + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + for i in range(shape_m_mix): + for k in range(shape_p_mix): + # inner difference + first_term = pulp.lpSum([map_V_vars[i,j] * bou_0[j][k] for j in range(shape_n_mix)]) + second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) + + # total expression + mixed_expression = (first_term - second_term) + + prob += mixed_expression == 0 + + # constraint 2: each column sums to 1 + for j in range(shape_n_mix): + prob += pulp.lpSum(map_V_vars[h,j] for h in range(shape_m_mix)) == 1 + + for k in range(shape_p_mix): + prob += pulp.lpSum(map_E_vars[l, k] for l in range(shape_o_mix)) == 1 + + + + for obj_type in ['V', 'E']: # deals with 3, 4, 5, 6, 7, 8, 9, 10 + if obj_type == 'E' and block == func_vals[-1]: # skip the last block for this type of diagrams. This is because we don't have an edge with the highest function value + continue + + # multiply inclusion matrices. Needed for the triangles + i_n_i_0 = myAssgn.I(starting_map, 'n', obj_type)[block].get_array() @ myAssgn.I(starting_map, '0', obj_type)[block].get_array() + + # write inclusion matrices for easier reference. Needed for the parallelograms + inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() + inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() + + + # set map matrices for easier reference + if starting_map == 'F': + map_0_para_vars = phi_vars[block]['0'][obj_type] + map_n_para_vars = phi_vars[block]['n'][obj_type] + + if starting_map == 'G': + map_0_para_vars = psi_vars[block]['0'][obj_type] + map_n_para_vars = psi_vars[block]['n'][obj_type] + + # set the dimensions + shape_m_tri = i_n_i_0.shape[0] # for triangles + + shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_n_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array().shape[0] # for parallelograms + + if starting_map == 'F': + shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles + + + + shape_m_para = myAssgn.phi('n', obj_type)[block].get_array().shape[0] # for parallelograms + shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms + else: + shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles + + # changed shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_m_para = myAssgn.psi('n', obj_type)[block].get_array().shape[0] # for parallelograms + shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms + + + + + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + # for triangles + for h in range(shape_m_tri): + + tri_expression = pulp.lpSum( (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for k in range(shape_o_tri)) + + # prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression + + + prob += tri_expression == 0 + + + + + # for parallelograms + for i in range(shape_m_para): + for k in range(shape_p_para): + # inner difference + first_term = pulp.lpSum([map_n_para_vars[i,j] * inc_0_para[j][k] for j in range(shape_n_para)]) + second_term = pulp.lpSum([inc_n_para[i][l] * map_0_para_vars[l,k] for l in range(shape_o_para)]) + + + # total expression + para_expression = first_term - second_term + + prob += para_expression == 0 + + + # constraint 2: map_multiplication and z relation. This is for triangles + for i in range(shape_m_tri): + for k in range(shape_o_tri): + prob += map_product_vars[block][starting_map][obj_type][i,k] == pulp.lpSum(z_vars[block][starting_map][obj_type][i,j,k] for j in range(shape_n_tri)) + + # constraint 3: z is less than either of the map values and greater than sum of the map values minus 1. This is for triangles + for i in range(shape_m_tri): + for j in range(shape_n_tri): + for k in range(shape_o_tri): + if starting_map == 'F': + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= psi_vars[block]['n'][obj_type][i,j] + phi_vars[block]['0'][obj_type][j,k] - 1 + else: + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= phi_vars[block]['n'][obj_type][i,j] + psi_vars[block]['0'][obj_type][j,k] - 1 + + # constraint 4: each column sums to 1 + if starting_map == 'F': + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + for k in range(shape_o_tri): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + else: + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + + for k in range(shape_o_tri): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + + # Set the objective function + prob += 0 + + # solve the problem + if pulp_solver == 'GUROBI': + prob.solve(pulp.GUROBI_CMD(msg=0)) + else: + prob.solve(pulp.PULP_CBC_CMD(msg=0)) + + status_str = pulp.LpStatus[prob.status] + + # prob.solve(pulp.GUROBI_CMD(msg=0)) + if prob.status != 1: + return None, status_str + + # create a dictionary to store the results + map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} + + # make the results a LabeledBlockMatrix + final_maps = build_map_matrices(myAssgn, map_results) + + if verbose: + print("Status:", pulp.LpStatus[prob.status]) + # prob.writeLP("model.lp") # Write the model in LP format + + + # return results + return final_maps, status_str + + + +##------------------- ILP Optimization with Distance Matrix -------------------## + +def solve_ilp_dist(myAssgn, pulp_solver = None, verbose=False): + + """ + Function to solve the ILP optimization problem for interleaving maps. The function creates a linear programming problem using the PuLP library and solves it to find the optimal interleaving maps. + + Parameters: + myAssgn (Assignment): the Assignment object containing the interleaving maps and other relevant data + pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. + verbose (bool): whether to print the optimization status and results + + Returns: + tuple: a tuple containing the final interleaving maps (as a dictionary of LabeledBlockMatrices) and the optimized loss value + """ + # function values + func_vals = myAssgn.all_func_vals() + + # create the ILP problem + prob = pulp.LpProblem("Interleave_Optimization_Problem", pulp.LpMinimize) + + # create empty dictionaries to store the decision variables + phi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + psi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + + # create the decision variables (NOTE: these are all the decision variables for all the diagrams) for thickening in ['0', 'n']: for obj_type in ['V', 'E']: @@ -216,7 +529,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_E_vars = psi_vars[block]['0']['E'] # set the dimensions - shape_m_mix = dist_n_other.shape[0] + shape_m_mix = map_V.shape[0] shape_n_mix = map_V.shape[1] shape_o_mix = bou_n.shape[1] shape_p_mix = map_E.shape[1] @@ -270,7 +583,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_n_para_vars = psi_vars[block]['n'][obj_type] # set the dimensions - shape_m_tri = dist_2n_starting.shape[0] # for triangles + shape_m_tri = i_n_i_0.shape[0] # for triangles shape_m_para = dist_2n_other.shape[0] # for parallelograms shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms @@ -393,12 +706,8 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): if verbose: print(f"The optimized loss is: {pulp.value(minmax_var)}") print("Status:", pulp.LpStatus[prob.status]) - prob.writeLP("model.lp") # Write the model in LP format + # prob.writeLP("model.lp") # Write the model in LP format - # if get_thickened_maps: - # final_maps = build_map_matrices(myAssgn, map_results, thickening = 'n') - - # return final_maps, pulp.value(minmax_var) # return results return final_maps, pulp.value(minmax_var) diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index fc4c4bc..37c3ce0 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -7,7 +7,7 @@ from ..compute.unionfind import UnionFind from .labeled_blocks import LabeledBlockMatrix as LBM from .labeled_blocks import LabeledMatrix as LM -from .ilp import solve_ilp +from .ilp import solve_ilp, solve_ilp_dist from ..compute.utils import HiddenPrints import sys, os @@ -35,118 +35,116 @@ def __init__(self, F, G): self.n = np.inf self.assignment = None - def old_fit(self, pulp_solver = None, verbose = False, max_n_for_error = 100, - ): + + + def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ - Compute the interleaving distance between the two Mapper graphs. + Finds the smallest feasible n using exponential + binary search and stores the optimal assignment in self.assignment and the n in self.n. Parameters: pulp_solver (pulp.LpSolver): - The solver to use for the ILP optimization. If None, the default solver is used. + solver to use for ILP. If None, the default solver is used. verbose (bool, optional): If True, print the progress of the optimization. Defaults to False. max_n_for_error (int, optional): - The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. - printOptimizerOutput (bool, optional): - If True, the output of the PULP optimizer is printed. Defaults to False. + The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. ####NOTE: this can be replaced by the bounding box. Returns: - Interleave : - The Interleave object with the computed interleaving distance. - """ - # -- Search through possible n values to find the smallest one that still allows for interleaving + int: smallest feasible n + """ + # catch the results to avoid recomputation + checked_results = {} + + + # Step 0: Check for smallest possible n (n=0 when they both have same function ranges) + + min_n = max(abs(self.F.min_f()-self.G.min_f()), abs(self.F.max_f()-self.G.max_f())) + + + if min_n not in checked_results: + myAssgn = Assignment(self.F, self.G, n=min_n) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + + checked_results[min_n] = (prob_status, myAssgn) + + prob_status, myAssgn = checked_results[min_n] + if verbose: + print(f"\n-\nTrying n = {min_n}...") + print(f"n = {min_n}, status = {prob_status}") + - # -- Dictionary to store the search data - # -- This will store the Loss for each value of n so we don't - # -- have to recompute it each time. - # -- The distance bound can be determined by search_data[key] = Loss implies d_I <= key + Loss - search_data = {} - # -- Set the initial value of the distance bound to be infinity - distance_bound = np.inf + if prob_status == 1: + self.n = min_n + self.assignment = myAssgn + return self.n + + # Step 1: Exponential search for first feasible n + low, high = min_n, min_n+1 + found_feasible_n = False - # -- Do an expoential search for the smallest n that allows for interleaving - # -- Start with n = 1 and double it until the Loss is less than or equal to 0 - N = 1 - found_max = False - - while not found_max: - # -- Compute the assignment for the current value of n - # -- If the Loss is less than or equal to 0, we have found a valid n - # -- Otherwise, double n and try again + while high <= max_n_for_error: + if high not in checked_results: + myAssgn = Assignment(self.F, self.G, n=high) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + checked_results[high] = (prob_status, myAssgn) - try: - if verbose: - print(f"\n-\nTrying n = {N}...") - myAssgn = Assignment(self.F, self.G, n = N) - - Loss = myAssgn.optimize() - search_data[N] = Loss - - if verbose: - print(f"n = {N}, Loss = {Loss}, distance_bound = {N + Loss}") - - except ValueError as e: - # -- If we get a ValueError, it means the current n is too small - # -- So we double n and try again - search_data[N] = np.inf - N *= 2 - continue + prob_status, myAssgn = checked_results[high] + + if verbose: + print(f"\n-\nTrying n = {high}...") + print(f"n = {high}, status = {prob_status}") - if Loss <= 0: - # -- If the Loss is less than or equal to 0, we have found a valid n - max_N = N - found_max = True - if verbose: - print(f"Found valid maximum n = {N} with Loss = {Loss}. Moving on to binary search.") - else: - # -- If the Loss is greater than 0, we need to try a larger n - N *= 2 - - if N > max_n_for_error: - # -- If we have exceeded the maximum value of n, raise an error. Useful while we're checking this while loop + + + if prob_status == 1: + found_feasible_n = True + break + low, high = high, high * 2 # double n + else: + if not found_feasible_n: raise ValueError(f"Interleaving distance not found for n <= {max_n_for_error}.") - - # now binary search in [0, max_N] to find the smallest n that gives a loss of 0 - - # -- Set the initial values for the binary search - low = 1 - high = max_N - while low < high: + high = min(high, max_n_for_error) # Clamp to max allowed + + # Step 2: Binary search for smallest feasible n in [low+1, high] + best_n = high + + while low <= high: mid = (low + high) // 2 + if mid not in checked_results: + myAssgn = Assignment(self.F, self.G, n=mid) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + checked_results[mid] = (prob_status, myAssgn) + + prob_status, myAssgn = checked_results[mid] + + if verbose: + print(f"\n-\nTrying n = {mid}...") + print(f"n = {mid}, status = {prob_status}") + - try: - myAssgn = Assignment(self.F, self.G, n = mid) - Loss = myAssgn.optimize() - search_data[mid] = Loss - except ValueError as e: - # -- If we get a ValueError, it means the current n is too small - search_data[mid] = np.inf - low = mid + 1 - continue - - if Loss <= 0: - # -- If the Loss is less than or equal to 0, we have found a valid n - high = mid + if prob_status == 1: + best_n = mid + high = mid - 1 # try smaller n + else: - low = mid + 1 - - # -- Set the final value of n to be the smallest one that gives a loss of 0 - self.n = low - # -- Set the assignment to be the one that minimizes the interleaving distance - self.assignment = Assignment(self.F, self.G, n = self.n) - Loss = self.assignment.optimize(pulp_solver = pulp_solver,) - - - # Raise error if the loss isn't 0 - if Loss > 0: - raise ValueError(f"Final fit object is not an interleaving. Loss = {Loss}. N = {self.n}.") + low = mid+1 # update to last infeasible n + + # Step 3: Final validation + self.n = best_n + prob_status, myAssgn = checked_results[self.n] + self.assignment = myAssgn + + if prob_status != 1: + raise ValueError(f"Final fit object is not an interleaving. Status = {prob_status} for n = {self.n}.") return self.n - - def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): + + + + def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ Compute the interleaving distance between the two Mapper graphs. @@ -163,28 +161,41 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): The interleaving distance, that is, the smallest n for which loss is zero. """ - # step 0: search for n=0 - myAssgn = Assignment(self.F, self.G, n = 0) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) + # catch the results to avoid recomputation + checked_results = {} + + # Step 0: Check for smallest possible n (n=0 when they both have same function ranges) + min_n = max(abs(self.F.min_f()-self.G.min_f()), abs(self.F.max_f()-self.G.max_f())) # minimum possible n based on function ranges + + + if min_n not in checked_results: + myAssgn = Assignment(self.F, self.G, n = min_n) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[min_n] = (Loss, myAssgn) + Loss, myAssgn = checked_results[min_n] if verbose: - print(f"\n-\nTrying n = 0...") - print(f"n = 0, Loss = {Loss}, distance_bound = {0 + Loss}") + print(f"\n-\nTrying n = {min_n}...") + print(f"n = {min_n}, Loss = {Loss}, distance_bound = {min_n + Loss}") # if loss is 0, we're done if Loss == 0: - self.n = 0 + self.n = min_n self.assignment = myAssgn return self.n # step 1: exponential search for the upperbound - low, high = 1, 1 + low, high = min_n, min_n+1 found_valid_n = False while high <= max_n_for_error: try: - myAssgn = Assignment(self.F, self.G, n = high) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) + if high not in checked_results: + myAssgn = Assignment(self.F, self.G, n = high) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[high] = (Loss, myAssgn) + + Loss, myAssgn = checked_results[high] if verbose: print(f"\n-\nTrying n = {high}...") @@ -209,9 +220,13 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): while low <= high: mid = (low + high) // 2 try: - myAssgn = Assignment(self.F, self.G, n = mid) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) - + if mid not in checked_results: + myAssgn = Assignment(self.F, self.G, n = mid) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[mid] = (Loss, myAssgn) + + Loss, myAssgn = checked_results[mid] + if verbose: print(f"\n-\nTrying n = {mid}...") print(f"n = {mid}, Loss = {Loss}, distance_bound = {mid + Loss}") @@ -226,12 +241,10 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # validate the final solution self.n = best_n - self.assignment = Assignment(self.F, self.G, n = self.n) - final_loss = self.assignment.optimize(pulp_solver = pulp_solver) - - if final_loss != 0: - raise ValueError(f"Unexpected non-zero loss (Loss={final_loss}) for n={self.n}") - + Loss, myAssgn = checked_results[self.n] + if Loss != 0: + raise ValueError(f"Unexpected non-zero loss (Loss={Loss}) for n={self.n}") + self.assignment = myAssgn return self.n @@ -383,8 +396,11 @@ def __init__(self, F, G, # Containers for matrices for later self.B_down_ = {'F':{}, 'G':{}} # boundary matrix + # self._B_down = None # don't build the boundary matrices unless needed self.B_up_ = {'F':{}, 'G':{}} # boundary matrix - self.D_ = {'F':{}, 'G':{}} # distance matrix + # self._B_up = None # don't build the boundary matrices unless needed + # self.D_ = {'F':{}, 'G':{}} # distance matrix + self._D = None # don't build the distance matrices unless needed self.I_ = {'F':{}, 'G':{}} # induced maps self.val_to_verts = {'F':{}, 'G':{}} # dictionaries from function values to vertices @@ -494,13 +510,17 @@ def __init__(self, F, G, # --- # --- - # Build the distance matrices - for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: - for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. - self.D_[name][key] = {'V':{}, 'E':{}} + # Build the distance matrices + # We don't actually need these for the ILP, so we won't build them unless needed. + # If you need them, you can access them using the property self.D_ which will build them if they haven't been built yet. - self.D_[name][key]['V'] = metagraph[key].thickening_distance_matrix() - self.D_[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') + + # for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: + # for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. + # self.D_[name][key] = {'V':{}, 'E':{}} + + # self.D_[name][key]['V'] = metagraph[key].thickening_distance_matrix() + # self.D_[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') # End distance matrices @@ -567,7 +587,91 @@ def __init__(self, F, G, # --- - + ### ---------------- + # Properties + ### ---------------- + + # @property + # def B_down_(self): + # if self._B_down is None: + # self._B_down = {'F':{}, 'G':{}} + # for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + + # B_down = LBM() + + # for i in self.val_to_verts[graph_name][key]: + # if i in self.val_to_edges[graph_name][key]: + # edges = self.val_to_edges[graph_name][key][i] + # verts_down = self.val_to_verts[graph_name][key][i] + # B_down[i] = LM(rows = verts_down, cols = edges) + + # for e in edges: + # B_down[i][e[0],e] = 1 + + # min_i = min(list(self.val_to_verts[graph_name][key].keys())) + # max_i = max(list(self.val_to_verts[graph_name][key].keys())) + + # min_verts = self.val_to_verts[graph_name][key][min_i] + # max_verts = self.val_to_verts[graph_name][key][max_i] + + # B_down[min_i-1] = LM(rows = min_verts, cols = []) + # B_down[max_i] = LM(rows = max_verts, cols = []) + + # self._B_down[graph_name][key] = B_down + # # self.B_up_[graph_name][key] = B_up + # return self._B_down + + # @property + # def B_up_(self): + # if self._B_up is None: + # self._B_up = {'F':{}, 'G':{}} + # for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + + # # B_down = LBM() + # B_up = LBM() + + # for i in self.val_to_verts[graph_name][key]: + # if i in self.val_to_edges[graph_name][key]: + # edges = self.val_to_edges[graph_name][key][i] + # # verts_down = self.val_to_verts[graph_name][key][i] + # verts_up = self.val_to_verts[graph_name][key][i+1] + # # B_down[i] = LM(rows = verts_down, cols = edges) + # B_up[i] = LM(rows = verts_up, cols = edges) + + # for e in edges: + # # B_down[i][e[0],e] = 1 + # B_up[i][e[1],e] = 1 + + # min_i = min(list(self.val_to_verts[graph_name][key].keys())) + # max_i = max(list(self.val_to_verts[graph_name][key].keys())) + + # min_verts = self.val_to_verts[graph_name][key][min_i] + # max_verts = self.val_to_verts[graph_name][key][max_i] + + # B_up[min_i-1] = LM(rows = min_verts, cols = []) + # B_up[max_i] = LM(rows = max_verts, cols = []) + + # # self.B_down_[graph_name][key] = B_down + # self._B_up[graph_name][key] = B_up + # return self._B_up + + @property + def D_(self): + if self._D is None: + self._D = {'F':{}, 'G':{}} + for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: + for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. + self._D[name][key] = {'V':{}, 'E':{}} + + self._D[name][key]['V'] = metagraph[key].thickening_distance_matrix() + self._D[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') + return self._D + + + + ### ---------------- # Functions for getting stuff out of all the dictionaries @@ -968,6 +1072,10 @@ def draw_all_graphs(self, figsize = (15,10), **kwargs): self.G('2n').draw(ax = axs[1,2], **kwargs) axs[1,2].set_title(r'$G_{2n}$') + # turn on grid for all axes + for ax in axs.flatten(): + ax.grid(True) + return fig, axs def draw_I(self, graph = 'F', key = '0', obj_type = 'V', ax = None, **kwargs): @@ -1334,6 +1442,7 @@ def parallelogram_Edge_Vert_matrix(self, Result_Dist = self.D(end_graph, 'n', 'V') @ Result else: if up_or_down == 'down': # tau_i \to \sigma_i + Top = self.get_interleaving_map(maptype, '0', 'V')[func_val] @ self.B_down(start_graph, '0')[func_val] Bottom = self.B_down(end_graph, 'n')[func_val] @ self.get_interleaving_map(maptype, '0', 'E')[func_val] Result = Top - Bottom @@ -1677,7 +1786,7 @@ def loss_by_block(self): return max(loss_list) - def all_func_vals(self): + def all_func_vals(self, map = None): """ Get all the function values that are in the graphs. @@ -1685,12 +1794,16 @@ def all_func_vals(self): list : A list of all the function values. """ + if map == 'F': + return self.F().get_function_values() + elif map == 'G': + return self.G().get_function_values() + else: + all_func_vals = set(self.F().get_function_values()) | set(self.G().get_function_values()) + all_func_vals = list(all_func_vals) + all_func_vals.sort() - all_func_vals = set(self.F().get_function_values()) | set(self.G().get_function_values()) - all_func_vals = list(all_func_vals) - all_func_vals.sort() - - return all_func_vals + return all_func_vals def optimize(self, pulp_solver = None): """Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively. @@ -1699,11 +1812,15 @@ def optimize(self, pulp_solver = None): Parameters: pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. Returns: - float : - The loss value found by the ILP solver. + int or None: + Returns 1 if an optimal solution was found and None otherwise. + """ - - map_dict, loss_val = solve_ilp(self, pulp_solver = pulp_solver) + + map_dict, prob_status = solve_ilp(self, pulp_solver = pulp_solver) + + if prob_status != 'Optimal': + return None self.phi_['0'] = {'V': map_dict['phi_0_V'], 'E': map_dict['phi_0_E']} self.phi_['n'] = {'V': map_dict['phi_n_V'], 'E': map_dict['phi_n_E']} @@ -1711,7 +1828,27 @@ def optimize(self, pulp_solver = None): self.psi_['n'] = {'V': map_dict['psi_n_V'], 'E': map_dict['psi_n_E']} - return loss_val - + return 1 - + def dist_optimize(self, pulp_solver = None): + """Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively. + ## NOTE: This version uses the D matrices in the constraints, which makes it slower. ## + This function requires the `pulp` package to be installed. + + Parameters: + pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. + Returns: + int or None: + Returns 1 if an optimal solution was found and None otherwise. + + """ + + map_dict, loss_val = solve_ilp_dist(self, pulp_solver = pulp_solver) + + self.phi_['0'] = {'V': map_dict['phi_0_V'], 'E': map_dict['phi_0_E']} + self.phi_['n'] = {'V': map_dict['phi_n_V'], 'E': map_dict['phi_n_E']} + self.psi_['0'] = {'V': map_dict['psi_0_V'], 'E': map_dict['psi_0_E']} + self.psi_['n'] = {'V': map_dict['psi_n_V'], 'E': map_dict['psi_n_E']} + + + return loss_val \ No newline at end of file