Skip to content

Instantly share code, notes, and snippets.

@xgao32

xgao32/computing-bottleneck-and-wasserstein-distance-exception.ipynb Secret

Created September 28, 2021 02:09
Show Gist options
  • Save xgao32/2a4ebff7061e8ea723bb2fa85fde4765 to your computer and use it in GitHub Desktop.
Save xgao32/2a4ebff7061e8ea723bb2fa85fde4765 to your computer and use it in GitHub Desktop.
Download ZIP
Computing Bottleneck and Wasserstein Distance exception.ipynb
Display the source blob
Display the rendered blob
Raw
computing-bottleneck-and-wasserstein-distance-exception.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Computing Bottleneck and Wasserstein Distance exception.ipynb",
"provenance": [],
"collapsed_sections": [],
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/xgao32/2a4ebff7061e8ea723bb2fa85fde4765/computing-bottleneck-and-wasserstein-distance-exception.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "Qp4wdUgjtZ_S",
"outputId": "fcdbb263-450c-4cda-d86a-ede54d81d94a"
},
"source": [
"import pandas as pd\n",
"import seaborn as sns\n",
"import networkx as nx\n",
"import numpy as np\n",
"from scipy.io import loadmat\n",
"from scipy import stats\n",
"from scipy import sparse\n",
"from numpy import array\n",
"from matplotlib import pyplot as plt\n",
"import time\n",
"\n",
"%pip install scikit-tda\n",
"from ripser import ripser\n",
"from persim import plot_diagrams\n",
"import persim\n",
"\n",
"np.random.seed(1)\n",
"\n",
"%config Completer.use_jedi = False\n",
"gdrive = False"
],
"execution_count": 9,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Requirement already satisfied: wrapt<2,>=1.10 in /usr/local/lib/python3.7/dist-packages (from deprecated->persim->scikit-tda) (1.12.1)\n"
]
},
{
"output_type": "stream",
"name": "stderr",
"text": [
"/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:19: UserWarning: Config option `use_jedi` not recognized by `IPCompleter`.\n"
]
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "a9ng6MgVMNcy"
},
"source": [
"#!/usr/bin/env python\n",
"\n",
"__author__ = 'David S. Campo, Ph.D.'\n",
"'''\n",
"NK_model_simulation.py is a program for creating fitness landscapes with the NK model. \n",
"\n",
"parameters\n",
"'-o', output folder \n",
"'-n', sequence length, N \n",
"'-k', Number of positions linked, K'\n",
"'-m', method desired: random, NK, NKP. In random, every sequence gets a random number. \n",
"'-p', parameter for NKP, fraction that are neutral. Those elements get fitness=1, all others get a random number \n",
"'-s', Number of symbols' \n",
"'-r', number of random trials, default 10\n",
"'-f', Make space network only, default is no. If you are not making the space, the file must be there on the same output folder \n",
"'-i', Measure mutual information matrix for each file, default is no. \n",
" If yes, it will create a csv file with the average mutual information of each pair of positions (over all randonm trials). \n",
" It assumes a proportional relationship between fitness and count number of each sequence. \n",
" The higher the fitness, the higher the counts. \n",
" The fitness is multiplied by 1000 and then the integer is taken to establish the counts.\n",
" \n",
"\n",
"'''\n",
"import os\n",
"import numpy as np\n",
"import optparse\n",
"import networkx as nx\n",
"from itertools import product\n",
"import random\n",
"from sklearn.metrics.cluster import adjusted_mutual_info_score, normalized_mutual_info_score, mutual_info_score\n",
"import pickle\n",
"\n",
"\n",
"def posVectorCalc(stateList, haploSize, reads, counts,stateNum):\n",
" ## states by position\n",
" freqCount = np.zeros((haploSize, stateNum))\n",
" readCounter = 0\n",
" for read in reads:\n",
" for pos in range(haploSize):\n",
" for state in range(stateNum):\n",
" if read[pos] == stateList[state]:\n",
" freqCount[pos, state] = freqCount[pos, state] + counts[readCounter]\n",
" readCounter = readCounter + 1\n",
" return freqCount\n",
"\n",
"def neiHamming(vector, stateList, notPosList, sizeVector):\n",
" allNei = []\n",
" for pos in range(sizeVector):\n",
" notPos = notPosList[vector[pos]]\n",
" for state in notPos:\n",
" temp = vector.copy()\n",
" temp[pos] = state\n",
" allNei.append(tuple(temp))\n",
" \n",
" return allNei\n",
"\n",
"def getFitParams(methodChoice, sizeVector,stateList, parameterMethod0, parameterMethod1): \n",
" if methodChoice == 'random':\n",
" fitParam = [methodChoice,sizeVector]\n",
"\n",
" elif methodChoice == 'NK':\n",
" kNum = parameterMethod0 \n",
" comboDict = {}\n",
" tempCombo = [stateList] * (kNum+1)\n",
" # print(\"tempCombo \" , tempCombo)\n",
" for combo in list(product(*tempCombo)):\n",
" comboDict[combo] = random.random() # assign random fitness to all possibele sequence\n",
" fitParam = [methodChoice, comboDict, kNum,sizeVector]\n",
" \n",
" # print(\"comboDict \" ,comboDict)\n",
"\n",
" elif methodChoice == 'NKP':\n",
" fracNeutral = parameterMethod1\n",
" kNum = parameterMethod0\n",
" comboDict = {}\n",
" tempCombo = [stateList] * (kNum+1)\n",
" comboList = []\n",
" for combo in list(product(*tempCombo)):\n",
" comboList.append(combo)\n",
" random.shuffle(comboList)\n",
" for combo in comboList:\n",
" coin1 = random.random()\n",
" if coin1 <= fracNeutral:\n",
" comboDict[combo] = 1\n",
" else:\n",
" comboDict[combo] = coin1\n",
" fitParam = [methodChoice, comboDict, kNum,sizeVector] \n",
"\n",
" \n",
" return fitParam # list of parameters\n",
"\n",
"# comment out to avoid duplicate fitness function\n",
"'''\n",
"def fitness(fitParam, vector):\n",
" # vector = vertex that is being assigned fitness\n",
" methodChoice = fitParam[0]\n",
" sizeVector = fitParam[-1]\n",
" minFit = 0#random.random()/100\n",
" if methodChoice == 'random':\n",
" fitVector = random.random()\n",
" fitnessValue = np.divide(minFit + fitVector, 1, dtype = float)\n",
" \n",
" elif methodChoice == 'NK' or methodChoice == 'NKP' : \n",
" comboDict = fitParam[1]\n",
" #print(\"comboDict \", comboDict)\n",
" kNum = fitParam[2]\n",
" vectorNew = vector.copy() # vector = vertex as a list of symbols\n",
" for pos in range(kNum):\n",
" vectorNew.append(vectorNew[pos])\n",
" # print(\"vectorNew \", vectorNew)\n",
" \n",
" fitVector = 0\n",
" for pos in range(sizeVector):\n",
" temp = []\n",
" for i in range(kNum+1):\n",
" temp.append(vectorNew[pos+i]) # \n",
" # print(\"temp \", temp)\n",
" fitVector = fitVector + comboDict[tuple(temp)]\n",
" fitnessValue = np.divide(fitVector+minFit, sizeVector, dtype = float)\n",
"\n",
" return fitnessValue\n",
"'''\n",
"\n",
"def infMatCalc1(seqs, useAdjusted, counts, stateList, freqVector):\n",
" haploNumPair = len(seqs)\n",
" haploSize = len(seqs[0])\n",
" posList = []\n",
" for p in range(haploSize):\n",
" tempList = []\n",
" for s in range(haploNumPair):\n",
" c = int(counts[s])\n",
" for i in range(c):\n",
" tempList.append(seqs[s][p])\n",
" posList.append(tempList) \n",
" \n",
" ## MI hamming\n",
" Matrix = np.zeros((haploSize,haploSize))\n",
" for pos1 in range(haploSize-1):\n",
" for pos2 in range(pos1+1, haploSize):\n",
" if useAdjusted == 0:\n",
" MI0 = mutual_info_score(posList[pos1],posList[pos2])\n",
" elif useAdjusted == 1:\n",
" MI0 = adjusted_mutual_info_score(posList[pos1],posList[pos2])\n",
" elif useAdjusted == 2:\n",
" MI0 = normalized_mutual_info_score(posList[pos1],posList[pos2])\n",
" MI = max(0, MI0)\n",
" Matrix[pos1][pos2] = MI\n",
" Matrix[pos2][pos1] = MI\n",
"\n",
" return Matrix\n",
"\n",
"#########################################################################\n",
"#argument_parser = optparse.OptionParser() \n",
"#argument_parser.add_option('-o', metavar='ODIR', action='store', type=str, dest='pathOutput', default='C:/Users/Perceval/Desktop/work/NK', help='Output directory') \n",
"#argument_parser.add_option('-n', action='store', dest='sizeVector', type= int, default= 10, help='sequence length, N') \n",
"#argument_parser.add_option('-k', action='store', dest='parameterMethod0', type= int, default= 1, help='Number of positions linked, K') \n",
"#argument_parser.add_option('-p', action='store', dest='parameterMethod1', type= int, default= 0.5, help='parameter for NKP') \n",
"#argument_parser.add_option('-s', action='store', dest='stateNum', type= int, default= 2, help='Number of symbols') \n",
"#argument_parser.add_option('-m', action='store', dest='methodChoice', type= str, default= 'NKP', help='methods: random, NK, NKP') \n",
"#argument_parser.add_option('-r', action='store', dest='randNum', type= int, default= 10, help='number of random trials, default 1') \n",
"#argument_parser.add_option('-f', action='store', dest='makeSpace', type= int, default= 0, help='Make space network only, default is no') \n",
"#argument_parser.add_option('-i', action='store', dest='mutualInf', type= int, default= 0, help='Measure mutual information matrix, default is no') \n",
"\n",
"#options, args = argument_parser.parse_args()\n",
"#pathOutput = options.pathOutput \n",
"\n",
"def make_skeleton_graph(N,K=0,symbols=2):\n",
" '''\n",
" Return networkx graph of genotypes connected to neighbors within 1 Hamming distance\n",
" '''\n",
"\n",
" sizeVector = N #options.sizeVector\n",
" parameterMethod0 = K #options.parameterMethod0\n",
" parameterMethod1 = 0 # options.parameterMethod1\n",
" randNum = 1 # options.randNum\n",
" methodChoice = 'NK' # options.methodChoice\n",
" makeSpace = 1 # options.makeSpace\n",
" stateNum = symbols # options.stateNum\n",
" mutualInf = 1 # options.mutualInf\n",
"\n",
" #if not os.path.exists(pathOutput):\n",
" # os.makedirs(pathOutput)\n",
"\n",
" ## initialize\n",
" figNum = 0\n",
" useAdjusted = 2 #0 raw mutual information, 1 adjusted, 2 normalized\n",
" spaceEdges = sizeVector * np.power(stateNum, sizeVector-1) # total number of edges\n",
" spaceNodes = stateNum**sizeVector # total number of vertices\n",
" stateList = [] # list of symbols represented as integers from 0 to number of symbols - 1\n",
" for i in range(stateNum):\n",
" stateList.append(i)\n",
" notPosList = []\n",
" for state in stateList:\n",
" temp = stateList.copy()\n",
" temp.remove(state)\n",
" notPosList.append(temp)\n",
"\n",
"\n",
" ## Define hamming space\n",
" vectorList = [] # list of symbols possible for each position in sequence\n",
" for pos in range(sizeVector):\n",
" vectorList.append(stateList)\n",
" # prevSpaceNet = pathOutput + '/allNet_N' + str(sizeVector)+ '_S' + str(stateNum) + '.pkl' \n",
" if makeSpace == 1:\n",
" # print('Making full hamming graph')\n",
" allNet = nx.Graph()\n",
" for nodeChosen in product(*vectorList): # add all sequence to graph\n",
" allNet.add_node(nodeChosen)\n",
" neiChosen = neiHamming(list(nodeChosen), stateList, notPosList, sizeVector) \n",
" for nei in neiChosen:\n",
" if nei in allNet:\n",
" allNet.add_edge(nodeChosen,nei)\n",
" #f = open(prevSpaceNet,\"wb\")\n",
" #pickle.dump(allNet,f)\n",
" #f.close()\n",
" \n",
" return allNet # return skeleton graph\n",
" else:\n",
" # allNet = pickle.load(open(prevSpaceNet, \"rb\")) \n",
" fitParam = getFitParams(methodChoice, sizeVector,stateList, parameterMethod0, parameterMethod1)\n",
" print('sizeVector=', sizeVector, ', stateNum=', stateNum,',methodChoice=',methodChoice, ',parameter0=', parameterMethod0, ',randNum=', randNum) \n",
" allMI = []\n",
" for sample in range(randNum):\n",
" #fitnessFileName = pathOutput + '/fitness_' + methodChoice + '_N_' + str(sizeVector)+ '_S_' + str(stateNum) + '_randNum_' + str(sample) + '.pkl'\n",
" fitnessDict = {} # dictionary of vertex and fitness of vertex\n",
" print('Sample',sample)\n",
" ## assign fitness\n",
" for nodeChosen in allNet.nodes():\n",
" fitnessDict[nodeChosen] = fitness(fitParam, list(nodeChosen)) # assigning fitness to vertex \n",
" #f = open(fitnessFileName,\"wb\")\n",
" #pickle.dump(fitnessDict,f)\n",
" #f.close() \n",
"\n",
" if mutualInf == 1:\n",
" reads = list(fitnessDict.keys())\n",
" counts = np.multiply(1000,list(fitnessDict.values())).astype(int)\n",
" freqVector = posVectorCalc(stateList, sizeVector, reads, counts, stateNum)\n",
" infoMatrix = infMatCalc1(reads, useAdjusted, counts, stateList, freqVector)\n",
" allMI.append(infoMatrix)\n",
" if mutualInf == 1: \n",
" meanMI = np.average(allMI, axis=0) # mean mutual information \n",
" #outputName0 = pathOutput + '/mutual_information_' + methodChoice + '_N_' + str(sizeVector)+ '_S_' + str(stateNum) + '_randNum_' + str(randNum) + '.csv'\n",
" #outputHandle0 = open(os.path.join(pathOutput, outputName0), \"w\")\n",
" #for i1 in range(sizeVector -1):\n",
" # for i2 in range(i1+1, sizeVector):\n",
" # outputHandle0.write(str(i1) + ',' + str(i2) + ',' + str(meanMI[i1][i2]) + '\\n') \n",
" #outputHandle0.close()\n",
" \n",
" return fitnessDict, meanMI\n",
" \n",
" \n"
],
"execution_count": 10,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "4yqmSdFqguE9"
},
"source": [
"def neighbors(genotype, genotypes):\n",
" return [g for g in genotypes if abs(g ^ genotype).sum() == 1]\n",
"\n",
"def int2bits(k, N):\n",
" x = list(map(int, bin(k)[2:]))\n",
" pad = N - len(x)\n",
" x = [0]*pad + x\n",
" return x\n",
"\n",
"def all_genotypes(N):\n",
" return np.array([int2bits(k, N) for k in range(2**N)], dtype=bool) # binary genotypes\n",
"\n",
"def fitness_i(genotype, i, contribs, mem):\n",
" key = tuple(zip(contribs[i], genotype[contribs[i]])) # assign fitness contribution for each group of interacting loci\n",
" if key not in mem:\n",
" mem[key] = np.random.uniform(0, 1)\n",
" return mem[key]\n",
"\n",
"def fitness(genotype, contribs, mem):\n",
" return np.mean([fitness_i(genotype, i, contribs, mem) for i in range(len(genotype))]) # mean of all fitness contributions\n",
"\n",
"def plot_NK_model(N, K, ax=None):\n",
" genotypes = all_genotypes(N)\n",
" contribs = {\n",
" i: sorted(np.random.choice(\n",
" [n for n in range(N) if n != i], \n",
" K, \n",
" replace=False\n",
" ).tolist() + [i])\n",
" for i in range(N)\n",
" }\n",
" fitness_mem = {}\n",
" ws = [fitness(g, contribs, fitness_mem) for g in genotypes]\n",
" max_w = max(ws)\n",
" min_w = min(ws)\n",
"\n",
" if ax is None:\n",
" fig, ax = plt.subplots(figsize=(12, 8))\n",
"\n",
" for i, genotype in enumerate(genotypes):\n",
" wi = (fitness(genotype, contribs, fitness_mem) - min_w) / (max_w - min_w)\n",
" maximum = True\n",
" minimum = True\n",
" for g in neighbors(genotype, genotypes): \n",
" w = (fitness(g, contribs, fitness_mem) - min_w) / (max_w - min_w)\n",
" if w > wi: \n",
" maximum = False\n",
" ax.plot(\n",
" [genotype.sum(), g.sum()], \n",
" [wi, w], \n",
" ls='-', color='k', marker='', alpha=0.2\n",
" )\n",
" if w < wi: \n",
" minimum = False\n",
" if maximum:\n",
" ax.plot(genotype.sum(), wi, '^', color=red, markersize=10)\n",
" elif minimum:\n",
" ax.plot(genotype.sum(), wi, 'v', color=green, markersize=10)\n",
" else:\n",
" ax.plot(genotype.sum(), wi, '.', color=blue, markersize=10)\n",
" ax.set(\n",
" xlabel='Distance from {}'.format('0'*N),\n",
" ylabel='Fitness',\n",
" )\n",
" sns.despine()\n",
" return ax\n",
"\n",
"#################################################################################################\n",
"# modified code based on Yoav Ram code\n",
"#################################################################################################\n",
"\n",
"def make_NK_model(N, K):\n",
" # for given N and K, make simulations\n",
" genotypes = all_genotypes(N)\n",
"\n",
" # dictionary of interacting loci\n",
" contribs = {i: sorted(np.random.choice([n for n in range(N) if n != i], K, replace=False).tolist() + [i]) for i in range(N)}\n",
" \n",
" fitness_mem = {}\n",
" \n",
" # fitness value of all binary genotypes of length N, consider casting to np.float32 to save space\n",
" ws = [fitness(g, contribs, fitness_mem) for g in genotypes]\n",
"\n",
" return ws, contribs, fitness_mem\n",
"\n",
"# simulate and save synthetic data \n",
"def make_NK_sim(N,nsim=10):\n",
" K = N-1\n",
" ff = np.zeros((nsim*N, 2+2**N)) # empty array to store fitness of genotypes\n",
" to_json = [] # list to store tuple of fitness, contribs and fitness_mem\n",
" col = [tuple(tt) for tt in all_genotypes(N)*1]\n",
" # col = [((\",\".join(map(str, tt)))) for tt in all_genotypes(N)*1]\n",
" col.insert(0, \"nsim\")\n",
" col.insert(1,\"K\")\n",
" # print(col)\n",
"\n",
" # for all K = 0 to N-1\n",
" for KK in range(N):\n",
" # print(\"K \", KK)\n",
" \n",
" for i in range(nsim*KK, nsim*KK + nsim):\n",
" # print(\" i \", i)\n",
" ws, contribs, fitness_mem = make_NK_model(N,KK)\n",
" ff[i,2:] = ws\n",
" ff[i,0] = np.mod(i,nsim)\n",
" ff[i,1] = KK\n",
"\n",
" df = pd.DataFrame(data = ff, columns = col )\n",
" # df.to_csv(\"N{}nsim{}.csv\".format(N,nsim))\n",
" return df\n",
"\n",
"def make_NK_graph(N,df,symbols = 2):\n",
" '''\n",
" Make list of networkx graphs with genotypes as vertices and edge weights as absolute fitness difference between vertices\n",
"\n",
" df - dataframe from make_NK_sim()\n",
"\n",
" '''\n",
"\n",
" NK_graph_list = [] # store graph objects\n",
"\n",
" # N = len(df.columns[2])\n",
"\n",
" for i in range(len(df)):\n",
" graph = make_skeleton_graph(N, symbols = symbols)\n",
"\n",
" # dictionary of sequences and fitness at row i of df\n",
" fff = dict(df.iloc[i,:])\n",
"\n",
" # assign fitness to genotype vertices\n",
" for v in graph.nodes():\n",
" graph.nodes[v]['fitness'] = fff[v]\n",
"\n",
" # reset edge weights\n",
" nx.set_edge_attributes(graph, np.Inf, 'weight')\n",
"\n",
" for e in graph.edges():\n",
" # print(test.nodes[e[0]])\n",
" fa = graph.nodes[e[0]]['fitness']\n",
" fb = graph.nodes[e[1]]['fitness']\n",
" # print(fa)\n",
"\n",
" if fa < 1E-10 or fb < 1E-10:\n",
" continue\n",
" # graph[e[0]][e[1]]['weight'] = np.Inf\n",
" else:\n",
" graph[e[0]][e[1]]['weight'] = abs(fa-fb)\n",
" # print(test[e[0]][e[1]]['weight'])\n",
"\n",
" NK_graph_list.append(graph)\n",
" return NK_graph_list\n",
"\n",
"def sanity_check(NK_graph_list, dim):\n",
" '''\n",
" Compute bottleneck, Wasserstein distances between persistence diagrams from list of graphs\n",
" '''\n",
" nn = len(NK_graph_list)\n",
" \n",
" a = time.time()\n",
"\n",
" # store pairwise distances for number of dimensions\n",
" bottleneck = np.zeros((nn,nn,len(dim)))\n",
" wasserstein = np.zeros((nn,nn,len(dim)))\n",
"\n",
" for i in range(nn):\n",
" g1 = NK_graph_list[i]\n",
" dm1 = nx.to_scipy_sparse_matrix(g1)\n",
" dgm1 = ripser(dm1, maxdim = len(dim)-1, distance_matrix=True)['dgms']\n",
" #plot_diagrams(dgm1, title = \"dgm1\" ,show=True)\n",
"\n",
" for j in range(i+1,nn):\n",
" # print(\" i j \", i , j)\n",
" g2 = NK_graph_list[j]\n",
" dm2 = nx.to_scipy_sparse_matrix(g2)\n",
" \n",
" dgm2 = ripser(dm2, maxdim = len(dim)-1, distance_matrix=True)['dgms']\n",
" #plot_diagrams(dgm2, title = \"dgm2\" ,show=True)\n",
"\n",
" # compute distances\n",
" for dd in dim:\n",
" print(\"dimension \", dd)\n",
" #print(dgm1[dim], dgm2[dim])\n",
" distance_bottleneck = persim.bottleneck(dgm1[dd], dgm2[dd], matching=False)\n",
" #print('bottle neck ', time.time() - a)\n",
" bottleneck[i,j,dim] = distance_bottleneck\n",
"\n",
" d_was = persim.wasserstein(dgm1[dd], dgm2[dd], matching=False)\n",
" #print(d_was)\n",
" wasserstein[i,j,dd] = d_was\n",
" #print('wasserstein ', time.time() - a)\n",
" \n",
" return bottleneck, wasserstein\n",
"\n"
],
"execution_count": 58,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "tqjIBxHGzypL"
},
"source": [
"N = 3\n",
"nsim = 5\n",
"df = make_NK_sim(N,nsim)"
],
"execution_count": 59,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "GRWMEWYNEiNf",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 630
},
"outputId": "356c84f9-441c-4bfb-9530-126cf4b45482"
},
"source": [
"NK_graph_list = make_NK_graph(N,df)\n",
"%time b,w = sanity_check(NK_graph_list, dim=[0,1,2])"
],
"execution_count": 60,
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"dimension 0\n",
"dimension 1\n"
]
},
{
"output_type": "stream",
"name": "stderr",
"text": [
"/usr/local/lib/python3.7/dist-packages/persim/bottleneck.py:57: UserWarning: dgm1 has points with non-finite death times;ignoring those points\n",
" \"dgm1 has points with non-finite death times;\"+\n",
"/usr/local/lib/python3.7/dist-packages/persim/bottleneck.py:67: UserWarning: dgm2 has points with non-finite death times;ignoring those points\n",
" \"dgm2 has points with non-finite death times;\"+\n",
"/usr/local/lib/python3.7/dist-packages/persim/wasserstein.py:52: UserWarning: dgm1 has points with non-finite death times;ignoring those points\n",
" \"dgm1 has points with non-finite death times;\"+\n",
"/usr/local/lib/python3.7/dist-packages/persim/wasserstein.py:62: UserWarning: dgm2 has points with non-finite death times;ignoring those points\n",
" \"dgm2 has points with non-finite death times;\"+\n",
"/usr/local/lib/python3.7/dist-packages/persim/bottleneck.py:57: UserWarning: dgm1 has points with non-finite death times;ignoring those points\n",
" \"dgm1 has points with non-finite death times;\"+\n",
"/usr/local/lib/python3.7/dist-packages/persim/bottleneck.py:67: UserWarning: dgm2 has points with non-finite death times;ignoring those points\n",
" \"dgm2 has points with non-finite death times;\"+\n"
]
},
{
"output_type": "error",
"ename": "IndexError",
"evalue": "ignored",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-60-c0f2279087cf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mNK_graph_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmake_NK_graph\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmagic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'time b,w = sanity_check(NK_graph_list, dim=[0,1,2])'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mmagic\u001b[0;34m(self, arg_s)\u001b[0m\n\u001b[1;32m 2158\u001b[0m \u001b[0mmagic_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marg_s\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpartition\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m' '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2159\u001b[0m \u001b[0mmagic_name\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmagic_name\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlstrip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprefilter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mESC_MAGIC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2160\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_name\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2161\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2162\u001b[0m \u001b[0;31m#-------------------------------------------------------------------------\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[0;34m(self, magic_name, line)\u001b[0m\n\u001b[1;32m 2079\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'local_ns'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getframe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_locals\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2080\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2081\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2082\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2083\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<decorator-gen-53>\u001b[0m in \u001b[0;36mtime\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m 186\u001b[0m \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 187\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 188\u001b[0;31m \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 190\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mtime\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[1;32m 1191\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1192\u001b[0m \u001b[0mst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mclock2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1193\u001b[0;31m \u001b[0mexec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mglob\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlocal_ns\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1194\u001b[0m \u001b[0mend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mclock2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1195\u001b[0m \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<timed exec>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n",
"\u001b[0;32m<ipython-input-58-c98a6ec3b03f>\u001b[0m in \u001b[0;36msanity_check\u001b[0;34m(NK_graph_list, dim)\u001b[0m\n\u001b[1;32m 181\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"dimension \"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdd\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 182\u001b[0m \u001b[0;31m#print(dgm1[dim], dgm2[dim])\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 183\u001b[0;31m \u001b[0mdistance_bottleneck\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpersim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbottleneck\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdgm1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdd\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdgm2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdd\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmatching\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 184\u001b[0m \u001b[0;31m#print('bottle neck ', time.time() - a)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 185\u001b[0m \u001b[0mbottleneck\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdistance_bottleneck\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/persim/bottleneck.py\u001b[0m in \u001b[0;36mbottleneck\u001b[0;34m(dgm1, dgm2, matching)\u001b[0m\n\u001b[1;32m 103\u001b[0m \u001b[0;31m# bottleneck distance\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 104\u001b[0m \u001b[0mds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mD\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m# Everything but np.inf\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 105\u001b[0;31m \u001b[0mbdist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 106\u001b[0m \u001b[0mmatching\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 107\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mIndexError\u001b[0m: index -1 is out of bounds for axis 0 with size 0"
]
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment