1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/Graph.h"
19 #include "polymake/PowerSet.h"
20 #include "polymake/topaz/Filtration.h"
21 #include "polymake/SparseMatrix.h"
22 #include "polymake/Rational.h"
23
24 namespace polymake { namespace topaz {
25
vietoris_rips_complex(const Matrix<Rational> & dist,Rational step)26 BigObject vietoris_rips_complex(const Matrix<Rational>& dist, Rational step)
27 {
28 BigObject NG = call_function("neighborhood_graph", dist, step);
29 BigObject vr_complex = call_function("clique_complex", NG);
30 vr_complex.set_description() << "Vietoris Rips complex of the input point set." << endl;
31 return vr_complex;
32 }
33
34 // this does not use the existing clique/flag complex functionality, as we need to calculate degrees for each face.
35 template <typename Coeff>
vietoris_rips_filtration(const Matrix<double> & dist,const Array<Int> & point_degs,double step,Int k)36 Filtration<SparseMatrix<Coeff>> vietoris_rips_filtration(const Matrix<double>& dist, const Array<Int>& point_degs, double step, Int k)
37 {
38 using pair = std::pair<Int, Int>;
39
40 Int n = dist.rows();
41
42 if (k>=n) k=n-1; //TODO this is ugly
43
44 Map<Set<Int>, pair> deg_idx; //first pair entry is degree in filtration, second is index in bd matrix
45
46 Int size = 0;
47 for (Int i = 1; i <= k+1; ++i)
48 size += Int(Integer::binom(n, i)); //number of all subsets of [n] with <=k members
49 Array<Cell> F(size); //TODO: IntType thingy
50
51 Int cell_index = 0; //index in filtration
52
53 Array<SparseMatrix<Coeff>> bd(k+1);
54 bd[0] = SparseMatrix<Coeff>(n,1);
55
56 for (Int index = 0; index < n; ++index) {
57 bd[0][index][0] = 1;
58 deg_idx[Set<Int>{index}] = pair(point_degs[index], index); //get point degrees from input array
59 Cell c(point_degs[index], 0, index);
60 F[cell_index++] = c;
61 }
62 Set<Int> facet = sequence(0, n); //the full-dimensional simplex
63
64 for (Int d = 1; d <= k; ++d) { //iterate over dimensions to max_dim
65 bd[d] = SparseMatrix<Coeff>(Int(Integer::binom(n,d+1)), Int(Integer::binom(n,d)));
66
67 Int index = 0; //index in bd matrix
68
69 for (auto i = entire(all_subsets_of_k(facet,d+1)); !i.at_end(); ++i) { //iterate over all simplices of dimension d
70
71 Int sgn = 1; //entry in bd matrix simply alternates as simplices are enumerated lexicographically
72 Int max = 0; //maximal degree of boundary simplices
73 Set<Int> simplex(*i);
74
75 for (auto s=entire(all_subsets_of_k(simplex,d)); !s.at_end(); ++s){ //iterate over boundary and find maximal degree
76
77 pair p = deg_idx[*s];
78
79 Int s_deg = p.first;
80 if (s_deg > max) max = s_deg;
81
82 bd[d][index][p.second] = sgn; // set bd matrix entry
83 sgn *= -1;
84 }
85 if (d == 1) { //for edges, calculate degree from distance matrix
86 const double l = double(dist(simplex.front(), simplex.back()));
87 assign_max(max, Int(ceil(l / step)));
88 }
89
90 deg_idx[simplex] = pair(max,index);
91 Cell c(max, d, index);
92 F[cell_index++] = c;
93 ++index;
94 }
95 }
96 return Filtration<SparseMatrix<Coeff> >(F,bd);
97 }
98
99
100 UserFunction4perl("# @category Producing a simplicial complex from other objects"
101 "# Computes the __Vietoris Rips complex__ of a point set. The set is passed as its so-called \"distance matrix\", whose (i,j)-entry is the distance between point i and j. This matrix can e.g. be computed using the distance_matrix function. The points corresponding to vertices of a common simplex will all have a distance less than //delta// from each other."
102 "# @param Matrix D the \"distance matrix\" of the point set (can be upper triangular)"
103 "# @param Rational delta"
104 "# @return SimplicialComplex"
105 "# @example The VR-complex from 3 points (vertices of a triangle with side lengths 3, 3, and 5) for different delta:"
106 "# > print vietoris_rips_complex(new Matrix([[0,3,3],[0,0,5],[0,0,0]]), 2)->FACETS;"
107 "# | {0}"
108 "# | {1}"
109 "# | {2}"
110 "# > print vietoris_rips_complex(new Matrix([[0,3,3],[0,0,5],[0,0,0]]), 4)->FACETS;"
111 "# | {0 1}"
112 "# | {0 2}"
113 "# > print vietoris_rips_complex(new Matrix([[0,3,3],[0,0,5],[0,0,0]]), 6)->FACETS;"
114 "# | {0 1 2}"
115 ,&vietoris_rips_complex, "vietoris_rips_complex($$)");
116
117
118 UserFunctionTemplate4perl("# @category Other"
119 "# Constructs the k-skeleton of the Vietrois Rips filtration of a point set. The set is passed as its so-called \"distance matrix\", whose (i,j)-entry is the distance between point i and j. This matrix can e.g. be computed using the distance_matrix function. The other inputs are an integer array containing the degree of each point, the desired distance step size between frames, and the dimension up to which to compute the skeleton. Redundant points will appear as separate vertices of the complex. Setting k to |S| will compute the entire VR-Complex for each frame."
120 "# @param Matrix D the \"distance matrix\" of the point set (can be upper triangular)"
121 "# @param Array<Int> deg the degrees of input points"
122 "# @param Float step_size "
123 "# @param Int k dimension of the resulting filtration"
124 "# @tparam Coeff desired coefficient type of the filtration"
125 "# @return Filtration<SparseMatrix<Coeff, NonSymmetric> >",
126 "vietoris_rips_filtration<Coeff>($$$$)");
127 } }
128
129 // Local Variables:
130 // mode:C++
131 // c-basic-offset:3
132 // indent-tabs-mode:nil
133 // End:
134