Skip to content
Snippets Groups Projects
Commit 4287baeb authored by Brad Nelson's avatar Brad Nelson
Browse files

modified homology alg

parent 208b6a22
Branches
No related tags found
No related merge requests found
......@@ -21,7 +21,7 @@ dgm, issublevel = layer1(y)
p = sum_finite(dgm[0])
p.backward()
algs = ['hom', 'cohom']
algs = ['hom', 'hom2', 'cohom']
tcs = {}
tfs = {}
......
......@@ -21,7 +21,7 @@ dgm, issublevel = layer1(y)
p = sum_finite(dgm[0])
p.backward()
algs = ['hom', 'cohom']
algs = ['hom', 'hom2', 'cohom']
tcs = {}
tfs = {}
......
from __future__ import print_function
from topologylayer.functional.cohom_cpp import SimplicialComplex, persistenceForward
from topologylayer.functional.persistence import SimplicialComplex, persistenceForwardCohom
from topologylayer.util.process import remove_zero_bars
import torch
......@@ -43,7 +43,7 @@ f = torch.Tensor([2., 0., 0., 0., 0.])
s.extendFloat(f)
# compute persistence with MAXDIM=1
ret = persistenceForward(s, 1)
ret = persistenceForwardCohom(s, 1)
for k in range(2):
print("dimension %d bars" % k)
......
import unittest
from topologylayer.functional.persistence import SimplicialComplex, persistenceForward
from topologylayer.functional.persistence import SimplicialComplex, persistenceForwardCohom
from topologylayer.util.process import remove_zero_bars
import torch
import numpy as np
......@@ -46,7 +46,7 @@ class BasicLevelset(unittest.TestCase):
s.extendFloat(f)
# compute persistence with MAXDIM=1
ret = persistenceForward(s, 1)
ret = persistenceForwardCohom(s, 1)
self.assertEqual(
torch.all(torch.eq(remove_zero_bars(ret[0]), torch.tensor([[0., np.inf]]))),
......
......@@ -10,7 +10,7 @@ class AlphaTest(unittest.TestCase):
from topologylayer.nn import AlphaLayer
# superlevel set
for alg in ['hom', 'cohom']:
for alg in ['hom', 'hom2', 'cohom']:
layer = AlphaLayer(maxdim=1, alg=alg)
x = torch.tensor([[1, 1], [1,-1], [-1,-1], [-1,1]], dtype=torch.float).requires_grad_(True)
......
......@@ -11,7 +11,7 @@ class Levelset1dsuper(unittest.TestCase):
from topologylayer.nn import LevelSetLayer1D
# superlevel set
for alg in ['hom', 'cohom']:
for alg in ['hom', 'hom2', 'cohom']:
layer = LevelSetLayer1D(size=3, sublevel=False, alg=alg)
y = torch.tensor([1,0,1], dtype=torch.float).requires_grad_(True)
......@@ -40,7 +40,7 @@ class Levelset1dsub(unittest.TestCase):
from topologylayer.nn import LevelSetLayer1D
# sublevel set
for alg in ['hom', 'cohom']:
for alg in ['hom', 'hom2', 'cohom']:
layer = LevelSetLayer1D(size=3, sublevel=True, alg=alg)
y = torch.tensor([1,0,1], dtype=torch.float).requires_grad_(True)
......@@ -69,7 +69,7 @@ class Levelset2dsuper(unittest.TestCase):
from topologylayer.nn import LevelSetLayer2D
# sublevel set
for alg in ['hom', 'cohom']:
for alg in ['hom', 'hom2', 'cohom']:
layer = LevelSetLayer2D(size=(3,3), maxdim=1, sublevel=False, alg=alg)
x = torch.tensor([[2, 1, 1],[1, 0.5, 1],[1, 1, 1]], dtype=torch.float).requires_grad_(True)
......
......@@ -9,7 +9,7 @@ class RipsTest(unittest.TestCase):
def test(self):
from topologylayer.nn import RipsLayer
for alg in ['hom', 'cohom']:
for alg in ['hom', 'hom2', 'cohom']:
# superlevel set
layer = RipsLayer(4, maxdim=1, alg=alg)
x = torch.tensor([[1, 1], [1,-1], [-1,-1], [-1,1]], dtype=torch.float).requires_grad_(True)
......
from __future__ import print_function
from torch.autograd import Variable, Function
from .persistence import SimplicialComplex, persistenceForward, persistenceBackwardFlag, persistenceForwardHom
from .persistence import SimplicialComplex, persistenceForwardCohom, persistenceBackwardFlag, persistenceForwardHom
class FlagDiagram(Function):
"""
......@@ -13,15 +13,18 @@ class FlagDiagram(Function):
maxdim - maximum homology dimension
alg - algorithm
'hom' = homology (default)
'hom2' = nz suppressing homology variant
'cohom' = cohomology
"""
@staticmethod
def forward(ctx, X, y, maxdim, alg='hom'):
X.extendFlag(y)
if alg == 'hom':
ret = persistenceForwardHom(X, maxdim)
ret = persistenceForwardHom(X, maxdim, 0)
elif alg == 'hom2':
ret = persistenceForwardHom(X, maxdim, 1)
elif alg == 'cohom':
ret = persistenceForward(X, maxdim)
ret = persistenceForwardCohom(X, maxdim)
ctx.X = X
ctx.save_for_backward(y)
return tuple(ret)
......
......@@ -19,13 +19,14 @@ INPUTS:
MAXDIM - maximum homology dimension
OUTPUTS: boundary matrix in List of Lists format.
*/
std::vector<SparseF2Vec<int>> sorted_boundary(SimplicialComplex &X) {
std::vector<SparseF2Vec<int>> sorted_boundary(SimplicialComplex &X, size_t MAXDIM) {
// empty boundary matrix
std::vector<SparseF2Vec<int>> B;
// build boundary in sorted order using filtration_perm
// should also use filtration_perm to permute nzs in rows of columns
std::vector<int> row_inds; // row indices for column
for (size_t j : X.filtration_perm ) {
if (X.dim(j) > MAXDIM+1) { continue; }
row_inds.clear(); // clear out column
// go through non-zeros in boundary
for (auto i : X.bdr[j].cochain.nzinds) {
......@@ -67,6 +68,40 @@ void homology_reduction_alg(std::vector<SparseF2Vec<int>> &B, std::map<int, int>
return;
}
/*
Function to reduce boundary matrix.
Tries to minimize nnz by continuing to reduce past pivot
INPUT:
B - boundary matrix
pivot_to_col - map from pivot to column
OUTPUT:
none - inputs are modified in-place.
*/
void homology_reduction_alg2(std::vector<SparseF2Vec<int>> &B, std::map<int, int> &pivot_to_col) {
// loop over columns of boundary matrix
for (size_t j = 0; j < B.size(); j++) {
// if nnz = 0, the reduction is complete
size_t offset = 0; // position from last
while (B[j].nnz() > offset) {
int piv = B[j].from_end(offset);
if (pivot_to_col.count(piv) > 0) {
int k = pivot_to_col[piv];
// there is a column with that pivot
B[j].add(B[k]);
} else {
// there is no column with that pivot
// see if we've found new pivot
if (offset == 0) {pivot_to_col[piv] = j;}
// increase offset
offset++;
}
} // end column reduction
} // end for loop
return;
}
/*
Standard reduction algorithm on simplicial complex.
......@@ -74,10 +109,12 @@ void homology_reduction_alg(std::vector<SparseF2Vec<int>> &B, std::map<int, int>
X - simplicial complex
IMPORTANT: assumes that X has been initialized, and filtration has been extended
MAXDIM - maximum homology dimension
alg_opt - 0 - standard reduction algorithm
1 - nz minimizing reduction algorithm
OUTPUTS: vector of tensors - t
t[k] is float32 tensor with barcode for dimension k
*/
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM) {
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM, size_t alg_opt) {
// produce sort permutation on X
X.sortedOrder();
......@@ -95,15 +132,18 @@ std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t
}
// produce boundary matrix
std::vector<SparseF2Vec<int>> B = sorted_boundary(X);
// for (size_t j = 0; j < B.size(); j++) {
// py::print(j);
// B[j].print();
// }
std::vector<SparseF2Vec<int>> B = sorted_boundary(X, MAXDIM);
// run standard reduction algorithm
std::map<int, int> pivot_to_col;
if (alg_opt == 0) {
// standard reduction algorithm
homology_reduction_alg(B, pivot_to_col);
} else {
// modified reduction algorithm
homology_reduction_alg2(B, pivot_to_col);
}
// keep track of how many pairs we've put in diagram
std::vector<size_t> nbars(MAXDIM+1);
......
......@@ -3,4 +3,4 @@
#include "complex.h"
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM);
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM, size_t alg_opt);
......@@ -24,7 +24,7 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
.def("printBoundary", &SimplicialComplex::printBoundary)
.def("numPairs", &SimplicialComplex::numPairs)
.def("printCells", &SimplicialComplex::printComplex);
m.def("persistenceForward", &persistence_forward);
m.def("persistenceForwardCohom", &persistence_forward);
m.def("persistenceBackward", &persistence_backward);
m.def("persistenceBackwardFlag", &persistence_backward_flag);
m.def("persistenceForwardHom", &persistence_forward_hom);
......
......@@ -104,14 +104,21 @@ class SparseF2Vec{
py::print(nzinds);
}
// return number of non-zeros
size_t nnz() {
return nzinds.size();
}
// return last nonzero index
T last() {
return nzinds.back();
}
// return offset element from last
T from_end(size_t offset) {
return nzinds.at(nzinds.size() - 1 - offset);
}
};
......
......@@ -3,7 +3,7 @@ from __future__ import print_function
import torch
from torch.autograd import Variable, Function
from .persistence import SimplicialComplex, persistenceForward, persistenceBackward, persistenceForwardHom
from .persistence import SimplicialComplex, persistenceForwardCohom, persistenceBackward, persistenceForwardHom
class SubLevelSetDiagram(Function):
"""
......@@ -14,6 +14,7 @@ class SubLevelSetDiagram(Function):
maxdim - maximum homology dimension
alg - algorithm
'hom' = homology (default)
'hom2' = nz suppressing homology variant
'cohom' = cohomology
"""
@staticmethod
......@@ -22,9 +23,11 @@ class SubLevelSetDiagram(Function):
f = f.view(-1)
X.extendFloat(f)
if alg == 'hom':
ret = persistenceForwardHom(X, maxdim)
ret = persistenceForwardHom(X, maxdim, 0)
elif alg == 'hom2':
ret = persistenceForwardHom(X, maxdim, 1)
elif alg == 'cohom':
ret = persistenceForward(X, maxdim)
ret = persistenceForwardCohom(X, maxdim)
ctx.X = X
return tuple(ret)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment