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

first complete homology reduction implementation

parent e41000e6
No related branches found
No related tags found
No related merge requests found
......@@ -41,7 +41,7 @@ void reduction_step(SimplicialComplex &X,\
size_t bindx = pivot->index;
size_t dindx = c.index;
// get birth dimension
size_t hdim = X.bdr[bindx].dim();
size_t hdim = X.dim(bindx);
// get location in diagram
size_t j = nbars[hdim]++;
......
......@@ -50,6 +50,12 @@ void SimplicialComplex::printFiltration() {
}
}
// return dimension of cell j
size_t SimplicialComplex::dim(size_t j) {
return cells[j].size() - 1;
//return X.bdr[bindx].dim()
}
void SimplicialComplex::initialize() {
// to call after complex is built.
......
......@@ -92,6 +92,9 @@ class SimplicialComplex{
// fill in filtration order
void sortedOrder();
// get dimension of cell j
size_t dim(size_t j);
// extend a filtration on vertices to whole complex
// template <typename T>
// std::vector<T> extend_lower_star(std::vector<T>);
......
......@@ -50,21 +50,18 @@ OUTPUT:
void homology_reduction_alg(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++) {
while (true) {
int piv = B[j].pivot();
if (piv == std::numeric_limits<int>::max()) {
// we have completely reduced column
break;
} else {
// if nnz = 0, the reduction is complete
while (B[j].nnz() > 0) {
int piv = B[j].last();
if (pivot_to_col.count(piv) > 0) {
int k = pivot_to_col[piv];
// there is a column with that pivot
B[j].add(B[piv]);
B[j].add(B[k]);
} else {
// there is no column with that pivot
pivot_to_col[piv] = j;
break;
}
}
} // end column reduction
} // end for loop
return;
......@@ -80,14 +77,14 @@ void homology_reduction_alg(std::vector<SparseF2Vec<int>> &B, std::map<int, int>
OUTPUTS: vector of tensors - t
t[k] is float32 tensor with barcode for dimension k
*/
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, int MAXDIM) {
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM) {
// produce sort permutation on X
X.sortedOrder();
// initialize reutrn diagram
std::vector<torch::Tensor> diagram(MAXDIM+1); // return array
for (int k = 0; k < MAXDIM+1; k++) {
for (size_t k = 0; k < MAXDIM+1; k++) {
int Nk = X.numPairs(k); // number of bars in dimension k
// allocate return tensor
diagram[k] = torch::empty({Nk,2},
......@@ -96,20 +93,56 @@ std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, int MAX
// TODO: do this in intialization since number of pairs is deterministic
X.backprop_lookup[k].resize(Nk);
}
// keep track of how many pairs we've put in diagram
std::vector<int> nbars(MAXDIM+1);
for (int k = 0; k < MAXDIM+1; k++) {
nbars[k] = 0;
}
// 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();
// }
// run standard reduction algorithm
std::map<int, int> pivot_to_col;
homology_reduction_alg(B, pivot_to_col);
// keep track of how many pairs we've put in diagram
std::vector<size_t> nbars(MAXDIM+1);
for (size_t k = 0; k < MAXDIM+1; k++) {
nbars[k] = 0;
}
// fill in diagram
for (size_t j = 0; j < B.size(); j++ ) {
// see if column j is completely reduced
if (B[j].nnz() == 0) {
// homology born in column j
// in filtration order:
// birth at j, death at k = pivot_to_col[j]
size_t bindx = X.filtration_perm[j];
size_t hdim = X.dim(bindx);
if (hdim > MAXDIM) { continue; }
// get location in diagram
size_t di = nbars[hdim]++;
// set birth time
(diagram[hdim][di].data<float>())[0] = (float) X.full_function[bindx].first;
if (pivot_to_col.count(j) > 0) {
// there is a finite death.
int k = pivot_to_col[j]; // column that has j as pivot
size_t dindx = X.filtration_perm[k];
// get dimension of cell with filtration poisition j
(diagram[hdim][di].data<float>())[1] = (float) X.full_function[dindx].first;
// need to fill in X.backprop_lookup as well
X.backprop_lookup[hdim][di] = {(int) bindx, (int) dindx};
} else {
// infinite death
(diagram[hdim][di].data<float>())[1] = std::numeric_limits<float>::infinity();
X.backprop_lookup[hdim][di] = {(int) bindx, -1};
}
}
}
return diagram;
}
......@@ -3,4 +3,4 @@
#include "complex.h"
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, int MAXDIM);
std::vector<torch::Tensor> persistence_forward_hom(SimplicialComplex &X, size_t MAXDIM);
......@@ -16,6 +16,7 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
.def("extendFloat", &SimplicialComplex::extend)
.def("extendFlag", &SimplicialComplex::extend_flag)
.def("sortedOrder", &SimplicialComplex::sortedOrder)
.def("dim", &SimplicialComplex::dim)
.def("printFiltration", &SimplicialComplex::printFiltration)
.def("printFunctionMap", &SimplicialComplex::printFunctionMap)
.def("printCritInds", &SimplicialComplex::printCritInds)
......
......@@ -44,8 +44,8 @@ class SparseF2Vec{
size_t i1 = 0;
size_t i2 = 0;
do {
size_t v1 = nzinds[i1];
size_t v2 = x.nzinds[i2];
T v1 = nzinds[i1];
T v2 = x.nzinds[i2];
if (v1 == v2) {
// F2 means sum is 0
i1++;
......@@ -83,8 +83,8 @@ class SparseF2Vec{
size_t i2 = 0;
size_t intersection = 0;
do {
auto v1 = nzinds[i1];
auto v2 = x.nzinds[i2];
T v1 = nzinds[i1];
T v2 = x.nzinds[i2];
if (v1 == v2) {
i1++;
i2++;
......@@ -104,11 +104,13 @@ class SparseF2Vec{
py::print(nzinds);
}
// returns last non-zero row in column
T pivot() {
return nzinds.size() == 0 ? std::numeric_limits<T>::max() : nzinds.back();
size_t nnz() {
return nzinds.size();
}
T last() {
return nzinds.back();
}
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment