From 46cef577ea816f9c6d75038f6376111d71d38804 Mon Sep 17 00:00:00 2001 From: Brenda Praggastis Date: Thu, 20 Feb 2020 09:48:10 -0500 Subject: [PATCH 1/5] improved performance of hom mod2 thanks to Andrew --- hypernetx/algorithms/homology_mod2.py | 67 +++++++++---------- .../algorithms/tests/test_homology_mod_2.py | 22 +++--- 2 files changed, 45 insertions(+), 44 deletions(-) diff --git a/hypernetx/algorithms/homology_mod2.py b/hypernetx/algorithms/homology_mod2.py index be7e51cb..f364fd0f 100644 --- a/hypernetx/algorithms/homology_mod2.py +++ b/hypernetx/algorithms/homology_mod2.py @@ -34,6 +34,9 @@ import warnings, copy from hypernetx import HyperNetXError +import datetime as dtm +tmstmp = lambda x: print(f'{x}: {dtm.datetime.now().strftime("%H.%M.%S")}') + def kchainbasis(h,k): """ Compute the set of k dimensional cells in the abstract simplicial @@ -60,6 +63,7 @@ def kchainbasis(h,k): - Hypergraph node uids must be sortable. """ + tmstmp('kchbasis-s') import itertools as it kchains = set() for e in h.edges(): @@ -70,6 +74,7 @@ def kchainbasis(h,k): return sorted(list(kchains)) def interpret(Ck,arr): + """ Returns the data as represented in Ck associated with the arr @@ -86,6 +91,7 @@ def interpret(Ck,arr): list of k-cells referenced by data in Ck """ + output = list() for vec in arr: if len(Ck) != len(vec): @@ -110,7 +116,7 @@ def bkMatrix(km1basis,kbasis): bk : np.array boundary matrix in $Z_2$ """ - + tmstmp('bkMatrix-s') bk = np.zeros((len(km1basis),len(kbasis)),dtype=int) for cell in kbasis: for idx in range(len(cell)): @@ -301,6 +307,7 @@ def matmulreduce(arr,reverse=False,mod=2): P : np.array Product of matrices in the list """ + tmstmp('matmulr-s') if reverse: items = range(len(arr)-1,-1,-1) else: @@ -333,6 +340,7 @@ def matsumreduce(arr,mod=2): ## Convenience methods for computing Smith Normal Form ## All of these operations have themselves as inverses + def _sr(i,j,M,L): return swap_rows(i,j,M,L) @@ -409,6 +417,8 @@ def smith_normal_form_mod2(M,track=False): where $L = L[s]L[s-1]...L[1]L[0]$ and $R = R[0]R[1]...R[t]$. """ + tmstmp('snf-start') + start = dtm.datetime.now() S = copy.copy(M) dimL,dimR = M.shape mod = 2 @@ -420,8 +430,7 @@ def smith_normal_form_mod2(M,track=False): L = np.eye(dimL,dtype=int) R = np.eye(dimR,dtype=int) - Linv = [IL] ## keep track of transformations to compute inverses at the end - Rinv = [IR] + Linv = np.eye(dimL,dtype=int) if track: print(L,'L\n') @@ -429,6 +438,7 @@ def smith_normal_form_mod2(M,track=False): print(R,'R\n') for s in range(min(dimL,dimR)): + print('.'), if track: print(f'\ns={s}\n') ## Find index pair (rdx,cdx) with value 1 in submatrix M[s:,s:] @@ -442,42 +452,28 @@ def smith_normal_form_mod2(M,track=False): ## Swap rows and columns as needed so that 1 is in the s,s position if rdx > s: S,L = _sr(s,rdx,S,L) -# This is equivalent to -# S = swap_rows(s,rdx,S)[0] -# L = swap_rows(s,rdx,L)[0] - Linv.append(swap_rows(s,rdx,IL)[0]) ## keep track of transformations on the left to reverse them for the inverse + Linv = swap_columns(s,rdx,Linv)[0] if track: print(L,f'L\n',S,'S\n',) - assert(np.all(S == matmulreduce([L,M,R],mod=mod))) if cdx > s: S,R= _sc(s,cdx,S,R) - Rinv.append(swap_columns(s,cdx,IR)[0]) if track: print(S,'S\n',R,f'R\n') - assert(np.all(S == matmulreduce([L,M,R],mod=mod))) # add sth row to every row with 1 in sth column & sth column to every column with 1 in sth row - row_indices = [idx for idx in range(dimL) if idx != s and S[idx][s] == 1] + row_indices = [idx for idx in range(s+1,dimL) if S[idx][s] == 1] for rdx in row_indices: S,L = _ar(rdx,s,S,L,mod=mod) - temp = add_to_row(IL,rdx,s,mod=mod) - Linv.append(temp) + Linv = add_to_column(Linv,s,rdx,mod=mod) if track: print(L,f'L\n',S,'S\n',) - assert(np.all(S == matmulreduce([L,M,R],mod=mod))) - column_indices = [jdx for jdx in range(dimR) if jdx != s and S[s][jdx] == 1] + column_indices = [jdx for jdx in range(s+1,dimR) if S[s][jdx] == 1] for jdx,cdx in enumerate(column_indices): S,R = _ac(cdx,s,S,R) - temp = add_to_column(IR,cdx,s,mod=mod) - Rinv.append(temp) if track: print(R,f'R\n',S,'S\n',) - assert(np.all(S == matmulreduce([L,M,R],mod=mod))) - Linv = matmulreduce(Linv,mod=mod) - Rinv = matmulreduce(Rinv,reverse=True,mod=mod) - assert(np.all(IL == matmulreduce([L,Linv],mod=mod))) - assert(np.all(IR == matmulreduce([R,Rinv],mod=mod))) - return L,R,S,Linv,Rinv + print(dtm.datetime.now() - start) + return L,R,S,Linv def reduced_row_echelon_form_mod2(M,track=False,mod=2): @@ -511,6 +507,7 @@ def reduced_row_echelon_form_mod2(M,track=False,mod=2): verify the product: L[s]L[s-1]...L[0] M = S . """ + tmstmp('ref-s') S = copy.copy(M) dimL,dimR = M.shape mod = 2 @@ -518,7 +515,7 @@ def reduced_row_echelon_form_mod2(M,track=False,mod=2): ## method with numpy IL = np.eye(dimL,dtype=int) - Linv = [np.eye(dimL,dtype=int)] + Linv = np.eye(dimL,dtype=int) L = np.eye(dimL,dtype=int) if track: @@ -540,22 +537,17 @@ def reduced_row_echelon_form_mod2(M,track=False,mod=2): continue if rdx > s1: S,L = _sr(s1,rdx,S,L) - Linv.append(swap_rows(s1,rdx,IL)[0]) ## keep track of transformations on the left to reverse them for the inverse + Linv = swap_columns(rdx,s1,Linv)[0] if track: print(L,f'L\n',S,'S\n',) - assert(np.all(S == matmulreduce([L,M],mod=mod))) - # add sth row to every nonzero row row_indices = [idx for idx in range(dimL) if idx != s1 and S[idx][cdx] == 1] for idx in row_indices: S,L = _ar(idx,s1,S,L,mod=mod) - temp = add_to_row(IL,idx,s1,mod=mod) - Linv.append(temp) + Linv = add_to_column(Linv,s1,idx,mod=mod) if track: print(L,f'L\n',S,'S\n',) - assert(np.all(S == matmulreduce([L,M],mod=mod))) - Linv = matmulreduce(Linv,mod=mod) - assert(np.all(IL == matmulreduce([L,Linv],mod=mod))) + tmstmp('ref-end') return L,S,Linv ## Private @@ -674,8 +666,8 @@ def homology_basis(bd,k,C=None,shortest=False): k-chains if shortest then returns a dictionary of shortest cycles for each coset. """ - L1,R1,S1,L1inv,R1inv = smith_normal_form_mod2(bd[k]) - L2,R2,S2,L2inv,R2inv = smith_normal_form_mod2(bd[k+1]) + L1,R1,S1,L1inv = smith_normal_form_mod2(bd[k]) + L2,R2,S2,L2inv = smith_normal_form_mod2(bd[k+1]) rank1 = np.sum(S1); print(f"rank{k} = {rank1}") rank2 = np.sum(S2); print(f"rank{k+1} = {rank2}") @@ -686,9 +678,12 @@ def homology_basis(bd,k,C=None,shortest=False): ker1 = R1[:,rank1:] im2 = L2inv[:,:rank2] cokernel2 = L2inv[:,rank2:] + cokproj2 = L2[rank2:,:] - proj = matmulreduce([L2,ker1])[rank2:,:] - proj = matmulreduce([L2inv[:,rank2:],proj]).transpose() + # proj = matmulreduce([L2,ker1])[rank2:,:] + # proj = matmulreduce([L2inv[:,rank2:],proj]).transpose() + tmstmp('proj-s') + proj = matmulreduce([cokernel2,cokproj2,ker1]).transpose() _,proj,_ = reduced_row_echelon_form_mod2(proj) proj = np.array([row for row in proj if np.sum(row)>0]) print('hom basis reps\n',proj) diff --git a/hypernetx/algorithms/tests/test_homology_mod_2.py b/hypernetx/algorithms/tests/test_homology_mod_2.py index 57f1a61f..703eb24b 100644 --- a/hypernetx/algorithms/tests/test_homology_mod_2.py +++ b/hypernetx/algorithms/tests/test_homology_mod_2.py @@ -24,22 +24,28 @@ def test_bkMatrix(triloop): def test_smith_normal_form_mod2(triloop): Ck = {k:hnx.kchainbasis(triloop.hypergraph,k) for k in range(0,2)} bd = bkMatrix(Ck[0],Ck[1]) - P1,Q1,S1,P1inv,Q1inv = smith_normal_form_mod2(bd) + P1,Q1,S1,P1inv = smith_normal_form_mod2(bd) assert np.array_equal(P1[0],np.array([1, 0, 0, 0])) assert np.array_equal(Q1[:,2],np.array([0, 1, 1, 0, 0])) + assert(np.all(S1 == matmulreduce([P1,bd,Q1],mod=2))) + r = len(P1) + assert(np.all(np.eye(r) == modmult(P1,P1inv,mod=2))) + assert(np.all(S1 == matmulreduce([P1,bd,Q1],mod=2))) def test_reduced_row_echelon_form_mod2(): m = np.array([[0,1,0,1,0,0,1,0,0,1], [0,0,0,0,0,0,0,0,0,1], [1,0,0,1,0,0,0,0,0,0]]) - L,rm,R = reduced_row_echelon_form_mod2(m) - assert np.array_equal(rm,np.array([[1, 0, 0, 1, 0, 0, 0, 0, 0, 0], - [0, 1, 0, 1, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])) + r = 3 + L,S,Linv = reduced_row_echelon_form_mod2(m) + assert np.array_equal(S,np.array([[1, 0, 0, 1, 0, 0, 0, 0, 0, 0], + [0, 1, 0, 1, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])) assert np.array_equal(L,np.array([[0, 0, 1], - [1, 1, 0], - [0, 1, 0]])) - + [1, 1, 0], + [0, 1, 0]])) + assert(np.all(S == modmult(L,m,mod=2))) + assert(np.all(np.eye(3) == modmult(L,Linv,mod=2))) def test_homology_basis(triloop): Ck = {k:hnx.kchainbasis(triloop.hypergraph,k) for k in range(0,3)} From 8d5cc2f0b3d52790eb9b3bc94d86d958aff47876 Mon Sep 17 00:00:00 2001 From: Brenda Praggastis Date: Thu, 20 Feb 2020 09:51:36 -0500 Subject: [PATCH 2/5] removed timings --- hypernetx/algorithms/homology_mod2.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/hypernetx/algorithms/homology_mod2.py b/hypernetx/algorithms/homology_mod2.py index f364fd0f..690d77aa 100644 --- a/hypernetx/algorithms/homology_mod2.py +++ b/hypernetx/algorithms/homology_mod2.py @@ -35,7 +35,7 @@ from hypernetx import HyperNetXError import datetime as dtm -tmstmp = lambda x: print(f'{x}: {dtm.datetime.now().strftime("%H.%M.%S")}') + def kchainbasis(h,k): """ @@ -63,7 +63,7 @@ def kchainbasis(h,k): - Hypergraph node uids must be sortable. """ - tmstmp('kchbasis-s') + import itertools as it kchains = set() for e in h.edges(): @@ -116,7 +116,7 @@ def bkMatrix(km1basis,kbasis): bk : np.array boundary matrix in $Z_2$ """ - tmstmp('bkMatrix-s') + bk = np.zeros((len(km1basis),len(kbasis)),dtype=int) for cell in kbasis: for idx in range(len(cell)): @@ -307,7 +307,7 @@ def matmulreduce(arr,reverse=False,mod=2): P : np.array Product of matrices in the list """ - tmstmp('matmulr-s') + if reverse: items = range(len(arr)-1,-1,-1) else: @@ -417,7 +417,7 @@ def smith_normal_form_mod2(M,track=False): where $L = L[s]L[s-1]...L[1]L[0]$ and $R = R[0]R[1]...R[t]$. """ - tmstmp('snf-start') + start = dtm.datetime.now() S = copy.copy(M) dimL,dimR = M.shape @@ -507,7 +507,7 @@ def reduced_row_echelon_form_mod2(M,track=False,mod=2): verify the product: L[s]L[s-1]...L[0] M = S . """ - tmstmp('ref-s') + S = copy.copy(M) dimL,dimR = M.shape mod = 2 @@ -547,7 +547,7 @@ def reduced_row_echelon_form_mod2(M,track=False,mod=2): Linv = add_to_column(Linv,s1,idx,mod=mod) if track: print(L,f'L\n',S,'S\n',) - tmstmp('ref-end') + return L,S,Linv ## Private @@ -680,9 +680,6 @@ def homology_basis(bd,k,C=None,shortest=False): cokernel2 = L2inv[:,rank2:] cokproj2 = L2[rank2:,:] - # proj = matmulreduce([L2,ker1])[rank2:,:] - # proj = matmulreduce([L2inv[:,rank2:],proj]).transpose() - tmstmp('proj-s') proj = matmulreduce([cokernel2,cokproj2,ker1]).transpose() _,proj,_ = reduced_row_echelon_form_mod2(proj) proj = np.array([row for row in proj if np.sum(row)>0]) From 6b6805e32acfda6504c3d3fec90fd6ac3cba8f10 Mon Sep 17 00:00:00 2001 From: Brenda Praggastis Date: Thu, 20 Feb 2020 10:05:30 -0500 Subject: [PATCH 3/5] removed unnecessary import --- hypernetx/algorithms/homology_mod2.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/hypernetx/algorithms/homology_mod2.py b/hypernetx/algorithms/homology_mod2.py index 690d77aa..12f73188 100644 --- a/hypernetx/algorithms/homology_mod2.py +++ b/hypernetx/algorithms/homology_mod2.py @@ -34,9 +34,6 @@ import warnings, copy from hypernetx import HyperNetXError -import datetime as dtm - - def kchainbasis(h,k): """ Compute the set of k dimensional cells in the abstract simplicial From 68068edb70d530e90cae82a2af22ca04a012cc29 Mon Sep 17 00:00:00 2001 From: Brenda Praggastis Date: Thu, 20 Feb 2020 10:07:30 -0500 Subject: [PATCH 4/5] removed calls to datetime --- hypernetx/algorithms/homology_mod2.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/hypernetx/algorithms/homology_mod2.py b/hypernetx/algorithms/homology_mod2.py index 12f73188..cb2c2c36 100644 --- a/hypernetx/algorithms/homology_mod2.py +++ b/hypernetx/algorithms/homology_mod2.py @@ -415,7 +415,6 @@ def smith_normal_form_mod2(M,track=False): """ - start = dtm.datetime.now() S = copy.copy(M) dimL,dimR = M.shape mod = 2 @@ -469,7 +468,6 @@ def smith_normal_form_mod2(M,track=False): S,R = _ac(cdx,s,S,R) if track: print(R,f'R\n',S,'S\n',) - print(dtm.datetime.now() - start) return L,R,S,Linv From c107a64489ef3cf704de291e97a560ff321d2349 Mon Sep 17 00:00:00 2001 From: Brenda Praggastis Date: Thu, 20 Feb 2020 10:09:55 -0500 Subject: [PATCH 5/5] updated version for hom optimization --- docs/source/conf.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index a2dba83a..d910b62d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -19,7 +19,7 @@ import os import shlex -__version__ = "0.3.0" +__version__ = "0.3.1" # If extensions (or modules to document with autodoc) are in another directory, diff --git a/setup.py b/setup.py index 273f3fce..6ffb1f62 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup import sys, os.path -__version__ = '0.3.0' +__version__ = '0.3.1' if sys.version_info < (3,6): sys.exit('HyperNetX requires Python 3.6 or later.')