diff --git a/experiments.ipynb b/experiments.ipynb index bb217a1..3d8263a 100644 --- a/experiments.ipynb +++ b/experiments.ipynb @@ -22,7 +22,7 @@ "from scipy import sparse\n", "from typing import Tuple, List\n", "import argparse\n", - "import zigzag.zz as zz\n", + "import gril.gril as gril\n", "import torch\n", "import torch.nn as nn\n", "from scipy.spatial import Delaunay" @@ -64,7 +64,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAicElEQVR4nO3da2xUdf7H8U872FY2pVKR4eKERrPeooAp0rRINm66kmhweWBoqAEkXnbYYgh1Ewa5VMXSrquGB1TnLytZH8gWlqhrVlJ27S4PjE1Yi03cXcCwyGXVKXQJbRcilZnzfzBOodiWOdOZ+Z0z5/1KJkcP58z8htMyn/n+LifPsixLAAAAhuSbbgAAAPA2wggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAo8aZbkAyYrGYvv76axUXFysvL890cwAAQBIsy1J/f7+mTZum/PyR6x+uCCNff/21AoGA6WYAAIAUnDp1SjfffPOIf+6KMFJcXCwp/mYmTJhguDUAACAZfX19CgQCg5/jI3FFGEl0zUyYMIEwAgCAy1xriAUDWAEAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAArhAOS2Vl8S1yC2EEAOAKzc3SiRPxLXILYQQA4AqhkDRjRnybaVRhsivPsizLdCOupa+vTyUlJert7dWECRNMNwcAkOPKyuJVmBkzpOPHTbfGvZL9/KYyAgDAVbJZhQFhBACQBrnWrREMxisiwaDplngDYQQAMGYMLsVYEEYAAGNGtwbGggGsAAAgIxjACgAAXIEwAgAAjCKMAAAAowgjAOBRuTYd1+28fD0YwAoAHhQOS6tWSdEoq4w6RS6u+soAVgDAiJqb40HE52M6rlN4eXo0YQQAPCjxwbdtG6uMOsVIq756ofuGbhoAABzMzd03dNMAAJADru6+ycVKCZURAMhR4XB8bEgoRFdMLnFTpYTKCAB4HDevy025ONCVMAIAOSoXP7Qw8kBXN6ObBgAAZATdNAAAwBUIIwAAwCjCCAAAMIowAgAAjCKMAIDL5OKiV0gvt/2MMJsGAFzGTYtewQyn/IwwmwYAchTrh+Ba3PYzQhgBAJdIlN6l3Fv0CunltoXRCCMA4BIs745UOX0MCWEEAFzCbaV3OIfTgyxhBABcwm2ldziH04MsYQQAHKy2Vho3Lr4FUuX0IEsYAQAH271bikbjW2CsnDp2hDACAA62eLHk88W3wFg5dewIYQQAHGznTunSpfgWGCunjh1hBVYAAJARrMAKAABcgTACAA7j1EGGQKYQRgDAYZw6yBDIFMIIADhIOCz190ulpc4bZAhkCmEEABykuVk6e1YqLnbuAlVAuhFGAMBBnDr1EsgkpvYCAICMYGovAABwBcIIAAAwijACAICHOeHO0IQRAAA8zAl3hiaMAADgYU64M3RKYaSlpUVlZWUqKipSRUWFDhw4MOrxW7du1e23367rr79egUBAa9as0bfffptSgwEAQPo44c7QtsPIrl27VF9fr4aGBh08eFCzZs3SggULdPr06WGP37lzp0KhkBoaGnTo0CG99dZb2rVrl5577rkxNx4AALif7TDy2muv6amnntKKFSt01113KRwOa/z48dqxY8ewx3/yySeaN2+eamtrVVZWpgcffFBLliy5ZjUFAAB4g60wMjAwoM7OTlVXV19+gvx8VVdXq6OjY9hzqqqq1NnZORg+jh07pr179+qhhx4a8XUuXryovr6+IQ8AAJCbxtk5uKenR9FoVH6/f8h+v9+vw4cPD3tObW2tenp6dP/998uyLF26dEnBYHDUbpqmpia98MILdpoGAABcKuOzafbv368tW7bo9ddf18GDB/Xuu+/qww8/1ObNm0c8Z926dert7R18nDp1KtPNBAAAhtiqjEyaNEk+n0/d3d1D9nd3d2vKlCnDnrNx40YtXbpUTz75pCTpnnvu0fnz5/X0009r/fr1ys//YR4qLCxUYWGhnaYBAACXslUZKSgoUHl5udrb2wf3xWIxtbe3q7KycthzLly48IPA4fP5JEkuuEcfAADIMFuVEUmqr6/X8uXLNWfOHM2dO1dbt27V+fPntWLFCknSsmXLNH36dDU1NUmSFi5cqNdee0333nuvKioqdPToUW3cuFELFy4cDCUAAMC7bIeRmpoanTlzRps2bVIkEtHs2bPV1tY2OKj15MmTQyohGzZsUF5enjZs2KCvvvpKN910kxYuXKjGxsb0vQsAAOBaeZYL+kr6+vpUUlKi3t5eTZgwwXRzAADIGeGw1NwshUJSMJje507285t70wCAQ4TDUllZfAtkS3OzdOJEfGsKYQQAHMIJHwrwnlBImjEjvjXF9pgRAEBmVFVJ//lPfAtkSzCY/u4Zu6iMAIBDfPKJFI3Gt4CXEEYAwCGcUC6HtzhlnBKzaQAA8Kiysvg4pRkzpOPH0//8zKYBAACjcko1jjACAA7jlNI5cl8wGK+IMIAVADBEYorvqlUEEngDYQQAHCYUkny++Mwa1hyBFxBGAMBhgkFp2zZn9OUj9zixG5DZNAAAeMiNN0pnz0qlpdJ//5vZ12I2DQAAcAXCCAAAHtLYGO8CbGw03ZLLCCMA4HBO7OOHezllOu+VGDMCAA6XzT5+IJ0YMwIAAFxRWSOMAIDDObGPH+6RWETPyWvWEEYAwOGc2McP93DK/WdGw5gRAACQEYwZAQAArkAYAQAARhFGAACAUYQRAABgFGEEAHKAG9aSAEZCGAGAHOCGtSSAkRBGACAHuGEtCaRfrlTEWGcEAHJMOByvkIRCLJSW68rK4hWxGTPiC+M5DeuMAIBHJbpsnn1WGjdOqq013SJkSq5UxKiMAECOSVRGTp2SYjHJ55MuXTLdKngRlREA8KjEvWxqauJBZPFi0y0CRkcYAYActXNnvCKyc6fplmAscmWQ6mgIIwAAOJgXpm0TRgDAo7zwjTsX5Mog1dEwgBUAPMrp00LhfgxgBQCMKvGNu6qKCgnMojICAB5HhQSZQmUEAJAUL4xJgLMRRgDA4xLrkoy2dDyDXVPH39210U0DALimRFdOaalUXMx9b+zwcjcY3TQAgLRJdOVIub/mRbrRDXZthBEAwDUlunIaG3/4wUo3xOiS6QbzOrppAABj4uVuCIyObhoAQFaM1g2RK1WTXHkfTkUYAQCMyWjdEMncV8WpH/RXtssL94cxiTACAMiYZAZvmvigTwSN2lrpxhvjj6vD0JXtYhBqZjFmBABgVKLykM3pwolxLj6fFI3G91095sVEu3JNsp/fhBEAgOckgkZVlbRvX3xfYyOhI90IIwAAwChm0wAAAFcgjAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMColMJIS0uLysrKVFRUpIqKCh04cGDU48+dO6e6ujpNnTpVhYWFuu2227R3796UGgwAAHLLOLsn7Nq1S/X19QqHw6qoqNDWrVu1YMECHTlyRJMnT/7B8QMDA/rZz36myZMna8+ePZo+fbpOnDihG264IR3tBwAALmf7rr0VFRW67777tG3bNklSLBZTIBDQM888o1Ao9IPjw+GwfvOb3+jw4cO67rrrUmokd+0FAMB9MnLX3oGBAXV2dqq6uvryE+Tnq7q6Wh0dHcOe88EHH6iyslJ1dXXy+/26++67tWXLFkWj0RFf5+LFi+rr6xvyAAAAuclWGOnp6VE0GpXf7x+y3+/3KxKJDHvOsWPHtGfPHkWjUe3du1cbN27Uq6++qpdeemnE12lqalJJScngIxAI2GkmAABwkYzPponFYpo8ebLefPNNlZeXq6amRuvXr1c4HB7xnHXr1qm3t3fwcerUqUw3EwAAGGJrAOukSZPk8/nU3d09ZH93d7emTJky7DlTp07VddddJ5/PN7jvzjvvVCQS0cDAgAoKCn5wTmFhoQoLC+00DQAAuJStykhBQYHKy8vV3t4+uC8Wi6m9vV2VlZXDnjNv3jwdPXpUsVhscN8XX3yhqVOnDhtEAACAt9jupqmvr9f27dv19ttv69ChQ1q5cqXOnz+vFStWSJKWLVumdevWDR6/cuVKnT17VqtXr9YXX3yhDz/8UFu2bFFdXV363gUAAHAt2+uM1NTU6MyZM9q0aZMikYhmz56ttra2wUGtJ0+eVH7+5YwTCAS0b98+rVmzRjNnztT06dO1evVqrV27Nn3vAgAAuJbtdUZMYJ0RAADcJyPrjAAAAKQbYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRKYWRlpYWlZWVqaioSBUVFTpw4EBS57W2tiovL0+LFi1K5WUBAEAOsh1Gdu3apfr6ejU0NOjgwYOaNWuWFixYoNOnT4963vHjx/WrX/1K8+fPT7mxAAAg99gOI6+99pqeeuoprVixQnfddZfC4bDGjx+vHTt2jHhONBrVY489phdeeEG33HLLmBoMAAByi60wMjAwoM7OTlVXV19+gvx8VVdXq6OjY8TzXnzxRU2ePFlPPPFEUq9z8eJF9fX1DXkgt4TD0o03xh/h8MjHlJUN/fPh9gEA3M1WGOnp6VE0GpXf7x+y3+/3KxKJDHvOxx9/rLfeekvbt29P+nWamppUUlIy+AgEAnaaCRdobpbOno0/mptHPubEiaF/Ptw+AIC7ZXQ2TX9/v5YuXart27dr0qRJSZ+3bt069fb2Dj5OnTqVwVYiVVdXKexULUIhqbQ0/giFRj5mxoyhfz7cvnRJploDAEi/PMuyrGQPHhgY0Pjx47Vnz54hM2KWL1+uc+fO6Y9//OOQ47u6unTvvffK5/MN7ovFYpLi3TtHjhzRrbfees3X7evrU0lJiXp7ezVhwoRkm4s0Cofj1YhQSAoG4/vKyuJVihkzpOPHf/j/bpNov+Te9wAATpLs57etykhBQYHKy8vV3t4+uC8Wi6m9vV2VlZU/OP6OO+7Q559/rq6ursHHI488ogceeEBdXV10v7jIcN0jV1cpMlm1yIZkqjWMWQGA9LNVGZHiU3uXL1+u//u//9PcuXO1detW7d69W4cPH5bf79eyZcs0ffp0NTU1DXv+448/rnPnzun9999P+jWpjJg3XGXEi9xe/QGAbMpIZUSSampq9Morr2jTpk2aPXu2urq61NbWNjio9eTJk/rmm29SbzkyKjEu4kc/sjc2IhiMf/h6OYhI167+UDkBAPtsV0ZMoDKSPleOi5D4hp9uVE4A4LKMVUbgbolxEePHjz42Aqlx+7gZADCByggAAMgIKiMAAMAVCCOAYSy2BsDrCCOAYcksjQ8AuYwwAhiWzGJrAJDLGMAKAAAyggGsQI5hQTUAuYow4kB86GA4w90fCAByAWHEQRIhZP16PnTwQ1cuqEZgBZBLGDPiIImlxEtLpeJibkqHkbHsPAA3YMyIS1z5DTfxzbexkZvSYXRUSQDkEiojhvENF2OV+Bny+aRt2wixAJyDyohLcGM1jFUoFA8i0SjjjAC4E2HEsGCQLhmMTTAYr4hcGWrpugHgJoSRLOHDAZl0dahNTANeuVKqrTXaNAC4JsJIlrBGBLLpym6/3bvNtQMAkkEYyRLGhiCbgkFpyZL4WJLFi023BgBGx2wawCPC4XhljvVrAGQLs2kADEFXIQCnIowAHsFCaQCcim4awINYbA9ANtBNA2BEiSpJVRUVEgDmURkBPOzGG6WzZ6X8fKmlhYGtANKLygiApMViDGwFYA5hBPCwxkaptDT+YA0cAKYQRlLATATkimBQ+u9/449gkJ9tAGYwZiQFiX720tL4P+JArmCWDYB0YswIANu4bQEAEwgjKWhsjP+D3dhouiVAel19918AyAbCSJKu7EvnH2wAANKHMJIk7usBAEBmEEaSRF86vIoZNgAyjdk0AEbFDBsAqWI2DYC0oCoIINMIIwBGlRiwLdFdAyAzCCMjuO8+KS8vvgXAIG4AmUMYGcGnnw7dAl4XCsVXHe7vpzoCIL0IIyOYM2foFvC6YFAqLo7fCoHqCIB0IoyM4O9/lywrvgUQx2BWAJnA1F4AAJARTO0FAACuQBgBAABGEUYAAIBRhJHvcf8NAADMIIx8jwWdAAAwgzCieDWkvz++oBNTFoHkUVEEkA5M7RV3JQVSxe8OgNEwtdeGqirJ54tvASSPRdAApIPnKyPhsLRqlRSN8u0OAIB0ojKSpObmeBDx+fh2BwCACZ4PI4ky87Zt8RuBAQCA7PJ8Nw0AAMgMumkAAIArEEYAjNl990l5efEtANhFGAEwZp9+OnQLAHYQRgCM2Zw5Q7cAYMc40w0A4H5//7vpFgBwMyojAADAKMIIAAAwKqUw0tLSorKyMhUVFamiokIHDhwY8djt27dr/vz5mjhxoiZOnKjq6upRjwcAAN5iO4zs2rVL9fX1amho0MGDBzVr1iwtWLBAp0+fHvb4/fv3a8mSJfrb3/6mjo4OBQIBPfjgg/rqq6/G3HgAzhEOx+/iGw6bbgkAt7G9AmtFRYXuu+8+bdu2TZIUi8UUCAT0zDPPKJTEzV2i0agmTpyobdu2admyZUm9JiuwAs5XViadOMENJwFclpEVWAcGBtTZ2anq6urLT5Cfr+rqanV0dCT1HBcuXNB3332n0tLSEY+5ePGi+vr6hjwAOFviPk/ccBKAXbbCSE9Pj6LRqPx+/5D9fr9fkUgkqedYu3atpk2bNiTQXK2pqUklJSWDj0AgYKeZAAwIBuMVEW44CcCurM6maW5uVmtrq9577z0VFRWNeNy6devU29s7+Dh16lQWWwkAALLJ1qJnkyZNks/nU3d395D93d3dmjJlyqjnvvLKK2pubtZHH32kmTNnjnpsYWGhCgsL7TQNAAC4lK3KSEFBgcrLy9Xe3j64LxaLqb29XZWVlSOe9/LLL2vz5s1qa2vTHNaLBgAAV7DdTVNfX6/t27fr7bff1qFDh7Ry5UqdP39eK1askCQtW7ZM69atGzz+17/+tTZu3KgdO3aorKxMkUhEkUhE//vf/9L3LgAYx9ReAKmyfW+ampoanTlzRps2bVIkEtHs2bPV1tY2OKj15MmTys+/nHHeeOMNDQwM6NFHHx3yPA0NDXr++efH1noAjtHcHJ/a29zMIFYA9theZ8QE1hkBnC8cjgeRUIgwAiAu2c9vwggAAMiIjCx6BgDDYbwIgLHwdBiprZXGjYtvAaTu2Wfj40WefdZ0SwC4kafDyO7dUjQa3wJI3bffDt0CgB2eDiOLF0s+X3wLIHU1NfHfpZoa0y0B4EaeDiM7d0rbtkmffEJfNzAWO3dKly7FtwBgl6fDiDR0bQQAAJB9ng8j3PYcAACzPB9GuO05kDqm9AJIB8+HEQCpo5sTQDoQRgCkjG5OAOnAcvAAACAjWA7eJvq+AQAwgzDyvUTf96pVBBIAALKJMPK9UCi+gmQ0ymA8AACyaZzpBjhFYmpvczOD8QAAyCYqI1dgzRHg2rjbNYB0I4wAsIW7XQNIN8IIAFu42zWAdGOdEQAAkBGsMwIAAFyBMDKKcFj60Y+k/HwG6wEAkCmEkVE0N0sXLkiWJbW2mm4NYAarEwPINMLIKEIhKS8v/t/XX2+2LYAp3JkXQKYRRkYRDEqvvx6/K+nPf863Q3hPOCz190ulpSwGCCBzmE2TpLKy+LfDGTPiC6MBXsDPPYCxYDZNmoVC8X+Q+XYIL+HnHkA2UBkBAAAZQWUkw5hhAABAehBGUsQMAwAA0oMwkiL60gEASA/GjAAAgIxgzAgAAHAFwggAADCKMJJmzLIBAMAewkiaJWbZrFpFIAEAIBmEkTQLhSSfT4pGmfYL56FyB8CJCCNpFgxK27Yx7RfOxPo4AJyIMJIBwWD8pmLBIN9E4SysjwPAiVhnJMMSdz0tLZWKi+MfAsGg6VYBAJB5rDPiEIlvot9+Gw8l69ebbhEAAM5CGMmwRJdNUZHplgAA4EyEkSxpbIxXSBobTbcEAABnIYxkyZWDWoF0Y6A0ADcjjAA5gCm7ANyMMALkAKbsAnAzwogDUXKHXXQDAnAzwogDUXIHAHgJYcSBriy5UyUBAOQ6VmB1uMQKrjNmxMvwAAC4BSuw5oirByZSKQEA5BoqIy7DvW4AAG5BZSRHJSolEoNccwXVLgBeRxhxmcQUzsTy8levK1FbK40bF9/Cua4MIMyeAuB1hBGXGmldid27pWhU2rWLb9tOM1IAYcEyAF7HmJEcU1sbDySFhdKFC8zCcZIrZ0aFQpeDCGN+AOQqxox41M6d0qVL0quvjvxtmzEK2XH13/OVFRBWTAWAy6iMeBBrl2QHf88AvI7KCEY00hgFKibpxVgQAEhOSmGkpaVFZWVlKioqUkVFhQ4cODDq8X/4wx90xx13qKioSPfcc4/27t2bUmORHiN1EVw5qJJgclmqfxd0xQBAcmyHkV27dqm+vl4NDQ06ePCgZs2apQULFuj06dPDHv/JJ59oyZIleuKJJ/TZZ59p0aJFWrRokf7xj3+MufFIryu/yXt1uulwwcOrfxcAkDWWTXPnzrXq6uoG/z8ajVrTpk2zmpqahj1+8eLF1sMPPzxkX0VFhfWLX/wi6dfs7e21JFm9vb12m4sUvfGGZc2YEd8mY8kSy/L54lunG+29zZhhWVJ8m8zxAICRJfv5basyMjAwoM7OTlVXVw/uy8/PV3V1tTo6OoY9p6OjY8jxkrRgwYIRj4cz2O1iSKxvsnt3aq+XyW6hxHPX1sa369ePXOkYbpwH3S0AkFm2wkhPT4+i0aj8fv+Q/X6/X5FIZNhzIpGIreMl6eLFi+rr6xvygLMtXiz5fPFtKjLZFZJ47t2741tp5IGlBA8AyD5HzqZpampSSUnJ4CMQCJhuEq4hsb7Jzp2pnZ/JmSeJ5168OL5tbCRwAICTjLNz8KRJk+Tz+dTd3T1kf3d3t6ZMmTLsOVOmTLF1vCStW7dO9fX1g//f19dHIMlxwWDmwkEmnxsAMHa2KiMFBQUqLy9Xe3v74L5YLKb29nZVVlYOe05lZeWQ4yXpL3/5y4jHS1JhYaEmTJgw5AEAAHKTrcqIJNXX12v58uWaM2eO5s6dq61bt+r8+fNasWKFJGnZsmWaPn26mpqaJEmrV6/WT37yE7366qt6+OGH1draqk8//VRvvvlmet8JAABwJdthpKamRmfOnNGmTZsUiUQ0e/ZstbW1DQ5SPXnypPLzLxdcqqqqtHPnTm3YsEHPPfecfvzjH+v999/X3Xffnb53AQAAXIt70wAAgIzg3jQAAMAVCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAo2wvB29CYpHYvr4+wy0BAADJSnxuX2uxd1eEkf7+fklSIBAw3BIAAGBXf3+/SkpKRvxzV9ybJhaL6euvv1ZxcbHy8vJMN8ex+vr6FAgEdOrUKe7h43BcK3fgOrkH18qZLMtSf3+/pk2bNuQmuldzRWUkPz9fN998s+lmuMaECRP4ZXQJrpU7cJ3cg2vlPKNVRBIYwAoAAIwijAAAAKMIIzmksLBQDQ0NKiwsNN0UXAPXyh24Tu7BtXI3VwxgBQAAuYvKCAAAMIowAgAAjCKMAAAAowgjAADAKMKIy7S0tKisrExFRUWqqKjQgQMHRjx2+/btmj9/viZOnKiJEyequrp61OORXnau1ZVaW1uVl5enRYsWZbaBkGT/Op07d051dXWaOnWqCgsLddttt2nv3r1Zaq232b1WW7du1e23367rr79egUBAa9as0bfffpul1sIWC67R2tpqFRQUWDt27LD++c9/Wk899ZR1ww03WN3d3cMeX1tba7W0tFifffaZdejQIevxxx+3SkpKrP/85z9Zbrn32L1WCV9++aU1ffp0a/78+dbPf/7z7DTWw+xep4sXL1pz5syxHnroIevjjz+2vvzyS2v//v1WV1dXllvuPXav1TvvvGMVFhZa77zzjvXll19a+/bts6ZOnWqtWbMmyy1HMggjLjJ37lyrrq5u8P+j0ag1bdo0q6mpKanzL126ZBUXF1tvv/12ppqI76VyrS5dumRVVVVZv/3tb63ly5cTRrLA7nV64403rFtuucUaGBjIVhPxPbvXqq6uzvrpT386ZF99fb01b968jLYTqaGbxiUGBgbU2dmp6urqwX35+fmqrq5WR0dHUs9x4cIFfffddyotLc1UM6HUr9WLL76oyZMn64knnshGMz0vlev0wQcfqLKyUnV1dfL7/br77ru1ZcsWRaPRbDXbk1K5VlVVVers7Bzsyjl27Jj27t2rhx56KCtthj2uuFEepJ6eHkWjUfn9/iH7/X6/Dh8+nNRzrF27VtOmTRvyC430S+Vaffzxx3rrrbfU1dWVhRZCSu06HTt2TH/961/12GOPae/evTp69Kh++ctf6rvvvlNDQ0M2mu1JqVyr2tpa9fT06P7775dlWbp06ZKCwaCee+65bDQZNlEZ8Yjm5ma1trbqvffeU1FRkenm4Ar9/f1aunSptm/frkmTJpluDkYRi8U0efJkvfnmmyovL1dNTY3Wr1+vcDhsumm4yv79+7Vlyxa9/vrrOnjwoN599119+OGH2rx5s+mmYRhURlxi0qRJ8vl86u7uHrK/u7tbU6ZMGfXcV155Rc3Nzfroo480c+bMTDYTsn+t/v3vf+v48eNauHDh4L5YLCZJGjdunI4cOaJbb701s432oFR+p6ZOnarrrrtOPp9vcN+dd96pSCSigYEBFRQUZLTNXpXKtdq4caOWLl2qJ598UpJ0zz336Pz583r66ae1fv165efzXdxJuBouUVBQoPLycrW3tw/ui8Viam9vV2Vl5Yjnvfzyy9q8ebPa2to0Z86cbDTV8+xeqzvuuEOff/65urq6Bh+PPPKIHnjgAXV1dSkQCGSz+Z6Ryu/UvHnzdPTo0cGwKElffPGFpk6dShDJoFSu1YULF34QOBIh0uKWbM5jegQtktfa2moVFhZav/vd76x//etf1tNPP23dcMMNViQSsSzLspYuXWqFQqHB45ubm62CggJrz5491jfffDP46O/vN/UWPMPutboas2myw+51OnnypFVcXGytWrXKOnLkiPWnP/3Jmjx5svXSSy+ZegueYfdaNTQ0WMXFxdbvf/9769ixY9af//xn69Zbb7UWL15s6i1gFHTTuEhNTY3OnDmjTZs2KRKJaPbs2Wpraxsc1HXy5Mkh3wTeeOMNDQwM6NFHHx3yPA0NDXr++eez2XTPsXutYIbd6xQIBLRv3z6tWbNGM2fO1PTp07V69WqtXbvW1FvwDLvXasOGDcrLy9OGDRv01Vdf6aabbtLChQvV2Nho6i1gFHmWRb0KAACYw1czAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUf8PCA3MRyIqlVIAAAAASUVORK5CYII=", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAi00lEQVR4nO3dbWxUZd7H8V87ta2wtG5Ba4HeQ3VBUaOEcoMMS8y62g0aDZvdtG43oi4mzu2yPKmbdrsRITatbjCuQqkrotmEsgVXN77oqn0FBfaJLiRGSDSCFJZW0hLbKlqW6blfHKd26ANzpjNznTPz/SQnxx7PdK7htJ3f/K+Hk2FZliUAAABDMk03AAAApDfCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjskw3IBqDg4M6c+aMpkyZooyMDNPNAQAAUbAsS/39/Zo+fboyM8euf3gijJw5c0bFxcWmmwEAAGJw6tQpzZw5c8z/74kwMmXKFEn2i8nLyzPcGgAAEI2+vj4VFxcPvY+PxRNhJNw1k5eXRxgBAMBjLjfEggGsAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAABPaGyUZs2y90gthBEAgCfU10snT9p7pBbCCADAE6qqJL/f3icaVZjkyrAsyzLdiMvp6+tTfn6+ent7lZeXZ7o5AIAUN2uWXYXx+6VPPzXdGu+K9v2byggAAJdIZhUGhBEAQBykWrdGMGhXRIJB0y1JD4QRAMCEMbgUE0EYAQBMGN0amAgGsAIAgIRgACsAAPAEwggAADCKMAIAAIwijABAmkq16bipIF2vCQNYASANNTZKq1ZJoRCrjLpJqq38ygBWAMCY6uvtIOLzMR3XTdJ1ijRhBADSUPhNb8sWVhl1k9FWfk2Hrhu6aQAAcDEvd93QTQMAQAq4tOsmFSslVEYAIEU1NtpjQ6qq6IpJJV6qlFAZAYA0x83rUlMqDnIljABAikrFNy2MPsjV6+imAQAACUE3DQAA8ATCCAAAMIowAgAAjCKMAAAAowgjAOAxqbjoFeLLaz8jzKYBAI/x0qJXMMMtPyPMpgGAFMX6Ibgcr/2MEEYAwCPCpXcp9Ra9Qnx5bWE0wggAeATLuyNWbh9DQhgBAI/wWukd7uH2IEsYAQCP8FrpHe7h9iBLGAEAF6uslLKy7D0QK7cHWcIIALjY7t1SKGTvgYly69gRwggAuFh5ueTz2Xtgotw6doQwAgAu1tQkXbxo74GJcuvYEVZgBQAACcEKrAAAwBMIIwDgMm4dZAgkCmEEAFzGrYMMgUQhjACAizQ2Sv39UkGB+wYZAolCGAEAF6mvl86dk6ZMce8CVUC8EUYAwEXcOvUSSCSm9gIAgIRgai8AAPAEwggAADCKMAIAQBpzw52hCSMAAKQxN9wZmjACAEAac8OdoWMKIw0NDSopKVFubq5KS0vV1tY27vk7d+7UbbfdpkmTJqmoqEiPPPKIenp6YmowAACIHzfcGdpxGGlubtbatWtVU1Ojw4cPa+nSpVq2bJk6OjpGPX///v1asWKFVq5cqQ8//FB79uzRv/71Lz366KMTbjwAAPA+x2HkhRde0MqVK/Xoo49q7ty5evHFF1VcXKxt27aNev7f//53zZo1S6tXr1ZJSYm+//3v67HHHtOhQ4cm3HgAAOB9jsLIhQsX1N7errKysojjZWVlOnjw4KiPCQQCOn36tFpaWmRZlj777DO9+eabuvfee8d8noGBAfX19UVsAAAgNTkKI93d3QqFQiosLIw4XlhYqK6urlEfEwgEtHPnTlVUVCg7O1vXXnutrrrqKr388stjPk9dXZ3y8/OHtuLiYifNBAAAHhLTANaMjIyIry3LGnEs7OjRo1q9erWefvpptbe3691339WJEycUHOcOUNXV1ert7R3aTp06FUszAQCAB2Q5OXnatGny+XwjqiBnz54dUS0Jq6ur05IlS/TUU09Jkm699VZNnjxZS5cu1bPPPquioqIRj8nJyVFOTo6TpgEAAI9yVBnJzs5WaWmpWltbI463trYqEAiM+pjz588rMzPyaXw+nyS7ogIAANKb426a9evXa/v27dqxY4eOHTumdevWqaOjY6jbpbq6WitWrBg6/7777tNbb72lbdu26fjx4zpw4IBWr16thQsXavr06fF7JQAAwJMcddNIUkVFhXp6erRp0yZ1dnbqlltuUUtLi/x+vySps7MzYs2Rhx9+WP39/dqyZYueeOIJXXXVVbrzzjv13HPPxe9VAAAAz8qwPNBX0tfXp/z8fPX29iovL890cwAASBmNjVJ9vVRVJY0ztyQm0b5/c28aAHCJxkZp1ix7DyRLfb108qS9N4UwAgAu4YY3BaSfqirJ77f3pjgeMwIASIxAQDp92t4DyRIMxr97xikqIwDgEgcPSqGQvQfSCWEEAFzCDeVypBe3jFNiNg0AAGlq1ix7nJLfL336afy/P7NpAADAuNxSjSOMAIDLuKV0jtQXDNoVEQawAgAihKf4rlpFIEF6IIwAgMtUVUk+nz2zhjVHkA4IIwDgMsGgtGWLO/rykXrc2A3IbBoAANLI1KnSuXNSQYHU05PY52I2DQAA8ATCCAAAaaS21u4CrK013ZJvEUYAwOXc2McP73LLdN7hGDMCAC6XzD5+IJ4YMwIAADxRWSOMAIDLubGPH94RXkTPzWvWEEYAwOXc2McP73DL/WfGw5gRAACQEIwZAQAAnkAYAQAARhFGAACAUYQRAABgFGEEAFKAF9aSAMZCGAGAFOCFtSSAsRBGACAFeGEtCcRfqlTEWGcEAFJMY6NdIamqYqG0VDdrll0R8/vthfHchnVGACBNhbtsnnhCysqSKitNtwiJkioVMSojAJBiwpWRU6ekwUHJ55MuXjTdKqQjKiMAkKbC97KpqLCDSHm56RYB4yOMAECKamqyKyJNTaZbgolIlUGq4yGMAADgYukwbZswAgBpKh0+caeCVBmkOh4GsAJAmnL7tFB4HwNYAQDjCn/iDgSokMAsKiMAkOaokCBRqIwAAKKSDmMS4G6EEQBIc+F1ScZbOp7BrhPDv9/46KYBAFxWuCunoECaMoX73jiVrl1hdNMAAOIm3JUjpf6aF4lAV9j4CCMAgMsKd+XU1o58U6UL4vKi6QpLZ3TTAAAmJF27IHB5dNMAAJJivC6IVKmapMrrcCvCCABgQsbrgojmvipufaMf3q50uD+MSYQRAEDCRDNw08QbfThoVFZKU6fa26VhaHi7GICaWIwZAQAYFa48JHO6cHici88nhUL2sUvHvJhoV6qJ9v2bMAIASDvhoBEISO+9Zx+rrSV0xBthBAAAGMVsGgAA4AmEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABgVUxhpaGhQSUmJcnNzVVpaqra2tnHPHxgYUE1Njfx+v3JycnT99ddrx44dMTUYAACkliynD2hubtbatWvV0NCgJUuW6JVXXtGyZct09OhR/c///M+ojykvL9dnn32m1157Td/73vd09uxZXbx4ccKNBwAA3uf4rr2LFi3S/PnztW3btqFjc+fO1fLly1VXVzfi/HfffVcPPPCAjh8/roKCgpgayV17AQDwnoTctffChQtqb29XWVlZxPGysjIdPHhw1Me88847WrBggZ5//nnNmDFDc+bM0ZNPPqmvvvpqzOcZGBhQX19fxAYAAFKTo26a7u5uhUIhFRYWRhwvLCxUV1fXqI85fvy49u/fr9zcXL399tvq7u7W448/rnPnzo05bqSurk4bN2500jQAAOBRMQ1gzcjIiPjasqwRx8IGBweVkZGhnTt3auHChbrnnnv0wgsv6I033hizOlJdXa3e3t6h7dSpU7E0EwAAeICjysi0adPk8/lGVEHOnj07oloSVlRUpBkzZig/P3/o2Ny5c2VZlk6fPq3Zs2ePeExOTo5ycnKcNA0AAHiUo8pIdna2SktL1draGnG8tbVVgUBg1McsWbJEZ86c0RdffDF07KOPPlJmZqZmzpwZQ5MBAEAqcdxNs379em3fvl07duzQsWPHtG7dOnV0dCgYDEqyu1hWrFgxdH5lZaWmTp2qRx55REePHtW+ffv01FNP6Re/+IWuvPLK+L0SAADgSY7XGamoqFBPT482bdqkzs5O3XLLLWppaZHf75ckdXZ2qqOjY+j873znO2ptbdWvfvUrLViwQFOnTlV5ebmeffbZ+L0KAADgWY7XGTGBdUYAAPCehKwzAgAAEG+EEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEbFFEYaGhpUUlKi3NxclZaWqq2tLarHHThwQFlZWZo3b14sTwsAAFKQ4zDS3NystWvXqqamRocPH9bSpUu1bNkydXR0jPu43t5erVixQj/84Q9jbiwAAEg9GZZlWU4esGjRIs2fP1/btm0bOjZ37lwtX75cdXV1Yz7ugQce0OzZs+Xz+fSXv/xFR44cifo5+/r6lJ+fr97eXuXl5TlpLgAAMCTa929HlZELFy6ovb1dZWVlEcfLysp08ODBMR/3+uuv65NPPtGGDRuiep6BgQH19fVFbEg9lZVSVpa9H01jozRrlr2P5jgAwJschZHu7m6FQiEVFhZGHC8sLFRXV9eoj/n4449VVVWlnTt3KisrK6rnqaurU35+/tBWXFzspJnwiN27pVDI3o+mvl46edLeR3McAOBNMQ1gzcjIiPjasqwRxyQpFAqpsrJSGzdu1Jw5c6L+/tXV1ert7R3aTp06FUszkWCXViicVizKyyWfz96PpqpK8vvtfTTHJ6qxUZo61d6ougBA8jgaM3LhwgVNmjRJe/bs0Y9//OOh42vWrNGRI0e0d+/eiPM///xzffe735XP5xs6Njg4KMuy5PP59P777+vOO++87PMyZsS8xka7ElFVJQWD9rFZs+wKhd8vffrpyK+9Jtx+ybuvAQDcJCFjRrKzs1VaWqrW1taI462trQoEAiPOz8vL0wcffKAjR44MbcFgUDfccIOOHDmiRYsWOXl6GDRa18ilFYpEVSySpapKKiiwt7FeA+NVACD+HM+maW5u1oMPPqjGxkYtXrxYf/jDH/Tqq6/qww8/lN/vV3V1tf7zn//oj3/846iPf+aZZ5hN40GjVUbSkderPwCQTAmpjEhSRUWFXnzxRW3atEnz5s3Tvn371NLSIr/fL0nq7Oy87JojMCc8LmLyZGdjI4JB+803nYOIdPnqD5UTAHDOcWXEBCoj8TN8XITEJ/x4o3ICAN9KWGUE3hYeFzFp0vhjIxAbr4+bAQATqIwAAICEoDICAAA8gTACGMZiawDSHWEEMKy+Xjp3zt5Y4h5AOiKMAIZFs9gaAKQyBrACAICEYAArkGJYUA1AqiKMuBBvOhjNaPcHAoBUQBhxkXAIqanhTQcjDV9QjcAKIJUwZsRFwkuJFxRIU6ZwUzqMjWXnAXgBY0Y8JPwpNxCw31xqa7kpHcZHlQRAKqEy4gJ8ysVEhH9+fD5pyxZCLAD3oDLiIdxcDRNRVWUHkVCIcUYAvIkw4gLBIN0yiF0waFdEhgdaum4AeAlhJEl4c0AiXRpow9OA/+//pMpKo00DgMsijCQJa0QgmYZ3+e3eba4dABANwkiSMC4EyRQMSj/7mT2WpLzcdGsAYHzMpgHSRGOjXZlj/RoAycJsGgAR6CoE4FaEESBNsFAaALeimwZIQyy0ByAZ6KYBMKZwlSQQoEICwDwqI0AamzpVOndOysyUtm5lYCuA+KIyAiBqg4MMbAVgDmEESGO1tVJBgb2xBg4AUwgjMWAmAlJFMCj19NhbMMjPNgAzGDMSg3A/e0GB/UccSBXMsgEQT4wZAeAYty0AYAJhJAa1tfYf7Npa0y0B4uvSu/8CQDIQRqI0vC+dP9gAAMQPYSRK3NcDAIDEIIxEib50pCtm2ABINGbTABgXM2wAxIrZNADigqoggEQjjAAYV3jAtkR3DYDEIIyMobJS8vmkyZP54wtIDOIGkDiEkTHs3m3fPOz8ef74ApLdTVNQIPX3E9ABxBdhZAzl5fZt1SdNoq8ckOzumilT7FshENABxBNhZAxNTVIoJH35JYubAWEMZgWQCEztBQAACcHUXgAA4AmEEQAAYBRhBAAAGEUY+Qb33wAAwAzCyDdY0AkAADMII7KrIf399oJOTFkEokdFEUA8MLVX3JUUiBW/OwDGw9ReBwIB+z40gYDplgDewiJoAOIh7SsjjY3SqlX2aqt8ugMAIH6ojESpvt4OIj4fn+4AADAh7cNIuMy8ZQv3oAEAwIS076YBAACJQTcNAADwBMIIgAmrrLTHXU2ezJojAJwjjACYsN27pcFB6fx5VjEG4BxhBMCElZdLmZnSpEnMSgPgHANYAQBAQjCAFQAAeAJhBAAAGBVTGGloaFBJSYlyc3NVWlqqtra2Mc996623dPfdd+vqq69WXl6eFi9erPfeey/mBgMAgNTiOIw0Nzdr7dq1qqmp0eHDh7V06VItW7ZMHR0do56/b98+3X333WppaVF7e7t+8IMf6L777tPhw4cn3HgA7tHYaN/Fl6m9AJxyPIB10aJFmj9/vrZt2zZ0bO7cuVq+fLnq6uqi+h4333yzKioq9PTTT0d1PgNYAfebNUs6eZIbTgL4VkIGsF64cEHt7e0qKyuLOF5WVqaDBw9G9T0GBwfV39+vgoKCMc8ZGBhQX19fxAbA3cL3eWJqLwCnHIWR7u5uhUIhFRYWRhwvLCxUV1dXVN9j8+bN+vLLL1VeXj7mOXV1dcrPzx/aiouLnTQTgAHBoF0R4YaTAJyKaQBrRkZGxNeWZY04Nppdu3bpmWeeUXNzs6655poxz6uurlZvb+/QdurUqViaCQAAPCDLycnTpk2Tz+cbUQU5e/bsiGrJpZqbm7Vy5Urt2bNHd91117jn5uTkKCcnx0nTAACARzmqjGRnZ6u0tFStra0Rx1tbWxUIBMZ83K5du/Twww+rqalJ9957b2wtBQAAKclxN8369eu1fft27dixQ8eOHdO6devU0dGh4DcdxdXV1VqxYsXQ+bt27dKKFSu0efNm3X777erq6lJXV5d6e3vj9yoAGMfUXgCxchxGKioq9OKLL2rTpk2aN2+e9u3bp5aWFvn9fklSZ2dnxJojr7zyii5evKhf/vKXKioqGtrWrFkTv1cBwLj6entqL3ftBeAUN8oDEBeNjXYQqapiRg0AW7Tv34QRAACQENy1F0DSMF4EwESkdRhpbJSmTrU3/ogCsXviCXu8yBNPmG4JAC9K6zBSXy+dO2dvDLoDYvf115F7AHAircNIVZVUUGBv3E8DiF1FheTz2XsAcMrRCqypJjzin6oIMDFNTfYGALFI68qIxNoIAACYlvZhhNueAwBgVtqHEW57DkwM03oBTFTahxEAE0NXJ4CJIowAmBC6OgFMFMvBAwCAhGA5+BjQ9w0AQPJRGRlm6lR7NdaCAqmnJ2FPAwBAWqAyAgAAPIEwMkxtrT0Qr7bWdEsAAEgfhJFhWHMEuLzKSikry94DQDwQRgA4snu3FArZewCIB8IIAEfKy+079JaXm24JgFTBbBoAAJAQzKYBAACeQBi5jP/9Xykjw94DAID4I4xcxqFD3+5ZmRXpiJWJASQaYeQyFiz49r+5KynSEXflBZBohJHL+Ne/pG3b7CXi+/v5dIj00tho/9wXFHBXXgCJw2yaKM2aZX869PvthdGAdMDPPYCJYDZNnFVV2X+Q+XSIdMLPPYBkoDICAAASgspIgjHDAACA+CCMxIgZBgAAxAdhJEb0pQMAEB+MGQEAAAnBmBEAAOAJhBEAAGAUYSTOmGUDAIAzhJE4C8+yWbWKQAIAQDQII3FWVSX5fFIoxLRfuA+VOwBuRBiJs2BQ2rKFab9wJ9bHAeBGhJEECAbtm4oFg3wShbuwPg4AN2KdkQQL3/W0oECaMsV+EwgGTbcKAIDEY50Rlwh/Ev36azuU1NSYbhEAAO5CGEmwcJdNbq7plgAA4E6EkSSprbUrJLW1plsCAIC7EEaSZPigViDeGCgNwMsII0AKYMouAC8jjAApgCm7ALyMMOJClNzhFN2AALyMMOJClNwBAOmEMOJCw0vuVEkAAKmOFVhdLryCq99vl+EBAPAKVmBNEZcOTKRSAgBINVRGPIZ73QAAvILKSIoKV0okBrmmCqpdANIdYcRjwlM4w8vLX7quRGWllJVl7+FewwMIs6cApDu6aVJMVpYUCkmZmVJxMd04bhMOH/390rlz3wbK+nquFYDUE+37N2EkxVRWSrt3Szk50vnzzMJxG8b8AEgnjBlJU01N0sWL0ubNYy8PzhiF5Bjt3zk85qe2lhVTASCMykgaYu2S5ODfGUC6ozKCMY11UzUqJvHFzesAIDoxhZGGhgaVlJQoNzdXpaWlamtrG/f8vXv3qrS0VLm5ubruuuvUyLudUWPdVG34rA6Cybdi/bfg5nUAEB3HYaS5uVlr165VTU2NDh8+rKVLl2rZsmXq6OgY9fwTJ07onnvu0dKlS3X48GH95je/0erVq/XnP/95wo1HfA3/JJ+u001HCx7p+m8BAEljObRw4UIrGAxGHLvxxhutqqqqUc//9a9/bd14440Rxx577DHr9ttvj/o5e3t7LUlWb2+v0+YiRtu2WZbfb++j8bOfWZbPZ+/dbrzX5vdblmTvozkfADC2aN+/HVVGLly4oPb2dpWVlUUcLysr08GDB0d9zN/+9rcR5//oRz/SoUOH9N///tdRcELyOO1i2L3bXt9k9+7Yni+R3ULDv3djo7Rq1diVjtHGedDdAgCJ5SiMdHd3KxQKqbCwMOJ4YWGhurq6Rn1MV1fXqOdfvHhR3d3doz5mYGBAfX19ERvcrbxc8vnsfSwS2RUy/HvX19uhyecbfWApwQMAki+mAawZGRkRX1uWNeLY5c4f7XhYXV2d8vPzh7bi4uJYmokkCq9v0tQU2+MTOfNk+PcO//eWLQQOAHCLLCcnT5s2TT6fb0QV5OzZsyOqH2HXXnvtqOdnZWVp6tSpoz6murpa69evH/q6r6+PQJLigsHEhYNLvzchBADcxVFlJDs7W6WlpWptbY043traqkAgMOpjFi9ePOL8999/XwsWLNAVV1wx6mNycnKUl5cXsQEAgNTkuJtm/fr12r59u3bs2KFjx45p3bp16ujoUPCbj5vV1dVasWLF0PnBYFAnT57U+vXrdezYMe3YsUOvvfaannzyyfi9CgAA4FmOumkkqaKiQj09Pdq0aZM6Ozt1yy23qKWlRX6/X5LU2dkZseZISUmJWlpatG7dOm3dulXTp0/XSy+9pJ/85CfxexUAAMCzuDcNAABICO5NAwAAPIEwAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDK8XLwJoQXie3r6zPcEgAAEK3w+/blFnv3RBjp7++XJBUXFxtuCQAAcKq/v1/5+flj/n9P3JtmcHBQZ86c0ZQpU5SRkWG6Oa7W19en4uJinTp1ivv4uBjXyTu4Vt7AdXIny7LU39+v6dOnKzNz7JEhnqiMZGZmaubMmaab4Sl5eXn8QnoA18k7uFbewHVyn/EqImEMYAUAAEYRRgAAgFGEkRSTk5OjDRs2KCcnx3RTMA6uk3dwrbyB6+RtnhjACgAAUheVEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRj2loaFBJSYlyc3NVWlqqtra2Mc996623dPfdd+vqq69WXl6eFi9erPfeey+JrU1vTq7VcAcOHFBWVpbmzZuX2AZCkvPrNDAwoJqaGvn9fuXk5Oj666/Xjh07ktTa9Ob0Wu3cuVO33XabJk2apKKiIj3yyCPq6elJUmvhiAXP+NOf/mRdccUV1quvvmodPXrUWrNmjTV58mTr5MmTo56/Zs0a67nnnrP++c9/Wh999JFVXV1tXXHFFda///3vJLc8/Ti9VmGff/65dd1111llZWXWbbfdlpzGprFYrtP9999vLVq0yGptbbVOnDhh/eMf/7AOHDiQxFanJ6fXqq2tzcrMzLR+//vfW8ePH7fa2tqsm2++2Vq+fHmSW45oEEY8ZOHChVYwGIw4duONN1pVVVVRf4+bbrrJ2rhxY7ybhkvEeq0qKiqs3/72t9aGDRsII0ng9Dr99a9/tfLz862enp5kNA/DOL1Wv/vd76zrrrsu4thLL71kzZw5M2FtROzopvGICxcuqL29XWVlZRHHy8rKdPDgwai+x+DgoPr7+1VQUJCIJuIbsV6r119/XZ988ok2bNiQ6CZCsV2nd955RwsWLNDzzz+vGTNmaM6cOXryySf11VdfJaPJaSuWaxUIBHT69Gm1tLTIsix99tlnevPNN3Xvvfcmo8lwyBM3yoPU3d2tUCikwsLCiOOFhYXq6uqK6nts3rxZX375pcrLyxPRRHwjlmv18ccfq6qqSm1tbcrK4tcyGWK5TsePH9f+/fuVm5urt99+W93d3Xr88cd17tw5xo0kUCzXKhAIaOfOnaqoqNDXX3+tixcv6v7779fLL7+cjCbDISojHpORkRHxtWVZI46NZteuXXrmmWfU3Nysa665JlHNwzDRXqtQKKTKykpt3LhRc+bMSVbz8A0nv1ODg4PKyMjQzp07tXDhQt1zzz164YUX9MYbb1AdSQIn1+ro0aNavXq1nn76abW3t+vdd9/ViRMnFAwGk9FUOMRHMI+YNm2afD7fiE8BZ8+eHfFp4VLNzc1auXKl9uzZo7vuuiuRzYScX6v+/n4dOnRIhw8f1qpVqyTZb3qWZSkrK0vvv/++7rzzzqS0PZ3E8jtVVFSkGTNmRNwSfe7cubIsS6dPn9bs2bMT2uZ0Fcu1qqur05IlS/TUU09Jkm699VZNnjxZS5cu1bPPPquioqKEtxvRozLiEdnZ2SotLVVra2vE8dbWVgUCgTEft2vXLj388MNqamqirzRJnF6rvLw8ffDBBzpy5MjQFgwGdcMNN+jIkSNatGhRspqeVmL5nVqyZInOnDmjL774YujYRx99pMzMTM2cOTOh7U1nsVyr8+fPKzMz8i3O5/NJsisqcBlzY2fhVHhq22uvvWYdPXrUWrt2rTV58mTr008/tSzLsqqqqqwHH3xw6PympiYrKyvL2rp1q9XZ2Tm0ff7556ZeQtpweq0uxWya5HB6nfr7+62ZM2daP/3pT60PP/zQ2rt3rzV79mzr0UcfNfUS0obTa/X6669bWVlZVkNDg/XJJ59Y+/fvtxYsWGAtXLjQ1EvAOAgjHrN161bL7/db2dnZ1vz58629e/cO/b+HHnrIuuOOO4a+vuOOOyxJI7aHHnoo+Q1PQ06u1aUII8nj9DodO3bMuuuuu6wrr7zSmjlzprV+/Xrr/PnzSW51enJ6rV566SXrpptusq688kqrqKjI+vnPf26dPn06ya1GNDIsi3oVAAAwhzEjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAo/4f3RH9TR4RiaIAAAAASUVORK5CYII=", "text/plain": [ "
" ] @@ -94,7 +94,7 @@ " self.sample_pts = self.sample_grid()\n", " self.hom_rank = hom_rank\n", " self.filt_layer = get_filtration\n", - " self.mpl = zz.MultiPers(hom_rank=hom_rank, l=l, res=res, ranks=list(range(1, 6)))\n", + " self.mpl = gril.MultiPers(hom_rank=hom_rank, l=l, step=2, res=res, ranks=list(range(1, 6)))\n", " self.mpl.set_max_jobs(40)\n", "\n", " def sample_grid(self):\n", @@ -123,8 +123,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "Took 58109 ms\n", - "Took 33645 ms\n" + "Took 19839 ms\n", + "Took 1842 ms\n" ] } ], @@ -134,11 +134,12 @@ "num_vertices = x.shape[0]\n", "\n", "layer1_0 = MultiPersLandscapeValLayer(res=0.01, hom_rank=0, step=2, l=2)\n", - "layer1_1 = MultiPersLandscapeValLayer(res=0.01, hom_rank=1, step=2, l=2)\n", + "# layer1_1 = MultiPersLandscapeValLayer(res=0.01, hom_rank=1, step=2, l=2)\n", "\n", "out1_0 = layer1_0(x, edges, triangles, tri_converted)\n", "lmbda1_0 = out1_0[0]\n", - "out1_1 = layer1_1(x, edges, triangles, tri_converted)\n", + "layer1_0.mpl.set_hom_rank(1)\n", + "out1_1 = layer1_0(x, edges, triangles, tri_converted)\n", "lmbda1_1 = out1_1[0]" ] }, @@ -149,7 +150,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAENCAYAAACSOWa6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAIU0lEQVR4nO3dT4jN/x7H8ff4lxp/JkX8yrAgl1koKxsLm7sj5YqkuDaUuxCl3N/v53fwy91IiZuNsrGwuKVsNHct9uqSsJC6SsiCjBzO3d3y58fQzHzPOa/HYzln6rxrmk/P+Zz3mTPQ6XQ6BQDEmtb0AABAs8QAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4M9JmxsbGaOXNmDQ0NVavVanocoEc4O7KJgT7Tbrfr4sWLtWzZsjpx4kQ9ffq06ZGAHuDsyCYG+sycOXNq9+7ddeTIkep0OnX79u2mRwJ6gLMjmxjoUyMjI1VVdffu3YYnAXqJsyOTGOhTY2NjVVV1586dhicBeomzI5MY6FOHDx+uKnUPfB9nRyYx0IcuX75cN2/erEWLFvmFBsbN2ZFroNPpdJoegonz6tWrWrVqVa1du7ZGRkbq9OnT9fz581qwYEHTowFdzNmRzc1Anzl16lQ9e/aszp49W2vWrKmqz1/7u3DhQq1bt65mzpzp/cRAVX377Hj79m3t3bu3hoeHa968ebV+/fq6detWU+MywcRAH3n48GGdOXOmDh8+XCtXrvzDreAlS5ZUq9WqrVu3NjEm0GXGc3a02+1avnx53bhxo16+fFkHDx6sTZs21atXr5oamwkkBvrIoUOHauHChfXzzz9XVdXq1aur6vObgS1bttTmzZtraGhoqkcEutB4zo7BwcE6duxYDQ8P17Rp02rHjh01a9asunfvXiMzM7FmND0AE2N0dLSuXbtWV65cqcHBwaqqmjt3bi1dutQiEPCHfvTsuH//fr148aJWrFgxVaMyidwM9IF3797VwYMHa+PGjbV9+/aPHluzZo0YAL7oR8+ON2/e1K5du+ro0aM1f/78qRiVSSYG+sD58+frwYMHde7cuc8eGxkZqcePH3tdD/jMj5wd7969q23bttWKFSvq2LFjUzUqk8xbC4Pt37+/Fi9e7B0FwLh8+PChdu7cWa9fv66rV6/WjBleae4XfpKB2u12tdvtev/+fbXb7f9/dOn06dObHg3oYvv27asnT57U6OioEOgzbgYCtVqtOn78+Edfu3TpUu3Zs6eZgYCu9+jRo1q+fHnNnj37oz8crl+/Xhs2bGhwMiaCGACAcBYIASCcGACAcGIAAMKJAQAIN+73hgwMtCZxDGA8Op1W0yN8N2cHNO9bZ4ebAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcGIAAMKJAQAIJwYAIJwYAIBwMyb9Gf7UmvSnAAB+nJsBAAgnBgAgnBgAgHBiAADCTf4C4af+NeXPCAB8hZsBAAgnBgAgnBgAgHBTvzPwBb+M/L3pEaBHnGp6AKAPuRkAgHBiAADCiQEACCcGACBcVywQfurkf//R9AjQnX6yQAhMPDcDABBODABAODEAAOG6cmfgi/7a9ADQBUabHqBBrVbTE0DfcjMAAOHEAACEEwMAEE4MAEC43lkg/ETr301PAFOv1fQAXeaX33ziKYzP1/9hmZsBAAgnBgAgnBgAgHA9uzMA8KmT//QhZ/BFB+wMAABfIQYAIJwYAIBwYgAAwvXsAmHrz01PAPSC1t+angCa1zrw9cfdDABAODEAAOHEAACE652dgUtNDwAA/cnNAACEEwMAEE4MAEA4MQAA4bpygfDXn442PQJ0pZNND9CDWuebngC6n5sBAAgnBgAgnBgAgHBdsTPw+39ONT0C9ISTI01P0N1+PWDfCL7kW/tGbgYAIJwYAIBwYgAAwokBAAg39QuEf5nyZ4T+cbfpAbrL78ctH8N4nPzt64+7GQCAcGIAAMKJAQAIN9DpdDrj+saB1iSPAnxLp9NqeoTv5uyA5n3r7HAzAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQb6HQ6naaHAACa42YAAMKJAQAIJwYAIJwYAIBwYgAAwokBAAgnBgAgnBgAgHBiAADC/Q+3jnnCsLimzQAAAABJRU5ErkJggg==", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAELCAYAAABEYIWnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAIGElEQVR4nO3dPWjVaRrG4ecYg0L8LCwcUQw24lSxEBstt/IjxWLnJ0TUrK2wknFOcJjFxsrFQoOFYGUf18JOO9FmWTFY2IhgIwhRScjZ2jGT6GDyPzn3dZU5gdyNLz9fX0yr0+l0CgCItarpAQBAs8QAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4M9JhPnz7V6tWra/369XXx4sWm5wArhLMjmxjoMa1Wqx49elT79u2rGzdu1MuXL5ueBKwAzo5sYqDHrFmzpg4ePFiXLl2qqqqnT582vAhYCZwd2cRAj9q9e3dVVT1//rzZIcCK4uzIJAZ61LVr16rKH2jg+zg7MomBHvTw4cO6efNmbd68uZ49e/bV5zdv3qy9e/dWf39/tdvt5R8IdKWFzo7Pnz/X6dOna/v27bVhw4bav39/PXnypKGl/GhioMe8f/++zpw5U0ePHq0LFy7Uu3fv6s2bN198z9atW2t8fLyGh4ebGQl0ncXOjtnZ2RocHKzHjx/X+/fv6/z583XkyJGanp5ucDU/ihjoMaOjozUzM1O3bt2qoaGhqvr6um94eLgOHz5cGzdubGAh0I0WOzsGBgbqypUrtWPHjlq1alWdPHmy5ubmampqqqHF/EhioIfcv3+/7t27VxMTE7Vly5bau3dvVfm3P2Bhf+XsePHiRX38+LF27dq1TCtZSmKgR7x9+7bOnTtXZ8+erUOHDlVV1eDgYG3atGnedwMAVX/t7Jienq7jx4/X2NhYrVu3bjnnskTEQI8YGRmpzZs31/Xr17/4+tDQkJsB4E9979kxMzNTx44dqz179tTly5eXaSVLTQz0gNu3b9fk5GTdvXu3BgYGvvhsaGioXr16VR8+fGhoHdCtvvfsmJubqxMnTlRfX19NTExUq9Va7skskVan0+k0PYLlNTs7W7Ozs3X+/Pnatm1bjY2NVX9/f/X19TU9DehiIyMjNTU1VQ8ePKi1a9c2PYcfSAwEarfbNT4+/sXX7ty5U6dOnWpmEND1Xr9+XTt37qy1a9d+8ReHycnJOnDgQIPL+BHEAACE82YAAMKJAQAIJwYAIJwYAIBwYgAAwq3+1m9stdpLOAP4Fp1Ou+kJ383ZAc1b7OxwMwAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQbvWy/8Td7WX/kQDAn3MzAADhxAAAhBMDABBODABAuOV/QDif+00PAIBcbgYAIJwYAIBwYgAAwnXHm4E/GPv5ctMToEv93vQAoAe5GQCAcGIAAMKJAQAIJwYAIFxXPiCcz9U3/2p6AjTvJw8IgR/PzQAAhBMDABBODABAuBXzZuArp5seAA34T9MDuky73fQC6AluBgAgnBgAgHBiAADCiQEACLdyHxDOo/2w6QWwtNpND+hyY7/6jacwv4X/wzI3AwAQTgwAQDgxAADheurNAMAfXf23X3IGNerNAACwADEAAOHEAACEEwMAEK6nHhC2/9b0AqDbtf/R9AJYfu3RhT93MwAA4cQAAIQTAwAQbuW+GbjT9AAA6A1uBgAgnBgAgHBiAADCiQEACLdiHhD+8tM/m54Ajbva9IAe0L7R9ALoPm4GACCcGACAcGIAAMJ15ZuB3/77e9MToCtd/bnpBSvPL6PeG8Fi743cDABAODEAAOHEAACEEwMAEK47HhD+vekBsEL8r+kB3e23cY+PYT5Xf134czcDABBODABAODEAAOFanU6n803f2Gov8RRgMZ1Ou+kJ383ZAc1b7OxwMwAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEEwMAEE4MAEA4MQAA4cQAAIQTAwAQTgwAQDgxAADhxAAAhBMDABBODABAODEAAOHEAACEa3U6nU7TIwCA5rgZAIBwYgAAwokBAAgnBgAgnBgAgHBiAADCiQEACCcGACCcGACAcP8HzKpHit+qcqkAAAAASUVORK5CYII=", "text/plain": [ "
" ] diff --git a/gril/__init__.py b/gril/__init__.py new file mode 100644 index 0000000..c5f4929 --- /dev/null +++ b/gril/__init__.py @@ -0,0 +1 @@ +from .gril import * diff --git a/gril/__pycache__/__init__.cpython-39.pyc b/gril/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000..b9ecfb0 Binary files /dev/null and b/gril/__pycache__/__init__.cpython-39.pyc differ diff --git a/gril/__pycache__/gril.cpython-39.pyc b/gril/__pycache__/gril.cpython-39.pyc new file mode 100644 index 0000000..69ffd3b Binary files /dev/null and b/gril/__pycache__/gril.cpython-39.pyc differ diff --git a/gril/__pycache__/zz.cpython-39.pyc b/gril/__pycache__/zz.cpython-39.pyc new file mode 100644 index 0000000..b90658f Binary files /dev/null and b/gril/__pycache__/zz.cpython-39.pyc differ diff --git a/gril/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so b/gril/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so new file mode 100755 index 0000000..731f4ca Binary files /dev/null and b/gril/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so differ diff --git a/gril/build/temp.linux-x86_64-cpython-39/.ninja_deps b/gril/build/temp.linux-x86_64-cpython-39/.ninja_deps new file mode 100644 index 0000000..d2679b6 Binary files /dev/null and b/gril/build/temp.linux-x86_64-cpython-39/.ninja_deps differ diff --git a/gril/build/temp.linux-x86_64-cpython-39/.ninja_log b/gril/build/temp.linux-x86_64-cpython-39/.ninja_log new file mode 100644 index 0000000..259d257 --- /dev/null +++ b/gril/build/temp.linux-x86_64-cpython-39/.ninja_log @@ -0,0 +1,15 @@ +# ninja log v5 +0 19238 1674245376000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +274 16567 1675726632000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +269 16888 1675728138000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +26 24273 1675902702086332969 /u/scratch1/ssamaga/mpml_altered_bdry/zigzag/build/temp.linux-x86_64-cpython-39/u/scratch1/ssamaga/mpml_altered_bdry/zigzag/multipers.o 72e776d2b1d272f3 +5160 33366 1675906707159817794 /u/scratch1/ssamaga/mpml_altered_bdry/zigzag/build/temp.linux-x86_64-cpython-39/u/scratch1/ssamaga/mpml_altered_bdry/zigzag/multipers.o 72e776d2b1d272f3 +0 17526 1675907214000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +31 17395 1675907889000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +257 17692 1675908048000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +308 16897 1675909436000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +231 16397 1675909668000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +257 16577 1675911626000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +290 17660 1675912935000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +44 17018 1679793102000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a +13 33166 1681688576281989022 /scratch1/mukher26/mpml_graph_repo/gril/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/gril/multipers.o 9d57bdffb2f12246 diff --git a/zigzag/build/temp.linux-x86_64-cpython-39/build.ninja b/gril/build/temp.linux-x86_64-cpython-39/build.ninja similarity index 85% rename from zigzag/build/temp.linux-x86_64-cpython-39/build.ninja rename to gril/build/temp.linux-x86_64-cpython-39/build.ninja index 35e5982..0bc18eb 100644 --- a/zigzag/build/temp.linux-x86_64-cpython-39/build.ninja +++ b/gril/build/temp.linux-x86_64-cpython-39/build.ninja @@ -12,7 +12,7 @@ rule compile -build /scratch1/mukher26/mpml_graph_repo/zigzag/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/zigzag/./multipers.o: compile /scratch1/mukher26/mpml_graph_repo/zigzag/multipers.cpp +build /scratch1/mukher26/mpml_graph_repo/gril/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/gril/./multipers.o: compile /scratch1/mukher26/mpml_graph_repo/gril/multipers.cpp diff --git a/gril/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o b/gril/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o new file mode 100644 index 0000000..9527c99 Binary files /dev/null and b/gril/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o differ diff --git a/gril/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/gril/multipers.o b/gril/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/gril/multipers.o new file mode 100644 index 0000000..ea85813 Binary files /dev/null and b/gril/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/gril/multipers.o differ diff --git a/gril/build/temp.linux-x86_64-cpython-39/u/scratch1/ssamaga/mpml_altered_bdry/zigzag/multipers.o b/gril/build/temp.linux-x86_64-cpython-39/u/scratch1/ssamaga/mpml_altered_bdry/zigzag/multipers.o new file mode 100644 index 0000000..0a2c5cc Binary files /dev/null and b/gril/build/temp.linux-x86_64-cpython-39/u/scratch1/ssamaga/mpml_altered_bdry/zigzag/multipers.o differ diff --git a/gril/dist/mpml-0.0.0-py3.9-linux-x86_64.egg b/gril/dist/mpml-0.0.0-py3.9-linux-x86_64.egg new file mode 100644 index 0000000..784851d Binary files /dev/null and b/gril/dist/mpml-0.0.0-py3.9-linux-x86_64.egg differ diff --git a/zigzag/zz.py b/gril/gril.py similarity index 85% rename from zigzag/zz.py rename to gril/gril.py index 84925a4..2c06543 100644 --- a/zigzag/zz.py +++ b/gril/gril.py @@ -35,7 +35,7 @@ class MultiPers: - def __init__(self, hom_rank: int, l: int, res: float, ranks: List[int]): + def __init__(self, hom_rank: int, l: int, res: float, step: int, ranks: List[int]): # try: # __M = load( # 'zigzag', @@ -45,8 +45,8 @@ def __init__(self, hom_rank: int, l: int, res: float, ranks: List[int]): # except Exception as ex: # print("Error was {}".format(ex)) - - self.mpl = mpml.Multipers(hom_rank, l, res, ranks) + # const int hom_rank, const int l, double res, int step, const std::vector ranks + self.mpl = mpml.Multipers(hom_rank, l, res, step, ranks) def compute_landscape(self, pts: List[Tuple[int]], batch: List[Tuple[Tensor, List[List[int]]]]): @@ -54,6 +54,12 @@ def compute_landscape(self, pts: List[Tuple[int]], batch: List[Tuple[Tensor, Lis def set_max_jobs(self, njobs: int): self.mpl.set_max_jobs(njobs) + + def set_hom_rank(self, hom_rank: int): + self.mpl.set_hom_rank(hom_rank) + + def refresh_rank_info(self): + self.mpl.refresh_rank_info() """ diff --git a/zigzag/mpml.egg-info/PKG-INFO b/gril/mpml.egg-info/PKG-INFO similarity index 100% rename from zigzag/mpml.egg-info/PKG-INFO rename to gril/mpml.egg-info/PKG-INFO diff --git a/zigzag/mpml.egg-info/SOURCES.txt b/gril/mpml.egg-info/SOURCES.txt similarity index 67% rename from zigzag/mpml.egg-info/SOURCES.txt rename to gril/mpml.egg-info/SOURCES.txt index 3f4f20c..7c7dc4a 100644 --- a/zigzag/mpml.egg-info/SOURCES.txt +++ b/gril/mpml.egg-info/SOURCES.txt @@ -1,5 +1,5 @@ setup.py -/scratch1/mukher26/mpml_graph_repo/zigzag/./multipers.cpp +/scratch1/mukher26/mpml_graph_repo/gril/./multipers.cpp mpml.egg-info/PKG-INFO mpml.egg-info/SOURCES.txt mpml.egg-info/dependency_links.txt diff --git a/zigzag/mpml.egg-info/dependency_links.txt b/gril/mpml.egg-info/dependency_links.txt similarity index 100% rename from zigzag/mpml.egg-info/dependency_links.txt rename to gril/mpml.egg-info/dependency_links.txt diff --git a/zigzag/mpml.egg-info/top_level.txt b/gril/mpml.egg-info/top_level.txt similarity index 100% rename from zigzag/mpml.egg-info/top_level.txt rename to gril/mpml.egg-info/top_level.txt diff --git a/zigzag/multipers.cpp b/gril/multipers.cpp similarity index 59% rename from zigzag/multipers.cpp rename to gril/multipers.cpp index 140da68..1e2a65d 100644 --- a/zigzag/multipers.cpp +++ b/gril/multipers.cpp @@ -17,64 +17,75 @@ Tensor Multipers::compute_l_worm(const int d){ assert (l > 0); auto tensopt = torch::TensorOptions().device(torch::kCPU); - // compute anchor points. There should be 4*l anchor points - - std::vector anchor_pts (4*l, std::make_pair(0, 0)); + // compute anchor points. There should be 4*(2l-1) anchor points + int num_anchor_pts = 8*l - 4; + std::vector anchor_pts (num_anchor_pts, std::make_pair(0, 0)); // bottom right: (4l-1)th anchor pt auto br = std::make_pair(l * d, - l * d); - anchor_pts[4*l-1] = br; - + anchor_pts[num_anchor_pts - 1] = br; - // Iteration over 2l-1 boxes to figure out 4l-2 anchor pts. After this + auto br_1 = std::make_pair(br.first - 2 * d, br.second); + auto br_2 = std::make_pair(br.first, br.second+ 2*d); + // Iteration over 2l-1 boxes to figure out 8l-6 anchor pts. After this // left will be 2l-1 th anchor pt - for (auto i = 0; i < 2*l - 1; i++){ - anchor_pts[i] = std::make_pair(br.first - 2 * d, br.second); - anchor_pts[4*l-2-i] = std::make_pair(br.first, br.second + 2 * d); - br = std::make_pair(br.first - d, br.second + d); + for (auto i = 0; i < 4*l - 3; i++){ + anchor_pts[i] = br_1; + anchor_pts[num_anchor_pts-2-i] = br_2; + br_1.first = br_1.first - (i % 2) * d; + br_1.second = br_1.second + ((i + 1) % 2) * d; + + br_2.first = br_2.first - ((i + 1) % 2) * d; + br_2.second = br_2.second + (i % 2) * d; + + // br_1 = std::make_pair(br_1.first - (i % 2) * d, br_1.second + ((i + 1) % 2) * d); + // br_2 = std::make_pair(br_2.first, br_2.second); } - anchor_pts[2*l-1] = std::make_pair(br.first - d, br.second + d); + anchor_pts[4*l-3] = std::make_pair(br_1.first, br_2.second); std::vector grid_pts; + for(auto i = 0; i < num_anchor_pts; i++){ + grid_pts.push_back(torch::tensor({anchor_pts[i].first, anchor_pts[i].second}, tensopt)); + } // Lower staircase - for (auto b = 0; b < 2*l-2; b++){ - auto cur_pt_x = anchor_pts[b].first; - auto cur_pt_y = anchor_pts[b].second; - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - for (auto i = 1; i < 2 * d; i++){ - cur_pt_x = cur_pt_x - ((i + 1) % 2); - cur_pt_y = cur_pt_y + (i % 2); - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - } - } - - // Box region - - for (auto b = 2*l-2; b < 2*l; b++){ - auto cur_pt_x = anchor_pts[b].first; - auto cur_pt_y = anchor_pts[b].second; - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - for (auto i = 1; i < 2 * d; i++){ - cur_pt_x = cur_pt_x + (b % 2); - cur_pt_y = cur_pt_y + ((b + 1) % 2); - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - } - } + // for (auto b = 0; b < 2*l-2; b++){ + // auto cur_pt_x = anchor_pts[b].first; + // auto cur_pt_y = anchor_pts[b].second; + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // for (auto i = 1; i < 2 * d; i++){ + // cur_pt_x = cur_pt_x - ((i + 1) % 2); + // cur_pt_y = cur_pt_y + (i % 2); + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // } + // } + + // // Box region + + // for (auto b = 2*l-2; b < 2*l; b++){ + // auto cur_pt_x = anchor_pts[b].first; + // auto cur_pt_y = anchor_pts[b].second; + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // for (auto i = 1; i < 2 * d; i++){ + // cur_pt_x = cur_pt_x + (b % 2); + // cur_pt_y = cur_pt_y + ((b + 1) % 2); + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // } + // } - // Upper staircase + // // Upper staircase - for (auto b = 2*l; b < 4*l -2; b++){ - auto cur_pt_x = anchor_pts[b].first; - auto cur_pt_y = anchor_pts[b].second; - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - for (auto i = 1; i < 2*d; i++){ - cur_pt_x = cur_pt_x + ((i + 1) % 2); - cur_pt_y = cur_pt_y - (i % 2); - grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); - } - } - - grid_pts.push_back(torch::tensor({anchor_pts[4*l-2].first, anchor_pts[4*l-2].second}, tensopt)); + // for (auto b = 2*l; b < 4*l -2; b++){ + // auto cur_pt_x = anchor_pts[b].first; + // auto cur_pt_y = anchor_pts[b].second; + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // for (auto i = 1; i < 2*d; i++){ + // cur_pt_x = cur_pt_x + ((i + 1) % 2); + // cur_pt_y = cur_pt_y - (i % 2); + // grid_pts.push_back(torch::tensor({cur_pt_x, cur_pt_y}, tensopt)); + // } + // } + + // grid_pts.push_back(torch::tensor({anchor_pts[4*l-2].first, anchor_pts[4*l-2].second}, tensopt)); auto grid_points_along_boundary_t = torch::stack(grid_pts, 0); return grid_points_along_boundary_t; @@ -83,6 +94,10 @@ Tensor Multipers::compute_l_worm(const int d){ std::vector> Multipers::compute_filtration_along_boundary_cap(const Tensor& grid_pts_along_boundary_t, const Tensor& f, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, int &manual_birth_pts, int &manual_death_pts){ //std::cout << " computing filtration along boundary " << std::endl; @@ -124,26 +139,47 @@ std::vector> Multipers::compute_filtration_along_bound // # Up arrow if (up_condition) { - temp = (y_new >= f.index({"...", 1})) & (x_old >= f.index({"...", 0})) & (y_old < f.index({"...", 1})); + // temp = (y_new >= f.index({"...", 1})) & (x_old >= f.index({"...", 0})) & (y_old < f.index({"...", 1})); + // ch = true; + + temp = (y_new >= f_y_sorted.index({"...", 1})) & (x_old >= f_y_sorted.index({"...", 0})) & (y_old < f_y_sorted.index({"...", 1})); + auto idx = torch::nonzero(temp); + temp = f_y_sorted_id.index({idx}); ch = true; } // # Right arrow else if (r_condition) { - temp = (x_new >= f.index({"...", 0})) & (y_old >= f.index({"...", 1})) & (x_old < f.index({"...", 0})); + // temp = (x_new >= f.index({"...", 0})) & (y_old >= f.index({"...", 1})) & (x_old < f.index({"...", 0})); + // ch = true; + + temp = (x_new >= f_x_sorted.index({"...", 0})) & (y_old >= f_x_sorted.index({"...", 1})) & (x_old < f_x_sorted.index({"...", 0})); + auto idx = torch::nonzero(temp); + temp = f_x_sorted_id.index({idx}); ch = true; } // # Left arrow else if (l_condition) { - temp = (x_new < f.index({"...", 0})) & (y_old >= f.index({"...", 1})) & (x_old >= f.index({"...", 0})); + // temp = (x_new < f.index({"...", 0})) & (y_old >= f.index({"...", 1})) & (x_old >= f.index({"...", 0})); + // ch = false; + + temp = (x_new < f_x_sorted.index({"...", 0})) & (y_old >= f_x_sorted.index({"...", 1})) & (x_old >= f_x_sorted.index({"...", 0})); + auto idx = torch::nonzero(temp); + temp = f_x_sorted_id.index({idx}); ch = false; } // # Down arrow else if (down_condition) { - temp = (y_new < f.index({"...", 1})) & (x_old >= f.index({"...", 0})) & (y_old >= f.index({"...", 1})); + // temp = (y_new < f.index({"...", 1})) & (x_old >= f.index({"...", 0})) & (y_old >= f.index({"...", 1})); + // ch = false; + + temp = (y_new < f_y_sorted.index({"...", 1})) & (x_old >= f_y_sorted.index({"...", 0})) & (y_old >= f_y_sorted.index({"...", 1})); + + auto idx = torch::nonzero(temp); + temp = f_y_sorted_id.index({idx}); ch = false; } - const Tensor ind = torch::where(temp)[0]; + const auto ind = temp; if (ch){ for (auto idx = 0; idx < ind.size(0); idx++){ auto simp = ind[idx].item(); @@ -177,10 +213,11 @@ std::vector> Multipers::compute_filtration_along_bound } -int Multipers::zigzag_pairs(std::vector> &simplices_birth_death, +void Multipers::zigzag_pairs(std::vector> &simplices_birth_death, const vector &simplices, const int manual_birth_pts, - const int manual_death_pts){ + const int manual_death_pts, + std::vector &num_full_bars){ std::vector orig_f_add_id; std::vector orig_f_del_id; @@ -217,6 +254,7 @@ int Multipers::zigzag_pairs(std::vector> &simplices_bi char op = std::get<0>(record); Integer simplex_id = std::get<1>(record); const auto& simp = simplices[simplex_id]; + // std::cout << "i: " << "simplex_id: " << simplex_id << simp << std::endl; if (op) { getBoundaryChainPhat(id_map, simp, bound_c); bound_chains.set_col(s_id, bound_c); @@ -271,7 +309,7 @@ int Multipers::zigzag_pairs(std::vector> &simplices_bi phat::persistence_pairs pairs; phat::compute_persistence_pairs< phat::twist_reduction >( pairs, bound_chains ); - int num_full_bars = 0; + // int num_full_bars = 0; for (phat::index idx = 0; idx < pairs.get_num_pairs(); idx++) { Integer b = pairs.get_pair(idx).first; Integer d = pairs.get_pair(idx).second - 1; @@ -283,16 +321,25 @@ int Multipers::zigzag_pairs(std::vector> &simplices_bi mapRelExtIntv(p, b, d, orig_f_add_id, orig_f_del_id, simp_num); } - if (b <= manual_birth_pts && d >= manual_death_pts && p == this->hom_rank) - num_full_bars++; + if (b <= manual_birth_pts && d >= manual_death_pts){ + if(p < 2) + num_full_bars[p] = num_full_bars[p] + 1; + } + } - return num_full_bars; - } -int Multipers::num_full_bars_for_specific_d(const Tensor& filtration, const vector& simplices, const Point& p, int d){ +void Multipers::num_full_bars_for_specific_d(const Tensor& filtration, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, + const vector& simplices, + const Point& p, + int d, + std::vector &num_full_bars){ int d1, d2; auto start = std::chrono::high_resolution_clock::now(); auto bd_cap = compute_l_worm(d); @@ -304,56 +351,91 @@ int Multipers::num_full_bars_for_specific_d(const Tensor& filtration, const vect // TODO: Needs bd_cap to be scaled and translated [DONE]. const auto shift = torch::tensor({p.first, p.second}, tensopt_real); - const auto scale = torch::tensor({step, step}, tensopt_real); + const auto scale = torch::tensor({res, res}, tensopt_real); const auto lower_left_corner = torch::tensor({ll_x, ll_y}, tensopt_real); bd_cap = ((bd_cap + shift) * scale) + lower_left_corner; start = std::chrono::high_resolution_clock::now(); - auto ret = compute_filtration_along_boundary_cap(bd_cap, filtration, d1, d2); + auto ret = compute_filtration_along_boundary_cap(bd_cap, + filtration, + f_x_sorted, + f_y_sorted, + f_x_sorted_id, + f_y_sorted_id, + d1, + d2); end = std::chrono::high_resolution_clock::now(); duration = std::chrono::duration_cast(end - start); if (ret.empty()){ - return 0; + // return 0; + return; } start = std::chrono::high_resolution_clock::now(); - int full_bars = zigzag_pairs(ret, simplices, d1, d2); + zigzag_pairs(ret, simplices, d1, d2, num_full_bars); end = std::chrono::high_resolution_clock::now(); - return full_bars; } -Tensor Multipers::find_maximal_worm_for_rank_k(const Tensor &filtration, const vector &simplices, const Point &p, const int rank) +Tensor Multipers::find_maximal_worm_for_rank_k(const Tensor &filtration, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, + const vector &simplices, + const Point &p, + const int rank, + std::vector*> rank_info) { auto tensopt_real = torch::TensorOptions().dtype(filtration.dtype()).device(filtration.device()); - int d_max = (int) (1.0 / step) + 1; + int d_max = (int) (1.0 / res) + 1; int d_min = 1; // std::cout << "Point: " << (p.first * grid_resolution_x) << " , " << (p.second * grid_resolution_x) << std::endl; int d = 1; - int res = -1; + int ans = -1; + // auto rank_info = (this->hom_rank == 0) ? rank_info_h0 : rank_info_h1; + int num_full_bars_for_this_k; while(d_min <= d_max){ - d = (d_min + d_max) / 2; - int num_full_bars = num_full_bars_for_specific_d(filtration, simplices, p, d); + d = (d_min + d_max) / 2; + if(rank_info[this->hom_rank]->count(d) == 0){ + std::vector num_full_bars = {0, 0}; + num_full_bars_for_specific_d(filtration, + f_x_sorted, + f_y_sorted, + f_x_sorted_id, + f_y_sorted_id, + simplices, + p, + d, + num_full_bars); + + (*rank_info[0])[d] = num_full_bars[0]; + (*rank_info[1])[d] = num_full_bars[1]; + num_full_bars_for_this_k = num_full_bars[this->hom_rank]; + } + else{ + num_full_bars_for_this_k = rank_info[this->hom_rank]->at(d); + } // std::cout <<" Rank = " << num_full_bars << " d = " << (d * grid_resolution_x) << " Query rank = "<< rank << std::endl; - if(num_full_bars == rank){ - res = d; + if(num_full_bars_for_this_k == rank){ + ans = d; d_min = d + 1; } - if (num_full_bars > rank) + if (num_full_bars_for_this_k > rank) d_min = d + 1; - else if (num_full_bars < rank) + else if (num_full_bars_for_this_k < rank) d_max = d - 1; } - if (res == -1) - res = d_max; + if (ans == -1) + ans = d_max; // auto scale = (grid_resolution_x * grid_resolution_x) + (grid_resolution_y * grid_resolution_y); // auto scaled_res = (double)res * std::sqrt(scale); - auto scaled_res = (double)res * step; + auto scaled_res = (double)ans * res; // std::cout <<"Landscape val " << scaled_res << " Query rank = "<< rank << std::endl; auto res_t = torch::tensor({scaled_res}, tensopt_real); @@ -396,11 +478,18 @@ std::vector Multipers::compute_landscape(const std::vector& pts, at::set_num_threads(this->max_threads); // int max_jobs = omp_get_max_threads(); // std::cout << "Max threads = " << max_jobs << std::endl; + const auto f_x_sorted_id = filtration.index({"...", 0}).argsort(); + const auto f_y_sorted_id = filtration.index({"...", 1}).argsort(); + const auto f_x_sorted = filtration.index({f_x_sorted_id}); + const auto f_y_sorted = filtration.index({f_y_sorted_id}); + // std::cout << "f_x: " << f_x_sorted << std::endl; + // std::cout << "f: " << filtration << std::endl; int j = 0; - #pragma omp parallel private(j) firstprivate(filtration) shared(num_points, num_ranks, simplices, pts, ranks, futures, hom_rank, std::cout) //default(none) shared(num_points, num_ranks, num_vertices, filtration, simplices, pts, ranks, futures, hom_rank, std::cout) + + #pragma omp parallel private(j) firstprivate(filtration, f_x_sorted, f_y_sorted, f_x_sorted_id, f_y_sorted_id) shared(num_points, num_ranks, simplices, pts, ranks, futures, hom_rank, std::cout) //default(none) shared(num_points, num_ranks, num_vertices, filtration, simplices, pts, ranks, futures, hom_rank, std::cout) { - #pragma omp for collapse(2) + #pragma omp for for(auto i = 0; i < num_points; i++) { for (j = 0; j < num_ranks; j++){ const auto &point = pts[i]; @@ -415,7 +504,8 @@ std::vector Multipers::compute_landscape(const std::vector& pts, // std::cout << "flattened k " << k << std::endl; // std::cout << "hom_rank " << hom_rank << std::endl; //} - auto ftr = find_maximal_worm_for_rank_k(filtration, simplices, point, rank); + std::vector*> rank_info = {this->rank_info_h0[i], this->rank_info_h1[i]}; + auto ftr = find_maximal_worm_for_rank_k(filtration, f_x_sorted, f_y_sorted, f_x_sorted_id, f_y_sorted_id, simplices, point, rank, rank_info); futures[k] = ftr; } @@ -466,7 +556,9 @@ void Multipers::set_grid_resolution_and_lower_left_corner(const Tensor& filtrati PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { py::class_(m, "Multipers") - .def(py::init>()) + .def(py::init>()) .def("compute_landscape", &Multipers::compute_landscape) + .def("set_hom_rank", &Multipers::set_hom_rank) + .def("refresh_rank_info", &Multipers::refresh_rank_info) .def("set_max_jobs", &Multipers::set_max_jobs); } diff --git a/gril/multipers.h b/gril/multipers.h new file mode 100644 index 0000000..4c1812a --- /dev/null +++ b/gril/multipers.h @@ -0,0 +1,120 @@ +#ifndef MULTIPERS_H +#define MULTIPERS_H +#include +#include +#include +#include +#include +#include +#include +#include +#include "utils.h" + +#include "./phat/compute_persistence_pairs.h" + +using torch::Tensor; +using namespace torch::indexing; +typedef std::pair Point; + +class Multipers{ + private: + int hom_rank; + std::vector ranks; + double step, ll_x, ll_y, res; + int px, py; + int l; + std::vector*> rank_info_h0; + std::vector*> rank_info_h1; + int num_points_guess; + void set_ranks(std::vector ranks_){ + this->ranks.insert(this->ranks.begin(), ranks_.begin(), ranks_.end()); + } + void set_step(double step){ + this->step = step; + } + void set_res(double res){ + this->res = res; + } + void set_l_for_worm(int l){ + this->l = l; + } + Tensor compute_l_worm(const int d); + std::vector> compute_filtration_along_boundary_cap(const Tensor& grid_pts_along_boundary_t, + const Tensor& f, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, + int &manual_birth_pts, + int &manual_death_pts); + + void zigzag_pairs(std::vector> &simplices_birth_death, + const vector &simplices, + const int manual_birth_pts, + const int manual_death_pts, + std::vector &num_full_bars); + + void num_full_bars_for_specific_d(const Tensor& filtration, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, + const vector& simplices, + const Point& p, + int d, + std::vector &num_full_bars); + + Tensor find_maximal_worm_for_rank_k(const Tensor &filtration, + const Tensor& f_x_sorted, + const Tensor& f_y_sorted, + const Tensor& f_x_sorted_id, + const Tensor& f_y_sorted_id, + const vector &simplices, + const Point &p, + const int rank, + std::vector*> rank_info); + + void set_grid_resolution_and_lower_left_corner(const Tensor& filtration); + + + + public: + int max_threads; + Multipers(const int hom_rank, const int l, double res, int step, const std::vector ranks){ + set_hom_rank(hom_rank); + set_l_for_worm(l); + // set_division_along_axes(px, py); + set_step(step); + set_res(res); + set_ranks(ranks); + this->max_threads = 1; + int num_cp = (int)(1.0 / res); + int num_div = num_cp / step; + this->num_points_guess = (num_div * num_div); + + for(auto i = 0; i < num_points_guess; i++){ + this->rank_info_h0.push_back(new std::map()); + this->rank_info_h1.push_back(new std::map()); + } + } + void set_max_jobs(int max_jobs){ + this->max_threads = max_jobs; + } + std::vector compute_landscape(const std::vector& pts, const std::vector>> &batch); + + void set_hom_rank(int hom_rank){ + this->hom_rank = hom_rank; + } + void refresh_rank_info(){ + for(auto i = 0; i < num_points_guess; i++){ + this->rank_info_h0[i] = new std::map(); + this->rank_info_h1[i] = new std::map(); + + } + } + + + +}; + +#endif \ No newline at end of file diff --git a/zigzag/phat/algorithms/twist_reduction.h b/gril/phat/algorithms/twist_reduction.h similarity index 100% rename from zigzag/phat/algorithms/twist_reduction.h rename to gril/phat/algorithms/twist_reduction.h diff --git a/zigzag/phat/boundary_matrix.h b/gril/phat/boundary_matrix.h similarity index 100% rename from zigzag/phat/boundary_matrix.h rename to gril/phat/boundary_matrix.h diff --git a/zigzag/phat/compute_persistence_pairs.h b/gril/phat/compute_persistence_pairs.h similarity index 100% rename from zigzag/phat/compute_persistence_pairs.h rename to gril/phat/compute_persistence_pairs.h diff --git a/zigzag/phat/helpers/dualize.h b/gril/phat/helpers/dualize.h similarity index 100% rename from zigzag/phat/helpers/dualize.h rename to gril/phat/helpers/dualize.h diff --git a/zigzag/phat/helpers/misc.h b/gril/phat/helpers/misc.h similarity index 100% rename from zigzag/phat/helpers/misc.h rename to gril/phat/helpers/misc.h diff --git a/zigzag/phat/helpers/thread_local_storage.h b/gril/phat/helpers/thread_local_storage.h similarity index 100% rename from zigzag/phat/helpers/thread_local_storage.h rename to gril/phat/helpers/thread_local_storage.h diff --git a/zigzag/phat/persistence_pairs.h b/gril/phat/persistence_pairs.h similarity index 100% rename from zigzag/phat/persistence_pairs.h rename to gril/phat/persistence_pairs.h diff --git a/zigzag/phat/representations/abstract_pivot_column.h b/gril/phat/representations/abstract_pivot_column.h similarity index 100% rename from zigzag/phat/representations/abstract_pivot_column.h rename to gril/phat/representations/abstract_pivot_column.h diff --git a/zigzag/phat/representations/bit_tree_pivot_column.h b/gril/phat/representations/bit_tree_pivot_column.h similarity index 100% rename from zigzag/phat/representations/bit_tree_pivot_column.h rename to gril/phat/representations/bit_tree_pivot_column.h diff --git a/zigzag/phat/representations/vector_vector.h b/gril/phat/representations/vector_vector.h similarity index 100% rename from zigzag/phat/representations/vector_vector.h rename to gril/phat/representations/vector_vector.h diff --git a/zigzag/setup.py b/gril/setup.py similarity index 100% rename from zigzag/setup.py rename to gril/setup.py diff --git a/gril/test.py b/gril/test.py new file mode 100644 index 0000000..f08acfa --- /dev/null +++ b/gril/test.py @@ -0,0 +1,6 @@ +from torch.utils.cpp_extension import load +from torch import Tensor + +test = load('test', sources=['test.cpp'], extra_cflags= ['-fopenmp'], verbose=True) + +test.test(32) \ No newline at end of file diff --git a/gril/torchph/.gitignore b/gril/torchph/.gitignore new file mode 100644 index 0000000..7798908 --- /dev/null +++ b/gril/torchph/.gitignore @@ -0,0 +1,9 @@ +develop.py +develop_2.py +.idea +.vscode +__pycache__ +*.pyc +.cache +.pytest_cache +pershom_dev/extensions_sandbox diff --git a/gril/torchph/chofer_torchex/__init__.py b/gril/torchph/chofer_torchex/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/chofer_torchex/nn/__init__.py b/gril/torchph/chofer_torchex/nn/__init__.py new file mode 100644 index 0000000..2607515 --- /dev/null +++ b/gril/torchph/chofer_torchex/nn/__init__.py @@ -0,0 +1,2 @@ +from .slayer import SLayerExponential, SLayerRational, SLayerRationalHat +from .modules import * diff --git a/gril/torchph/chofer_torchex/nn/fuctional.py b/gril/torchph/chofer_torchex/nn/fuctional.py new file mode 100644 index 0000000..e169be3 --- /dev/null +++ b/gril/torchph/chofer_torchex/nn/fuctional.py @@ -0,0 +1,49 @@ +import torch + + +def histogram_intersection_loss(input: torch.Tensor, + target: torch.Tensor, + size_average: bool=True, + reduce: bool=True, + symetric_version: bool=True)->torch.Tensor: + r""" + This loss function is based on the `Histogram Intersection` score. + + The output is the *negative* Histogram Intersection Score. + + Args: + input (Tensor): :math:`(N, B)` where `N = batch size` and `B = number of classes` + target (Tensor): :math:`(N, B)` where `N = batch size` and `B = number of classes` + size_average (bool, optional): By default, the losses are averaged + over observations for each minibatch. However, if the field + :attr:`size_average` is set to ``False``, the losses are instead summed + for each minibatch. Ignored if :attr:`reduce` is ``False``. Default: ``True`` + reduce (bool, optional): + symetric_version (bool, optional): By default, the symetric version of histogram intersection + is used. If false the asymetric version is used. Default: ``True`` + + Returns: Tensor. + + """ + assert input.size() == target.size(), \ + "input.size() != target.size(): {} != {}!".format(input.size(), target.size()) + assert input.dim() == target.dim() == 2, \ + "input, target must be 2 dimensional. Got dim {} resp. {}".format(input.dim(), target.dim()) + + minima = input.min(target) + summed_minima = minima.sum(dim=1) + + if symetric_version: + normalization_factor = (input.sum(dim=1)).max(target.sum(dim=1)) + else: + normalization_factor = target.sum(dim=1) + + loss = summed_minima / normalization_factor + + if reduce: + loss = sum(loss) + + if size_average: + loss = loss / input.size(0) + + return -loss diff --git a/gril/torchph/chofer_torchex/nn/modules.py b/gril/torchph/chofer_torchex/nn/modules.py new file mode 100644 index 0000000..395741f --- /dev/null +++ b/gril/torchph/chofer_torchex/nn/modules.py @@ -0,0 +1,45 @@ +import torch +import torch.nn as nn + + +class LinearView(nn.Module): + def __init__(self): + super().__init__() + + def forward(self, x): + return x.view(x.size()[0], -1) + + +class Apply(nn.Module): + def __init__(self, function): + super().__init__() + self.function = function + + def forward(self, *args, **kwargs): + return self.function(*args, **kwargs) + + +class IndependentBranchesLinear(nn.Linear): + def __init__(self, in_features, out_features_branch, n_branches, bias=True): + assert in_features % n_branches == 0 + in_features_branch = int(in_features/n_branches) + super().__init__(in_features, out_features_branch*n_branches, bias) + + mask = torch.zeros_like(self.weight) + for i in range(n_branches): + mask[i*out_features_branch:(i+1)*out_features_branch, + i*in_features_branch:(i+1)*in_features_branch] = 1 + + self.register_buffer('mask', mask) + + def forward(self, inputs): + return torch.nn.functional.linear(inputs, self.weight * self.mask, self.bias) + + +class View(nn.Module): + def __init__(self, view_args): + super().__init__() + self.view_args = view_args + + def forward(self, input): + return input.view(*self.view_args) diff --git a/gril/torchph/chofer_torchex/nn/slayer.py b/gril/torchph/chofer_torchex/nn/slayer.py new file mode 100644 index 0000000..aaf6e87 --- /dev/null +++ b/gril/torchph/chofer_torchex/nn/slayer.py @@ -0,0 +1,492 @@ +r""" + +Implementation of **differentiable vectorization layers** for persistent +homology barcodes. + +For a basic tutorial click `here `_. +""" +import torch +import numpy as np +from torch import Tensor +from torch.nn.parameter import Parameter +from torch.nn.modules.module import Module + +from typing import List, Tuple + +import warnings + + +# region helper functions + + +def prepare_batch( + batch: List[Tensor], + point_dim: int=None + )->Tuple[Tensor, Tensor, int, int]: + """ + This method 'vectorizes' the multiset in order to take advances of GPU + processing. The policy is to embed all multisets in batch to the highest + dimensionality occurring in batch, i.e., ``max(t.size()[0]`` for ``t`` in batch). + + Args: + batch: + The input batch to process as a list of tensors. + + point_dim: + The dimension of the points the inputs consist of. + + Returns: + A four-tuple consisting of (1) the constructed ``batch``, i.e., a + tensor with size + ``batch_size`` x ``n_max_points`` x ``point_dim``; (2) a tensor + ``not_dummy`` of size ``batch_size`` x ``n_max_points``, where + ``1`` at position (i,j) indicates if the point is a dummy point, + whereas ``0`` indicates a dummy point used for padding; (3) + the max. number of points and (4) the batch size. + + Example:: + + >>> from chofer_torchex.nn.slayer import prepare_batch + >>> import torch + >>> x = [torch.rand(10,2), torch.rand(20,2)] + >>> batch, not_dummy, max_pts, batch_size = prepare_batch(x) + """ + if point_dim is None: + point_dim = batch[0].size(1) + assert (all(x.size(1) == point_dim for x in batch if len(x) != 0)) + + batch_size = len(batch) + batch_max_points = max([t.size(0) for t in batch]) + input_device = batch[0].device + + if batch_max_points == 0: + # if we are here, batch consists only of empty diagrams. + batch_max_points = 1 + + # This will later be used to set the dummy points to zero in the output. + not_dummy_points = torch.zeros( + batch_size, + batch_max_points, + device=input_device) + + prepared_batch = [] + + for i, multi_set in enumerate(batch): + n_points = multi_set.size(0) + + prepared_dgm = torch.zeros( + batch_max_points, + point_dim, + device=input_device) + + if n_points > 0: + index_selection = torch.tensor(range(n_points), + device=input_device) + + prepared_dgm.index_add_(0, index_selection, multi_set) + + not_dummy_points[i, :n_points] = 1 + + prepared_batch.append(prepared_dgm) + + prepared_batch = torch.stack(prepared_batch) + + return prepared_batch, not_dummy_points, batch_max_points, batch_size + + +def is_prepared_batch(input): + if not (isinstance(input, tuple) and len(input) == 4): + return False + else: + batch, not_dummy_points, max_points, batch_size = input + return isinstance(batch, Tensor) and isinstance(not_dummy_points, Tensor) and max_points > 0 and batch_size > 0 + + +def is_list_of_tensors(input): + try: + return all([isinstance(x, Tensor) for x in input]) + + except TypeError: + return False + + +def prepare_batch_if_necessary(input, point_dimension=None): + batch, not_dummy_points, max_points, batch_size = None, None, None, None + + if is_prepared_batch(input): + batch, not_dummy_points, max_points, batch_size = input + elif is_list_of_tensors(input): + if point_dimension is None: + point_dimension = input[0].size(1) + + batch, not_dummy_points, max_points, batch_size = prepare_batch( + input, + point_dimension) + + else: + raise ValueError( + 'SLayer does not recognize input format! Expecting [Tensor] or \ + prepared batch. Not {}'.format(input)) + + return batch, not_dummy_points, max_points, batch_size + + +def parameter_init_from_arg(arg, size, default, scalar_is_valid=False): + if isinstance(arg, (int, float)): + if not scalar_is_valid: + raise ValueError('Scalar initialization values are not valid. \ + Got {} expected Tensor of size {}.' + .format(arg, size)) + return torch.Tensor(*size).fill_(arg) + elif isinstance(arg, torch.Tensor): + assert(arg.size() == size) + return arg + elif arg is None: + if default in [torch.rand, torch.randn, torch.ones, torch.ones_like]: + return default(*size) + else: + return default(size) + else: + raise ValueError('Cannot handle parameter initialization. \ + Got "{}" '.format(arg)) + +# endregion + + +class SLayerExponential(Module): + """ + Proposed input layer for multisets [1]. + """ + def __init__(self, n_elements: int, + point_dimension: int=2, + centers_init: Tensor=None, + sharpness_init: Tensor=None): + """ + Args: + n_elements: + Number of structure elements used. + + point_dimension: D + Dimensionality of the points of which the + input multi set consists of. + + centers_init: + The initialization for the centers of the structure elements. + + sharpness_init: + Initialization for the sharpness of the structure elements. + """ + super().__init__() + + self.n_elements = n_elements + self.point_dimension = point_dimension + + expected_init_size = (self.n_elements, self.point_dimension) + + centers_init = parameter_init_from_arg( + centers_init, + expected_init_size, + torch.rand, scalar_is_valid=False) + sharpness_init = parameter_init_from_arg( + sharpness_init, + expected_init_size, + lambda size: torch.ones(*size)*3) + + self.centers = Parameter(centers_init) + self.sharpness = Parameter(sharpness_init) + + def forward(self, input)->Tensor: + batch, not_dummy_points, max_points, batch_size = prepare_batch_if_necessary( + input, + point_dimension=self.point_dimension) + + batch = torch.cat([batch] * self.n_elements, 1) + + not_dummy_points = torch.cat([not_dummy_points] * self.n_elements, 1) + + centers = torch.cat([self.centers] * max_points, 1) + centers = centers.view(-1, self.point_dimension) + centers = torch.stack([centers] * batch_size, 0) + + sharpness = torch.pow(self.sharpness, 2) + sharpness = torch.cat([sharpness] * max_points, 1) + sharpness = sharpness.view(-1, self.point_dimension) + sharpness = torch.stack([sharpness] * batch_size, 0) + + x = centers - batch + x = x.pow(2) + x = torch.mul(x, sharpness) + x = torch.sum(x, 2) + x = torch.exp(-x) + x = torch.mul(x, not_dummy_points) + x = x.view(batch_size, self.n_elements, -1) + x = torch.sum(x, 2) + x = x.squeeze() + + return x + + def __repr__(self): + return 'SLayerExponential (... -> {} )'.format(self.n_elements) + + +class SLayerRational(Module): + """ + """ + def __init__(self, n_elements: int, + point_dimension: int=2, + centers_init: Tensor=None, + sharpness_init: Tensor=None, + exponent_init: Tensor=None, + pointwise_activation_threshold=None, + share_sharpness=False, + share_exponent=False, + freeze_exponent=True): + """ + Args: + n_elements: + Number of structure elements used. + + point_dimension: + Dimensionality of the points of which the input multi set + consists of. + + centers_init: + The initialization for the centers of the structure elements. + + sharpness_init: + Initialization for the sharpness of the structure elements. + """ + super().__init__() + + self.n_elements = int(n_elements) + self.point_dimension = int(point_dimension) + self.pointwise_activation_threshold = float(pointwise_activation_threshold) \ + if pointwise_activation_threshold is not None else None + self.share_sharpness = bool(share_sharpness) + self.share_exponent = bool(share_exponent) + self.freeze_exponent = freeze_exponent + + if self.pointwise_activation_threshold is not None: + self.pointwise_activation_threshold = float(self.pointwise_activation_threshold) + + centers_init = parameter_init_from_arg( + arg=centers_init, + size=(self.n_elements, self.point_dimension), + default=torch.rand) + + sharpness_init = parameter_init_from_arg( + arg=sharpness_init, + size=(1,) if self.share_sharpness else (self.n_elements, self.point_dimension), + default=torch.ones, + scalar_is_valid=True) + + exponent_init = parameter_init_from_arg( + arg=exponent_init, + size=(1,) if self.share_exponent else (self.n_elements,), + default=torch.ones, + scalar_is_valid=True) + + self.centers = Parameter(centers_init) + self.sharpness = Parameter(sharpness_init) + + if self.freeze_exponent: + self.register_buffer('exponent', exponent_init) + else: + self.exponent = Parameter(exponent_init) + + def forward(self, input)->Tensor: + batch, not_dummy_points, max_points, batch_size = prepare_batch_if_necessary( + input, + point_dimension=self.point_dimension) + + batch = batch.unsqueeze(1).expand( + batch_size, + self.n_elements, + max_points, + self.point_dimension) + + not_dummy_points = not_dummy_points.unsqueeze(1).expand(-1, self.n_elements, -1) + + centers = self.centers.unsqueeze(1).expand( + self.n_elements, + max_points, + self.point_dimension) + centers = centers.unsqueeze(0).expand(batch_size, *centers.size()) + + sharpness = self.sharpness + if self.share_sharpness: + sharpness = sharpness.expand(self.n_elements, self.point_dimension) + sharpness = sharpness.unsqueeze(1).expand(-1, max_points, -1) + sharpness = sharpness.unsqueeze(0).expand(batch_size, *sharpness.size()) + + exponent = self.exponent + if self.share_exponent: + exponent = exponent.expand(self.n_elements) + + exponent = exponent.unsqueeze(1).expand(-1, max_points) + exponent = exponent.unsqueeze(0).expand(batch_size, *exponent.size()) + + x = centers - batch + x = x.abs() + x = torch.mul(x, sharpness.abs()) + x = torch.sum(x, 3) + x = 1/(1+x).pow(exponent.abs()) + + if self.pointwise_activation_threshold is not None: + x[(x < self.pointwise_activation_threshold).data] = 0 + + x = torch.mul(x, not_dummy_points) + x = torch.sum(x, 2) + + return x + + def __repr__(self): + return 'SLayerRational (... -> {} )'.format(self.n_elements) + + +class SLayerRationalHat(Module): + """ + """ + def __init__(self, n_elements: int, + point_dimension: int=2, + centers_init: Tensor=None, + radius_init: float=1, + exponent: int=1 + ): + """ + Args: + n_elements: + Number of structure elements used. + + point_dimension: + Dimensionality of the points of which the input multi + set consists of. + + centers_init: + The initialization for the centers of the structure elements. + + radius_init: + Initialization for radius of zero level-set of the hat. + + exponent: + Exponent of the rationals forming the hat. + + """ + super().__init__() + + self.n_elements = int(n_elements) + self.point_dimension = int(point_dimension) + self.exponent = int(exponent) + + centers_init = parameter_init_from_arg(arg=centers_init, + size=(self.n_elements, self.point_dimension), + default=torch.rand) + + radius_init = parameter_init_from_arg(arg=radius_init, + size=(self.n_elements,), + default=torch.ones, + scalar_is_valid=True) + + self.centers = Parameter(centers_init) + self.radius = Parameter(radius_init) + self.norm_p = 1 + + def forward(self, input)->Tensor: + batch, not_dummy_points, max_points, batch_size = prepare_batch_if_necessary( + input, + point_dimension=self.point_dimension) + + batch = batch.unsqueeze(1).expand( + batch_size, + self.n_elements, + max_points, + self.point_dimension) + + not_dummy_points = not_dummy_points.unsqueeze(1).expand(-1, self.n_elements, -1) + + centers = self.centers.unsqueeze(1).expand( + self.n_elements, + max_points, + self.point_dimension) + centers = centers.unsqueeze(0).expand(batch_size, *centers.size()) + + radius = self.radius + radius = radius.unsqueeze(1).expand(-1, max_points) + radius = radius.unsqueeze(0).expand(batch_size, *radius.size()) + radius = radius.abs() + + norm_to_center = centers - batch + norm_to_center = torch.norm(norm_to_center, p=self.norm_p, dim=3) + + positive_part = 1/(1+norm_to_center).pow_(self.exponent) + + negative_part = 1/(1 + (radius - norm_to_center).abs_()).pow_(self.exponent) + + x = positive_part - negative_part + + x = torch.mul(x, not_dummy_points) + x = torch.sum(x, 2) + + return x + + def __repr__(self): + return 'SLayerRationalHat (... -> {} )'.format(self.n_elements) + + +class LinearRationalStretchedBirthLifeTimeCoordinateTransform: + def __init__(self, nu): + self._nu = nu + self._nu_squared = nu**2 + self._2_nu = 2*nu + + def __call__(self, dgm): + if len(dgm) == 0: + return dgm + + x, y = dgm[:, 0], dgm[:, 1] + y = y - x + + i = (y <= self._nu) + y[i] = - self._nu_squared/y[i] + self._2_nu + + return torch.stack([x, y], dim=1) + + +class LogStretchedBirthLifeTimeCoordinateTransform: + def __init__(self, nu): + self.nu = nu + + def __call__(self, dgm): + if len(dgm) == 0: + return dgm + + x, y = dgm[:, 0], dgm[:, 1] + y = y - x + + i = (y <= self.nu) + y[i] = torch.log(y[i] / self.nu)*self.nu + self.nu + + return torch.stack([x, y], dim=1) + + +class UpperDiagonalThresholdedLogTransform: + def __init__(self, nu): + self.b_1 = (torch.Tensor([1, 1]) / np.sqrt(2)) + self.b_2 = (torch.Tensor([-1, 1]) / np.sqrt(2)) + self.nu = nu + + def __call__(self, dgm): + if len(dgm) == 0: + return dgm + + self.b_1 = self.b_1.to(dgm.device) + self.b_2 = self.b_2.to(dgm.device) + + x = torch.mul(dgm, self.b_1.repeat(dgm.size(0), 1)) + x = torch.sum(x, 1).squeeze() + y = torch.mul(dgm, self.b_2.repeat(dgm.size(0), 1)) + y = torch.sum(y, 1).squeeze() + i = (y <= self.nu) + y[i] = torch.log(y[i] / self.nu)*self.nu + self.nu + ret = torch.stack([x, y], 1) + return ret diff --git a/gril/torchph/chofer_torchex/pershom/__init__.py b/gril/torchph/chofer_torchex/pershom/__init__.py new file mode 100644 index 0000000..0cb7d2c --- /dev/null +++ b/gril/torchph/chofer_torchex/pershom/__init__.py @@ -0,0 +1 @@ +from .pershom_backend import * diff --git a/gril/torchph/chofer_torchex/pershom/pershom_backend.py b/gril/torchph/chofer_torchex/pershom/pershom_backend.py new file mode 100644 index 0000000..d3deace --- /dev/null +++ b/gril/torchph/chofer_torchex/pershom/pershom_backend.py @@ -0,0 +1,266 @@ +r""" +This module exposes the C++/CUDA backend functionality for Python. + +Terminology +----------- + +Descending sorted boundary array: + Boundary array which encodes the boundary matrix (BM) for a given + filtration in column first order. + Let BA be the descending_sorted_boundary of BM, then + ``BA[i, :]`` is the i-th column of BM. + Content encoded as decreasingly sorted list, embedded into the array + with -1 padding from the right. + + Example : + ``BA[3, :] = [2, 0, -1, -1]`` + then :math:`\partial(v_3) = v_0 + v_2` + + ``BA[6, :] = [5, 4, 3, -1]`` + then :math:`\partial(v_6) = v_3 + v_4 + v_5` + + +Compressed descending sorted boundary array: + Same as *descending sorted boundary array* but rows consisting only of -1 + are omitted. + This is sometimes used for efficiency purposes and is usually accompanied + by a vector, ``v``, telling which row of the reduced BA corresponds to + which row of the uncompressed BA, i.e., ``v[3] = 5`` means that the 3rd + row of the reduced BA corresponds to the 5th row in the uncompressed + version. + +Birth/Death-time: + Index of the coresponding birth/death event in the filtration. + This is always an *integer*. + +Birth/Death-value: + If a filtration is induced by a real-valued function, this corresponds + to the value of this function corresponding to the birth/death event. + This is always *real*-valued. + +Limitations +----------- + +Currently all ``cuda`` backend functionality **only** supports ``int64_t`` and +``float32_t`` typing. + +""" +import os.path as pth +from typing import List +from torch import Tensor +from glob import glob + + +from torch.utils.cpp_extension import load + + +__module_file_dir = pth.dirname(pth.realpath(__file__)) +__cpp_src_dir = pth.join(__module_file_dir, 'pershom_cpp_src') +src_files = [] + +for extension in ['*.cpp', '*.cu']: + src_files += glob(pth.join(__cpp_src_dir, extension)) + +# jit compiling the c++ extension + +try: + __C = load( + 'pershom_cuda_ext', + src_files, + verbose=False) + +except Exception as ex: + print("Error in {}. Failed jit compilation. Maybe your CUDA environment \ + is messed up?".format(__file__)) + print("Error was {}".format(ex)) + + +def find_merge_pairings( + pivots: Tensor, + max_pairs: int = -1 + )->Tensor: + """Finds the pairs which have to be merged in the current iteration of the + matrix reduction. + + Args: + pivots: + The pivots of a descending sorted boundary array. + Expected size is ``Nx1``, where N is the number of columns of the + underlying descending sorted boundary array. + + max_pairs: + The output is at most a ``max_pairs x 2`` Tensor. If set to + default all possible merge pairs are returned. + + Returns: + The merge pairs, ``p``, for the current iteration of the reduction. + ``p[i]`` is a merge pair. + In boundary matrix notation this would mean column ``p[i][0]`` has to + be merged into column ``p[i][1]``. + """ + return __C.CalcPersCuda__find_merge_pairings(pivots, max_pairs) + + +def merge_columns_( + compr_desc_sort_ba: Tensor, + merge_pairs: Tensor + )->None: + r"""Executes the given merging operations inplace on the descending + sorted boundary array. + + Args: + compr_desc_sort_ba: + see module description top. + + merge_pairs: + output of a ``find_merge_pairings`` call. + + Returns: + None + """ + __C.CalcPersCuda__merge_columns_(compr_desc_sort_ba, merge_pairs) + + +def read_barcodes( + pivots: Tensor, + simplex_dimension: Tensor, + max_dim_to_read_of_reduced_ba: int + )->List[List[Tensor]]: + """Reads the barcodes using the pivot of a reduced boundary array. + + Arguments: + pivots: + pivots is the first column of a compr_desc_sort_ba + + simplex_dimension: + Vector whose i-th entry is the dimension if the i-th simplex in + the given filtration. + + max_dim_to_read_of_reduced_ba: + features up to max_dim_to_read_of_reduced_ba are read from the + reduced boundary array + + Returns: + List of birth/death times. + + ``ret[0][n]`` are non essential birth/death-times of dimension ``n``. + + ``ret[1][n]`` are birth-times of essential classes of dimension ``n``. + """ + return __C.CalcPersCuda__read_barcodes( + pivots, + simplex_dimension, + max_dim_to_read_of_reduced_ba) + + +def calculate_persistence( + compr_desc_sort_ba: Tensor, + ba_row_i_to_bm_col_i: Tensor, + simplex_dimension: Tensor, + max_dim_to_read_of_reduced_ba: int, + max_pairs: int = -1 + )->List[List[Tensor]]: + """Returns the barcodes of the given encoded boundary array. + + Arguments: + compr_desc_sort_ba: + A `compressed descending sorted boundary array`, + see readme section top. + + ba_row_i_to_bm_col_i: + Vector whose i-th entry is the column index of the boundary + matrix the i-th row in ``compr_desc_sort_ba corresponds`` to. + + simplex_dimension: + Vector whose i-th entry is the dimension if the i-th simplex in + the given filtration + + max_pairs: see ``find_merge_pairings``. + + max_dim_to_read_of_reduced_ba: + features up to max_dim_to_read_of_reduced_ba are read from the + reduced boundary array. + + Returns: + List of birth/death times. + + ``ret[0][n]`` are non essential birth/death-times of dimension ``n``. + + ``ret[1][n]`` are birth-times of essential classes of dimension ``n``. + """ + return __C.CalcPersCuda__calculate_persistence( + compr_desc_sort_ba, + ba_row_i_to_bm_col_i, + simplex_dimension, + max_dim_to_read_of_reduced_ba, + max_pairs) + + +def vr_persistence_l1( + point_cloud: Tensor, + max_dimension: int, + max_ball_diameter: float=0.0 + )->List[List[Tensor]]: + """Returns the barcodes of the Vietoris-Rips complex of a given point cloud + w.r.t. the l1 (manhatten) distance. + + Args: + point_cloud: + Point cloud from which the Vietoris-Rips complex is generated. + + max_dimension: + The dimension of the used Vietoris-Rips complex. + + max_ball_diameter: + If not 0, edges whose two defining vertices' distance is greater + than ``max_ball_diameter`` are ignored. + + Returns: + List of birth/death times. + + ``ret[0][n]`` are non essential birth/death-*values* of dimension ``n``. + + ``ret[1][n]`` are birth-*values* of essential classes of + dimension ``n``. + """ + return __C.VRCompCuda__vr_persistence_l1( + point_cloud, + max_dimension, + max_ball_diameter) + + +def vr_persistence( + distance_matrix: Tensor, + max_dimension: int, + max_ball_diameter: float=0.0 + )->List[List[Tensor]]: + """Returns the barcodes of the Vietoris-Rips complex of a given distance + matrix. + + **Note**: ``distance_matrix`` is assumed to be a square matrix. + Practically, symmetry is *not* required and the upper diagonal part is + *ignored*. For the computation, just the *lower* diagonal part is used. + + Args: + distance_matrix: + Distance matrix the Vietoris-Rips complex is based on. + + max_dimension: + The dimension of the used Vietoris-Rips complex. + + max_ball_diameter: + If not 0, edges whose two defining vertices' distance is greater + than ``max_ball_diameter`` are ignored. + + Returns: + List of birth/death times. + + ``ret[0][n]`` are non essential birth/death-*values* of dimension ``n``. + + ``ret[1][n]`` are birth-*values* of essential classes of + dimension ``n``. + """ + return __C.VRCompCuda__vr_persistence( + distance_matrix, + max_dimension, + max_ball_diameter) diff --git a/gril/torchph/chofer_torchex/pershom/pershom_cpp_src/__init__.py b/gril/torchph/chofer_torchex/pershom/pershom_cpp_src/__init__.py new file mode 100644 index 0000000..0cb7d2c --- /dev/null +++ b/gril/torchph/chofer_torchex/pershom/pershom_cpp_src/__init__.py @@ -0,0 +1 @@ +from .pershom_backend import * diff --git a/gril/torchph/chofer_torchex/utils/__init__.py b/gril/torchph/chofer_torchex/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/chofer_torchex/utils/boiler_plate.py b/gril/torchph/chofer_torchex/utils/boiler_plate.py new file mode 100644 index 0000000..ebe7314 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/boiler_plate.py @@ -0,0 +1,45 @@ +import torch + + +def apply_model(model, dataset, batch_size=100, device='cpu'): + dl = torch.utils.data.DataLoader( + dataset, + batch_size=batch_size, + num_workers=0 + ) + X, Y = [], [] + + model.eval() + model.to(device) + + with torch.no_grad(): + for x, y in dl: + + x = x.to(device) + x = model(x) + + # If the model returns more than one tensor, e.g., encoder of an + # autoencoder model, take the first one as output... + if isinstance(x, (tuple, list)): + x = x[0] + + X += x.cpu().tolist() + Y += (y.tolist()) + + return X, Y + + +def argmax_and_accuracy(X, Y, binary=False): + + if not binary: + decision = torch.tensor(X).max(1)[1] + else: + X = torch.tensor(X).squeeze() + assert X.ndimension() == 1 + decision = (X > 0.0).long() + + Y = torch.tensor(Y) + + correct = (decision == Y).sum().float().item() + num_samples = float(Y.size(0)) + return (correct/num_samples)*100.0 diff --git a/gril/torchph/chofer_torchex/utils/data/__init__.py b/gril/torchph/chofer_torchex/utils/data/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/chofer_torchex/utils/data/collate.py b/gril/torchph/chofer_torchex/utils/data/collate.py new file mode 100644 index 0000000..7ba1d66 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/data/collate.py @@ -0,0 +1,39 @@ +from collections import defaultdict + + +def dict_sample_target_iter_concat(sample_target_iter: iter): + """ + Gets an sample target iterator of dict samples. Returns + a concatenation of the samples based on each key and the + target list. + + Example: + ``` + sample_target_iter = [({'a': 'a1', 'b': 'b1'}, 0), ({'a': 'a2', 'b': 'b2'}, 1)] + x = dict_sample_iter_concat([({'a': 'a1', 'b': 'b1'}, 0), ({'a': 'a2', 'b': 'b2'}, 1)]) + print(x) + ({'a': ['a1', 'a2'], 'b': ['b1', 'b2']}, [0, 1]) + ``` + + :param sample_target_iter: + :return: + """ + + samples = defaultdict(list) + targets = [] + + for sample_dict, y in sample_target_iter: + for k, v in sample_dict.items(): + samples[k].append(v) + + targets.append(y) + + samples = dict(samples) + + length = len(samples[next(iter(samples))]) + assert all(len(samples[k]) == length for k in samples) + + return samples, targets + + + diff --git a/gril/torchph/chofer_torchex/utils/data/dataset.py b/gril/torchph/chofer_torchex/utils/data/dataset.py new file mode 100644 index 0000000..adbca73 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/data/dataset.py @@ -0,0 +1,67 @@ +import requests +import torch + + +class Dataset(object): + def __getitem__(self, index): + raise NotImplementedError + + def __len__(self): + raise NotImplementedError + + @property + def labels(self): + raise NotImplementedError + + @property + def sample_labels(self): + raise NotImplementedError + + +class SimpleDataset(torch.utils.data.dataset.Dataset): + def __init__(self, X, Y): + super().__init__() + self.X = X + self.Y = Y + + def __len__(self): + return len(self.Y) + + def __getitem__(self, item): + return self.X[item], self.Y[item] + + def __iter__(self): + return zip(self.X, self.Y) + + +def _download_file_from_google_drive(id, destination): + def get_confirm_token(response): + for key, value in response.cookies.items(): + if key.startswith('download_warning'): + return value + + return None + + def save_response_content(response, destination): + CHUNK_SIZE = 32768 + content_iter = response.iter_content(CHUNK_SIZE) + with open(destination, "wb") as f: + + for i, chunk in enumerate(content_iter): + if chunk: # filter out keep-alive new chunks + f.write(chunk) + print(i, end='\r') + + URL = "https://docs.google.com/uc?export=download" + + session = requests.Session() + + response = session.get(URL, params={'id': id}, stream=True) + + token = get_confirm_token(response) + + if token: + params = {'id': id, 'confirm': token} + response = session.get(URL, params=params, stream=True) + + save_response_content(response, destination) diff --git a/gril/torchph/chofer_torchex/utils/data/ds_operations.py b/gril/torchph/chofer_torchex/utils/data/ds_operations.py new file mode 100644 index 0000000..8561b0c --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/data/ds_operations.py @@ -0,0 +1,69 @@ +import torch +import numpy as np +from typing import Union +import torch.utils.data + + +def ds_random_subset( + dataset: torch.utils.data.dataset.Dataset, + percentage: float=None, + absolute_size: int=None, + replace: bool=False): + r""" + Represents a fixed random subset of the given dataset. + + Args: + dataset (torch.utils.data.dataset.Dataset): Target dataset. + percentage (float): Percentage of target dataset to use (within [0,1]). + absolute_size (int): Absolute size of the subset to use. + replace (bool): Draw with replacement. + + Returns: + A ``torch.utils.data.dataset.Dataset`` with randomly selected samples. + + .. note:: + ``percentage`` and ``absolute_size`` are mutally exclusive. So only + one of them can be specified. + """ + assert isinstance(dataset, torch.utils.data.dataset.Dataset) + assert percentage is not None or absolute_size is not None + assert not (percentage is None and absolute_size is None) + if percentage is not None: + assert 0 < percentage and percentage < 1, "percentage assumed to be > 0 and < 1" + if absolute_size is not None: + assert absolute_size <= len(dataset) + + n_samples = int(percentage*len(dataset)) if percentage is not None else absolute_size + indices = np.random.choice( + list(range(len(dataset))), + n_samples, + replace=replace) + + indices = [int(i) for i in indices] + + return torch.utils.data.dataset.Subset(dataset, indices) + + +def ds_label_filter( + dataset: torch.utils.data.dataset.Dataset, + labels: Union[tuple, list]): + """ + Returns a dataset with samples having selected labels. + + Args: + dataset (torch.utils.data.dataset.Dataset): Target dataset. + labels (tuple or list): White list of labels to use. + + Returns: + A ``torch.utils.data.dataset.Dataset`` only containing samples having + the selected labels. + """ + assert isinstance(dataset, torch.utils.data.dataset.Dataset) + assert isinstance(labels, (tuple, list)), "labels is expected to be list or tuple." + assert len(set(labels)) == len(labels), "labels is expected to have unique elements." + assert hasattr(dataset, 'targets'), "dataset is expected to have 'targets' attribute" + assert set(labels) <= set(dataset.targets), "labels is expected to contain only valid labels of dataset" + + indices = [i for i in range(len(dataset)) if dataset.targets[i] in labels] + + return torch.utils.data.dataset.Subset(dataset, indices) diff --git a/gril/torchph/chofer_torchex/utils/data/iterable_decorators.py b/gril/torchph/chofer_torchex/utils/data/iterable_decorators.py new file mode 100644 index 0000000..aeb8bea --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/data/iterable_decorators.py @@ -0,0 +1,22 @@ +class ConditionalIterable: + def __init__(self, iterable, pred): + self.decoratee = iterable + self.pred = pred + + def __iter__(self): + return (x for x in self.decoratee if self.pred(x)) + + def __len__(self): + return len(self.decoratee) + + +class TransformingIterable: + def __init__(self, iterable, transform): + self.decoratee = iterable + self.transform = transform + + def __iter__(self): + return (self.transform(x) for x in self.decoratee) + + def __len__(self): + return len(self.decoratee) diff --git a/gril/torchph/chofer_torchex/utils/functional.py b/gril/torchph/chofer_torchex/utils/functional.py new file mode 100644 index 0000000..f51ced7 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/functional.py @@ -0,0 +1,28 @@ +# TODO: Possibly delete + +import torch + + +def collection_cascade(input, stop_predicate: callable, function_to_apply: callable): + if stop_predicate(input): + return function_to_apply(input) + elif isinstance(input, list or tuple): + return [collection_cascade( + x, + stop_predicate=stop_predicate, + function_to_apply=function_to_apply) for x in input] + elif isinstance(input, dict): + return {k: collection_cascade( + v, + stop_predicate=stop_predicate, + function_to_apply=function_to_apply) for k, v in input.items()} + else: + raise ValueError('Unknown type collection type. Expected list, \ + tuple, dict but got {}'.format(type(input))) + + +def cuda_cascade(input, **kwargs): + return collection_cascade( + input, + stop_predicate=lambda x: isinstance(x, torch.Tensor), + function_to_apply=lambda x: x.cuda(**kwargs)) diff --git a/gril/torchph/chofer_torchex/utils/ipynb/plot.py b/gril/torchph/chofer_torchex/utils/ipynb/plot.py new file mode 100644 index 0000000..99439dc --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/ipynb/plot.py @@ -0,0 +1,16 @@ +import matplotlib.pyplot as plt +import numpy as np +from torchvision.utils import make_grid + + +def plt_img_grid(images, nrow=8, figure=None, scale_each=True): + """ + Plot images on a grid. + """ + if figure is None: + plt.figure(figsize=(16, 16)) + + ax = plt.gca() + img = make_grid(images, normalize=True, nrow=nrow, scale_each=scale_each) + npimg = img.cpu().detach().numpy() + ax.imshow(np.transpose(npimg, (1, 2, 0)), interpolation='nearest') diff --git a/gril/torchph/chofer_torchex/utils/run_experiments.py b/gril/torchph/chofer_torchex/utils/run_experiments.py new file mode 100644 index 0000000..5674431 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/run_experiments.py @@ -0,0 +1,220 @@ +from typing import Callable +import torch +import multiprocessing as mp +import itertools +import contextlib +import sys +import traceback +import time + + +from pathlib import Path + + +class DummyFile(object): + def write(self, x): pass + + def flush(self, *args, **kwargs): pass + + +@contextlib.contextmanager +def nostdout(): + save_stdout = sys.stdout + sys.stdout = DummyFile() + yield + sys.stdout = save_stdout + + +class _ErrorCatchingProcess(mp.Process): + def __init__(self, *args, **kwargs): + mp.Process.__init__(self, *args, **kwargs) + self._pconn, self._cconn = mp.Pipe() + self._exception = None + + def run(self): + try: + mp.Process.run(self) + self._cconn.send(None) + except Exception as e: + tb = traceback.format_exc() + self._cconn.send((e, tb)) + + @property + def exception(self): + if self._pconn.poll(): + self._exception = self._pconn.recv() + return self._exception + + +class _GPUTask(): + def __init__( + self, + fn, + args_kwargs, + gpu_id): + + self.fn = fn + self.args_kwargs = args_kwargs + self.gpu_id = gpu_id + + def __call__(self): + with torch.cuda.device(self.gpu_id): + with nostdout(): + args, kwargs = self.args_kwargs + self.fn(*args, **kwargs) + + +class Distributor: + + def __init__(self, + fn: Callable, + args_kwargs: list, + visible_devices: list, + num_max_jobs_on_device: int, + log_file_path=None): + + self.visible_devices = visible_devices + self.counter = {k: 0 for k in visible_devices} + + self.fn = fn + self.args_kwargs = args_kwargs + self.num_max_jobs = num_max_jobs_on_device + self.log_file_path = None if log_file_path is None else Path( + log_file_path) + + if self.log_file_path is not None: + assert self.log_file_path.parent.is_dir() + + self.progress_wheel = itertools.cycle(['|', '/', '-', '\\']) + + def _log(self, *args): + print(*args) + + if self.log_file_path is not None: + with open(self.log_file_path, 'w') as fid: + print(*args, file=fid) + + def _iter_free_devices(self): + for k, v in self.counter.items(): + if v < self.num_max_jobs: + yield k + + def run(self): + + q = list(enumerate(self.args_kwargs)) + + processes = [] + + finished_jobs = 0 + num_jobs = len(q) + + error_args_kwargs = [] + + s = '=== Starting {} jobs on devices {} with {} jobs per device ===' + self._log(s.format( + num_jobs, + ', '.join((str(i) for i in self.visible_devices)), + self.num_max_jobs)) + + while True: + + for gpu_id, (id, args_kwargs) in zip(self._iter_free_devices(), q): + + task = _GPUTask(self.fn, args_kwargs, gpu_id) + proc = _ErrorCatchingProcess(group=None, target=task) + proc.gpu_id = gpu_id + proc.args_kwargs_id = id + + self.counter[gpu_id] += 1 + q.remove((id, args_kwargs)) + processes.append(proc) + + proc.start() + time.sleep(1.0) # to make sure startup phases of process do not + # interleave + + for p in processes: + + if p.exitcode is not None: + processes.remove(p) + self.counter[p.gpu_id] -= 1 + finished_jobs += 1 + print('Finished job {}/{}'.format(finished_jobs, num_jobs)) + + if p.exception is not None: + self._log('=== ERROR (job {}) ==='.format( + p.args_kwargs_id)) + self._log(p.exception) + self._log('=============') + error_args_kwargs.append( + (p.exception, p.args_kwargs_id)) + + if len(processes) == 0: + break + + time.sleep(0.25) + print(next(self.progress_wheel), end='\r') + + return error_args_kwargs + + +def scatter_fn_on_devices( + fn: Callable, + fn_args_kwargs: list, + visible_devices: list, + max_process_on_device: int): + + assert isinstance(fn_args_kwargs, list) + assert isinstance(visible_devices, list) + assert all((i < torch.cuda.device_count() for i in visible_devices)) + + d = Distributor( + fn, + fn_args_kwargs, + visible_devices, + max_process_on_device + ) + + return d.run() + + +def keychain_value_iter(d, key_chain=None, allowed_values=None): + key_chain = [] if key_chain is None else list(key_chain).copy() + + if not isinstance(d, dict): + if allowed_values is not None: + assert isinstance(d, allowed_values), 'Value needs to be of type {}!'.format( + allowed_values) + yield key_chain, d + else: + for k, v in d.items(): + yield from keychain_value_iter( + v, + key_chain + [k], + allowed_values=allowed_values) + + +def configs_from_grid(grid): + tmp = list(keychain_value_iter(grid, allowed_values=(list, tuple))) + values = [x[1] for x in tmp] + key_chains = [x[0] for x in tmp] + + ret = [] + + for v in itertools.product(*values): + + ret_i = {} + + for kc, kc_v in zip(key_chains, v): + tmp = ret_i + for k in kc[:-1]: + if k not in tmp: + tmp[k] = {} + + tmp = tmp[k] + + tmp[kc[-1]] = kc_v + + ret.append(ret_i) + + return ret diff --git a/gril/torchph/chofer_torchex/utils/trainer/__init__.py b/gril/torchph/chofer_torchex/utils/trainer/__init__.py new file mode 100644 index 0000000..260e4c8 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/__init__.py @@ -0,0 +1 @@ +from .trainer import Trainer diff --git a/gril/torchph/chofer_torchex/utils/trainer/plugins/__init__.py b/gril/torchph/chofer_torchex/utils/trainer/plugins/__init__.py new file mode 100644 index 0000000..7962bd7 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/plugins/__init__.py @@ -0,0 +1,5 @@ +# TODO: Possibly delete + +from .training_schedule import LearningRateScheduler +from .scores import PredictionMonitor, LossMonitor +from .progress import ConsoleBatchProgress diff --git a/gril/torchph/chofer_torchex/utils/trainer/plugins/plugin.py b/gril/torchph/chofer_torchex/utils/trainer/plugins/plugin.py new file mode 100644 index 0000000..3b0c88d --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/plugins/plugin.py @@ -0,0 +1,3 @@ +class Plugin(object): + def register(self, trainer): + raise NotImplementedError() diff --git a/gril/torchph/chofer_torchex/utils/trainer/plugins/progress.py b/gril/torchph/chofer_torchex/utils/trainer/plugins/progress.py new file mode 100644 index 0000000..e9fad7e --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/plugins/progress.py @@ -0,0 +1,40 @@ + +from .plugin import Plugin + + +class ConsoleBatchProgress(Plugin): + def __init__(self, print_to_std_out=True): + super(ConsoleBatchProgress, self).__init__() + self._print_to_std_out = print_to_std_out + self.n_batches = None + self.ith_batch = None + self.max_epochs = None + self.current_epoch_number = None + self.max_batches = None + + self._blanks = 50*' ' + + def register(self, trainer): + trainer.events.pre_epoch.append(self.pre_epoch_handler) + trainer.events.post_batch.append(self.post_batch_handler) + trainer.events.post_epoch.append(self.post_epoch_handler) + + def pre_epoch_handler(self, **kwargs): + self.max_batches = len(kwargs['train_data']) + self.max_epochs = kwargs['max_epochs'] + self.current_epoch_number = kwargs['current_epoch_number'] + + def post_batch_handler(self, **kwargs): + current_batch_number = kwargs['current_batch_number'] + print(self._blanks, end='\r') + + str = """Epoch {}/{} Batch {}/{} ({:.2f} %)""".format(self.current_epoch_number, + self.max_epochs, + current_batch_number, + self.max_batches, + 100*current_batch_number/self.max_batches) + + print(str, end='\r') + + def post_epoch_handler(self, **kwargs): + print('') diff --git a/gril/torchph/chofer_torchex/utils/trainer/plugins/scores.py b/gril/torchph/chofer_torchex/utils/trainer/plugins/scores.py new file mode 100644 index 0000000..0e9cf4a --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/plugins/scores.py @@ -0,0 +1,97 @@ +from .plugin import Plugin +from torch.utils.data.dataloader import DataLoader +from sklearn.metrics import accuracy_score, confusion_matrix +from collections import OrderedDict +from torch.autograd import Variable + + +class PredictionMonitor(Plugin): + def __init__(self, test_data: DataLoader, eval_every_n_epochs=1, verbose=False, variable_created_by_model=False): + super(PredictionMonitor, self).__init__() + + self.eval_ever_n_epochs = eval_every_n_epochs + self._test_data = test_data + self.verbose = verbose + self.variable_created_by_model = variable_created_by_model + + self.accuracies = OrderedDict() + self.confusion_matrices = OrderedDict() + + self.evaluated_this_epoch = False + + def register(self, trainer): + trainer.events.post_epoch.append(self.post_epoch_handler) + + def post_epoch_handler(self, **kwargs): + epoch_count = kwargs['epoch_count'] + + if epoch_count % self.eval_ever_n_epochs != 0: + self.evaluated_this_epoch = False + return + + else: + trainer = kwargs['trainer'] + model = kwargs['model'] + model.eval() + + target_list = [] + predictions_list = [] + + if self.verbose: + print('testing...', end=' ') + + for batch_input, target in self._test_data: + + batch_input, target = trainer.data_typing(batch_input, target) + + if not self.variable_created_by_model: + batch_input = Variable(batch_input) + + output = model(batch_input).data + predictions = output.max(1)[1].type_as(target) + + target_list += target.view(-1).tolist() + predictions_list += predictions.view(-1).tolist() + + self.accuracies[epoch_count] = accuracy_score(target_list, predictions_list) + self.confusion_matrices[epoch_count] = confusion_matrix(target_list, predictions_list) + self.evaluated_this_epoch = True + + if self.verbose: + print(str(self)) + + @staticmethod + def last_of_orddict(ordered_dict): + return ordered_dict[next(reversed(ordered_dict))] + + @property + def accuracy(self): + return self.last_of_orddict(self.accuracies) + + @property + def confusion_matrix(self): + return self.last_of_orddict(self.confusion_matrices) + + def __str__(self): + if self.evaluated_this_epoch: + return """Accuracy: {:.2f} %""".format(100*self.accuracy) + else: + return "Accuracy: na" + + +class LossMonitor(Plugin): + def __init__(self): + super(LossMonitor, self).__init__() + self.losses_by_epoch = [] + self._losses_current_epoch = [] + + def register(self, trainer): + trainer.events.post_batch_backward.append(self.post_batch_backward_handler) + trainer.events.post_epoch.append(self.post_epoch_handler) + + def post_batch_backward_handler(self, **kwargs): + self._losses_current_epoch.append(kwargs['loss']) + + def post_epoch_handler(self, **kwargs): + self.losses_by_epoch.append(self._losses_current_epoch) + self._losses_current_epoch = [] diff --git a/gril/torchph/chofer_torchex/utils/trainer/plugins/training_schedule.py b/gril/torchph/chofer_torchex/utils/trainer/plugins/training_schedule.py new file mode 100644 index 0000000..f263ce6 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/plugins/training_schedule.py @@ -0,0 +1,38 @@ +from .plugin import Plugin +from collections import OrderedDict + + +class LearningRateScheduler(Plugin): + def __init__(self, learning_rate_shedule: callable, verbose=False): + super(LearningRateScheduler, self).__init__() + self.learning_rate_schedule = learning_rate_shedule + self.learning_rate_by_epoch = OrderedDict() + self.optimizer = None + self.verbose = verbose + + def register(self, trainer): + self.optimizer = trainer.optimizer + trainer.events.pre_epoch.append(self.pre_epoch_handler) + + def set_optimizer_to_new_lr(self, lr, optimizer): + + for param_group in optimizer.param_groups: + param_group['lr'] = lr + + def pre_epoch_handler(self, **kwargs): + if kwargs['epoch_count'] == 1: + return + + optimizer = kwargs['optimizer'] + + new_learning_rate = self.learning_rate_schedule(self, **kwargs) + + if new_learning_rate is not None: + self.set_optimizer_to_new_lr(new_learning_rate, optimizer) + self.learning_rate_by_epoch[kwargs['epoch_count']] = new_learning_rate + + if self.verbose: + print('lr -> {}'.format(new_learning_rate)) + + def __str__(self): + return 'lr: {}'.format(self.learning_rate_by_epoch[reversed(self.learning_rate_by_epoch).__next__()]) \ No newline at end of file diff --git a/gril/torchph/chofer_torchex/utils/trainer/trainer.py b/gril/torchph/chofer_torchex/utils/trainer/trainer.py new file mode 100644 index 0000000..6d5b230 --- /dev/null +++ b/gril/torchph/chofer_torchex/utils/trainer/trainer.py @@ -0,0 +1,163 @@ +import torch +from torch.autograd import Variable +from torch.nn import Module +from torch.optim import Optimizer +from torch.utils.data import DataLoader +from typing import Callable +from torch import Tensor + + + +class Event: + def __init__(self): + self._callbacks = [] + + def __call__(self, default_kwargs: {}, **kwargs): + kwargs.update(default_kwargs) + for callback in self._callbacks: + callback(**kwargs) + + def append(self, callback): + if not callable(callback): + raise ValueError('Expected callable.') + + self._callbacks.append(callback) + + +class TrainerEvents: + def __init__(self): + self.pre_run = Event() + + self.pre_epoch = Event() + self.post_epoch = Event() + self.post_batch_backward = Event() + self.pre_batch = Event() + self.post_batch = Event() + + self.post_epoch_gui = Event() + + self.post_run = Event() + + +class Trainer(object): + def __init__(self, model: Module, loss: Callable, optimizer: Optimizer, train_data: DataLoader, n_epochs, + cuda=False, + cuda_device_id=None, + variable_created_by_model=False): + self.n_epochs = n_epochs + self.model = model + self.criterion = loss + self.optimizer = optimizer + self.train_data = train_data + self.epoch_count = 1 + self.cuda = cuda + self.cuda_device_id = cuda_device_id + self.variable_created_by_model = variable_created_by_model + + self.return_value = {} + + self.events = TrainerEvents() + + def _get_default_event_kwargs(self): + return { + 'trainer': self, + 'model': self.model, + 'epoch_count': self.epoch_count, + 'cuda': self.cuda + } + + @property + def iteration_count(self): + return self.batch_count * self.epoch_count + + def register_plugin(self, plugin): + plugin.register(self) + + def run(self): + self.events.pre_run(self._get_default_event_kwargs()) + + if self.cuda: + self.model.cuda(self.cuda_device_id) + + for i in range(1, self.n_epochs + 1): + self.epoch_count = i + + self.events.pre_epoch(self._get_default_event_kwargs(), + optimizer=self.optimizer, + train_data=self.train_data, + max_epochs=self.n_epochs, + current_epoch_number=i) + self._train_epoch() + self.events.post_epoch(self._get_default_event_kwargs(), trainer=self) + self.events.post_epoch_gui(self._get_default_event_kwargs()) + + self.events.post_run(self._get_default_event_kwargs()) + return self.return_value + + def _train_epoch(self): + self.model.train() + + for i, (batch_input, batch_target) in enumerate(self.train_data, start=1): + + self.events.pre_batch(self._get_default_event_kwargs(), + batch_input=batch_input, + batch_target=batch_target) + + batch_input, batch_target = self.data_typing(batch_input, batch_target) + + target_var = Variable(batch_target) + + if not self.variable_created_by_model: + batch_input = Variable(batch_input) + + def closure(): + batch_output = self.model(batch_input) + loss = self.criterion(batch_output, target_var) + loss.backward() + + assert len(loss.data) == 1 + self.events.post_batch_backward(self._get_default_event_kwargs(), + batch_output=batch_output, + loss=float(loss.data[0])) + + return loss + + self.optimizer.zero_grad() + self.optimizer.step(closure) + + self.events.post_batch(self._get_default_event_kwargs(), + batch_input=batch_input, + batch_target=batch_target, + current_batch_number=i) + + @staticmethod + def before_data_typing_hook(batch_input, batch_targets): + return batch_input, batch_targets + + @staticmethod + def after_data_typing_hook(batch_input, batch_targets): + return batch_input, batch_targets + + def data_typing(self, batch_input, batch_targets): + batch_input, batch_targets = self.before_data_typing_hook(batch_input, batch_targets) + + tensor_cast = Tensor.cuda if self.cuda else Tensor.cpu + + def cast(x): + if isinstance(x, torch.tensor._TensorBase): + return tensor_cast(x) + elif isinstance(x, list): + return [cast(v) for v in x] + elif isinstance(x, dict): + return {k: cast(v) for k, v in x.items()} + elif isinstance(x, tuple): + return tuple(cast(v) for v in x) + else: + return x + + batch_input = cast(batch_input) + batch_targets = cast(batch_targets) + + batch_input, batch_targets = self.after_data_typing_hook(batch_input, batch_targets) + + return batch_input, batch_targets diff --git a/gril/torchph/docs_src/source/_static/.gitignore b/gril/torchph/docs_src/source/_static/.gitignore new file mode 100644 index 0000000..5e7d273 --- /dev/null +++ b/gril/torchph/docs_src/source/_static/.gitignore @@ -0,0 +1,4 @@ +# Ignore everything in this directory +* +# Except this file +!.gitignore diff --git a/gril/torchph/docs_src/source/conf.py b/gril/torchph/docs_src/source/conf.py new file mode 100644 index 0000000..81d92c3 --- /dev/null +++ b/gril/torchph/docs_src/source/conf.py @@ -0,0 +1,201 @@ +# -*- coding: utf-8 -*- +# +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +# import os +# import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- Project information ----------------------------------------------------- + +project = 'chofer_torchex' +copyright = '2019, Christoph D. Hofer, Roland Kwitt' +author = 'Christoph D. Hofer, Roland Kwitt' + +# The short X.Y version +version = '' +# The full version, including alpha/beta/rc tags +release = '0.0.1' + + +# -- General configuration --------------------------------------------------- + +# If your documentation needs a minimal Sphinx version, state it here. +# +# needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'nbsphinx', + 'sphinx.ext.autosummary', + 'sphinx.ext.coverage', + 'sphinx.ext.napoleon', + 'sphinx.ext.autodoc', + 'sphinx.ext.intersphinx', + 'sphinx.ext.mathjax', + 'sphinx.ext.viewcode', + 'sphinx.ext.githubpages', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', '**.ipynb_checkpoints'] + + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = None + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' # 'alabaster' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Custom sidebar templates, must be a dictionary that maps document names +# to template names. +# +# The default sidebars (for documents that don't match any pattern) are +# defined by theme itself. Builtin themes are using these templates by +# default: ``['localtoc.html', 'relations.html', 'sourcelink.html', +# 'searchbox.html']``. +# +# html_sidebars = {} + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = 'chofer_torchexdoc' + + +# -- Options for LaTeX output ------------------------------------------------ + +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + # + # 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + # + # 'pointsize': '10pt', + + # Additional stuff for the LaTeX preamble. + # + # 'preamble': '', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'chofer_torchex.tex', 'chofer\\_torchex Documentation', + 'Christoph D. Hofer', 'manual'), +] + + +# -- Options for manual page output ------------------------------------------ + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'chofer_torchex', 'chofer_torchex Documentation', + [author], 1) +] + + +# -- Options for Texinfo output ---------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'chofer_torchex', 'chofer_torchex Documentation', + author, 'chofer_torchex', 'One line description of project.', + 'Miscellaneous'), +] + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ['search.html'] + + +# -- Extension configuration ------------------------------------------------- + +# -- Options for intersphinx extension --------------------------------------- + +# Example configuration for intersphinx: refer to the Python standard library. +intersphinx_mapping = {'https://docs.python.org/': None} + + +# -- Customization chofer ---------------------------------------------------- + +def skip(app, what, name, obj, would_skip, options): + if name == "__init__": + return False + return would_skip + +def setup(app): + app.connect("autodoc-skip-member", skip) diff --git a/gril/torchph/docs_src/source/tutorials/.gitignore b/gril/torchph/docs_src/source/tutorials/.gitignore new file mode 100644 index 0000000..763513e --- /dev/null +++ b/gril/torchph/docs_src/source/tutorials/.gitignore @@ -0,0 +1 @@ +.ipynb_checkpoints diff --git a/gril/torchph/docs_src/source/tutorials/shared_code.py b/gril/torchph/docs_src/source/tutorials/shared_code.py new file mode 100644 index 0000000..47a92c4 --- /dev/null +++ b/gril/torchph/docs_src/source/tutorials/shared_code.py @@ -0,0 +1,22 @@ +import sys +import os + + +def check_chofer_torchex_availability(): + try: + import chofer_torchex + + except ImportError: + sys.path.append(os.path.dirname(os.getcwd())) + + try: + import chofer_torchex + + except ImportError as ex: + raise ImportError( + """ + Could not import chofer_torchex. Running your python \ + interpreter in the 'tutorials' sub folder could resolve \ + this issue. + """ + ) from ex diff --git a/gril/torchph/pershom_dev/bug_1.py b/gril/torchph/pershom_dev/bug_1.py new file mode 100644 index 0000000..1c5bbaf --- /dev/null +++ b/gril/torchph/pershom_dev/bug_1.py @@ -0,0 +1,64 @@ +import os +import os.path as pth +import pickle +from simplicial_complex import * +from cpu_sorted_boundary_array_implementation import SortedListBoundaryMatrix +import torch +from time import time +from pershombox import toplex_persistence_diagrams +import chofer_torchex.pershom.pershom_backend as pershom_backend +import yep +from chofer_torchex.pershom.calculate_persistence import calculate_persistence +import cProfile +from collections import Counter + +os.environ['CUDA_LAUNCH_BLOCKING'] = str(1) + + +def test(): + c = None + + with open('bug_1_sp.pickle', 'br') as f: + c = pickle.load(f) + print('|C| = ', len(c)) + max_red_by_iteration = 10000 + + # cpu_impl = SortedListBoundaryMatrix(c) + # cpu_impl.max_pairs = max_red_by_iteration + bm, col_dim = descending_sorted_boundary_array_from_filtrated_sp(c) + bm, col_dim = bm.to('cuda'), col_dim.to('cuda') + + + # barcodes_true = toplex_persistence_diagrams(c, list(range(len(c)))) + # dgm_true = [Counter(((float(b), float(d)) for b, d in dgm )) for dgm in barcodes_true] + + + def my_output_to_dgms(input): + ret = [] + b, b_e = input + + for dim, (b_dim, b_dim_e) in enumerate(zip(b, b_e)): + b_dim, b_dim_e = b_dim.float(), b_dim_e.float() + + tmp = torch.empty_like(b_dim_e) + tmp.fill_(float('inf')) + b_dim_e = torch.cat([b_dim_e, tmp], dim=1) + + + dgm = torch.cat([b_dim, b_dim_e], dim=0) + dgm = dgm.tolist() + dgm = Counter(((float(b), float(d)) for b, d in dgm )) + + ret.append(dgm) + + return ret + + output = calculate_persistence(bm, col_dim, max(col_dim), max_red_by_iteration) + + + + + + + +test() diff --git a/gril/torchph/pershom_dev/calc_pers.py b/gril/torchph/pershom_dev/calc_pers.py new file mode 100644 index 0000000..4facfc6 --- /dev/null +++ b/gril/torchph/pershom_dev/calc_pers.py @@ -0,0 +1,15 @@ +import torch +import chofer_torchex.pershom.pershom_backend as pershom_backend + +device = torch.device('cuda') +dtype = torch.int64 + +ba = torch.empty([0, 2], device=device, dtype=dtype) +ba_row_i_to_bm_col_i = ba +simplex_dimension = torch.zeros(10, device=device, dtype=dtype) +print(simplex_dimension) + +# ret = pershom_backend.__C.CalcPersCuda__calculate_persistence(ba, ba_row_i_to_bm_col_i, simplex_dimension, 3, -1) +ret = pershom_backend.__C.CalcPersCuda__my_test_f(ba); + +print(ret) \ No newline at end of file diff --git a/gril/torchph/pershom_dev/cpu_sorted_boundary_array_implementation.py b/gril/torchph/pershom_dev/cpu_sorted_boundary_array_implementation.py new file mode 100644 index 0000000..b4d87f4 --- /dev/null +++ b/gril/torchph/pershom_dev/cpu_sorted_boundary_array_implementation.py @@ -0,0 +1,114 @@ +from collections import Counter, defaultdict +import numpy as np +from simplicial_complex import boundary_operator + +class SortedListBoundaryMatrix: + + def __init__(self, filtrated_sc): + self._data = None + self._simplex_dims = None + self.stats = defaultdict(list) + self._init_from_filtrated_sp(filtrated_sc) + self.max_pairs = 100 + + def _init_from_filtrated_sp(self, filtrated_sp): + simplex_to_ordering_position = {s: i for i, s in enumerate(filtrated_sp)} + + self._data = [[]]*len(filtrated_sp) + self._simplex_dims = {} + + for col_i, s in enumerate(filtrated_sp): + boundary = boundary_operator(s) + orderings_of_boundaries = sorted((simplex_to_ordering_position[b] for b in boundary)) + self._data[col_i] = orderings_of_boundaries + self._simplex_dims[col_i] = len(s) - 1 + + + + def add_column_i_to_j(self, i, j): + col_i = self._data[i] + col_j = self._data[j] + + self._data[j] = sorted((k for k, v in Counter(col_i + col_j).items() if v == 1)) + + def get_pivot(self): + return {col_i: col_i_entries[-1] for col_i, col_i_entries in enumerate(self._data) + if len(col_i_entries) > 0} + + def pairs_for_reduction(self): + tmp = defaultdict(list) + for column_i, pivot_i in self.get_pivot().items(): + tmp[pivot_i].append(column_i) + + ret = [] + for k, v in tmp.items(): + if len(v) > 1: + for j in v[1:]: + ret.append((v[0], j)) + + if self.max_pairs is not None: + ret = ret[:self.max_pairs] + return ret + + def reduction_step(self): + pairs = self.pairs_for_reduction() + + if len(pairs) == 0: + raise StopIteration() + + self.stats['n_pairs'].append(len(pairs)) + + for i, j in pairs: + self.add_column_i_to_j(i, j) + + self.stats['longest_column'].append(max(len(c) for c in self._data)) + + def reduce(self): + iterations = 0 + try: + while True: + self.reduction_step() + iterations += 1 + + except StopIteration: + print('Reached end of reduction after ', iterations, ' iterations') + + def row_i_contains_lowest_one(self): + pass + + def read_barcodes(self): + assert len(self.pairs_for_reduction()) == 0, 'Matrix is not reduced' + + barcodes = defaultdict(list) + pivot = self.get_pivot() + + for death, birth in pivot.items(): + assert birth < death + barcodes[self._simplex_dims[birth]].append((birth, death)) + + rows_with_lowest_one = set() + for k, v in pivot.items(): + rows_with_lowest_one.add(v) + + + for i in range(len(self)): + if (i not in pivot) and i not in rows_with_lowest_one: + barcodes[self._simplex_dims[i]].append((i, float('inf'))) + + + return barcodes + + + def __repr__(self): + n = len(self._data) + matrix = np.zeros((n, n)) + + for i, col_i in enumerate(self._data): + for j in col_i: + matrix[j, i] = 1 + + + return str(matrix) + + def __len__(self): + return len(self._data) diff --git a/gril/torchph/pershom_dev/generate_profile_data.py b/gril/torchph/pershom_dev/generate_profile_data.py new file mode 100644 index 0000000..c3f344e --- /dev/null +++ b/gril/torchph/pershom_dev/generate_profile_data.py @@ -0,0 +1,50 @@ +import sys +import os.path as pth + +import torch +import numpy as np + +from simplicial_complex import * +from cpu_sorted_boundary_array_implementation import SortedListBoundaryMatrix + +def generate_profile_input_dump(): + c = random_simplicial_complex(1000, 2000, 4000) + max_dimension = 2 + + bm, simplex_dim = descending_sorted_boundary_array_from_filtrated_sp(c) + np.savetxt('profiling_pershom/data/boundary_array.np_txt', bm) + np.savetxt('profiling_pershom/data/boundary_array.np_txt', simplex_dim) + + ind_not_reduced = torch.tensor(list(range(simplex_dim.size(0)))) + ind_not_reduced = ind_not_reduced.masked_select(bm[:, 0] >= 0) + + bm = bm.index_select(0, ind_not_reduced) + + with open('profiling_pershom/data/boundary_array.txt', 'w') as f: + for row in bm: + for v in row: + f.write(str(int(v))) + f.write(" ") + + with open ('profiling_pershom/data/boundary_array_size.txt', 'w') as f: + for v in bm.size(): + f.write(str(int(v))) + f.write(" ") + + with open ('profiling_pershom/data/ind_not_reduced.txt', 'w') as f: + for v in ind_not_reduced: + f.write(str(int(v))) + f.write(" ") + + with open ('profiling_pershom/data/column_dimension.txt', 'w') as f: + for v in simplex_dim: + f.write(str(int(v))) + f.write(" ") + + with open ('profiling_pershom/data/max_dimension.txt', 'w') as f: + f.write(str(int(max_dimension))) + + + +if __name__ == "__main__": + generate_profile_input_dump() \ No newline at end of file diff --git a/gril/torchph/pershom_dev/generate_test_data.py b/gril/torchph/pershom_dev/generate_test_data.py new file mode 100644 index 0000000..3f986d4 --- /dev/null +++ b/gril/torchph/pershom_dev/generate_test_data.py @@ -0,0 +1,46 @@ +import os +import os.path as pth +import pickle +from simplicial_complex import * +from cpu_sorted_boundary_array_implementation import SortedListBoundaryMatrix +import torch +from pershombox import toplex_persistence_diagrams +from collections import Counter + + +def generate(): + + output_path = "/tmp/random_simplicial_complexes" + try: + os.mkdir(output_path) + except: + pass + + args = [(100, 200, 300), + (100, 200, 300, 400), + (100, 100, 100, 100, 100)] + + + for arg in args: + c = random_simplicial_complex(*arg) + file_name_stub = os.path.join(output_path, "random_sp__args__" + "_".join(str(x) for x in arg)) + + bm, col_dim = descending_sorted_boundary_array_from_filtrated_sp(c) + bm, col_dim = bm.to('cuda'), col_dim.to('cuda') + + ind_not_reduced = torch.tensor(list(range(col_dim.size(0)))).to('cuda') + ind_not_reduced = ind_not_reduced.masked_select(bm[:, 0] >= 0).long() + bm = bm.index_select(0, ind_not_reduced) + + barcodes_true = toplex_persistence_diagrams(c, list(range(len(c)))) + dgm_true = [Counter(((float(b), float(d)) for b, d in dgm )) for dgm in barcodes_true] + + with open(file_name_stub + ".pickle", 'wb') as f: + pickle.dump({'calculate_persistence_args': (bm, ind_not_reduced, col_dim, max(col_dim)), + 'expected_result': dgm_true}, + f) + + + +if __name__ == "__main__": + generate() diff --git a/gril/torchph/pershom_dev/pershom_script.py b/gril/torchph/pershom_dev/pershom_script.py new file mode 100644 index 0000000..c2df6a0 --- /dev/null +++ b/gril/torchph/pershom_dev/pershom_script.py @@ -0,0 +1,118 @@ +import os +import os.path as pth +import pickle +from simplicial_complex import * +from cpu_sorted_boundary_array_implementation import SortedListBoundaryMatrix +import torch +from time import time +# from pershombox import toplex_persistence_diagrams +import chofer_torchex.pershom.pershom_backend as pershom_backend +import yep +from chofer_torchex.pershom.calculate_persistence import calculate_persistence +import cProfile +from collections import Counter + +# os.environ['CUDA_LAUNCH_BLOCKING'] = str(1) + + +def test(): + c = None + + use_cache = False + + if use_cache: + random_simplicial_complex_path = './random_simplicial_complex.pickle' + if pth.exists(random_simplicial_complex_path): + with open(random_simplicial_complex_path, 'br') as f: + c = pickle.load(f) + else: + c = random_simplicial_complex(100, 100, 100, 100, 100, 100) + with open(random_simplicial_complex_path, 'bw') as f: + pickle.dump(c, f) + else: + c = random_simplicial_complex(100, 100, 100, 100, 100) + + print('|C| = ', len(c)) + max_red_by_iteration = -1 + + # cpu_impl = SortedListBoundaryMatrix(c) + # cpu_impl.max_pairs = max_red_by_iteration + bm, col_dim = descending_sorted_boundary_array_from_filtrated_sp(c) + + print(bm[-1]) + + bm, col_dim = bm.to('cuda'), col_dim.to('cuda') + + + barcodes_true = toplex_persistence_diagrams(c, list(range(len(c)))) + dgm_true = [Counter(((float(b), float(d)) for b, d in dgm )) for dgm in barcodes_true] + + + def my_output_to_dgms(input): + ret = [] + b, b_e = input + + for dim, (b_dim, b_dim_e) in enumerate(zip(b, b_e)): + b_dim, b_dim_e = b_dim.float(), b_dim_e.float() + + tmp = torch.empty_like(b_dim_e) + tmp.fill_(float('inf')) + b_dim_e = torch.cat([b_dim_e, tmp], dim=1) + + + dgm = torch.cat([b_dim, b_dim_e], dim=0) + dgm = dgm.tolist() + dgm = Counter(((float(b), float(d)) for b, d in dgm )) + + ret.append(dgm) + + return ret + + + # pr = cProfile.Profile() + # pr.enable() + + ind_not_reduced = torch.tensor(list(range(col_dim.size(0)))).to('cuda').detach() + ind_not_reduced = ind_not_reduced.masked_select(bm[:, 0] >= 0).long().detach() + bm = bm.index_select(0, ind_not_reduced).detach() + + yep.start('profiling_pershom/profile.google-pprof') + + for i in range(10): + time_start = time() + output = pershom_backend.calculate_persistence(bm.clone(), ind_not_reduced.clone(), col_dim.clone(), max(col_dim)) + print(time() - time_start) + yep.stop() + + # pr.disable() + # pr.dump_stats('high_level_profile.cProfile') + + print([[len(x) for x in y] for y in output ]) + + dgm_test = my_output_to_dgms(output) + + print('dgm_true lengths:', [len(dgm) for dgm in dgm_true]) + print('dgm_test lengths:', [len(dgm) for dgm in dgm_test]) + + for dgm_test, dgm_true in zip(dgm_test, dgm_true): + assert(dgm_test == dgm_true) + + +def vr_l1_persistence_performance_test(): + pc = torch.randn(20,3) + pc = torch.tensor(pc, device='cuda', dtype=torch.float) + + max_dimension = 2 + max_ball_radius = 0 + + yep.start('profiling_pershom/profile.google-pprof') + time_start = time() + res = pershom_backend.__C.VRCompCuda__vr_persistence(pc, max_dimension, max_ball_radius, 'l1') + print(time() - time_start) + yep.stop() + print(res[0][0]) + + + + +vr_l1_persistence_performance_test() diff --git a/gril/torchph/pershom_dev/profiling_pershom/.gitignore b/gril/torchph/pershom_dev/profiling_pershom/.gitignore new file mode 100644 index 0000000..92946be --- /dev/null +++ b/gril/torchph/pershom_dev/profiling_pershom/.gitignore @@ -0,0 +1,2 @@ +*.google-pprof +profile diff --git a/gril/torchph/pershom_dev/refactor_vc_comp.py b/gril/torchph/pershom_dev/refactor_vc_comp.py new file mode 100644 index 0000000..b0c66eb --- /dev/null +++ b/gril/torchph/pershom_dev/refactor_vc_comp.py @@ -0,0 +1,40 @@ +import chofer_torchex.pershom.pershom_backend as pershom_backend +import torch +import time +from scipy.special import binom +from itertools import combinations + + +from collections import Counter + + +point_cloud = [(0, 0), (1, 0), (0, 0.5), (1, 1.5)] +point_cloud = torch.tensor(point_cloud, device='cuda', dtype=torch.float, requires_grad=True) + +def l1_norm(x, y): + return float((x-y).abs().sum()) + +testee = pershom_backend.__C.VRCompCuda__PointCloud2VR_factory("l1") + +args = testee(point_cloud, 2, 2) +print(args[3]) +args[3].sum().backward() +print(point_cloud.grad) + + + + + + + + + + + + + + + + +# print(c(torch.rand((10, 3), device='cuda'), 2, 0)) +# print(c.filtration_values_by_dim) \ No newline at end of file diff --git a/gril/torchph/pershom_dev/setup.py b/gril/torchph/pershom_dev/setup.py new file mode 100644 index 0000000..1e7c547 --- /dev/null +++ b/gril/torchph/pershom_dev/setup.py @@ -0,0 +1,14 @@ +from setuptools import setup +from torch.utils.cpp_extension import BuildExtension, CUDAExtension + +setup( + name='pershom_backend_C', + ext_modules=[ + CUDAExtension('pershom_backend_C_cuda', [ + 'pershom_cpp_src/pershom_cuda.cu', + 'pershom_cpp_src/pershom.cpp', + ]) + ], + cmdclass={ + 'build_ext': BuildExtension + }) \ No newline at end of file diff --git a/gril/torchph/pershom_dev/simplicial_complex.py b/gril/torchph/pershom_dev/simplicial_complex.py new file mode 100644 index 0000000..fe09ed0 --- /dev/null +++ b/gril/torchph/pershom_dev/simplicial_complex.py @@ -0,0 +1,82 @@ +import numpy as np +import torch +from collections import defaultdict +from itertools import combinations +from scipy.special import binom + + +def boundary_operator(simplex): + s = tuple(simplex) + + if len(simplex) == 1: + return () + + else: + return [s[:i] + s[(i + 1):] for i in range(len(s))] + + +def random_simplicial_complex(*args): + simplices_by_dim = defaultdict(set) + + vertex_ids = np.array((range(args[0]))) + + vertices = [(i,) for i in vertex_ids] + simplices_by_dim[0] = set(vertices) + + for dim_i, n_simplices in enumerate(args[1:], start=1): + + n_simplices = min(n_simplices, int(binom(len(vertices), dim_i+1))) + while len(simplices_by_dim[dim_i]) != n_simplices: + chosen_vertices = np.random.choice(vertex_ids, replace=False, size=dim_i+1) + simplex = tuple(sorted(chosen_vertices)) + + if simplex not in simplices_by_dim[dim_i]: + simplices_by_dim[dim_i].add(simplex) + + for dim_i in sorted(simplices_by_dim, reverse=True): + if dim_i == 0 : + break + + for s in simplices_by_dim[dim_i]: + for boundary_s in boundary_operator(s): + if boundary_s not in simplices_by_dim[dim_i-1]: + simplices_by_dim[dim_i-1].add(boundary_s) + + sp = [] + for dim_i in range(len(args)): + sp += list(simplices_by_dim[dim_i]) + + return sp + + +def descending_sorted_boundary_array_from_filtrated_sp(filtrated_sp, + dtype=torch.int32, + resize_factor=2): + simplex_to_ordering_position = {s: i for i, s in enumerate(filtrated_sp)} + + max_boundary_size = max(len(s) for s in filtrated_sp) + n_cols = len(filtrated_sp) + n_rows = resize_factor*max_boundary_size + + bm = torch.empty(size=(n_rows, n_cols), + dtype=dtype) + bm.fill_(-1) + + col_to_dim = torch.empty(size=(n_cols,), + dtype=dtype) + + for col_i, s in enumerate(filtrated_sp): + boundary = boundary_operator(s) + orderings_of_boundaries = sorted((simplex_to_ordering_position[b] for b in boundary), + reverse=True) + + col_to_dim[col_i] = len(s) - 1 + + for row_i, entry in enumerate(orderings_of_boundaries): + bm[row_i, col_i] = entry + + # boundary array is delivered in column first order for efficency when merging + bm = bm.transpose_(0, 1) + + + return bm, col_to_dim diff --git a/gril/torchph/pershom_dev/vc_comp_script.py b/gril/torchph/pershom_dev/vc_comp_script.py new file mode 100644 index 0000000..7c1e69e --- /dev/null +++ b/gril/torchph/pershom_dev/vc_comp_script.py @@ -0,0 +1,87 @@ +import chofer_torchex.pershom.pershom_backend as pershom_backend +import torch +import time +from scipy.special import binom + + + #torch.tensor([[-0.6690, 1.5059], [ 0.4220, 1.2434], [-0.3436, -0.0053], [-0.1569, 0.0627]], device='cuda', requires_grad=True).float() + + +point_cloud = torch.randn(5,3, device='cuda', requires_grad=True).float() +# point_cloud = torch.tensor([(0, 0), (1, 0), (0, 0.5), (1, 1.5)], device='cuda', requires_grad=True) + +# loss = x.sum() +# loss.backward() +# print(point_cloud.grad) + + +print(point_cloud) +# pershom_backend.__C.CalcPersCuda__my_test_f(point_cloud).sum().backward(); + +# time_start = time.time() +try: + r = pershom_backend.__C.VRCompCuda__vr_l1_persistence(point_cloud, 0, 0) + +except Exception as ex: + print("=== Error ===") + print(ex) + exit() + +non_essentials = r[0] +essentials = r[1] + +print(len(non_essentials), len(essentials)) + +print("=== non-essentials ===") +for x in non_essentials: print(x) + +print("=== essentials ===") +for x in essentials: print(x) + +print("=== grad ===") +loss = non_essentials[0].sum() +loss.backward() +print(point_cloud.grad) + +# ba = r[0] +# simp_dim = r[2] +# filt_val = r[3] + +# dim = 3 +# for simp_id, boundaries in enumerate(ba): +# simp_filt_val = filt_val[simp_id + point_cloud.size(0)] + +# for boundary in boundaries.tolist(): +# if boundary == -1: continue +# cond = True +# cond = cond and (simp_dim[simp_id + point_cloud.size(0)] - 1 == simp_dim[boundary]) +# cond = cond and filt_val[boundary] <= simp_filt_val + +# if not cond: +# print("{}, {}".format(simp_id, boundary)) +# print(ba[simp_id]) +# print(simp_dim[simp_id + point_cloud.size(0)], simp_dim[boundary]) + +# raise Exception() + + + + + +# print(print(r[2][int(5+binom(5,2)):int(5+binom(5,2))+1000])) + + +# loss = t.sum() +# loss.backward() +# +# for i, r in enumerate(ba): +# x = point_cloud[r[0]-point_cloud.size(0)] +# y = point_cloud[r[1]-point_cloud.size(0)] + +# assert torch.equal((x-y).abs().sum(), t[i]) + + +# for row_i in range(t.size(0)): +# for col_i in range(t.size(1)): +# assert int(t[row_i, col_i]) == binom(col_i, row_i+1) +# print(time.time() - time_start) diff --git a/gril/torchph/tests/.gitignore b/gril/torchph/tests/.gitignore new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/tests/__init__.py b/gril/torchph/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/tests/pershom/__init__.py b/gril/torchph/tests/pershom/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/gril/torchph/tests/pershom/test_VRCompCuda.py b/gril/torchph/tests/pershom/test_VRCompCuda.py new file mode 100644 index 0000000..4627715 --- /dev/null +++ b/gril/torchph/tests/pershom/test_VRCompCuda.py @@ -0,0 +1,444 @@ +import pytest +import torch +import chofer_torchex.pershom.pershom_backend as pershom_backend +import glob +import os +import numpy + +from itertools import combinations +from scipy.special import binom +from itertools import combinations +from collections import Counter + + +__C = pershom_backend.__C +EPSILON = 0.0001 + + +def test_l2_distance_matrix(): + point_cloud = [(0, 0), (1, 0), (0, 0.5), (1, 1.5)] + point_cloud = torch.tensor(point_cloud, device='cuda', dtype=torch.float) + + D = __C.VRCompCuda__l2_norm_distance_matrix(point_cloud) + + for i, j in combinations(range(point_cloud.size(0)), 2): + l2_dist = (point_cloud[i] - point_cloud[j]).pow(2).sum().sqrt() + l2_dist = float(l2_dist) + + assert float(D[i, j]) == l2_dist + + +def test_l1_distance_matrix(): + point_cloud = [(0, 0), (1, 0), (0, 0.5), (1, 1.5)] + point_cloud = torch.tensor(point_cloud, device='cuda', dtype=torch.float) + + D = __C.VRCompCuda__l1_norm_distance_matrix(point_cloud) + + for i, j in combinations(range(point_cloud.size(0)), 2): + l1_dist = (point_cloud[i] - point_cloud[j]).abs().sum() + l1_dist = float(l1_dist) + + assert float(D[i, j]) == l1_dist + + +def test_co_faces_from_combinations__ground_truth(): + faces = [(1, 0), (2, 0), (2, 1), (3, 0), (3, 1)] + comb = [sorted(cf, reverse=True) for cf in combinations(range(len(faces)), 3)] + + faces_t = torch.tensor(faces, device='cuda', dtype=torch.int64) + comb_t = torch.tensor(comb, device='cuda', dtype=torch.int64) + + result = __C.VRCompCuda__co_faces_from_combinations(comb_t, faces_t) + + expected_result = [[2, 1, 0], [4, 3, 0]] + + assert result.tolist() == expected_result + + +def test_co_faces_from_combinations__correct_size_if_just_one_co_face(): + faces = [(1, 0), (2, 0), (2, 1), (3, 0)] + comb = [sorted(cf, reverse=True) for cf in combinations(range(len(faces)), 3)] + + faces_t = torch.tensor(faces, device='cuda', dtype=torch.int64) + comb_t = torch.tensor(comb, device='cuda', dtype=torch.int64) + + result = __C.VRCompCuda__co_faces_from_combinations(comb_t, faces_t) + + assert result.dim() == 2 + assert result.size(0) == 1 + + expected_result = [[2, 1, 0]] + + assert result.tolist() == expected_result + + +def test_co_faces_from_combinations__no_co_face(): + faces = [(1, 0), (2, 0), (3, 0), (4, 0), (5, 0)] + comb = [sorted(cf, reverse=True) for cf in combinations(range(len(faces)), 3)] + + faces_t = torch.tensor(faces, device='cuda', dtype=torch.int64) + comb_t = torch.tensor(comb, device='cuda', dtype=torch.int64) + + result = __C.VRCompCuda__co_faces_from_combinations(comb_t, faces_t) + + expected_result = [] + + assert result.tolist() == expected_result + + +@pytest.mark.parametrize("n_vertices, dim_faces", [(10, 2), (8, 3)]) +def test_co_faces_from_combinations__result(n_vertices, dim_faces): + faces = [sorted(f) for f in combinations(range(n_vertices), dim_faces)] + comb = [sorted(cf) for cf in combinations(range(len(faces)), dim_faces + 1)] + + faces_t = torch.tensor(faces, device='cuda', dtype=torch.int64) + comb_t = torch.tensor(comb, device='cuda', dtype=torch.int64) + + expected_result = [] + for c in comb: + tmp = sum([faces[i] for i in c], []) + tmp = Counter(tmp) + + if all(v == 2 for v in tmp.values()): + expected_result.append(c) + + expected_result = {tuple(x) for x in expected_result} + + result = __C.VRCompCuda__co_faces_from_combinations(comb_t, faces_t) + + result = {tuple(x) for x in result.tolist()} + + assert expected_result == result + + +@pytest.mark.parametrize( + "max_dimension, max_ball_radius", + [(2, 0.0), (2, 1.0), (2, 2.0), (3, 0.0), (3, 1.0), (3, 2.0)]) +def test_VietorisRipsArgsGenerator__l1_excessive_state_testing_of_all_intermediate_steps(max_dimension, max_ball_radius): + # ASSERTS that the number of edges is > 0!!! + point_cloud = [(0, 0), (1, 0), (0, 0.5), (1, 1.5)] + point_cloud = torch.tensor(point_cloud, device='cuda', dtype=torch.float) + distance_matrix = pershom_backend.__C.VRCompCuda__l1_norm_distance_matrix(point_cloud) + + def l1_norm(x, y): + return float((x-y).abs().sum()) + + testee = pershom_backend.__C.VRCompCuda__VietorisRipsArgsGenerator() + # + # Test: init_state + # + testee.init_state(distance_matrix, max_dimension, max_ball_radius) + + assert len(testee.filtration_values_by_dim) == 1 + assert testee.n_simplices_by_dim[0] == point_cloud.size(0) + assert testee.filtration_values_by_dim[0].size(0) == point_cloud.size(0) + assert all(testee.filtration_values_by_dim[0] == 0) + + # + # Test: make_boundary_info_edges + # + testee.make_boundary_info_edges() + + assert len(testee.filtration_values_by_dim) == 2 + + expected_n_dim_1_simplices = 0 + expected_filtration_values = [] + expected_boundaries = [] + threshold = max_ball_radius if max_ball_radius > 0 else float('inf') + for i, j in combinations(range(point_cloud.size(0)), 2): + + norm = l1_norm(point_cloud[i], point_cloud[j]) + + if norm <= threshold: + expected_n_dim_1_simplices += 1 + expected_boundaries.append(sorted((i, j), reverse=True)) + + expected_boundaries = sorted(expected_boundaries) + + for i, j in expected_boundaries: + norm = l1_norm(point_cloud[i], point_cloud[j]) + expected_filtration_values.append(norm) + + expected_boundaries = torch.tensor( + expected_boundaries, + device='cuda', + dtype=torch.int64) + if expected_boundaries.numel() == 0: + expected_boundaries = torch.empty( + (0, 2), + device='cuda', + dtype=torch.int64) + + expected_filtration_values = torch.tensor( + expected_filtration_values, + device='cuda') + + assert testee.n_simplices_by_dim[1] == expected_n_dim_1_simplices + assert testee.boundary_info_non_vertices[0].size(0) == testee.filtration_values_by_dim[1].size(0) + assert testee.boundary_info_non_vertices[0].equal(expected_boundaries) + assert testee.filtration_values_by_dim[1].equal(expected_filtration_values) + assert testee.filtration_values_by_dim[1].dim() == 1 + + # + # Test: make_boundary_info_non_edges + # + testee.make_boundary_info_non_edges() + + for dim in range(2, max_dimension + 1): + n_dim_min_one_simplices = testee.n_simplices_by_dim[dim-1] + dim_min_one_filtration_values = testee.filtration_values_by_dim[dim-1] + dim_min_one_boundary_info = testee.boundary_info_non_vertices[dim-2].tolist() + expected_n_simplices = binom(n_dim_min_one_simplices, dim + 1) + + expected_boundaries = [] + expected_filtration_values = [] + + for boundaries in combinations(range(n_dim_min_one_simplices), dim + 1): + + boundaries = list(boundaries) + boundaries_of_boundaries = sum((dim_min_one_boundary_info[i] for i in boundaries), []) + c = Counter(boundaries_of_boundaries) + if all(v == 2 for v in c.values()): + expected_boundaries.append(tuple(sorted(boundaries, reverse=True))) + + expected_boundaries = sorted(expected_boundaries) + + for boundaries in expected_boundaries: + expected_filtration_values.append( + max(dim_min_one_filtration_values[b] for b in boundaries) + ) + + expected_filtration_values = torch.tensor( + expected_filtration_values, + device='cuda') + + assert expected_boundaries == [tuple(x) for x in testee.boundary_info_non_vertices[dim-1].tolist()] + assert expected_filtration_values.equal(testee.filtration_values_by_dim[dim]) + assert testee.n_simplices_by_dim[dim] == len(expected_boundaries) + + # + # Test: make_simplex_ids_compatible_within_dimensions + # + testee.make_simplex_ids_compatible_within_dimensions() + for dim, boundary_info in enumerate(testee.boundary_info_non_vertices[1:], start=2): + if boundary_info.numel() > 0: + min_id = sum(testee.n_simplices_by_dim[:dim-1]) + max_id = min_id + testee.n_simplices_by_dim[dim-1] + + simp_ids = boundary_info.view(-1).tolist() + + assert min_id <= min(simp_ids) + assert max_id >= max(simp_ids) + else: + assert testee.n_simplices_by_dim[dim] == 0 + + # + # Test: make_simplex_dimension_vector + # + testee.make_simplex_dimension_vector() + + tmp = testee.n_simplices_by_dim + for i in range(max_dimension+1): + a = sum(tmp[:i]) + b = sum(tmp[:i+1]) + + slice_i = testee.simplex_dimension_vector[a:b] + assert slice_i.numel() == 0 or (slice_i.equal(torch.empty_like(slice_i).fill_(i))) + + assert testee.simplex_dimension_vector.dim() == 1 + assert testee.simplex_dimension_vector.size(0) == sum(tmp) + + # + # Test: make_filtration_values_vector_without_vertices + # + testee.make_filtration_values_vector_without_vertices() + assert testee.filtration_values_vector_without_vertices.dim() == 1 + assert testee.filtration_values_vector_without_vertices.size(0) == sum(testee.n_simplices_by_dim[1:]) + + # + # Test: do_filtration_add_eps_hack + # + testee.do_filtration_add_eps_hack() + + tmp = testee.n_simplices_by_dim + for i, boundary_info in enumerate(testee.boundary_info_non_vertices): + if i == 0: continue + + id_offset = sum(tmp[1:i+1]) + + for j, boundaries in enumerate(boundary_info.tolist()): + id = id_offset + j + + for b in boundaries: + assert testee.filtration_values_vector_without_vertices[int(b) - point_cloud.size(0)] <= \ + testee.filtration_values_vector_without_vertices[id] + + # + # Test: make_sorting_infrastructure + # + testee.make_sorting_infrastructure() + assert testee.sort_indices_without_vertices.dim() == 1 + assert testee.sort_indices_without_vertices.size(0) == sum(testee.n_simplices_by_dim[1:]) + + tmp = testee.sort_indices_without_vertices_inverse.index_select( + 0, testee.sort_indices_without_vertices) + expected_tmp = torch.tensor(list(range(sum(testee.n_simplices_by_dim[1:]))), device='cuda') + assert tmp.equal(expected_tmp) + + # Test + # undo_filtration_add_eps_hack + # + testee.undo_filtration_add_eps_hack() + + tmp = testee.n_simplices_by_dim + for i, boundary_info in enumerate(testee.boundary_info_non_vertices): + if i == 0: continue + + id_offset = sum(tmp[1:i+1]) + + for j, boundaries in enumerate(boundary_info.tolist()): + id = id_offset + j + + filt_b = [] + + for b in boundaries: + filt_b.append( + testee.filtration_values_vector_without_vertices[int(b) - point_cloud.size(0)]) + assert testee.filtration_values_vector_without_vertices[id] == max(filt_b) + + # + # Test: make_sorted_filtration_values_vector + # + testee.make_sorted_filtration_values_vector() + assert testee.sorted_filtration_values_vector.equal( + testee.sorted_filtration_values_vector.sort(0)[0] + ) + + # + # Test: make_boundary_array_rows_unsorted + # + testee.make_boundary_array_rows_unsorted() + + assert testee.boundary_array.size() \ + == (sum(testee.n_simplices_by_dim[1:]), 4 if max_dimension == 0 else 2*(max_dimension +1)) + + # + # Test: apply_sorting_to_rows + # + testee.apply_sorting_to_rows() + for i, row in enumerate(testee.boundary_array.tolist()): + simplex_id = i + testee.n_simplices_by_dim[0] + row = [x for x in row if x != -1] + simplex_dim = len(row) - 1 + assert simplex_dim == testee.simplex_dimension_vector[simplex_id] + b_filt_vals = [] + for b in row: + assert testee.simplex_dimension_vector[b] == simplex_dim - 1 + b_filt_vals.append(float(testee.sorted_filtration_values_vector[b])) + + + assert sorted(b_filt_vals, reverse=True) == b_filt_vals + + if simplex_dim > 1: + assert float(testee.sorted_filtration_values_vector[simplex_id]) == max(b_filt_vals) + + # + # Test: make_ba_row_i_to_bm_col_i_vector + # + testee.make_ba_row_i_to_bm_col_i_vector() + assert testee.ba_row_i_to_bm_col_i_vector.dim() == 1 + assert testee.ba_row_i_to_bm_col_i_vector.size(0) == sum(testee.n_simplices_by_dim[1:]) + + expected_ba_row_i_to_bm_col_i_vector = range(testee.n_simplices_by_dim[0], sum(testee.n_simplices_by_dim)) + expected_ba_row_i_to_bm_col_i_vector = list(expected_ba_row_i_to_bm_col_i_vector) + expected_ba_row_i_to_bm_col_i_vector = torch.tensor(expected_ba_row_i_to_bm_col_i_vector, device='cuda') + assert testee.ba_row_i_to_bm_col_i_vector.equal(expected_ba_row_i_to_bm_col_i_vector) + + +random_point_clouds = [] +for pth in glob.glob(os.path.join(os.path.dirname(__file__), "test_pershom_backend_data/random_point_clouds/*.txt")): + X = numpy.loadtxt(pth).tolist() + X = torch.tensor(X, device='cuda', dtype=torch.float) + random_point_clouds.append(X) + + +@pytest.mark.parametrize("test_args", random_point_clouds) +def test_VietorisRipsArgsGenerator__l1_valid_calculate_persistence_args(test_args): + point_cloud = test_args + point_cloud = point_cloud.float().to('cuda') + distance_matrix = pershom_backend.__C.VRCompCuda__l1_norm_distance_matrix(point_cloud) + max_dimension = 2 + max_ball_radius = 0 + + testee = pershom_backend.__C.VRCompCuda__VietorisRipsArgsGenerator() + ba, ba_row_i_to_bm_col_i, simplex_dimension, sorted_filtration_values_vector = testee( + distance_matrix, + max_dimension, + max_ball_radius) + + if ba.size(0) > 0: + + assert ba.dim() == 2 + + assert ba_row_i_to_bm_col_i.dim() == 1 + assert ba.size(0) == ba_row_i_to_bm_col_i.size(0) + n_simplices = point_cloud.size(0) + ba.size(0) + + assert simplex_dimension.dim() == 1 + assert simplex_dimension.size(0) == n_simplices + + assert sorted_filtration_values_vector.dim() == 1 + assert sorted_filtration_values_vector.size(0) == n_simplices + + # discard right half of columns, which are filled with -1 ... + content_column_border = (max_dimension if max_dimension != 0 else 1) + 1 + ba_right = ba[:, content_column_border:].contiguous() + ba_right = ba_right.view(ba_right.numel()) + ba = ba[:, :content_column_border].tolist() + # ba = [[x for x in row if x != -1] for row in ba] + + assert all(x == -1 for x in ba_right) + + simplex_dimension = simplex_dimension.tolist() + assert max(simplex_dimension) == max_dimension + + ba_row_i_to_bm_col_i = ba_row_i_to_bm_col_i.tolist() + sorted_filtration_values_vector = sorted_filtration_values_vector.tolist() + + for i, row_i in enumerate(ba): + simplex_id = i + point_cloud.size(0) + simplex_dim = simplex_dimension[simplex_id] + simplex_filt_val = sorted_filtration_values_vector[simplex_id] + + # sorted descendingly by id? + assert row_i == sorted(row_i, reverse=True) + boundary_ids = [x for x in row_i if x != -1] + + # dim n simplex has n+1 boundaries? + assert len(boundary_ids) == simplex_dim + 1 + boundary_filt_vals = [sorted_filtration_values_vector[x] for x in boundary_ids] + + # boundary filtration values are coherent with ids? + assert boundary_filt_vals == sorted(boundary_filt_vals, reverse=True) + + # simplex_filt_val is maximum of boundary filtration values? + if simplex_dim > 1 : + assert abs(max(boundary_filt_vals) - simplex_filt_val) < EPSILON + + # boundary of simplex has dim - 1 ? + boundary_dims = [simplex_dimension[b_id] for b_id in boundary_ids] + assert all(list(dim == simplex_dim - 1 for dim in boundary_dims)) + + if simplex_dim > 1: + boundary_of_boundaries = sum([ba[b_id - point_cloud.size(0)] for b_id in boundary_ids], []) + boundary_of_boundaries = [x for x in boundary_of_boundaries if x != -1] + + c = Counter(boundary_of_boundaries) + + assert all(list(v == 2 for v in c.values())) + # ba.size(0) == 0 + else: + + assert ba_row_i_to_bm_col_i.numel() == 0 + assert simplex_dimension.tolist() == [0]*point_cloud.size(0) + assert sorted_filtration_values_vector.tolist() == [0]*point_cloud.size(0) diff --git a/gril/torchph/tests/pershom/test_pershom_backend.py b/gril/torchph/tests/pershom/test_pershom_backend.py new file mode 100644 index 0000000..87f01e1 --- /dev/null +++ b/gril/torchph/tests/pershom/test_pershom_backend.py @@ -0,0 +1,169 @@ +import pytest +import torch +import glob +import pickle +import chofer_torchex.pershom.pershom_backend as pershom_backend +from collections import Counter + + +class Test_find_merge_pairings: + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_return_value_dtype(self, device, dtype): + pivots = torch.tensor([1, 1], device=device, dtype=dtype) + + result = pershom_backend.find_merge_pairings(pivots) + + assert result.dtype == torch.int64 + + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_parameter_max_pairs(self, device, dtype): + pivots = torch.tensor([1]*1000, device=device, dtype=dtype).unsqueeze(1) + + result = pershom_backend.find_merge_pairings(pivots, max_pairs=100) + + assert result.size(0) == 100 + + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_no_merge_pairs(self, device, dtype): + pivots = torch.tensor(list(range(100)), device=device, dtype=dtype).unsqueeze(1) + assert pershom_backend.find_merge_pairings(pivots).numel() == 0 + + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_result_1(self, device, dtype): + pivots = [6, 3, 3, 3, 5, 6, 6, 0, 5, 5] + pivots = torch.tensor(pivots, device=device, dtype=dtype).unsqueeze(1) + + result = pershom_backend.find_merge_pairings(pivots) + + assert result.dtype == torch.int64 + + expected_result = set([(0, 5), (0, 6), + (1, 2), (1, 3), + (4, 8), (4, 9)]) + + result = set(tuple(x) for x in result.tolist()) + + assert result == (expected_result) + + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_result_2(self, device, dtype): + pivots = sum([100*[i] for i in range(100)], []) + pivots = torch.tensor(pivots, device=device, dtype=dtype).unsqueeze(1) + + expected_result = torch.tensor([(int(i/100) * 100, i) for i in range(100 * 100) if i % 100 != 0]) + expected_result = expected_result.long().to('cuda') + + result = pershom_backend.find_merge_pairings(pivots) + + assert expected_result.equal(result) + + +class Test_calculate_persistence: + + @staticmethod + def calculate_persistence_output_to_barcode_list(input): + ret = [] + b, b_e = input + + for dim, (b_dim, b_dim_e) in enumerate(zip(b, b_e)): + b_dim, b_dim_e = b_dim.float(), b_dim_e.float() + + tmp = torch.empty_like(b_dim_e) + tmp.fill_(float('inf')) + b_dim_e = torch.cat([b_dim_e, tmp], dim=1) + + dgm = torch.cat([b_dim, b_dim_e], dim=0) + dgm = dgm.tolist() + dgm = Counter(((float(b), float(d)) for b, d in dgm)) + + ret.append(dgm) + + return ret + + @pytest.mark.parametrize("device, dtype", [ + (torch.device('cuda'), torch.int64) + ]) + def test_empty_input(self, device, dtype): + ba = torch.empty([0], device=device, dtype=dtype) + ba_row_i_to_bm_col_i = ba + simplex_dimension = torch.zeros(10, device=device, dtype=dtype) + + not_ess, ess = pershom_backend.calculate_persistence( + ba, + ba_row_i_to_bm_col_i, + simplex_dimension, + 2, + -1) + + for pairings in not_ess: + assert len(pairings) == 0 + + assert len(ess[0]) == 10 + + for birth_i in ess[1:]: + assert len(birth_i) == 0 + + def test_simple_1(self): + device = torch.device('cuda') + dtype = torch.int64 + + ba = torch.empty((1, 4)) + ba.fill_(-1) + ba[0, 0:2] = torch.tensor([1, 0]) + ind_not_reduced = torch.tensor([2]) + + row_dim = torch.tensor([0, 0, 1]) + + max_dim_to_read_of_reduced_ba = 1 + + ba = ba.to(device).type(dtype) + ind_not_reduced = ind_not_reduced.to(device).long() + row_dim = row_dim.to(device).type(dtype) + + out = pershom_backend.calculate_persistence( + ba, + ind_not_reduced, + row_dim, + max_dim_to_read_of_reduced_ba, + 100) + + barcodes = Test_calculate_persistence.calculate_persistence_output_to_barcode_list(out) + + assert barcodes[0] == Counter([(1.0, 2.0), (0.0, float('inf'))]) + + def test_random_simplicial_complexes(self): + device = torch.device('cuda') + dtype = torch.int64 + + for sp_path in glob.glob('test_pershom_backend_data/random_simplicial_complexes/*'): + + with open(sp_path, 'br') as f: + data = pickle.load(f) + + assert len(data) == 2 + + ba, ind_not_reduced, row_dim, max_dim_to_read_of_reduced_ba = data['calculate_persistence_args'] + expected_result = data['expected_result'] + + ba, ind_not_reduced, row_dim = ba.to(device).type(dtype), ind_not_reduced.to(device).long(), row_dim.to(device).type(dtype) + + result = pershom_backend.calculate_persistence( + ba, + ind_not_reduced, + row_dim, + max_dim_to_read_of_reduced_ba, + 10000) + result = Test_calculate_persistence.calculate_persistence_output_to_barcode_list(result) + + for dgm, dgm_exp in zip(result, expected_result): + assert dgm == dgm_exp diff --git a/zigzag/utils.h b/gril/utils.h similarity index 100% rename from zigzag/utils.h rename to gril/utils.h diff --git a/zigzag/__init__.py b/zigzag/__init__.py deleted file mode 100644 index 158fb9c..0000000 --- a/zigzag/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .zz import * diff --git a/zigzag/__pycache__/__init__.cpython-39.pyc b/zigzag/__pycache__/__init__.cpython-39.pyc deleted file mode 100644 index 09659ba..0000000 Binary files a/zigzag/__pycache__/__init__.cpython-39.pyc and /dev/null differ diff --git a/zigzag/__pycache__/zz.cpython-39.pyc b/zigzag/__pycache__/zz.cpython-39.pyc deleted file mode 100644 index d1f8d33..0000000 Binary files a/zigzag/__pycache__/zz.cpython-39.pyc and /dev/null differ diff --git a/zigzag/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so b/zigzag/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so deleted file mode 100755 index e12375f..0000000 Binary files a/zigzag/build/lib.linux-x86_64-cpython-39/mpml.cpython-39-x86_64-linux-gnu.so and /dev/null differ diff --git a/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_deps b/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_deps deleted file mode 100644 index cc0b40a..0000000 Binary files a/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_deps and /dev/null differ diff --git a/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_log b/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_log deleted file mode 100644 index daf6e36..0000000 --- a/zigzag/build/temp.linux-x86_64-cpython-39/.ninja_log +++ /dev/null @@ -1,3 +0,0 @@ -# ninja log v5 -0 19238 1674245376000000000 /scratch/bell/mukher26/mpml_graph_static/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o 9a2b7b576f7f16a -6 21939 1674668331865933337 /scratch1/mukher26/mpml_graph_repo/zigzag/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/zigzag/multipers.o a2a98d9a7402a0f5 diff --git a/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o b/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o deleted file mode 100644 index 82259ce..0000000 Binary files a/zigzag/build/temp.linux-x86_64-cpython-39/scratch/bell/mukher26/mpml_graph_static/zigzag/multipers.o and /dev/null differ diff --git a/zigzag/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/zigzag/multipers.o b/zigzag/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/zigzag/multipers.o deleted file mode 100644 index d7a68de..0000000 Binary files a/zigzag/build/temp.linux-x86_64-cpython-39/scratch1/mukher26/mpml_graph_repo/zigzag/multipers.o and /dev/null differ diff --git a/zigzag/dist/mpml-0.0.0-py3.9-linux-x86_64.egg b/zigzag/dist/mpml-0.0.0-py3.9-linux-x86_64.egg deleted file mode 100644 index f960199..0000000 Binary files a/zigzag/dist/mpml-0.0.0-py3.9-linux-x86_64.egg and /dev/null differ diff --git a/zigzag/multipers.h b/zigzag/multipers.h deleted file mode 100644 index 4c860cb..0000000 --- a/zigzag/multipers.h +++ /dev/null @@ -1,77 +0,0 @@ -#ifndef MULTIPERS_H -#define MULTIPERS_H -#include -#include -#include -#include -#include -#include -#include -#include "utils.h" - -#include "./phat/compute_persistence_pairs.h" - -using torch::Tensor; -using namespace torch::indexing; -typedef std::pair Point; - -class Multipers{ - private: - int hom_rank; - std::vector ranks; - double step, ll_x, ll_y; - int px, py; - int l; - - - void set_hom_rank(int hom_rank){ - this->hom_rank = hom_rank; - } - void set_ranks(std::vector ranks_){ - this->ranks.insert(this->ranks.begin(), ranks_.begin(), ranks_.end()); - } - void set_step(double step){ - this->step = step; - } - void set_l_for_worm(int l){ - this->l = l; - } - Tensor compute_l_worm(const int d); - std::vector> compute_filtration_along_boundary_cap(const Tensor& grid_pts_along_boundary_t, - const Tensor& f, - int &manual_birth_pts, - int &manual_death_pts); - - int zigzag_pairs(std::vector> &simplices_birth_death, - const vector &simplices, - const int manual_birth_pts, - const int manual_death_pts); - - int num_full_bars_for_specific_d(const Tensor& filtration, const vector& simplices, const Point& p, int d); - - Tensor find_maximal_worm_for_rank_k(const Tensor& filtration, const vector& simplices, const Point& p, const int rank); - - void set_grid_resolution_and_lower_left_corner(const Tensor& filtration); - - - - public: - int max_threads; - Multipers(const int hom_rank, const int l, double step, const std::vector ranks){ - set_hom_rank(hom_rank); - set_l_for_worm(l); - // set_division_along_axes(px, py); - set_step(step); - set_ranks(ranks); - this->max_threads = 1; - } - void set_max_jobs(int max_jobs){ - this->max_threads = max_jobs; - } - std::vector compute_landscape(const std::vector& pts, const std::vector>> &batch); - - - -}; - -#endif \ No newline at end of file diff --git a/zigzag/phat/algorithms/chunk_reduction.h b/zigzag/phat/algorithms/chunk_reduction.h deleted file mode 100644 index ed864f2..0000000 --- a/zigzag/phat/algorithms/chunk_reduction.h +++ /dev/null @@ -1,226 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - template class chunk_reduction_impl { - public: - enum column_type { GLOBAL - , LOCAL_POSITIVE - , LOCAL_NEGATIVE }; - - public: - template< typename Representation > - void operator() ( boundary_matrix< Representation >& boundary_matrix ) { - - - const index nr_columns = boundary_matrix.get_num_cols(); - if( omp_get_max_threads( ) > nr_columns ) - omp_set_num_threads( 1 ); - - const dimension max_dim = boundary_matrix.get_max_dim(); - - std::vector< index > lowest_one_lookup( nr_columns, -1 ); - std::vector < column_type > column_type( nr_columns, GLOBAL ); - std::vector< char > is_active( nr_columns, false ); - - const index chunk_size = use_sqrt ? (index)sqrt( (double)nr_columns ) : nr_columns / omp_get_max_threads(); - - std::vector< index > chunk_boundaries; - for( index cur_boundary = 0; cur_boundary < nr_columns; cur_boundary += chunk_size ) - chunk_boundaries.push_back( cur_boundary ); - chunk_boundaries.push_back( nr_columns ); - - for( dimension cur_dim = max_dim; cur_dim >= 1; cur_dim-- ) { - // Phase 1: Reduce chunks locally -- 1st pass - #pragma omp parallel for schedule( guided, 1 ) - for( index chunk_id = 0; chunk_id < (index)chunk_boundaries.size() - 1; chunk_id++ ) - _local_chunk_reduction( boundary_matrix, lowest_one_lookup, column_type, cur_dim, - chunk_boundaries[ chunk_id ], chunk_boundaries[ chunk_id + 1 ], chunk_boundaries[ chunk_id ] ); - boundary_matrix.sync(); - - // Phase 1: Reduce chunks locally -- 2nd pass - #pragma omp parallel for schedule( guided, 1 ) - for( index chunk_id = 1; chunk_id < (index)chunk_boundaries.size( ) - 1; chunk_id++ ) - _local_chunk_reduction( boundary_matrix, lowest_one_lookup, column_type, cur_dim, - chunk_boundaries[ chunk_id ], chunk_boundaries[ chunk_id + 1 ], chunk_boundaries[ chunk_id - 1 ] ); - boundary_matrix.sync( ); - } - - // get global columns - std::vector< index > global_columns; - for( index cur_col_idx = 0; cur_col_idx < nr_columns; cur_col_idx++ ) - if( column_type[ cur_col_idx ] == GLOBAL ) - global_columns.push_back( cur_col_idx ); - - // get active columns - #pragma omp parallel for - for( index idx = 0; idx < (index)global_columns.size(); idx++ ) - is_active[ global_columns[ idx ] ] = true; - _get_active_columns( boundary_matrix, lowest_one_lookup, column_type, global_columns, is_active ); - - // Phase 2+3: Simplify columns and reduce them - for( dimension cur_dim = max_dim; cur_dim >= 1; cur_dim-- ) { - // Phase 2: Simplify columns - std::vector< index > temp_col; - #pragma omp parallel for schedule( guided, 1 ), private( temp_col ) - for( index idx = 0; idx < (index)global_columns.size(); idx++ ) - if( boundary_matrix.get_dim( global_columns[ idx ] ) == cur_dim ) - _global_column_simplification( global_columns[ idx ], boundary_matrix, lowest_one_lookup, column_type, is_active, temp_col ); - boundary_matrix.sync(); - - // Phase 3: Reduce columns - for( index idx = 0; idx < (index)global_columns.size(); idx++ ) { - index cur_col = global_columns[ idx ]; - if( boundary_matrix.get_dim( cur_col ) == cur_dim && column_type[ cur_col ] == GLOBAL ) { - index lowest_one = boundary_matrix.get_max_index( cur_col ); - while( lowest_one != -1 && lowest_one_lookup[ lowest_one ] != -1 ) { - boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); - lowest_one = boundary_matrix.get_max_index( cur_col ); - } - if( lowest_one != -1 ) { - lowest_one_lookup[ lowest_one ] = cur_col; - boundary_matrix.clear( lowest_one ); - } - boundary_matrix.finalize( cur_col ); - } - } - } - - boundary_matrix.sync(); - } - - protected: - template< typename Representation > - void _local_chunk_reduction( boundary_matrix< Representation >& boundary_matrix - , std::vector& lowest_one_lookup - , std::vector< column_type >& column_type - , const dimension cur_dim - , const index chunk_begin - , const index chunk_end - , const index row_begin ) { - - for( index cur_col = chunk_begin; cur_col < chunk_end; cur_col++ ) { - if( column_type[ cur_col ] == GLOBAL && boundary_matrix.get_dim( cur_col ) == cur_dim ) { - index lowest_one = boundary_matrix.get_max_index( cur_col ); - while( lowest_one != -1 && lowest_one >= row_begin && lowest_one_lookup[ lowest_one ] != -1 ) { - boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); - lowest_one = boundary_matrix.get_max_index( cur_col ); - } - if( lowest_one >= row_begin ) { - lowest_one_lookup[ lowest_one ] = cur_col; - column_type[ cur_col ] = LOCAL_NEGATIVE; - column_type[ lowest_one ] = LOCAL_POSITIVE; - boundary_matrix.clear( lowest_one ); - boundary_matrix.finalize( cur_col ); - } - } - } - } - - template< typename Representation > - void _get_active_columns( const boundary_matrix< Representation >& boundary_matrix - , const std::vector< index >& lowest_one_lookup - , const std::vector< column_type >& column_type - , const std::vector< index >& global_columns - , std::vector< char >& is_active ) { - - const index nr_columns = boundary_matrix.get_num_cols(); - std::vector< char > finished( nr_columns, false ); - - std::vector< std::pair < index, index > > stack; - std::vector< index > cur_col_values; - #pragma omp parallel for schedule( guided, 1 ), private( stack, cur_col_values ) - for( index idx = 0; idx < (index)global_columns.size(); idx++ ) { - bool pop_next = false; - index start_col = global_columns[ idx ]; - stack.push_back( std::pair< index, index >( start_col, -1 ) ); - while( !stack.empty() ) { - index cur_col = stack.back().first; - index prev_col = stack.back().second; - if( pop_next ) { - stack.pop_back(); - pop_next = false; - if( prev_col != -1 ) { - if( is_active[ cur_col ] ) { - is_active[ prev_col ] = true; - } - if( prev_col == stack.back().first ) { - finished[ prev_col ] = true; - pop_next = true; - } - } - } else { - pop_next = true; - boundary_matrix.get_col( cur_col, cur_col_values ); - for( index idx = 0; idx < (index) cur_col_values.size(); idx++ ) { - index cur_row = cur_col_values[ idx ]; - if( ( column_type[ cur_row ] == GLOBAL ) ) { - is_active[ cur_col ] = true; - } else if( column_type[ cur_row ] == LOCAL_POSITIVE ) { - index next_col = lowest_one_lookup[ cur_row ]; - if( next_col != cur_col && !finished[ cur_col ] ) { - stack.push_back( std::make_pair( next_col, cur_col ) ); - pop_next = false; - } - } - } - } - } - } - } - - template< typename Representation > - void _global_column_simplification( const index col_idx - , boundary_matrix< Representation >& boundary_matrix - , const std::vector< index >& lowest_one_lookup - , const std::vector< column_type >& column_type - , const std::vector< char >& is_active - , std::vector< index >& temp_col ) - { - temp_col.clear(); - while( !boundary_matrix.is_empty( col_idx ) ) { - index cur_row = boundary_matrix.get_max_index( col_idx ); - switch( column_type[ cur_row ] ) { - case GLOBAL: - temp_col.push_back( cur_row ); - boundary_matrix.remove_max( col_idx ); - break; - case LOCAL_NEGATIVE: - boundary_matrix.remove_max( col_idx ); - break; - case LOCAL_POSITIVE: - if( is_active[ lowest_one_lookup[ cur_row ] ] ) - boundary_matrix.add_to( lowest_one_lookup[ cur_row ], col_idx ); - else - boundary_matrix.remove_max( col_idx ); - break; - } - } - std::reverse( temp_col.begin(), temp_col.end() ); - boundary_matrix.set_col( col_idx, temp_col ); - } - }; - - class chunk_reduction : public chunk_reduction_impl {}; - class chunk_reduction_sqrt : public chunk_reduction_impl {}; -} diff --git a/zigzag/phat/algorithms/row_reduction.h b/zigzag/phat/algorithms/row_reduction.h deleted file mode 100644 index cdd1a8f..0000000 --- a/zigzag/phat/algorithms/row_reduction.h +++ /dev/null @@ -1,56 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class row_reduction { - public: - template< typename Representation > - void operator() ( boundary_matrix< Representation >& boundary_matrix ) { - - const index nr_columns = boundary_matrix.get_num_cols(); - std::vector< std::vector< index > > lowest_one_lookup( nr_columns ); - - for( index cur_col = nr_columns - 1; cur_col >= 0; cur_col-- ) { - if( !boundary_matrix.is_empty( cur_col ) ) - lowest_one_lookup[ boundary_matrix.get_max_index( cur_col ) ].push_back( cur_col ); - - if( !lowest_one_lookup[ cur_col ].empty() ) { - boundary_matrix.clear( cur_col ); - boundary_matrix.finalize( cur_col ); - std::vector< index >& cols_with_cur_lowest = lowest_one_lookup[ cur_col ]; - index source = *min_element( cols_with_cur_lowest.begin(), cols_with_cur_lowest.end() ); - for( index idx = 0; idx < (index)cols_with_cur_lowest.size(); idx++ ) { - index target = cols_with_cur_lowest[ idx ]; - if( target != source && !boundary_matrix.is_empty( target ) ) { - boundary_matrix.add_to( source, target ); - if( !boundary_matrix.is_empty( target ) ) { - index lowest_one_of_target = boundary_matrix.get_max_index( target ); - lowest_one_lookup[ lowest_one_of_target ].push_back( target ); - } - } - } - } - } - } - }; -} \ No newline at end of file diff --git a/zigzag/phat/algorithms/spectral_sequence_reduction.h b/zigzag/phat/algorithms/spectral_sequence_reduction.h deleted file mode 100644 index bf442e6..0000000 --- a/zigzag/phat/algorithms/spectral_sequence_reduction.h +++ /dev/null @@ -1,80 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class spectral_sequence_reduction { - public: - template< typename Representation > - void operator () ( boundary_matrix< Representation >& boundary_matrix ) { - - const index nr_columns = boundary_matrix.get_num_cols(); - std::vector< index > lowest_one_lookup( nr_columns, -1 ); - - //const index num_stripes = (index) sqrt( (double)nr_columns ); - const index num_stripes = omp_get_max_threads(); - - index block_size = ( nr_columns % num_stripes == 0 ) ? nr_columns / num_stripes : block_size = nr_columns / num_stripes + 1; - - std::vector< std::vector< index > > unreduced_cols_cur_pass( num_stripes ); - std::vector< std::vector< index > > unreduced_cols_next_pass( num_stripes ); - - for( index cur_dim = boundary_matrix.get_max_dim(); cur_dim >= 1 ; cur_dim-- ) { - #pragma omp parallel for schedule( guided, 1 ) - for( index cur_stripe = 0; cur_stripe < num_stripes; cur_stripe++ ) { - index col_begin = cur_stripe * block_size; - index col_end = std::min( (cur_stripe+1) * block_size, nr_columns ); - for( index cur_col = col_begin; cur_col < col_end; cur_col++ ) - if( boundary_matrix.get_dim( cur_col ) == cur_dim && boundary_matrix.get_max_index( cur_col ) != -1 ) - unreduced_cols_cur_pass[ cur_stripe ].push_back( cur_col ); - } - for( index cur_pass = 0; cur_pass < num_stripes; cur_pass++ ) { - boundary_matrix.sync(); - #pragma omp parallel for schedule( guided, 1 ) - for( int cur_stripe = 0; cur_stripe < num_stripes; cur_stripe++ ) { - index row_begin = (cur_stripe - cur_pass) * block_size; - index row_end = row_begin + block_size; - unreduced_cols_next_pass[ cur_stripe ].clear(); - for( index idx = 0; idx < (index)unreduced_cols_cur_pass[ cur_stripe ].size(); idx++ ) { - index cur_col = unreduced_cols_cur_pass[ cur_stripe ][ idx ]; - index lowest_one = boundary_matrix.get_max_index( cur_col ); - while( lowest_one != -1 && lowest_one >= row_begin && lowest_one < row_end && lowest_one_lookup[ lowest_one ] != -1 ) { - boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); - lowest_one = boundary_matrix.get_max_index( cur_col ); - } - if( lowest_one != -1 ) { - if( lowest_one >= row_begin && lowest_one < row_end ) { - lowest_one_lookup[ lowest_one ] = cur_col; - boundary_matrix.clear( lowest_one ); - boundary_matrix.finalize( cur_col ); - } else { - unreduced_cols_next_pass[ cur_stripe ].push_back( cur_col ); - } - } - } - unreduced_cols_next_pass[ cur_stripe ].swap( unreduced_cols_cur_pass[ cur_stripe ] ); - } - } - } - } - }; -} diff --git a/zigzag/phat/algorithms/standard_reduction.h b/zigzag/phat/algorithms/standard_reduction.h deleted file mode 100644 index d8762fa..0000000 --- a/zigzag/phat/algorithms/standard_reduction.h +++ /dev/null @@ -1,46 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class standard_reduction { - public: - template< typename Representation > - void operator() ( boundary_matrix< Representation >& boundary_matrix ) { - - const index nr_columns = boundary_matrix.get_num_cols(); - std::vector< index > lowest_one_lookup( nr_columns, -1 ); - - for( index cur_col = 0; cur_col < nr_columns; cur_col++ ) { - index lowest_one = boundary_matrix.get_max_index( cur_col ); - while( lowest_one != -1 && lowest_one_lookup[ lowest_one ] != -1 ) { - boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); - lowest_one = boundary_matrix.get_max_index( cur_col ); - } - if( lowest_one != -1 ) { - lowest_one_lookup[ lowest_one ] = cur_col; - } - boundary_matrix.finalize( cur_col ); - } - } - }; -} diff --git a/zigzag/phat/representations/full_pivot_column.h b/zigzag/phat/representations/full_pivot_column.h deleted file mode 100644 index c2e9e3c..0000000 --- a/zigzag/phat/representations/full_pivot_column.h +++ /dev/null @@ -1,100 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class full_column { - - protected: - std::priority_queue< index > history; - std::vector< char > is_in_history; - std::vector< char > col_bit_field; - - public: - void init( const index total_size ) { - col_bit_field.resize( total_size, false ); - is_in_history.resize( total_size, false ); - } - - void add_col( const column& col ) { - for( index idx = 0; idx < (index) col.size(); idx++ ) { - add_index( col[ idx ] ); - } - } - - void add_index( const index idx ) { - if( !is_in_history[ idx ] ) { - history.push( idx ); - is_in_history[ idx ] = true; - } - - col_bit_field[ idx ] = !col_bit_field[ idx ]; - } - - index get_max_index() { - while( history.size() > 0 ) { - index topIndex = history.top(); - if( col_bit_field[ topIndex ] ) { - return topIndex; - } else { - history.pop(); - is_in_history[ topIndex ] = false; - } - } - - return -1; - } - - void get_col_and_clear( column& col ) { - while( !is_empty() ) { - col.push_back( get_max_index() ); - add_index( get_max_index() ); - } - std::reverse( col.begin(), col.end() ); - } - - bool is_empty() { - return (get_max_index() == -1); - } - - void clear() { - while( !is_empty() ) - add_index( get_max_index() ); - } - - void remove_max() { - add_index( get_max_index() ); - } - - void set_col( const column& col ) { - clear(); - add_col( col ); - } - - void get_col( column& col ) { - get_col_and_clear( col ); - add_col( col ); - } - }; - - typedef abstract_pivot_column< full_column > full_pivot_column; -} diff --git a/zigzag/phat/representations/heap_pivot_column.h b/zigzag/phat/representations/heap_pivot_column.h deleted file mode 100644 index 33cd07b..0000000 --- a/zigzag/phat/representations/heap_pivot_column.h +++ /dev/null @@ -1,126 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class heap_column { - - protected: - std::priority_queue< index > data; - - column temp_col; - index inserts_since_last_prune; - - void prune() - { - temp_col.clear( ); - index max_index = pop_max_index( ); - while( max_index != -1 ) { - temp_col.push_back( max_index ); - max_index = pop_max_index( ); - } - - for( index idx = 0; idx < (index)temp_col.size( ); idx++ ) - data.push( temp_col[ idx ] ); - - inserts_since_last_prune = 0; - } - - index pop_max_index() - { - if( data.empty( ) ) - return -1; - else { - index max_element = data.top( ); - data.pop(); - while( !data.empty( ) && data.top( ) == max_element ) { - data.pop( ); - if( data.empty( ) ) - return -1; - else { - max_element = data.top( ); - data.pop( ); - } - } - return max_element; - } - } - - public: - void init( const index total_size ) { - inserts_since_last_prune = 0; - clear(); - } - - void add_col( const column& col ) { - for( index idx = 0; idx < (index) col.size(); idx++ ) - data.push( col[ idx ] ); - inserts_since_last_prune += col.size( ); - if( 2 * inserts_since_last_prune >( index ) data.size( ) ) - prune(); - } - - index get_max_index() { - index max_element = pop_max_index( ); - if( max_element == -1 ) - return -1; - else { - data.push( max_element ); - return max_element; - } - } - - void get_col_and_clear( column& col ) { - col.clear(); - index max_index = pop_max_index( ); - while( max_index != -1 ) { - col.push_back( max_index ); - max_index = pop_max_index( ); - } - std::reverse( col.begin(), col.end() ); - } - - bool is_empty() { - return get_max_index() == -1; - } - - void clear() { - data = std::priority_queue< index >(); - } - - void remove_max() { - pop_max_index(); - } - - void set_col( const column& col ) { - clear(); - add_col( col ); - } - - void get_col( column& col ) { - get_col_and_clear( col ); - add_col( col ); - } - }; - - typedef abstract_pivot_column< heap_column > heap_pivot_column; -} diff --git a/zigzag/phat/representations/sparse_pivot_column.h b/zigzag/phat/representations/sparse_pivot_column.h deleted file mode 100644 index 390fd91..0000000 --- a/zigzag/phat/representations/sparse_pivot_column.h +++ /dev/null @@ -1,79 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include -#include - -namespace phat { - class sparse_column { - - protected: - std::set< index > data; - - void add_index( const index idx ) { - std::pair< std::set< index >::iterator, bool > result = data.insert( idx ); - if( result.second == false ) - data.erase( result.first ); - } - - public: - void init( const index total_size ) { - data.clear(); - } - - void add_col( const column& col ) { - for( index idx = 0; idx < (index) col.size(); idx++ ) - add_index( col[ idx ] ); - } - - index get_max_index() { - return data.empty() ? -1 : *data.rbegin(); - } - - void get_col_and_clear( column& col ) { - col.assign( data.begin(), data.end() ); - data.clear(); - } - - bool is_empty() { - return data.empty(); - } - - void clear() { - data.clear(); - } - - void remove_max() { - add_index( get_max_index() ); - } - - void set_col( const column& col ) { - clear(); - add_col( col ); - } - - void get_col( column& col ) { - get_col_and_clear( col ); - add_col( col ); - } - }; - - typedef abstract_pivot_column< sparse_column > sparse_pivot_column; -} diff --git a/zigzag/phat/representations/vector_heap.h b/zigzag/phat/representations/vector_heap.h deleted file mode 100644 index db0420f..0000000 --- a/zigzag/phat/representations/vector_heap.h +++ /dev/null @@ -1,170 +0,0 @@ -/* Copyright 2013 IST Austria -Contributed by: Jan Reininghaus - -This file is part of PHAT. - -PHAT is free software: you can redistribute it and/or modify -it under the terms of the GNU Lesser General Public License as published by -the Free Software Foundation, either version 3 of the License, or -(at your option) any later version. - -PHAT is distributed in the hope that it will be useful, -but WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU Lesser General Public License for more details. - -You should have received a copy of the GNU Lesser General Public License -along with PHAT. If not, see . */ - -#pragma once - -#include - -namespace phat { - class vector_heap { - - protected: - std::vector< dimension > dims; - std::vector< column > matrix; - - std::vector< index > inserts_since_last_prune; - - mutable thread_local_storage< column > temp_column_buffer; - - protected: - void _prune( index idx ) - { - column& col = matrix[ idx ]; - column& temp_col = temp_column_buffer(); - temp_col.clear(); - index max_index = _pop_max_index( col ); - while( max_index != -1 ) { - temp_col.push_back( max_index ); - max_index = _pop_max_index( col ); - } - col = temp_col; - std::reverse( col.begin( ), col.end( ) ); - std::make_heap( col.begin( ), col.end( ) ); - inserts_since_last_prune[ idx ] = 0; - } - - index _pop_max_index( index idx ) - { - return _pop_max_index( matrix[ idx ] ); - } - - index _pop_max_index( column& col ) const - { - if( col.empty( ) ) - return -1; - else { - index max_element = col.front( ); - std::pop_heap( col.begin( ), col.end( ) ); - col.pop_back( ); - while( !col.empty( ) && col.front( ) == max_element ) { - std::pop_heap( col.begin( ), col.end( ) ); - col.pop_back( ); - if( col.empty( ) ) - return -1; - else { - max_element = col.front( ); - std::pop_heap( col.begin( ), col.end( ) ); - col.pop_back( ); - } - } - return max_element; - } - } - - public: - // overall number of cells in boundary_matrix - index _get_num_cols( ) const - { - return (index)matrix.size( ); - } - void _set_num_cols( index nr_of_columns ) - { - dims.resize( nr_of_columns ); - matrix.resize( nr_of_columns ); - inserts_since_last_prune.assign( nr_of_columns, 0 ); - } - - // dimension of given index - dimension _get_dim( index idx ) const - { - return dims[ idx ]; - } - void _set_dim( index idx, dimension dim ) - { - dims[ idx ] = dim; - } - - // replaces(!) content of 'col' with boundary of given index - void _get_col( index idx, column& col ) const - { - temp_column_buffer( ) = matrix[ idx ]; - - index max_index = _pop_max_index( temp_column_buffer() ); - while( max_index != -1 ) { - col.push_back( max_index ); - max_index = _pop_max_index( temp_column_buffer( ) ); - } - std::reverse( col.begin( ), col.end( ) ); - } - void _set_col( index idx, const column& col ) - { - matrix[ idx ] = col; - std::make_heap( matrix[ idx ].begin( ), matrix[ idx ].end( ) ); - } - - // true iff boundary of given idx is empty - bool _is_empty( index idx ) const - { - return _get_max_index( idx ) == -1; - } - - // largest row index of given column idx (new name for lowestOne()) - index _get_max_index( index idx ) const - { - column& col = const_cast< column& >( matrix[ idx ] ); - index max_element = _pop_max_index( col ); - col.push_back( max_element ); - std::push_heap( col.begin( ), col.end( ) ); - return max_element; - } - - // removes the maximal index of a column - void _remove_max( index idx ) - { - _pop_max_index( idx ); - } - - // clears given column - void _clear( index idx ) - { - matrix[ idx ].clear( ); - } - - // syncronizes all data structures (essential for openmp stuff) - void _sync( ) {} - - // adds column 'source' to column 'target' - void _add_to( index source, index target ) - { - for( index idx = 0; idx < (index)matrix[ source ].size( ); idx++ ) { - matrix[ target ].push_back( matrix[ source ][ idx ] ); - std::push_heap( matrix[ target ].begin(), matrix[ target ].end() ); - } - inserts_since_last_prune[ target ] += matrix[ source ].size(); - - if( 2 * inserts_since_last_prune[ target ] > ( index )matrix[ target ].size() ) - _prune( target ); - } - - // finalizes given column - void _finalize( index idx ) { - _prune( idx ); - } - - }; -} diff --git a/zigzag/phat/representations/vector_list.h b/zigzag/phat/representations/vector_list.h deleted file mode 100644 index ca0b5b8..0000000 --- a/zigzag/phat/representations/vector_list.h +++ /dev/null @@ -1,101 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include - -namespace phat { - class vector_list { - - protected: - std::vector< dimension > dims; - std::vector< std::list< index > > matrix; - - public: - // overall number of cells in boundary_matrix - index _get_num_cols() const { - return (index)matrix.size(); - } - void _set_num_cols( index nr_of_columns ) { - dims.resize( nr_of_columns ); - matrix.resize( nr_of_columns ); - } - - // dimension of given index - dimension _get_dim( index idx ) const { - return dims[ idx ]; - } - void _set_dim( index idx, dimension dim ) { - dims[ idx ] = dim; - } - - // replaces(!) content of 'col' with boundary of given index - void _get_col( index idx, column& col ) const { - col.clear(); - col.reserve( matrix[idx].size() ); - std::copy (matrix[idx].begin(), matrix[idx].end(), std::back_inserter(col) ); - } - - void _set_col( index idx, const column& col ) { - matrix[ idx ].clear(); - matrix[ idx ].resize( col.size() ); - std::copy (col.begin(), col.end(), matrix[ idx ].begin() ); - } - - // true iff boundary of given idx is empty - bool _is_empty( index idx ) const { - return matrix[ idx ].empty(); - } - - // largest row index of given column idx (new name for lowestOne()) - index _get_max_index( index idx ) const { - return matrix[ idx ].empty() ? -1 : *matrix[ idx ].rbegin(); - } - - // removes the maximal index of a column - void _remove_max( index idx ) { - std::list< index >::iterator it = matrix[ idx ].end(); - it--; - matrix[ idx ].erase( it ); - } - - // clears given column - void _clear( index idx ) { - matrix[ idx ].clear(); - } - - // syncronizes all data structures (essential for openmp stuff) - void _sync() {} - - // adds column 'source' to column 'target' - void _add_to( index source, index target ) { - std::list< index >& source_col = matrix[ source ]; - std::list< index >& target_col = matrix[ target ]; - std::list< index > temp_col; - target_col.swap( temp_col ); - std::set_symmetric_difference( temp_col.begin(), temp_col.end(), - source_col.begin(), source_col.end(), - std::back_inserter( target_col ) ); - } - - // finalizes given column - void _finalize( index idx ) { - } - }; -} diff --git a/zigzag/phat/representations/vector_set.h b/zigzag/phat/representations/vector_set.h deleted file mode 100644 index 6878a27..0000000 --- a/zigzag/phat/representations/vector_set.h +++ /dev/null @@ -1,99 +0,0 @@ -/* Copyright 2013 IST Austria - Contributed by: Ulrich Bauer, Michael Kerber, Jan Reininghaus - - This file is part of PHAT. - - PHAT is free software: you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - PHAT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with PHAT. If not, see . */ - -#pragma once - -#include - -namespace phat { - class vector_set { - - protected: - std::vector< dimension > dims; - std::vector< std::set< index > > matrix; - - public: - // overall number of cells in boundary_matrix - index _get_num_cols() const { - return (index)matrix.size(); - } - void _set_num_cols( index nr_of_columns ) { - dims.resize( nr_of_columns ); - matrix.resize( nr_of_columns ); - } - - // dimension of given index - dimension _get_dim( index idx ) const { - return dims[ idx ]; - } - void _set_dim( index idx, dimension dim ) { - dims[ idx ] = dim; - } - - // replaces(!) content of 'col' with boundary of given index - void _get_col( index idx, column& col ) const { - col.clear(); - col.reserve( matrix[idx].size() ); - std::copy (matrix[idx].begin(), matrix[idx].end(), std::back_inserter(col) ); - } - void _set_col( index idx, const column& col ) { - matrix[ idx ].clear(); - matrix[ idx ].insert( col.begin(), col.end() ); - } - - // true iff boundary of given idx is empty - bool _is_empty( index idx ) const { - return matrix[ idx ].empty(); - } - - // largest row index of given column idx (new name for lowestOne()) - index _get_max_index( index idx ) const { - return matrix[ idx ].empty() ? -1 : *matrix[ idx ].rbegin(); - } - - // removes the maximal index of a column - void _remove_max( index idx ) { - std::set< index >::iterator it = matrix[ idx ].end(); - it--; - matrix[ idx ].erase( it ); - } - - // clears given column - void _clear( index idx ) { - matrix[ idx ].clear(); - } - - // syncronizes all data structures (essential for openmp stuff) - void _sync() {} - - // adds column 'source' to column 'target' - void _add_to( index source, index target ) { - for( std::set< index >::iterator it = matrix[ source ].begin(); it != matrix[ source ].end(); it++ ) { - std::set< index >& col = matrix[ target ]; - std::pair< std::set< index >::iterator, bool > result = col.insert( *it ); - if( !result.second ) - col.erase( result.first ); - } - } - - // finalizes given column - void _finalize( index idx ) { - } - - }; -} diff --git a/zigzag/utils.cpp b/zigzag/utils.cpp deleted file mode 100644 index 35915d0..0000000 --- a/zigzag/utils.cpp +++ /dev/null @@ -1,102 +0,0 @@ -#include "utils.h" -#include -#include -#include - -size_t IntVecHash::operator()(const vector &v) const { - std::size_t seed = 0; - for (auto i = 0; i < v.size(); i ++) { - boost::hash_combine(seed, v[i]); - } - - return seed; -} - -bool operator==(const vector &v1, const vector &v2) { - if (v1.size() != v2.size()) { - return false; - } - - for (auto i = 0; i < v1.size(); i ++) { - if (v1[i] != v2[i]) { - return false; - } - } - - return true; -} - -ostream& operator<<(ostream& os, const array& a) { - if (a.size() == 0) { - os << "()"; - return os; - } - - os << '(' << a[0]; - for (auto i = 1; i < a.size(); i ++) { - os << ',' << a[i]; - } - os << ')'; - - return os; -} - -ostream& operator<<(ostream& os, const array& a) { - if (a.size() == 0) { - os << "()"; - return os; - } - - os << '(' << a[0]; - for (auto i = 1; i < a.size(); i ++) { - os << ',' << a[i]; - } - os << ')'; - - return os; -} - -void splitStrLast( - const std::string& str, const std::string splitter, - std::string* before, std::string* after) { - - auto pos = str.find_last_of(splitter); - - if (before != NULL) { - *before = str.substr(0, pos); - } - - if (after != NULL) { - if (pos != std::string::npos) { - *after = str.substr(pos + 1); - } else { - *after = ""; - } - } -} - -void getFilePurename(const std::string& filename, std::string *purename) { - auto pos = filename.find_last_of("/\\"); - if (pos != std::string::npos) { - *purename = filename.substr(pos + 1); - } else { - *purename = filename; - } - - pos = purename->find_last_of("."); - *purename = purename->substr(0, pos); -} - -void printPers(const std::string outfname, - const std::vector< std::tuple > &pers) { - - std::ofstream fout(outfname); - if (!fout) { - ERR() << "open '" << outfname << "' for write failed!" << std::endl; - exit(-1); - } - - for (auto& e : pers) { - fout << std::get<0>(e) << " " << std::get<1>(e) << " " << std::get<2>(e) << std::endl; - } -}