diff options
32 files changed, 2742 insertions, 1079 deletions
diff --git a/data/points/grid_10_10_10_in_0_1.weights b/data/points/grid_10_10_10_in_0_1.weights new file mode 100644 index 00000000..48926e09 --- /dev/null +++ b/data/points/grid_10_10_10_in_0_1.weights @@ -0,0 +1,1000 @@ +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 +1e-6 diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp index 67735ba1..9d1bc08c 100644 --- a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp +++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp @@ -20,7 +20,6 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ - #include <gudhi/reader_utils.h> #include <gudhi/Bitmap_cubical_complex.h> #include <gudhi/Persistent_cohomology.h> @@ -32,16 +31,18 @@ #include <cstddef> int main(int argc, char** argv) { - std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " << - "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a" << - "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let " << - "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are " << - "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See " << - "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n" << std::endl; + std::cout + << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " + << "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a" + << "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let " + << "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are " + << "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See " + << "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n" + << std::endl; if (argc != 2) { - std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " << - "the input. The program will now terminate.\n"; + std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " + << "the input. The program will now terminate.\n"; return 1; } @@ -54,7 +55,7 @@ int main(int argc, char** argv) { // Compute the persistence diagram of the complex Persistent_cohomology pcoh(b); - int p = 2; + int p = 11; double min_persistence = 0; pcoh.init_coefficients(p); // initializes the coefficient field for homology @@ -66,7 +67,7 @@ int main(int argc, char** argv) { std::size_t last_in_path = output_file_name.find_last_of("/\\"); if (last_in_path != std::string::npos) { - output_file_name = output_file_name.substr(last_in_path+1); + output_file_name = output_file_name.substr(last_in_path + 1); } std::ofstream out(output_file_name.c_str()); @@ -77,4 +78,3 @@ int main(int argc, char** argv) { return 0; } - diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp index 122160a2..c812cb3a 100644 --- a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp +++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp @@ -20,7 +20,6 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ - #include <gudhi/reader_utils.h> #include <gudhi/Bitmap_cubical_complex.h> #include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h> @@ -33,22 +32,24 @@ #include <string> int main(int argc, char** argv) { - std::cout << "This program computes persistent homology, by using " << - "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in " << - "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 " << - "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers " << - "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We " << - "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for" << - " example.\n" << std::endl; + std::cout + << "This program computes persistent homology, by using " + << "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in " + << "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 " + << "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers " + << "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We " + << "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for" + << " example.\n" + << std::endl; if (argc != 2) { - std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " << - "the input. The program will now terminate.\n"; + std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " + << "the input. The program will now terminate.\n"; return 1; } typedef Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> Bitmap_base; - typedef Gudhi::cubical_complex::Bitmap_cubical_complex< Bitmap_base > Bitmap_cubical_complex; + typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_base> Bitmap_cubical_complex; Bitmap_cubical_complex b(argv[1]); @@ -57,7 +58,7 @@ int main(int argc, char** argv) { // Compute the persistence diagram of the complex Persistent_cohomology pcoh(b, true); - int p = 2; + int p = 11; double min_persistence = 0; pcoh.init_coefficients(p); // initializes the coefficient field for homology pcoh.compute_persistent_cohomology(min_persistence); @@ -68,7 +69,7 @@ int main(int argc, char** argv) { std::size_t last_in_path = output_file_name.find_last_of("/\\"); if (last_in_path != std::string::npos) { - output_file_name = output_file_name.substr(last_in_path+1); + output_file_name = output_file_name.substr(last_in_path + 1); } std::ofstream out(output_file_name.c_str()); @@ -79,4 +80,3 @@ int main(int argc, char** argv) { return 0; } - diff --git a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp index 16ad65a0..f70558f2 100644 --- a/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp +++ b/src/Bitmap_cubical_complex/example/Random_bitmap_cubical_complex.cpp @@ -20,7 +20,6 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ - // for persistence algorithm #include <gudhi/reader_utils.h> #include <gudhi/Bitmap_cubical_complex.h> @@ -34,10 +33,11 @@ int main(int argc, char** argv) { srand(time(0)); - std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " << - "complexes. The first parameter of the program is the dimension D of the bitmap. The next D parameters are " << - "number of top dimensional cubes in each dimension of the bitmap. The program will create random cubical " << - "complex of that sizes and compute persistent homology of it." << std::endl; + std::cout + << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " + << "complexes. The first parameter of the program is the dimension D of the bitmap. The next D parameters are " + << "number of top dimensional cubes in each dimension of the bitmap. The program will create random cubical " + << "complex of that sizes and compute persistent homology of it." << std::endl; int p = 2; double min_persistence = 0; @@ -47,16 +47,16 @@ int main(int argc, char** argv) { return 1; } - size_t dimensionOfBitmap = (size_t) atoi(argv[1]); - std::vector< unsigned > sizes; + size_t dimensionOfBitmap = (size_t)atoi(argv[1]); + std::vector<unsigned> sizes; size_t multipliers = 1; for (size_t dim = 0; dim != dimensionOfBitmap; ++dim) { - unsigned sizeInThisDimension = (unsigned) atoi(argv[2 + dim]); + unsigned sizeInThisDimension = (unsigned)atoi(argv[2 + dim]); sizes.push_back(sizeInThisDimension); multipliers *= sizeInThisDimension; } - std::vector< double > data; + std::vector<double> data; for (size_t i = 0; i != multipliers; ++i) { data.push_back(rand() / static_cast<double>(RAND_MAX)); } @@ -80,4 +80,3 @@ int main(int argc, char** argv) { return 0; } - diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index f395de65..f82e8ce9 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -31,7 +31,7 @@ #endif #include <limits> -#include <utility> // for pair<> +#include <utility> // for pair<> #include <algorithm> // for sort #include <vector> #include <numeric> // for iota @@ -43,7 +43,8 @@ namespace cubical_complex { // global variable, was used just for debugging. const bool globalDbg = false; -template <typename T> class is_before_in_filtration; +template <typename T> +class is_before_in_filtration; /** * @brief Cubical complex represented as a bitmap. @@ -64,7 +65,6 @@ class Bitmap_cubical_complex : public T { typedef typename T::filtration_type Filtration_value; typedef Simplex_key Simplex_handle; - //*********************************************// // Constructors //*********************************************// @@ -77,8 +77,8 @@ class Bitmap_cubical_complex : public T { /** * Constructor form a Perseus-style file. **/ - Bitmap_cubical_complex(const char* perseus_style_file) : - T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) { + Bitmap_cubical_complex(const char* perseus_style_file) + : T(perseus_style_file), key_associated_to_simplex(this->total_number_of_cells + 1) { if (globalDbg) { std::cerr << "Bitmap_cubical_complex( const char* perseus_style_file )\n"; } @@ -97,9 +97,8 @@ class Bitmap_cubical_complex : public T { * with filtration on top dimensional cells. **/ Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, - const std::vector<Filtration_value>& top_dimensional_cells) : - T(dimensions, top_dimensional_cells), - key_associated_to_simplex(this->total_number_of_cells + 1) { + const std::vector<Filtration_value>& top_dimensional_cells) + : T(dimensions, top_dimensional_cells), key_associated_to_simplex(this->total_number_of_cells + 1) { for (size_t i = 0; i != this->total_number_of_cells; ++i) { this->key_associated_to_simplex[i] = i; } @@ -118,9 +117,9 @@ class Bitmap_cubical_complex : public T { **/ Bitmap_cubical_complex(const std::vector<unsigned>& dimensions, const std::vector<Filtration_value>& top_dimensional_cells, - std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) : - T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed), - key_associated_to_simplex(this->total_number_of_cells + 1) { + std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed) + : T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed), + key_associated_to_simplex(this->total_number_of_cells + 1) { for (size_t i = 0; i != this->total_number_of_cells; ++i) { this->key_associated_to_simplex[i] = i; } @@ -142,9 +141,7 @@ class Bitmap_cubical_complex : public T { /** * Returns number of all cubes in the complex. **/ - size_t num_simplices()const { - return this->total_number_of_cells; - } + size_t num_simplices() const { return this->total_number_of_cells; } /** * Returns a Simplex_handle to a cube that do not exist in this complex. @@ -159,14 +156,12 @@ class Bitmap_cubical_complex : public T { /** * Returns dimension of the complex. **/ - inline size_t dimension()const { - return this->sizes.size(); - } + inline size_t dimension() const { return this->sizes.size(); } /** * Return dimension of a cell pointed by the Simplex_handle. **/ - inline unsigned dimension(Simplex_handle sh)const { + inline unsigned dimension(Simplex_handle sh) const { if (globalDbg) { std::cerr << "unsigned dimension(const Simplex_handle& sh)\n"; } @@ -199,7 +194,7 @@ class Bitmap_cubical_complex : public T { /** * Return the key of a cube pointed by the Simplex_handle. **/ - Simplex_key key(Simplex_handle sh)const { + Simplex_key key(Simplex_handle sh) const { if (globalDbg) { std::cerr << "Simplex_key key(const Simplex_handle& sh)\n"; } @@ -217,7 +212,7 @@ class Bitmap_cubical_complex : public T { std::cerr << "Simplex_handle simplex(Simplex_key key)\n"; } if (key != null_key()) { - return this->simplex_associated_to_key[ key ]; + return this->simplex_associated_to_key[key]; } return null_simplex(); } @@ -246,8 +241,8 @@ class Bitmap_cubical_complex : public T { /** * Boundary_simplex_range class provides ranges for boundary iterators. **/ - typedef typename std::vector< Simplex_handle >::iterator Boundary_simplex_iterator; - typedef typename std::vector< Simplex_handle > Boundary_simplex_range; + typedef typename std::vector<Simplex_handle>::iterator Boundary_simplex_iterator; + typedef typename std::vector<Simplex_handle> Boundary_simplex_range; /** * Filtration_simplex_iterator class provides an iterator though the whole structure in the order of filtration. @@ -257,13 +252,13 @@ class Bitmap_cubical_complex : public T { **/ class Filtration_simplex_range; - class Filtration_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + class Filtration_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> { // Iterator over all simplices of the complex in the order of the indexing scheme. // 'value_type' must be 'Simplex_handle'. public: - Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) { } + Filtration_simplex_iterator(Bitmap_cubical_complex* b) : b(b), position(0) {} - Filtration_simplex_iterator() : b(NULL), position(0) { } + Filtration_simplex_iterator() : b(NULL), position(0) {} Filtration_simplex_iterator operator++() { if (globalDbg) { @@ -288,14 +283,14 @@ class Bitmap_cubical_complex : public T { return (*this); } - bool operator==(const Filtration_simplex_iterator& rhs)const { + bool operator==(const Filtration_simplex_iterator& rhs) const { if (globalDbg) { std::cerr << "bool operator == ( const Filtration_simplex_iterator& rhs )\n"; } - return ( this->position == rhs.position); + return (this->position == rhs.position); } - bool operator!=(const Filtration_simplex_iterator& rhs)const { + bool operator!=(const Filtration_simplex_iterator& rhs) const { if (globalDbg) { std::cerr << "bool operator != ( const Filtration_simplex_iterator& rhs )\n"; } @@ -306,7 +301,7 @@ class Bitmap_cubical_complex : public T { if (globalDbg) { std::cerr << "Simplex_handle operator*()\n"; } - return this->b->simplex_associated_to_key[ this->position ]; + return this->b->simplex_associated_to_key[this->position]; } friend class Filtration_simplex_range; @@ -326,7 +321,7 @@ class Bitmap_cubical_complex : public T { typedef Filtration_simplex_iterator const_iterator; typedef Filtration_simplex_iterator iterator; - Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) { } + Filtration_simplex_range(Bitmap_cubical_complex<T>* b) : b(b) {} Filtration_simplex_iterator begin() { if (globalDbg) { @@ -348,8 +343,6 @@ class Bitmap_cubical_complex : public T { Bitmap_cubical_complex<T>* b; }; - - //*********************************************// // Methods to access iterators from the container: @@ -357,9 +350,7 @@ class Bitmap_cubical_complex : public T { * boundary_simplex_range creates an object of a Boundary_simplex_range class * that provides ranges for the Boundary_simplex_iterator. **/ - Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { - return this->get_boundary_of_a_cell(sh); - } + Boundary_simplex_range boundary_simplex_range(Simplex_handle sh) { return this->get_boundary_of_a_cell(sh); } /** * filtration_simplex_range creates an object of a Filtration_simplex_range class @@ -374,8 +365,6 @@ class Bitmap_cubical_complex : public T { } //*********************************************// - - //*********************************************// // Elements which are in Gudhi now, but I (and in all the cases I asked also Marc) do not understand why they are // there. @@ -390,25 +379,25 @@ class Bitmap_cubical_complex : public T { * Function needed for compatibility with Gudhi. Not useful for other purposes. **/ std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) { - std::vector< size_t > bdry = this->get_boundary_of_a_cell(sh); + std::vector<size_t> bdry = this->get_boundary_of_a_cell(sh); if (globalDbg) { std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; std::cerr << "bdry.size() : " << bdry.size() << std::endl; } // this method returns two first elements from the boundary of sh. if (bdry.size() < 2) - throw("Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the " - "boundary."); + throw( + "Error in endpoints in Bitmap_cubical_complex class. The cell have less than two elements in the " + "boundary."); return std::make_pair(bdry[0], bdry[1]); } - /** * Class needed for compatibility with Gudhi. Not useful for other purposes. **/ class Skeleton_simplex_range; - class Skeleton_simplex_iterator : std::iterator< std::input_iterator_tag, Simplex_handle > { + class Skeleton_simplex_iterator : std::iterator<std::input_iterator_tag, Simplex_handle> { // Iterator over all simplices of the complex in the order of the indexing scheme. // 'value_type' must be 'Simplex_handle'. public: @@ -418,15 +407,13 @@ class Bitmap_cubical_complex : public T { } // find the position of the first simplex of a dimension d this->position = 0; - while ( - (this->position != b->data.size()) && - (this->b->get_dimension_of_a_cell(this->position) != this->dimension) - ) { + while ((this->position != b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) { ++this->position; } } - Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) { } + Skeleton_simplex_iterator() : b(NULL), position(0), dimension(0) {} Skeleton_simplex_iterator operator++() { if (globalDbg) { @@ -434,10 +421,8 @@ class Bitmap_cubical_complex : public T { } // increment the position as long as you did not get to the next element of the dimension dimension. ++this->position; - while ( - (this->position != this->b->data.size()) && - (this->b->get_dimension_of_a_cell(this->position) != this->dimension) - ) { + while ((this->position != this->b->data.size()) && + (this->b->get_dimension_of_a_cell(this->position) != this->dimension)) { ++this->position; } return (*this); @@ -459,14 +444,14 @@ class Bitmap_cubical_complex : public T { return (*this); } - bool operator==(const Skeleton_simplex_iterator& rhs)const { + bool operator==(const Skeleton_simplex_iterator& rhs) const { if (globalDbg) { std::cerr << "bool operator ==\n"; } - return ( this->position == rhs.position); + return (this->position == rhs.position); } - bool operator!=(const Skeleton_simplex_iterator& rhs)const { + bool operator!=(const Skeleton_simplex_iterator& rhs) const { if (globalDbg) { std::cerr << "bool operator != ( const Skeleton_simplex_iterator& rhs )\n"; } @@ -481,6 +466,7 @@ class Bitmap_cubical_complex : public T { } friend class Skeleton_simplex_range; + private: Bitmap_cubical_complex<T>* b; size_t position; @@ -497,7 +483,7 @@ class Bitmap_cubical_complex : public T { typedef Skeleton_simplex_iterator const_iterator; typedef Skeleton_simplex_iterator iterator; - Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) { } + Skeleton_simplex_range(Bitmap_cubical_complex<T>* b, unsigned dimension) : b(b), dimension(dimension) {} Skeleton_simplex_iterator begin() { if (globalDbg) { @@ -533,8 +519,8 @@ class Bitmap_cubical_complex : public T { friend class is_before_in_filtration<T>; protected: - std::vector< size_t > key_associated_to_simplex; - std::vector< size_t > simplex_associated_to_key; + std::vector<size_t> key_associated_to_simplex; + std::vector<size_t> simplex_associated_to_key; }; // Bitmap_cubical_complex template <typename T> @@ -552,16 +538,15 @@ void Bitmap_cubical_complex<T>::initialize_simplex_associated_to_key() { #endif // we still need to deal here with a key_associated_to_simplex: - for ( size_t i = 0 ; i != simplex_associated_to_key.size() ; ++i ) { - this->key_associated_to_simplex[ simplex_associated_to_key[i] ] = i; + for (size_t i = 0; i != simplex_associated_to_key.size(); ++i) { + this->key_associated_to_simplex[simplex_associated_to_key[i]] = i; } } template <typename T> class is_before_in_filtration { public: - explicit is_before_in_filtration(Bitmap_cubical_complex<T> * CC) - : CC_(CC) { } + explicit is_before_in_filtration(Bitmap_cubical_complex<T>* CC) : CC_(CC) {} bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1, const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const { diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h index 0442ac34..4adadce6 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h @@ -32,7 +32,8 @@ #include <algorithm> #include <iterator> #include <limits> -#include <utility> // for pair<> +#include <utility> +#include <stdexcept> namespace Gudhi { @@ -65,8 +66,7 @@ class Bitmap_cubical_complex_base { /** *Default constructor **/ - Bitmap_cubical_complex_base() : - total_number_of_cells(0) { } + Bitmap_cubical_complex_base() : total_number_of_cells(0) {} /** * There are a few constructors of a Bitmap_cubical_complex_base class. * First one, that takes vector<unsigned>, creates an empty bitmap of a dimension equal @@ -90,7 +90,7 @@ class Bitmap_cubical_complex_base { /** * Destructor of the Bitmap_cubical_complex_base class. **/ - virtual ~Bitmap_cubical_complex_base() { } + virtual ~Bitmap_cubical_complex_base() {} /** * The functions get_boundary_of_a_cell, get_coboundary_of_a_cell, get_dimension_of_a_cell @@ -100,8 +100,10 @@ class Bitmap_cubical_complex_base { * non-negative integer, indicating a position of a cube in the data structure. * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to * the positions of (co)boundary element of the input cell. + * The boundary elements are guaranteed to be returned so that the + * incidence coefficients of boundary elements are alternating. */ - virtual inline std::vector< size_t > get_boundary_of_a_cell(size_t cell)const; + virtual inline std::vector<size_t> get_boundary_of_a_cell(size_t cell) const; /** * The functions get_coboundary_of_a_cell, get_coboundary_of_a_cell, * get_dimension_of_a_cell and get_cell_data are the basic @@ -112,13 +114,74 @@ class Bitmap_cubical_complex_base { * In the case of functions that compute (co)boundary, the output is a vector if * non-negative integers pointing to the * positions of (co)boundary element of the input cell. + * Note that unlike in the case of boundary, over here the elements are + * not guaranteed to be returned with alternating incidence numbers. + * **/ - virtual inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell)const; + virtual inline std::vector<size_t> get_coboundary_of_a_cell(size_t cell) const; + /** - * In the case of get_dimension_of_a_cell function, the output is a non-negative integer - * indicating the dimension of a cell. - **/ - inline unsigned get_dimension_of_a_cell(size_t cell)const; + * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of + * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence + * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$. + * Note that first parameter is a cube of dimension n, + * and the second parameter is an adjusted cube in dimension n-1. + * Given \f$A = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [b_{j},e_{j}] \times [b_{j+1},e_{j+1}] \times \ldots + *\times [b_{n},e_{n}] \f$ + * such that \f$ b_{j} \neq e_{j} \f$ + * and \f$B = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [a,a] \times [b_{j+1},e_{j+1}] \times \ldots \times + *[b_{n},e_{n}] \f$ + * where \f$ a = b_{j}\f$ or \f$ a = e_{j}\f$, the incidence between \f$A\f$ and \f$B\f$ + * computed by this procedure is given by formula: + * \f$ c\ (-1)^{\sum_{i=1}^{j-1} dim [b_{i},e_{i}]} \f$ + * Where \f$ dim [b_{i},e_{i}] = 0 \f$ if \f$ b_{i}=e_{i} \f$ and 1 in other case. + * c is -1 if \f$ a = b_{j}\f$ and 1 if \f$ a = e_{j}\f$. + * @exception std::logic_error In case when the cube \f$B\f$ is not n-1 + * dimensional face of a cube \f$A\f$. + **/ + virtual int compute_incidence_between_cells(size_t coface, size_t face) const { + // first get the counters for coface and face: + std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface); + std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face); + + // coface_counter and face_counter should agree at all positions except from one: + int number_of_position_in_which_counters_do_not_agree = -1; + size_t number_of_full_faces_that_comes_before = 0; + for (size_t i = 0; i != coface_counter.size(); ++i) { + if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) { + ++number_of_full_faces_that_comes_before; + } + if (coface_counter[i] != face_counter[i]) { + if (number_of_position_in_which_counters_do_not_agree != -1) { + std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n"; + throw std::logic_error( + "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face."); + } + number_of_position_in_which_counters_do_not_agree = i; + } + } + + int incidence = 1; + if (number_of_full_faces_that_comes_before % 2) incidence = -1; + // if the face cell is on the right from coface cell: + if (coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 == + face_counter[number_of_position_in_which_counters_do_not_agree]) { + incidence *= -1; + } + + return incidence; + } + + /** +* In the case of get_dimension_of_a_cell function, the output is a non-negative integer +* indicating the dimension of a cell. +* Note that unlike in the case of boundary, over here the elements are +* not guaranteed to be returned with alternating incidence numbers. +* To compute incidence between cells use compute_incidence_between_cells +* procedure +**/ + inline unsigned get_dimension_of_a_cell(size_t cell) const; + /** * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position. * This allows reading and changing the value of filtration. Note that if the value of a filtration is changed, the @@ -127,7 +190,6 @@ class Bitmap_cubical_complex_base { **/ inline T& get_cell_data(size_t cell); - /** * Typical input used to construct a baseBitmap class is a filtration given at the top dimensional cells. * Then, there are a few ways one can pick the filtration of lower dimensional @@ -141,23 +203,19 @@ class Bitmap_cubical_complex_base { /** * Returns dimension of a complex. **/ - inline unsigned dimension()const { - return sizes.size(); - } + inline unsigned dimension() const { return sizes.size(); } /** * Returns number of all cubes in the data structure. **/ - inline unsigned size()const { - return this->data.size(); - } + inline unsigned size() const { return this->data.size(); } /** * Writing to stream operator. By using it we get the values T of cells in order in which they are stored in the * structure. This procedure is used for debugging purposes. **/ template <typename K> - friend std::ostream& operator<<(std::ostream & os, const Bitmap_cubical_complex_base<K>& b); + friend std::ostream& operator<<(std::ostream& os, const Bitmap_cubical_complex_base<K>& b); /** * Function that put the input data to bins. By putting data to bins we mean rounding them to a sequence of values @@ -184,7 +242,7 @@ class Bitmap_cubical_complex_base { /** * Functions to find min and max values of filtration. **/ - std::pair< T, T > min_max_filtration(); + std::pair<T, T> min_max_filtration(); // ITERATORS @@ -192,11 +250,9 @@ class Bitmap_cubical_complex_base { * @brief Iterator through all cells in the complex (in order they appear in the structure -- i.e. * in lexicographical order). **/ - class All_cells_iterator : std::iterator< std::input_iterator_tag, T > { + class All_cells_iterator : std::iterator<std::input_iterator_tag, T> { public: - All_cells_iterator() { - this->counter = 0; - } + All_cells_iterator() { this->counter = 0; } All_cells_iterator operator++() { // first find first element of the counter that can be increased: @@ -215,14 +271,12 @@ class Bitmap_cubical_complex_base { return *this; } - bool operator==(const All_cells_iterator& rhs)const { - if (this->counter != rhs.counter)return false; + bool operator==(const All_cells_iterator& rhs) const { + if (this->counter != rhs.counter) return false; return true; } - bool operator!=(const All_cells_iterator& rhs)const { - return !(*this == rhs); - } + bool operator!=(const All_cells_iterator& rhs) const { return !(*this == rhs); } /* * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as @@ -231,10 +285,9 @@ class Bitmap_cubical_complex_base { * boundary and coboundary and dimension * and in function get_cell_data to get a filtration of a cell. */ - size_t operator*() { - return this->counter; - } + size_t operator*() { return this->counter; } friend class Bitmap_cubical_complex_base; + protected: size_t counter; }; @@ -261,57 +314,47 @@ class Bitmap_cubical_complex_base { **/ class All_cells_range { public: - All_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + All_cells_range(Bitmap_cubical_complex_base* b) : b(b) {} - All_cells_iterator begin() { - return b->all_cells_iterator_begin(); - } + All_cells_iterator begin() { return b->all_cells_iterator_begin(); } + + All_cells_iterator end() { return b->all_cells_iterator_end(); } - All_cells_iterator end() { - return b->all_cells_iterator_end(); - } private: Bitmap_cubical_complex_base<T>* b; }; - All_cells_range all_cells_range() { - return All_cells_range(this); - } - + All_cells_range all_cells_range() { return All_cells_range(this); } /** * Boundary_range class provides ranges for boundary iterators. **/ - typedef typename std::vector< size_t >::const_iterator Boundary_iterator; - typedef typename std::vector< size_t > Boundary_range; + typedef typename std::vector<size_t>::const_iterator Boundary_iterator; + typedef typename std::vector<size_t> Boundary_range; /** * boundary_simplex_range creates an object of a Boundary_simplex_range class * that provides ranges for the Boundary_simplex_iterator. **/ - Boundary_range boundary_range(size_t sh) { - return this->get_boundary_of_a_cell(sh); - } + Boundary_range boundary_range(size_t sh) { return this->get_boundary_of_a_cell(sh); } /** * Coboundary_range class provides ranges for boundary iterators. **/ - typedef typename std::vector< size_t >::const_iterator Coboundary_iterator; - typedef typename std::vector< size_t > Coboundary_range; + typedef typename std::vector<size_t>::const_iterator Coboundary_iterator; + typedef typename std::vector<size_t> Coboundary_range; /** * boundary_simplex_range creates an object of a Boundary_simplex_range class * that provides ranges for the Boundary_simplex_iterator. **/ - Coboundary_range coboundary_range(size_t sh) { - return this->get_coboundary_of_a_cell(sh); - } + Coboundary_range coboundary_range(size_t sh) { return this->get_coboundary_of_a_cell(sh); } /** * @brief Iterator through top dimensional cells of the complex. The cells appear in order they are stored * in the structure (i.e. in lexicographical order) **/ - class Top_dimensional_cells_iterator : std::iterator< std::input_iterator_tag, T > { + class Top_dimensional_cells_iterator : std::iterator<std::input_iterator_tag, T> { public: Top_dimensional_cells_iterator(Bitmap_cubical_complex_base& b) : b(b) { this->counter = std::vector<size_t>(b.dimension()); @@ -321,7 +364,7 @@ class Bitmap_cubical_complex_base { Top_dimensional_cells_iterator operator++() { // first find first element of the counter that can be increased: size_t dim = 0; - while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1))++dim; + while ((dim != this->b.dimension()) && (this->counter[dim] == this->b.sizes[dim] - 1)) ++dim; if (dim != this->b.dimension()) { ++this->counter[dim]; @@ -346,18 +389,16 @@ class Bitmap_cubical_complex_base { return *this; } - bool operator==(const Top_dimensional_cells_iterator& rhs)const { - if (&this->b != &rhs.b)return false; - if (this->counter.size() != rhs.counter.size())return false; + bool operator==(const Top_dimensional_cells_iterator& rhs) const { + if (&this->b != &rhs.b) return false; + if (this->counter.size() != rhs.counter.size()) return false; for (size_t i = 0; i != this->counter.size(); ++i) { - if (this->counter[i] != rhs.counter[i])return false; + if (this->counter[i] != rhs.counter[i]) return false; } return true; } - bool operator!=(const Top_dimensional_cells_iterator& rhs)const { - return !(*this == rhs); - } + bool operator!=(const Top_dimensional_cells_iterator& rhs) const { return !(*this == rhs); } /* * The operator * returns position of a cube in the structure of cubical complex. This position can be then used as @@ -366,11 +407,9 @@ class Bitmap_cubical_complex_base { * boundary and coboundary and dimension * and in function get_cell_data to get a filtration of a cell. */ - size_t operator*() { - return this->compute_index_in_bitmap(); - } + size_t operator*() { return this->compute_index_in_bitmap(); } - size_t compute_index_in_bitmap()const { + size_t compute_index_in_bitmap() const { size_t index = 0; for (size_t i = 0; i != this->counter.size(); ++i) { index += (2 * this->counter[i] + 1) * this->b.multipliers[i]; @@ -378,14 +417,15 @@ class Bitmap_cubical_complex_base { return index; } - void print_counter()const { + void print_counter() const { for (size_t i = 0; i != this->counter.size(); ++i) { std::cout << this->counter[i] << " "; } } friend class Bitmap_cubical_complex_base; + protected: - std::vector< size_t > counter; + std::vector<size_t> counter; Bitmap_cubical_complex_base& b; }; @@ -414,32 +454,24 @@ class Bitmap_cubical_complex_base { **/ class Top_dimensional_cells_range { public: - Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) { } + Top_dimensional_cells_range(Bitmap_cubical_complex_base* b) : b(b) {} - Top_dimensional_cells_iterator begin() { - return b->top_dimensional_cells_iterator_begin(); - } + Top_dimensional_cells_iterator begin() { return b->top_dimensional_cells_iterator_begin(); } + + Top_dimensional_cells_iterator end() { return b->top_dimensional_cells_iterator_end(); } - Top_dimensional_cells_iterator end() { - return b->top_dimensional_cells_iterator_end(); - } private: Bitmap_cubical_complex_base<T>* b; }; - Top_dimensional_cells_range top_dimensional_cells_range() { - return Top_dimensional_cells_range(this); - } - + Top_dimensional_cells_range top_dimensional_cells_range() { return Top_dimensional_cells_range(this); } //****************************************************************************************************************// //****************************************************************************************************************// //****************************************************************************************************************// //****************************************************************************************************************// - inline size_t number_cells()const { - return this->total_number_of_cells; - } + inline size_t number_cells() const { return this->total_number_of_cells; } //****************************************************************************************************************// //****************************************************************************************************************// @@ -463,7 +495,7 @@ class Bitmap_cubical_complex_base { this->total_number_of_cells = multiplier; } - size_t compute_position_in_bitmap(const std::vector< unsigned >& counter) { + size_t compute_position_in_bitmap(const std::vector<unsigned>& counter) { size_t position = 0; for (size_t i = 0; i != this->multipliers.size(); ++i) { position += this->multipliers[i] * counter[i]; @@ -471,7 +503,7 @@ class Bitmap_cubical_complex_base { return position; } - std::vector<unsigned> compute_counter_for_given_cell(size_t cell)const { + std::vector<unsigned> compute_counter_for_given_cell(size_t cell) const { std::vector<unsigned> counter; counter.reserve(this->sizes.size()); for (size_t dim = this->sizes.size(); dim != 0; --dim) { @@ -486,8 +518,7 @@ class Bitmap_cubical_complex_base { const std::vector<T>& top_dimensional_cells); Bitmap_cubical_complex_base(const char* perseus_style_file, std::vector<bool> directions); Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes, std::vector<bool> directions); - Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, - const std::vector<T>& top_dimensional_cells, + Bitmap_cubical_complex_base(const std::vector<unsigned>& dimensions, const std::vector<T>& top_dimensional_cells, std::vector<bool> directions); }; @@ -495,8 +526,8 @@ template <typename T> void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) { bool bdg = false; - std::pair< T, T > min_max = this->min_max_filtration(); - T dx = (min_max.second - min_max.first) / (T) number_of_bins; + std::pair<T, T> min_max = this->min_max_filtration(); + T dx = (min_max.second - min_max.first) / (T)number_of_bins; // now put the data into the appropriate bins: for (size_t i = 0; i != this->data.size(); ++i) { @@ -514,7 +545,7 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(size_t number_of_bins) { template <typename T> void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { bool bdg = false; - std::pair< T, T > min_max = this->min_max_filtration(); + std::pair<T, T> min_max = this->min_max_filtration(); size_t number_of_bins = (min_max.second - min_max.first) / diameter_of_bin; // now put the data into the appropriate bins: @@ -531,33 +562,32 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { } template <typename T> -std::pair< T, T > Bitmap_cubical_complex_base<T>::min_max_filtration() { - std::pair< T, T > min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min()); +std::pair<T, T> Bitmap_cubical_complex_base<T>::min_max_filtration() { + std::pair<T, T> min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min()); for (size_t i = 0; i != this->data.size(); ++i) { - if (this->data[i] < min_max.first)min_max.first = this->data[i]; - if (this->data[i] > min_max.second)min_max.second = this->data[i]; + if (this->data[i] < min_max.first) min_max.first = this->data[i]; + if (this->data[i] > min_max.second) min_max.second = this->data[i]; } return min_max; } template <typename K> -std::ostream& operator<<(std::ostream & out, const Bitmap_cubical_complex_base<K>& b) { - for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator - it = b.all_cells_const_begin(); it != b.all_cells_const_end(); ++it) { +std::ostream& operator<<(std::ostream& out, const Bitmap_cubical_complex_base<K>& b) { + for (typename Bitmap_cubical_complex_base<K>::all_cells_const_iterator it = b.all_cells_const_begin(); + it != b.all_cells_const_end(); ++it) { out << *it << " "; } return out; } template <typename T> -Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base -(const std::vector<unsigned>& sizes) { +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes) { this->set_up_containers(sizes); } template <typename T> -void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list(const std::vector<unsigned>& sizes_in_following_directions, - const std::vector<T>& top_dimensional_cells) { +void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells_list( + const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { this->set_up_containers(sizes_in_following_directions); size_t number_of_top_dimensional_elements = 1; @@ -566,12 +596,13 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells } if (number_of_top_dimensional_elements != top_dimensional_cells.size()) { std::cerr << "Error in constructor Bitmap_cubical_complex_base ( std::vector<size_t> sizes_in_following_directions" - << ", std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " - << "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector." - << std::endl; - throw("Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions," - "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " - "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector."); + << ", std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " + << "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector." + << std::endl; + throw( + "Error in constructor Bitmap_cubical_complex_base( std::vector<size_t> sizes_in_following_directions," + "std::vector<T> top_dimensional_cells ). Number of top dimensional elements that follow from " + "sizes_in_following_directions vector is different than the size of top_dimensional_cells vector."); } Bitmap_cubical_complex_base<T>::Top_dimensional_cells_iterator it(*this); @@ -584,8 +615,8 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells } template <typename T> -Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base -(const std::vector<unsigned>& sizes_in_following_directions, const std::vector<T>& top_dimensional_cells) { +Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions, + const std::vector<T>& top_dimensional_cells) { this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells); } @@ -621,11 +652,9 @@ void Bitmap_cubical_complex_base<T>::read_perseus_style_file(const char* perseus T filtrationLevel; inFiltration >> filtrationLevel; if (dbg) { - std::cerr << "Cell of an index : " - << it.compute_index_in_bitmap() - << " and dimension: " - << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) - << " get the value : " << filtrationLevel << std::endl; + std::cerr << "Cell of an index : " << it.compute_index_in_bitmap() + << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) + << " get the value : " << filtrationLevel << std::endl; } this->get_cell_data(*it) = filtrationLevel; ++it; @@ -668,28 +697,36 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_ } template <typename T> -std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const { - std::vector< size_t > boundary_elements; +std::vector<size_t> Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell) const { + std::vector<size_t> boundary_elements; // Speed traded of for memory. Check if it is better in practice. - boundary_elements.reserve(this->dimension()*2); + boundary_elements.reserve(this->dimension() * 2); + size_t sum_of_dimensions = 0; size_t cell1 = cell; for (size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; if (position % 2 == 1) { - boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); - boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (sum_of_dimensions % 2) { + boundary_elements.push_back(cell + this->multipliers[i - 1]); + boundary_elements.push_back(cell - this->multipliers[i - 1]); + } else { + boundary_elements.push_back(cell - this->multipliers[i - 1]); + boundary_elements.push_back(cell + this->multipliers[i - 1]); + } + ++sum_of_dimensions; } cell1 = cell1 % this->multipliers[i - 1]; } + return boundary_elements; } template <typename T> -std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell)const { +std::vector<size_t> Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell) const { std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); - std::vector< size_t > coboundary_elements; + std::vector<size_t> coboundary_elements; size_t cell1 = cell; for (size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; @@ -697,8 +734,7 @@ std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(s if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) { coboundary_elements.push_back(cell - this->multipliers[i - 1]); } - if ( - (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) { + if ((cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) { coboundary_elements.push_back(cell + this->multipliers[i - 1]); } } @@ -708,7 +744,7 @@ std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(s } template <typename T> -unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell)const { +unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell) const { bool dbg = false; if (dbg) std::cerr << "\n\n\n Computing position o a cell of an index : " << cell << std::endl; unsigned dimension = 0; @@ -746,7 +782,7 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { size_t size_to_reserve = 1; for (size_t i = 0; i != this->multipliers.size(); ++i) { - size_to_reserve *= (size_t) ((this->multipliers[i] - 1) / 2); + size_to_reserve *= (size_t)((this->multipliers[i] - 1) / 2); } std::vector<size_t> indices_to_consider; @@ -771,22 +807,22 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { std::vector<size_t> bd = this->get_boundary_of_a_cell(indices_to_consider[i]); for (size_t boundaryIt = 0; boundaryIt != bd.size(); ++boundaryIt) { if (dbg) { - std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[ bd[boundaryIt] ] - << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[ indices_to_consider[i] ] - << std::endl; + std::cerr << "filtration of a cell : " << bd[boundaryIt] << " is : " << this->data[bd[boundaryIt]] + << " while of a cell: " << indices_to_consider[i] << " is: " << this->data[indices_to_consider[i]] + << std::endl; getchar(); } - if (this->data[ bd[boundaryIt] ] > this->data[ indices_to_consider[i] ]) { - this->data[ bd[boundaryIt] ] = this->data[ indices_to_consider[i] ]; + if (this->data[bd[boundaryIt]] > this->data[indices_to_consider[i]]) { + this->data[bd[boundaryIt]] = this->data[indices_to_consider[i]]; if (dbg) { - std::cerr << "Setting the value of a cell : " << bd[boundaryIt] << " to : " - << this->data[ indices_to_consider[i] ] << std::endl; + std::cerr << "Setting the value of a cell : " << bd[boundaryIt] + << " to : " << this->data[indices_to_consider[i]] << std::endl; getchar(); } } - if (is_this_cell_considered[ bd[boundaryIt] ] == false) { + if (is_this_cell_considered[bd[boundaryIt]] == false) { new_indices_to_consider.push_back(bd[boundaryIt]); - is_this_cell_considered[ bd[boundaryIt] ] = true; + is_this_cell_considered[bd[boundaryIt]] = true; } } } @@ -795,8 +831,8 @@ void Bitmap_cubical_complex_base<T>::impose_lower_star_filtration() { } template <typename T> -bool compareFirstElementsOfTuples(const std::pair< std::pair< T, size_t >, char >& first, - const std::pair< std::pair< T, size_t >, char >& second) { +bool compareFirstElementsOfTuples(const std::pair<std::pair<T, size_t>, char>& first, + const std::pair<std::pair<T, size_t>, char>& second) { if (first.first.first < second.first.first) { return true; } else { diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h index c3cc93dd..e2f86f3b 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -28,6 +28,7 @@ #include <cmath> #include <limits> // for numeric_limits<> #include <vector> +#include <stdexcept> namespace Gudhi { @@ -41,7 +42,8 @@ namespace cubical_complex { /** * @brief Cubical complex with periodic boundary conditions represented as a bitmap. * @ingroup cubical_complex - * @details This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the functions are + * @details This is a class implementing a bitmap data structure with periodic boundary conditions. Most of the + * functions are * identical to the functions from Bitmap_cubical_complex_base. * The ones that needed to be updated are the constructors and get_boundary_of_a_cell and get_coboundary_of_a_cell. */ @@ -53,7 +55,7 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c /** * Default constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class. */ - Bitmap_cubical_complex_periodic_boundary_conditions_base() { } + Bitmap_cubical_complex_periodic_boundary_conditions_base() {} /** * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the following * parameters: (1) vector with numbers of top dimensional cells in all dimensions and (2) vector of booleans. If @@ -61,8 +63,9 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c * imposed in this direction. In case of false, the periodic boundary conditions will not be imposed in the direction * i. */ - Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, - const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); + Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); /** * A constructor of Bitmap_cubical_complex_periodic_boundary_conditions_base class that takes the name of Perseus * style file as an input. Please consult the documentation about the specification of the file. @@ -75,9 +78,9 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c * value, that means that periodic boundary conditions are to be imposed in this direction. In case of false, the * periodic boundary conditions will not be imposed in the direction i. */ - Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, - const std::vector<T>& topDimensionalCells, - const std::vector< bool >& directions_in_which_periodic_b_cond_are_to_be_imposed); + Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); /** * Destructor of the Bitmap_cubical_complex_periodic_boundary_conditions_base class. @@ -88,17 +91,77 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c /** * A version of a function that return boundary of a given cell for an object of * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + * The boundary elements are guaranteed to be returned so that the + * incidence coefficients are alternating. */ - virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell) const; + virtual std::vector<size_t> get_boundary_of_a_cell(size_t cell) const; /** * A version of a function that return coboundary of a given cell for an object of * Bitmap_cubical_complex_periodic_boundary_conditions_base class. + * Note that unlike in the case of boundary, over here the elements are + * not guaranteed to be returned with alternating incidence numbers. + * To compute incidence between cells use compute_incidence_between_cells + * procedure */ - virtual std::vector< size_t > get_coboundary_of_a_cell(size_t cell) const; + virtual std::vector<size_t> get_coboundary_of_a_cell(size_t cell) const; + + /** + * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of + * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence + * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$. + * Note that first parameter is a cube of dimension n, + * and the second parameter is an adjusted cube in dimension n-1. + * Given \f$A = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [b_{j},e_{j}] \times [b_{j+1},e_{j+1}] \times \ldots + *\times [b_{n},e_{n}] \f$ + * such that \f$ b_{j} \neq e_{j} \f$ + * and \f$B = [b_1,e_1] \times \ldots \ [b_{j-1},e_{j-1}] \times [a,a] \times [b_{j+1},e_{j+1}] \times \ldots \times + *[b_{n},e_{n}]s \f$ + * where \f$ a = b_{j}\f$ or \f$ a = e_{j}\f$, the incidence between \f$A\f$ and \f$B\f$ + * computed by this procedure is given by formula: + * \f$ c\ (-1)^{\sum_{i=1}^{j-1} dim [b_{i},e_{i}]} \f$ + * Where \f$ dim [b_{i},e_{i}] = 0 \f$ if \f$ b_{i}=e_{i} \f$ and 1 in other case. + * c is -1 if \f$ a = b_{j}\f$ and 1 if \f$ a = e_{j}\f$. + * @exception std::logic_error In case when the cube \f$B\f$ is not n-1 + * dimensional face of a cube \f$A\f$. + **/ + virtual int compute_incidence_between_cells(size_t coface, size_t face) { + // first get the counters for coface and face: + std::vector<unsigned> coface_counter = this->compute_counter_for_given_cell(coface); + std::vector<unsigned> face_counter = this->compute_counter_for_given_cell(face); + + // coface_counter and face_counter should agree at all positions except from one: + int number_of_position_in_which_counters_do_not_agree = -1; + size_t number_of_full_faces_that_comes_before = 0; + for (size_t i = 0; i != coface_counter.size(); ++i) { + if ((coface_counter[i] % 2 == 1) && (number_of_position_in_which_counters_do_not_agree == -1)) { + ++number_of_full_faces_that_comes_before; + } + if (coface_counter[i] != face_counter[i]) { + if (number_of_position_in_which_counters_do_not_agree != -1) { + std::cout << "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face.\n"; + throw std::logic_error( + "Cells given to compute_incidence_between_cells procedure do not form a pair of coface-face."); + } + number_of_position_in_which_counters_do_not_agree = i; + } + } + + int incidence = 1; + if (number_of_full_faces_that_comes_before % 2) incidence = -1; + // if the face cell is on the right from coface cell: + if ((coface_counter[number_of_position_in_which_counters_do_not_agree] + 1 == + face_counter[number_of_position_in_which_counters_do_not_agree]) || + ((coface_counter[number_of_position_in_which_counters_do_not_agree] != 1) && + (face_counter[number_of_position_in_which_counters_do_not_agree] == 0))) { + incidence *= -1; + } + + return incidence; + } protected: - std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed; + std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed; void set_up_containers(const std::vector<unsigned>& sizes) { unsigned multiplier = 1; @@ -119,15 +182,19 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes); Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells); - void construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, - const std::vector<T>& topDimensionalCells, - const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); + + /** + * A procedure used to construct the data structures in the class. + **/ + void construct_complex_based_on_top_dimensional_cells( + const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed); }; template <typename T> -void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells(const std::vector<unsigned>& dimensions, - const std::vector<T>& topDimensionalCells, - const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { +void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_complex_based_on_top_dimensional_cells( + const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { this->directions_in_which_periodic_b_cond_are_to_be_imposed = directions_in_which_periodic_b_cond_are_to_be_imposed; this->set_up_containers(dimensions); @@ -140,14 +207,16 @@ void Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::construct_comp } template <typename T> -Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes, - const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& sizes, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { this->directions_in_which_periodic_b_cond_are_to_be_imposed(directions_in_which_periodic_b_cond_are_to_be_imposed); this->set_up_containers(sizes); } template <typename T> -Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const char* perseus_style_file) { +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( + const char* perseus_style_file) { // for Perseus style files: bool dbg = false; @@ -176,14 +245,12 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_comp while (!inFiltration.eof()) { double filtrationLevel; inFiltration >> filtrationLevel; - if (inFiltration.eof())break; + if (inFiltration.eof()) break; if (dbg) { - std::cerr << "Cell of an index : " - << it.compute_index_in_bitmap() - << " and dimension: " - << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) - << " get the value : " << filtrationLevel << std::endl; + std::cerr << "Cell of an index : " << it.compute_index_in_bitmap() + << " and dimension: " << this->get_dimension_of_a_cell(it.compute_index_in_bitmap()) + << " get the value : " << filtrationLevel << std::endl; } this->get_cell_data(*it) = filtrationLevel; ++it; @@ -193,24 +260,24 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_comp } template <typename T> -Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes) { +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& sizes) { this->directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(sizes.size(), false); this->set_up_containers(sizes); } template <typename T> -Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, - const std::vector<T>& topDimensionalCells) { +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells) { std::vector<bool> directions_in_which_periodic_b_cond_are_to_be_imposed = std::vector<bool>(dimensions.size(), false); this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, directions_in_which_periodic_b_cond_are_to_be_imposed); } template <typename T> -Bitmap_cubical_complex_periodic_boundary_conditions_base<T>:: -Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& dimensions, - const std::vector<T>& topDimensionalCells, - const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { +Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::Bitmap_cubical_complex_periodic_boundary_conditions_base( + const std::vector<unsigned>& dimensions, const std::vector<T>& topDimensionalCells, + const std::vector<bool>& directions_in_which_periodic_b_cond_are_to_be_imposed) { this->construct_complex_based_on_top_dimensional_cells(dimensions, topDimensionalCells, directions_in_which_periodic_b_cond_are_to_be_imposed); } @@ -218,46 +285,65 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsig // ***********************Methods************************ // template <typename T> -std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell(size_t cell) const { +std::vector<size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_boundary_of_a_cell( + size_t cell) const { bool dbg = false; if (dbg) { std::cerr << "Computations of boundary of a cell : " << cell << std::endl; } - std::vector< size_t > boundary_elements; + std::vector<size_t> boundary_elements; + boundary_elements.reserve(this->dimension() * 2); size_t cell1 = cell; + size_t sum_of_dimensions = 0; + for (size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction. - if (position % 2 == 1) { // if there are no periodic boundary conditions in this direction, we do not have to do anything. if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) { // std::cerr << "A\n"; - boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); - boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (sum_of_dimensions % 2) { + boundary_elements.push_back(cell - this->multipliers[i - 1]); + boundary_elements.push_back(cell + this->multipliers[i - 1]); + } else { + boundary_elements.push_back(cell + this->multipliers[i - 1]); + boundary_elements.push_back(cell - this->multipliers[i - 1]); + } if (dbg) { - std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + std::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " "; } } else { // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end. - if (position != 2 * this->sizes[ i - 1 ] - 1) { + if (position != 2 * this->sizes[i - 1] - 1) { // std::cerr << "B\n"; - boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); - boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); + if (sum_of_dimensions % 2) { + boundary_elements.push_back(cell - this->multipliers[i - 1]); + boundary_elements.push_back(cell + this->multipliers[i - 1]); + } else { + boundary_elements.push_back(cell + this->multipliers[i - 1]); + boundary_elements.push_back(cell - this->multipliers[i - 1]); + } if (dbg) { - std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " "; + std::cerr << cell - this->multipliers[i - 1] << " " << cell + this->multipliers[i - 1] << " "; } } else { // std::cerr << "C\n"; - boundary_elements.push_back(cell - this->multipliers[ i - 1 ]); - boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]); + if (sum_of_dimensions % 2) { + boundary_elements.push_back(cell - this->multipliers[i - 1]); + boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]); + } else { + boundary_elements.push_back(cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]); + boundary_elements.push_back(cell - this->multipliers[i - 1]); + } if (dbg) { - std::cerr << cell - this->multipliers[ i - 1 ] << " " << - cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ] << " "; + std::cerr << cell - this->multipliers[i - 1] << " " + << cell - (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1] << " "; } } } + ++sum_of_dimensions; } cell1 = cell1 % this->multipliers[i - 1]; } @@ -265,9 +351,10 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T } template <typename T> -std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell(size_t cell) const { +std::vector<size_t> Bitmap_cubical_complex_periodic_boundary_conditions_base<T>::get_coboundary_of_a_cell( + size_t cell) const { std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell); - std::vector< size_t > coboundary_elements; + std::vector<size_t> coboundary_elements; size_t cell1 = cell; for (size_t i = this->multipliers.size(); i != 0; --i) { unsigned position = cell1 / this->multipliers[i - 1]; @@ -289,7 +376,7 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T } else { // in this case counter[i-1] == 0. coboundary_elements.push_back(cell + this->multipliers[i - 1]); - coboundary_elements.push_back(cell + (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[i - 1]); + coboundary_elements.push_back(cell + (2 * this->sizes[i - 1] - 1) * this->multipliers[i - 1]); } } } diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp index db90eb94..4af699e9 100644 --- a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp +++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp @@ -33,17 +33,16 @@ #include <sstream> #include <vector> - typedef Gudhi::cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base; typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex; typedef Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> -Bitmap_cubical_complex_periodic_boundary_conditions_base; + Bitmap_cubical_complex_periodic_boundary_conditions_base; typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_periodic_boundary_conditions_base> -Bitmap_cubical_complex_periodic_boundary_conditions; + Bitmap_cubical_complex_periodic_boundary_conditions; BOOST_AUTO_TEST_CASE(check_dimension) { - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); std::vector<unsigned> dimensions({3, 3}); @@ -52,29 +51,28 @@ BOOST_AUTO_TEST_CASE(check_dimension) { } BOOST_AUTO_TEST_CASE(topDimensionalCellsIterator_test) { - std::vector< double > expectedFiltrationValues1({0, 0, 0, 0, 100, 0, 0, 0, 0}); + std::vector<double> expectedFiltrationValues1({0, 0, 0, 0, 100, 0, 0, 0, 0}); - std::vector< double > expectedFiltrationValues2({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> expectedFiltrationValues2({1, 2, 3, 4, 5, 6, 7, 8, 9}); - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); - std::vector< double > oneDimensionalCycle({0, 0, 0, 0, 100, 0, 0, 0, 0}); + std::vector<double> oneDimensionalCycle({0, 0, 0, 0, 100, 0, 0, 0, 0}); std::vector<unsigned> dimensions({3, 3}); Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); Bitmap_cubical_complex hole(dimensions, oneDimensionalCycle); - int i = 0; - for (Bitmap_cubical_complex::Top_dimensional_cells_iterator - it = increasing.top_dimensional_cells_iterator_begin(); it != increasing.top_dimensional_cells_iterator_end(); ++it) { + for (Bitmap_cubical_complex::Top_dimensional_cells_iterator it = increasing.top_dimensional_cells_iterator_begin(); + it != increasing.top_dimensional_cells_iterator_end(); ++it) { BOOST_CHECK(increasing.get_cell_data(*it) == expectedFiltrationValues2[i]); ++i; } i = 0; - for (Bitmap_cubical_complex::Top_dimensional_cells_iterator - it = hole.top_dimensional_cells_iterator_begin(); it != hole.top_dimensional_cells_iterator_end(); ++it) { + for (Bitmap_cubical_complex::Top_dimensional_cells_iterator it = hole.top_dimensional_cells_iterator_begin(); + it != hole.top_dimensional_cells_iterator_end(); ++it) { BOOST_CHECK(hole.get_cell_data(*it) == expectedFiltrationValues1[i]); ++i; } @@ -100,24 +98,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { std::vector<double> boundary8; boundary8.push_back(1); boundary8.push_back(15); - boundary8.push_back(7); boundary8.push_back(9); + boundary8.push_back(7); std::vector<double> boundary9; boundary9.push_back(2); boundary9.push_back(16); std::vector<double> boundary10; boundary10.push_back(3); boundary10.push_back(17); - boundary10.push_back(9); boundary10.push_back(11); + boundary10.push_back(9); std::vector<double> boundary11; boundary11.push_back(4); boundary11.push_back(18); std::vector<double> boundary12; boundary12.push_back(5); boundary12.push_back(19); - boundary12.push_back(11); boundary12.push_back(13); + boundary12.push_back(11); std::vector<double> boundary13; boundary13.push_back(6); boundary13.push_back(20); @@ -140,24 +138,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { std::vector<double> boundary22; boundary22.push_back(15); boundary22.push_back(29); - boundary22.push_back(21); boundary22.push_back(23); + boundary22.push_back(21); std::vector<double> boundary23; boundary23.push_back(16); boundary23.push_back(30); std::vector<double> boundary24; boundary24.push_back(17); boundary24.push_back(31); - boundary24.push_back(23); boundary24.push_back(25); + boundary24.push_back(23); std::vector<double> boundary25; boundary25.push_back(18); boundary25.push_back(32); std::vector<double> boundary26; boundary26.push_back(19); boundary26.push_back(33); - boundary26.push_back(25); boundary26.push_back(27); + boundary26.push_back(25); std::vector<double> boundary27; boundary27.push_back(20); boundary27.push_back(34); @@ -180,24 +178,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { std::vector<double> boundary36; boundary36.push_back(29); boundary36.push_back(43); - boundary36.push_back(35); boundary36.push_back(37); + boundary36.push_back(35); std::vector<double> boundary37; boundary37.push_back(30); boundary37.push_back(44); std::vector<double> boundary38; boundary38.push_back(31); boundary38.push_back(45); - boundary38.push_back(37); boundary38.push_back(39); + boundary38.push_back(37); std::vector<double> boundary39; boundary39.push_back(32); boundary39.push_back(46); std::vector<double> boundary40; boundary40.push_back(33); boundary40.push_back(47); - boundary40.push_back(39); boundary40.push_back(41); + boundary40.push_back(39); std::vector<double> boundary41; boundary41.push_back(34); boundary41.push_back(48); @@ -214,7 +212,7 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { boundary47.push_back(46); boundary47.push_back(48); std::vector<double> boundary48; - std::vector< std::vector<double> > boundaries; + std::vector<std::vector<double> > boundaries; boundaries.push_back(boundary0); boundaries.push_back(boundary1); boundaries.push_back(boundary2); @@ -265,15 +263,13 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { boundaries.push_back(boundary47); boundaries.push_back(boundary48); - - - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); std::vector<unsigned> dimensions({3, 3}); Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); for (size_t i = 0; i != increasing.size(); ++i) { - std::vector< size_t > bd = increasing.get_boundary_of_a_cell(i); + std::vector<size_t> bd = increasing.get_boundary_of_a_cell(i); for (size_t j = 0; j != bd.size(); ++j) { BOOST_CHECK(boundaries[i][j] == bd[j]); } @@ -281,13 +277,12 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) { } BOOST_AUTO_TEST_CASE(compute_boundary_test_2) { - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); std::vector<unsigned> dimensions({3, 3}); Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); - std::vector<double> coboundaryElements; coboundaryElements.push_back(7); coboundaryElements.push_back(1); @@ -373,9 +368,10 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_2) { coboundaryElements.push_back(40); coboundaryElements.push_back(41); coboundaryElements.push_back(47); + size_t number = 0; for (size_t i = 0; i != increasing.size(); ++i) { - std::vector< size_t > bd = increasing.get_coboundary_of_a_cell(i); + std::vector<size_t> bd = increasing.get_coboundary_of_a_cell(i); for (size_t j = 0; j != bd.size(); ++j) { BOOST_CHECK(coboundaryElements[number] == bd[j]); ++number; @@ -384,7 +380,7 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_2) { } BOOST_AUTO_TEST_CASE(compute_boundary_test_3) { - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); std::vector<unsigned> dimensions({3, 3}); @@ -447,13 +443,13 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_3) { } BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) { - std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); + std::vector<double> increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9}); std::vector<unsigned> dimensions({3, 3}); Bitmap_cubical_complex increasing(dimensions, increasingFiltrationOfTopDimensionalCells); - std::vector< unsigned > dim; + std::vector<unsigned> dim; dim.push_back(0); dim.push_back(0); dim.push_back(0); @@ -555,7 +551,6 @@ BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) { fil.push_back(9); fil.push_back(9); - Bitmap_cubical_complex::Filtration_simplex_range range = increasing.filtration_simplex_range(); size_t position = 0; for (Bitmap_cubical_complex::Filtration_simplex_iterator it = range.begin(); it != range.end(); ++it) { @@ -566,7 +561,7 @@ BOOST_AUTO_TEST_CASE(Filtration_simplex_iterator_test) { } BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { - std::vector< double > filtration({0, 0, 0, 0}); + std::vector<double> filtration({0, 0, 0, 0}); std::vector<unsigned> dimensions({2, 2}); @@ -575,57 +570,56 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions); BOOST_CHECK(cmplx.dimension() == 2); - std::vector<double> boundary0; std::vector<double> boundary1; - boundary1.push_back(0); boundary1.push_back(2); + boundary1.push_back(0); std::vector<double> boundary2; std::vector<double> boundary3; - boundary3.push_back(2); boundary3.push_back(0); + boundary3.push_back(2); std::vector<double> boundary4; - boundary4.push_back(0); boundary4.push_back(8); + boundary4.push_back(0); std::vector<double> boundary5; - boundary5.push_back(1); boundary5.push_back(9); + boundary5.push_back(1); boundary5.push_back(4); boundary5.push_back(6); std::vector<double> boundary6; - boundary6.push_back(2); boundary6.push_back(10); + boundary6.push_back(2); std::vector<double> boundary7; - boundary7.push_back(3); boundary7.push_back(11); + boundary7.push_back(3); boundary7.push_back(6); boundary7.push_back(4); std::vector<double> boundary8; std::vector<double> boundary9; - boundary9.push_back(8); boundary9.push_back(10); + boundary9.push_back(8); std::vector<double> boundary10; std::vector<double> boundary11; - boundary11.push_back(10); boundary11.push_back(8); + boundary11.push_back(10); std::vector<double> boundary12; - boundary12.push_back(8); boundary12.push_back(0); + boundary12.push_back(8); std::vector<double> boundary13; - boundary13.push_back(9); boundary13.push_back(1); + boundary13.push_back(9); boundary13.push_back(12); boundary13.push_back(14); std::vector<double> boundary14; - boundary14.push_back(10); boundary14.push_back(2); + boundary14.push_back(10); std::vector<double> boundary15; - boundary15.push_back(11); boundary15.push_back(3); + boundary15.push_back(11); boundary15.push_back(14); boundary15.push_back(12); - std::vector< std::vector<double> > boundaries; + std::vector<std::vector<double> > boundaries; boundaries.push_back(boundary0); boundaries.push_back(boundary1); boundaries.push_back(boundary2); @@ -644,7 +638,7 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { boundaries.push_back(boundary15); for (size_t i = 0; i != cmplx.size(); ++i) { - std::vector< size_t > bd = cmplx.get_boundary_of_a_cell(i); + std::vector<size_t> bd = cmplx.get_boundary_of_a_cell(i); for (size_t j = 0; j != bd.size(); ++j) { BOOST_CHECK(boundaries[i][j] == bd[j]); } @@ -652,7 +646,7 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) { } BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { - std::vector< double > filtration({0, 0, 0, 0}); + std::vector<double> filtration({0, 0, 0, 0}); std::vector<unsigned> dimensions({2, 2}); @@ -661,7 +655,6 @@ BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtration, periodic_directions); BOOST_CHECK(cmplx.dimension() == 2); - std::vector<double> coboundary0; coboundary0.push_back(4); coboundary0.push_back(12); @@ -711,7 +704,7 @@ BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { coboundary14.push_back(15); std::vector<double> coboundary15; - std::vector< std::vector<double> > coboundaries; + std::vector<std::vector<double> > coboundaries; coboundaries.push_back(coboundary0); coboundaries.push_back(coboundary1); coboundaries.push_back(coboundary2); @@ -730,7 +723,7 @@ BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { coboundaries.push_back(coboundary15); for (size_t i = 0; i != cmplx.size(); ++i) { - std::vector< size_t > cbd = cmplx.get_coboundary_of_a_cell(i); + std::vector<size_t> cbd = cmplx.get_coboundary_of_a_cell(i); for (size_t j = 0; j != cbd.size(); ++j) { BOOST_CHECK(coboundaries[i][j] == cbd[j]); } @@ -738,7 +731,7 @@ BOOST_AUTO_TEST_CASE(coboudary_operator_2d_bitmap_with_periodic_bcond) { } BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) { - std::vector< double > filtrationOrg({0, 1, 2, 3}); + std::vector<double> filtrationOrg({0, 1, 2, 3}); std::vector<unsigned> dimensions({2, 2}); @@ -747,7 +740,6 @@ BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) { Bitmap_cubical_complex_periodic_boundary_conditions cmplx(dimensions, filtrationOrg, periodic_directions); BOOST_CHECK(cmplx.dimension() == 2); - std::vector<double> filtration; filtration.push_back(0); // 0 filtration.push_back(0); // 1 @@ -766,613 +758,821 @@ BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) { filtration.push_back(2); // 14 filtration.push_back(3); // 15 - for (size_t i = 0; i != cmplx.size(); ++i) { BOOST_CHECK(filtration[i] == cmplx.get_cell_data(i)); } } -BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check) -{ - std::vector< double > expected_filtration; - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - - std::vector<unsigned> expected_dimension; - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - - std::vector<size_t> expected_boundary; - expected_boundary.push_back(0); - expected_boundary.push_back(2); - expected_boundary.push_back(2); - expected_boundary.push_back(4); - expected_boundary.push_back(0); - expected_boundary.push_back(10); - expected_boundary.push_back(1); - expected_boundary.push_back(11); - expected_boundary.push_back(5); - expected_boundary.push_back(7); - expected_boundary.push_back(2); - expected_boundary.push_back(12); - expected_boundary.push_back(3); - expected_boundary.push_back(13); - expected_boundary.push_back(7); - expected_boundary.push_back(9); - expected_boundary.push_back(4); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(12); - expected_boundary.push_back(12); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(20); - expected_boundary.push_back(11); - expected_boundary.push_back(21); - expected_boundary.push_back(15); - expected_boundary.push_back(17); - expected_boundary.push_back(12); - expected_boundary.push_back(22); - expected_boundary.push_back(13); - expected_boundary.push_back(23); - expected_boundary.push_back(17); - expected_boundary.push_back(19); - expected_boundary.push_back(14); - expected_boundary.push_back(24); - expected_boundary.push_back(20); - expected_boundary.push_back(22); - expected_boundary.push_back(22); - expected_boundary.push_back(24); - - - std::vector<size_t> expected_coboundary; - expected_coboundary.push_back(5); - expected_coboundary.push_back(1); - expected_coboundary.push_back(6); - expected_coboundary.push_back(7); - expected_coboundary.push_back(1); - expected_coboundary.push_back(3); - expected_coboundary.push_back(8); - expected_coboundary.push_back(9); - expected_coboundary.push_back(3); - expected_coboundary.push_back(6); - expected_coboundary.push_back(6); - expected_coboundary.push_back(8); - expected_coboundary.push_back(8); - expected_coboundary.push_back(5); - expected_coboundary.push_back(15); - expected_coboundary.push_back(11); - expected_coboundary.push_back(6); - expected_coboundary.push_back(16); - expected_coboundary.push_back(7); - expected_coboundary.push_back(17); - expected_coboundary.push_back(11); - expected_coboundary.push_back(13); - expected_coboundary.push_back(8); - expected_coboundary.push_back(18); - expected_coboundary.push_back(9); - expected_coboundary.push_back(19); - expected_coboundary.push_back(13); - expected_coboundary.push_back(16); - expected_coboundary.push_back(16); - expected_coboundary.push_back(18); - expected_coboundary.push_back(18); - expected_coboundary.push_back(15); - expected_coboundary.push_back(21); - expected_coboundary.push_back(16); - expected_coboundary.push_back(17); - expected_coboundary.push_back(21); - expected_coboundary.push_back(23); - expected_coboundary.push_back(18); - expected_coboundary.push_back(19); - expected_coboundary.push_back(23); - - - - std::vector< unsigned > sizes(2); - sizes[0] = 2; - sizes[1] = 2; - - std::vector< double > data(4); - data[0] = 0; - data[1] = 1; - data[2] = 2; - data[3] = 3; - - Bitmap_cubical_complex_base ba( sizes , data ); - int i = 0; - int bd_it = 0; - int cbd_it = 0; - for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it ) - { - BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); - BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); - - Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); - for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) - { - BOOST_CHECK( expected_boundary[bd_it] == *bd ); - ++bd_it; - } - - Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); - for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) - { - BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); - ++cbd_it; - } - ++i; + +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check) { + std::vector<double> expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(7); + expected_boundary.push_back(5); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(9); + expected_boundary.push_back(7); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(17); + expected_boundary.push_back(15); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(19); + expected_boundary.push_back(17); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + std::vector<unsigned> sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector<double> data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba(sizes, data); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + for (Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + BOOST_CHECK(expected_filtration[i] == ba.get_cell_data(*it)); + BOOST_CHECK(expected_dimension[i] == ba.get_dimension_of_a_cell(*it)); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + BOOST_CHECK(expected_boundary[bd_it] == *bd); + ++bd_it; } + + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for (Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin(); cbd != cbdrange.end(); ++cbd) { + BOOST_CHECK(expected_coboundary[cbd_it] == *cbd); + ++cbd_it; + } + ++i; + } } +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check_2) { + std::vector<double> expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(7); + expected_boundary.push_back(5); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(9); + expected_boundary.push_back(7); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(17); + expected_boundary.push_back(15); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(19); + expected_boundary.push_back(17); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + std::vector<unsigned> sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector<double> data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba(sizes, data); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + + Bitmap_cubical_complex_base::All_cells_range range(&ba); + for (Bitmap_cubical_complex_base::All_cells_iterator it = range.begin(); it != range.end(); ++it) { + BOOST_CHECK(expected_filtration[i] == ba.get_cell_data(*it)); + BOOST_CHECK(expected_dimension[i] == ba.get_dimension_of_a_cell(*it)); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + BOOST_CHECK(expected_boundary[bd_it] == *bd); + ++bd_it; + } + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for (Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin(); cbd != cbdrange.end(); ++cbd) { + BOOST_CHECK(expected_coboundary[cbd_it] == *cbd); + ++cbd_it; + } + ++i; + } +} +BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check) { + std::vector<double> expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> expected_dimension; + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(2); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + expected_dimension.push_back(1); + expected_dimension.push_back(0); + + std::vector<size_t> expected_boundary; + expected_boundary.push_back(0); + expected_boundary.push_back(2); + expected_boundary.push_back(2); + expected_boundary.push_back(4); + expected_boundary.push_back(0); + expected_boundary.push_back(10); + expected_boundary.push_back(1); + expected_boundary.push_back(11); + expected_boundary.push_back(7); + expected_boundary.push_back(5); + expected_boundary.push_back(2); + expected_boundary.push_back(12); + expected_boundary.push_back(3); + expected_boundary.push_back(13); + expected_boundary.push_back(9); + expected_boundary.push_back(7); + expected_boundary.push_back(4); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(12); + expected_boundary.push_back(12); + expected_boundary.push_back(14); + expected_boundary.push_back(10); + expected_boundary.push_back(20); + expected_boundary.push_back(11); + expected_boundary.push_back(21); + expected_boundary.push_back(17); + expected_boundary.push_back(15); + expected_boundary.push_back(12); + expected_boundary.push_back(22); + expected_boundary.push_back(13); + expected_boundary.push_back(23); + expected_boundary.push_back(19); + expected_boundary.push_back(17); + expected_boundary.push_back(14); + expected_boundary.push_back(24); + expected_boundary.push_back(20); + expected_boundary.push_back(22); + expected_boundary.push_back(22); + expected_boundary.push_back(24); + + std::vector<size_t> expected_coboundary; + expected_coboundary.push_back(5); + expected_coboundary.push_back(1); + expected_coboundary.push_back(6); + expected_coboundary.push_back(7); + expected_coboundary.push_back(1); + expected_coboundary.push_back(3); + expected_coboundary.push_back(8); + expected_coboundary.push_back(9); + expected_coboundary.push_back(3); + expected_coboundary.push_back(6); + expected_coboundary.push_back(6); + expected_coboundary.push_back(8); + expected_coboundary.push_back(8); + expected_coboundary.push_back(5); + expected_coboundary.push_back(15); + expected_coboundary.push_back(11); + expected_coboundary.push_back(6); + expected_coboundary.push_back(16); + expected_coboundary.push_back(7); + expected_coboundary.push_back(17); + expected_coboundary.push_back(11); + expected_coboundary.push_back(13); + expected_coboundary.push_back(8); + expected_coboundary.push_back(18); + expected_coboundary.push_back(9); + expected_coboundary.push_back(19); + expected_coboundary.push_back(13); + expected_coboundary.push_back(16); + expected_coboundary.push_back(16); + expected_coboundary.push_back(18); + expected_coboundary.push_back(18); + expected_coboundary.push_back(15); + expected_coboundary.push_back(21); + expected_coboundary.push_back(16); + expected_coboundary.push_back(17); + expected_coboundary.push_back(21); + expected_coboundary.push_back(23); + expected_coboundary.push_back(18); + expected_coboundary.push_back(19); + expected_coboundary.push_back(23); + + std::vector<unsigned> sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector<double> data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba(sizes, data); + int i = 0; + int bd_it = 0; + int cbd_it = 0; + + Bitmap_cubical_complex_base::All_cells_range range = ba.all_cells_range(); + for (Bitmap_cubical_complex_base::All_cells_iterator it = range.begin(); it != range.end(); ++it) { + BOOST_CHECK(expected_filtration[i] == ba.get_cell_data(*it)); + BOOST_CHECK(expected_dimension[i] == ba.get_dimension_of_a_cell(*it)); + + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + BOOST_CHECK(expected_boundary[bd_it] == *bd); + ++bd_it; + } + Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); + for (Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin(); cbd != cbdrange.end(); ++cbd) { + BOOST_CHECK(expected_coboundary[cbd_it] == *cbd); + ++cbd_it; + } + ++i; + } +} +BOOST_AUTO_TEST_CASE(Top_dimensional_cells_iterator_range_check) { + std::vector<double> expected_filtration; + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(0); + expected_filtration.push_back(1); + expected_filtration.push_back(1); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(2); + expected_filtration.push_back(3); + expected_filtration.push_back(3); + + std::vector<unsigned> sizes(2); + sizes[0] = 2; + sizes[1] = 2; + + std::vector<double> data(4); + data[0] = 0; + data[1] = 1; + data[2] = 2; + data[3] = 3; + + Bitmap_cubical_complex_base ba(sizes, data); + int i = 0; + Bitmap_cubical_complex_base::Top_dimensional_cells_range range = ba.top_dimensional_cells_range(); + for (Bitmap_cubical_complex_base::Top_dimensional_cells_iterator it = range.begin(); it != range.end(); ++it) { + BOOST_CHECK(data[i] == ba.get_cell_data(*it)); + BOOST_CHECK(ba.get_dimension_of_a_cell(*it) == 2); + ++i; + } +} +BOOST_AUTO_TEST_CASE(check_if_boundary_of_boundary_is_zero_non_periodic_case_3_d) { + std::vector<unsigned> sizes(3); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + + std::vector<double> data(27, 0); + + int number_of_all_elements = (2 * sizes[0] + 1) * (2 * sizes[1] + 1) * (2 * sizes[2] + 1); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + Bitmap_cubical_complex_base ba(sizes, data); + for (Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + int i = 1; + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd); + int j = 1; + for (Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin(); bd2 != second_bdrange.end(); + ++bd2) { + elems_in_boundary[*bd2] += i * j; + j *= -1; + } + i *= -1; + } + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } +} -BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check_2) -{ - std::vector< double > expected_filtration; - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - - std::vector<unsigned> expected_dimension; - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - - std::vector<size_t> expected_boundary; - expected_boundary.push_back(0); - expected_boundary.push_back(2); - expected_boundary.push_back(2); - expected_boundary.push_back(4); - expected_boundary.push_back(0); - expected_boundary.push_back(10); - expected_boundary.push_back(1); - expected_boundary.push_back(11); - expected_boundary.push_back(5); - expected_boundary.push_back(7); - expected_boundary.push_back(2); - expected_boundary.push_back(12); - expected_boundary.push_back(3); - expected_boundary.push_back(13); - expected_boundary.push_back(7); - expected_boundary.push_back(9); - expected_boundary.push_back(4); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(12); - expected_boundary.push_back(12); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(20); - expected_boundary.push_back(11); - expected_boundary.push_back(21); - expected_boundary.push_back(15); - expected_boundary.push_back(17); - expected_boundary.push_back(12); - expected_boundary.push_back(22); - expected_boundary.push_back(13); - expected_boundary.push_back(23); - expected_boundary.push_back(17); - expected_boundary.push_back(19); - expected_boundary.push_back(14); - expected_boundary.push_back(24); - expected_boundary.push_back(20); - expected_boundary.push_back(22); - expected_boundary.push_back(22); - expected_boundary.push_back(24); - - - std::vector<size_t> expected_coboundary; - expected_coboundary.push_back(5); - expected_coboundary.push_back(1); - expected_coboundary.push_back(6); - expected_coboundary.push_back(7); - expected_coboundary.push_back(1); - expected_coboundary.push_back(3); - expected_coboundary.push_back(8); - expected_coboundary.push_back(9); - expected_coboundary.push_back(3); - expected_coboundary.push_back(6); - expected_coboundary.push_back(6); - expected_coboundary.push_back(8); - expected_coboundary.push_back(8); - expected_coboundary.push_back(5); - expected_coboundary.push_back(15); - expected_coboundary.push_back(11); - expected_coboundary.push_back(6); - expected_coboundary.push_back(16); - expected_coboundary.push_back(7); - expected_coboundary.push_back(17); - expected_coboundary.push_back(11); - expected_coboundary.push_back(13); - expected_coboundary.push_back(8); - expected_coboundary.push_back(18); - expected_coboundary.push_back(9); - expected_coboundary.push_back(19); - expected_coboundary.push_back(13); - expected_coboundary.push_back(16); - expected_coboundary.push_back(16); - expected_coboundary.push_back(18); - expected_coboundary.push_back(18); - expected_coboundary.push_back(15); - expected_coboundary.push_back(21); - expected_coboundary.push_back(16); - expected_coboundary.push_back(17); - expected_coboundary.push_back(21); - expected_coboundary.push_back(23); - expected_coboundary.push_back(18); - expected_coboundary.push_back(19); - expected_coboundary.push_back(23); - - - - std::vector< unsigned > sizes(2); - sizes[0] = 2; - sizes[1] = 2; - - std::vector< double > data(4); - data[0] = 0; - data[1] = 1; - data[2] = 2; - data[3] = 3; - - Bitmap_cubical_complex_base ba( sizes , data ); - int i = 0; - int bd_it = 0; - int cbd_it = 0; - - Bitmap_cubical_complex_base::All_cells_range range(&ba); - for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it ) - { - BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); - BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); - - Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); - for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) - { - BOOST_CHECK( expected_boundary[bd_it] == *bd ); - ++bd_it; - } - - Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); - for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) - { - BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); - ++cbd_it; - } - ++i; +BOOST_AUTO_TEST_CASE(check_if_boundary_of_boundary_is_zero_non_periodic_case_4_d) { + std::vector<unsigned> sizes(4); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + sizes[3] = 3; + + std::vector<double> data(81, 0); + + int number_of_all_elements = (2 * sizes[0] + 1) * (2 * sizes[1] + 1) * (2 * sizes[2] + 1) * (2 * sizes[3] + 1); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + Bitmap_cubical_complex_base ba(sizes, data); + for (Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + int i = 1; + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd); + int j = 1; + for (Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin(); bd2 != second_bdrange.end(); + ++bd2) { + elems_in_boundary[*bd2] += i * j; + j *= -1; + } + i *= -1; } + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } } +BOOST_AUTO_TEST_CASE(check_if_boundary_of_boundary_is_zero_periodic_case_2d) { + std::vector<unsigned> sizes(2); + sizes[0] = 3; + sizes[1] = 3; + + std::vector<bool> directions_of_periodicity(2, true); + std::vector<double> data(9, 0); + + int number_of_all_elements = (2 * sizes[0]) * (2 * sizes[1]); // *(2*sizes[2]); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + Bitmap_cubical_complex_periodic_boundary_conditions ba(sizes, data, directions_of_periodicity); + for (Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + int i = 1; + + // std::cout << "Element : " << *it << std::endl; + + Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin(); + bd != bdrange.end(); ++bd) { + // std::cout << *bd << " "; + Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd); + int j = 1; + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin(); + bd2 != second_bdrange.end(); ++bd2) { + elems_in_boundary[*bd2] += i * j; + j *= -1; + } + i *= -1; + } + // getchar(); + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } +} +BOOST_AUTO_TEST_CASE(check_if_boundary_of_boundary_is_zero_periodic_case_3d) { + std::vector<unsigned> sizes(3); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + + std::vector<bool> directions_of_periodicity(3, true); + + std::vector<double> data(27, 0); + + int number_of_all_elements = (2 * sizes[0]) * (2 * sizes[1]) * (2 * sizes[2]); + Bitmap_cubical_complex_periodic_boundary_conditions ba(sizes, data, directions_of_periodicity); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + for (Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + // std::cout << "Element : " << *it << std::endl; + + int i = 1; + + Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin(); + bd != bdrange.end(); ++bd) { + Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd); + // std::cout << *bd << " "; + int j = 1; + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin(); + bd2 != second_bdrange.end(); ++bd2) { + elems_in_boundary[*bd2] += i * j; + j *= -1; + } + i *= -1; + } + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } +} +BOOST_AUTO_TEST_CASE(check_if_boundary_of_boundary_is_zero_periodic_case_4d) { + std::vector<unsigned> sizes(4); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + sizes[3] = 3; + + std::vector<bool> directions_of_periodicity(4, true); + + std::vector<double> data(81, 0); + + int number_of_all_elements = (2 * sizes[0]) * (2 * sizes[1]) * (2 * sizes[2]) * (2 * sizes[3]); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + Bitmap_cubical_complex_periodic_boundary_conditions ba(sizes, data, directions_of_periodicity); + for (Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + int i = 1; + + Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin(); + bd != bdrange.end(); ++bd) { + Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd); + int j = 1; + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin(); + bd2 != second_bdrange.end(); ++bd2) { + elems_in_boundary[*bd2] += i * j; + j *= -1; + } + i *= -1; + } -BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check_range_check) -{ - std::vector< double > expected_filtration; - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - - std::vector<unsigned> expected_dimension; - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(2); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - expected_dimension.push_back(1); - expected_dimension.push_back(0); - - std::vector<size_t> expected_boundary; - expected_boundary.push_back(0); - expected_boundary.push_back(2); - expected_boundary.push_back(2); - expected_boundary.push_back(4); - expected_boundary.push_back(0); - expected_boundary.push_back(10); - expected_boundary.push_back(1); - expected_boundary.push_back(11); - expected_boundary.push_back(5); - expected_boundary.push_back(7); - expected_boundary.push_back(2); - expected_boundary.push_back(12); - expected_boundary.push_back(3); - expected_boundary.push_back(13); - expected_boundary.push_back(7); - expected_boundary.push_back(9); - expected_boundary.push_back(4); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(12); - expected_boundary.push_back(12); - expected_boundary.push_back(14); - expected_boundary.push_back(10); - expected_boundary.push_back(20); - expected_boundary.push_back(11); - expected_boundary.push_back(21); - expected_boundary.push_back(15); - expected_boundary.push_back(17); - expected_boundary.push_back(12); - expected_boundary.push_back(22); - expected_boundary.push_back(13); - expected_boundary.push_back(23); - expected_boundary.push_back(17); - expected_boundary.push_back(19); - expected_boundary.push_back(14); - expected_boundary.push_back(24); - expected_boundary.push_back(20); - expected_boundary.push_back(22); - expected_boundary.push_back(22); - expected_boundary.push_back(24); - - - std::vector<size_t> expected_coboundary; - expected_coboundary.push_back(5); - expected_coboundary.push_back(1); - expected_coboundary.push_back(6); - expected_coboundary.push_back(7); - expected_coboundary.push_back(1); - expected_coboundary.push_back(3); - expected_coboundary.push_back(8); - expected_coboundary.push_back(9); - expected_coboundary.push_back(3); - expected_coboundary.push_back(6); - expected_coboundary.push_back(6); - expected_coboundary.push_back(8); - expected_coboundary.push_back(8); - expected_coboundary.push_back(5); - expected_coboundary.push_back(15); - expected_coboundary.push_back(11); - expected_coboundary.push_back(6); - expected_coboundary.push_back(16); - expected_coboundary.push_back(7); - expected_coboundary.push_back(17); - expected_coboundary.push_back(11); - expected_coboundary.push_back(13); - expected_coboundary.push_back(8); - expected_coboundary.push_back(18); - expected_coboundary.push_back(9); - expected_coboundary.push_back(19); - expected_coboundary.push_back(13); - expected_coboundary.push_back(16); - expected_coboundary.push_back(16); - expected_coboundary.push_back(18); - expected_coboundary.push_back(18); - expected_coboundary.push_back(15); - expected_coboundary.push_back(21); - expected_coboundary.push_back(16); - expected_coboundary.push_back(17); - expected_coboundary.push_back(21); - expected_coboundary.push_back(23); - expected_coboundary.push_back(18); - expected_coboundary.push_back(19); - expected_coboundary.push_back(23); - - - - std::vector< unsigned > sizes(2); - sizes[0] = 2; - sizes[1] = 2; - - std::vector< double > data(4); - data[0] = 0; - data[1] = 1; - data[2] = 2; - data[3] = 3; - - Bitmap_cubical_complex_base ba( sizes , data ); - int i = 0; - int bd_it = 0; - int cbd_it = 0; - - Bitmap_cubical_complex_base::All_cells_range range = ba.all_cells_range(); - for ( Bitmap_cubical_complex_base::All_cells_iterator it = range.begin() ; it != range.end() ; ++it ) - { - BOOST_CHECK( expected_filtration[i] == ba.get_cell_data( *it ) ); - BOOST_CHECK( expected_dimension[i] == ba.get_dimension_of_a_cell( *it ) ); - - Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); - for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd ) - { - BOOST_CHECK( expected_boundary[bd_it] == *bd ); - ++bd_it; - } - - Bitmap_cubical_complex_base::Coboundary_range cbdrange = ba.coboundary_range(*it); - for ( Bitmap_cubical_complex_base::Coboundary_iterator cbd = cbdrange.begin() ; cbd != cbdrange.end() ; ++cbd ) - { - BOOST_CHECK( expected_coboundary[cbd_it] == *cbd ); - ++cbd_it; - } - ++i; + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); } + } } -BOOST_AUTO_TEST_CASE(Top_dimensional_cells_iterator_range_check) -{ - std::vector< double > expected_filtration; - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(0); - expected_filtration.push_back(1); - expected_filtration.push_back(1); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(2); - expected_filtration.push_back(3); - expected_filtration.push_back(3); - - - std::vector< unsigned > sizes(2); - sizes[0] = 2; - sizes[1] = 2; - - std::vector< double > data(4); - data[0] = 0; - data[1] = 1; - data[2] = 2; - data[3] = 3; - - Bitmap_cubical_complex_base ba( sizes , data ); - int i = 0; - - Bitmap_cubical_complex_base::Top_dimensional_cells_range range = ba.top_dimensional_cells_range(); - for ( Bitmap_cubical_complex_base::Top_dimensional_cells_iterator it = range.begin() ; it != range.end() ; ++it ) - { - BOOST_CHECK( data[i] == ba.get_cell_data( *it ) ); - BOOST_CHECK( ba.get_dimension_of_a_cell( *it ) == 2 ); - ++i; +BOOST_AUTO_TEST_CASE(compute_incidence_between_cells_test) { + std::vector<unsigned> sizes(3); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + + std::vector<double> data(27, 0); + + int number_of_all_elements = (2 * sizes[0] + 1) * (2 * sizes[1] + 1) * (2 * sizes[1] + 1); + Bitmap_cubical_complex_base ba(sizes, data); + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + for (Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin(); bd != bdrange.end(); ++bd) { + Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd); + for (Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin(); bd2 != second_bdrange.end(); + ++bd2) { + elems_in_boundary[*bd2] += + ba.compute_incidence_between_cells(*it, *bd) * ba.compute_incidence_between_cells(*bd, *bd2); + } } + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } } +BOOST_AUTO_TEST_CASE(compute_incidence_between_cells_test_periodic_boundary_conditions) { + std::vector<unsigned> sizes(3); + sizes[0] = 3; + sizes[1] = 3; + sizes[2] = 3; + + std::vector<bool> directions_of_periodicity(3, true); + std::vector<double> data(27, 0); + + int number_of_all_elements = (2 * sizes[0]) * (2 * sizes[1]) * (2 * sizes[2]); + Bitmap_cubical_complex_periodic_boundary_conditions ba(sizes, data, directions_of_periodicity); + + std::vector<int> elems_in_boundary(number_of_all_elements, 0); + for (Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin(); + it != ba.all_cells_iterator_end(); ++it) { + Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it); + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin(); + bd != bdrange.end(); ++bd) { + // std::cout << *bd << " "; + Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd); + for (Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin(); + bd2 != second_bdrange.end(); ++bd2) { + elems_in_boundary[*bd2] += + ba.compute_incidence_between_cells(*it, *bd) * ba.compute_incidence_between_cells(*bd, *bd2); + } + } + // check if there is anything nonzero in elems_in_boundary + for (size_t i = 0; i != elems_in_boundary.size(); ++i) { + BOOST_CHECK(elems_in_boundary[i] == 0); + } + } +} diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt index 871e46a6..461b6db2 100644 --- a/src/Nerve_GIC/example/CMakeLists.txt +++ b/src/Nerve_GIC/example/CMakeLists.txt @@ -2,22 +2,28 @@ cmake_minimum_required(VERSION 2.6) project(Nerve_GIC_examples) add_executable ( Nerve Nerve.cpp ) -target_link_libraries(Nerve ${Boost_SYSTEM_LIBRARY}) - add_executable ( CoordGIC CoordGIC.cpp ) -target_link_libraries(CoordGIC ${Boost_SYSTEM_LIBRARY}) - add_executable ( FuncGIC FuncGIC.cpp ) -target_link_libraries(FuncGIC ${Boost_SYSTEM_LIBRARY}) - add_executable ( VoronoiGIC VoronoiGIC.cpp ) -target_link_libraries(VoronoiGIC ${Boost_SYSTEM_LIBRARY}) - -#[[file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)]] if (TBB_FOUND) target_link_libraries(Nerve ${TBB_LIBRARIES}) target_link_libraries(CoordGIC ${TBB_LIBRARIES}) target_link_libraries(FuncGIC ${TBB_LIBRARIES}) target_link_libraries(VoronoiGIC ${TBB_LIBRARIES}) -endif()
\ No newline at end of file +endif() + +file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +add_test(NAME Nerve_GIC_example_nerve COMMAND $<TARGET_FILE:Nerve> + "${CMAKE_SOURCE_DIR}/data/points/human.off" "2" "10" "0.3") + +add_test(NAME Nerve_GIC_example_VoronoiGIC COMMAND $<TARGET_FILE:VoronoiGIC> + "${CMAKE_SOURCE_DIR}/data/points/human.off" "100") + +add_test(NAME Nerve_GIC_example_CoordGIC COMMAND $<TARGET_FILE:CoordGIC> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "0") + +add_test(NAME Nerve_GIC_example_FuncGIC COMMAND $<TARGET_FILE:FuncGIC> + "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off" + "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat_PCA1") diff --git a/src/Nerve_GIC/example/CoordGIC.cpp b/src/Nerve_GIC/example/CoordGIC.cpp index e24943df..f23e3f55 100644 --- a/src/Nerve_GIC/example/CoordGIC.cpp +++ b/src/Nerve_GIC/example/CoordGIC.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp index 99e3fb78..1601ddab 100644 --- a/src/Nerve_GIC/example/FuncGIC.cpp +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp index b8b5cd05..cdfe73a7 100644 --- a/src/Nerve_GIC/example/GIC.cpp +++ b/src/Nerve_GIC/example/GIC.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py b/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py index 06d22abd..406264ba 100755 --- a/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py +++ b/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py @@ -2,6 +2,32 @@ import km import numpy as np from collections import defaultdict +"""This file is part of the Gudhi Library. The Gudhi library + (Geometric Understanding in Higher Dimensions) is a generic C++ + library for computational topology. + + Author(s): Mathieu Carriere + + Copyright (C) 2017 INRIA + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see <http://www.gnu.org/licenses/>. +""" + +__author__ = "Mathieu Carriere" +__copyright__ = "Copyright (C) 2017 INRIA" +__license__ = "GPL v3" + network = {} mapper = km.KeplerMapper(verbose=0) data = np.zeros((3,3)) @@ -41,4 +67,4 @@ network["meta"] = lens mapper.visualize(network, color_function = color, path_html="SC.html", title=dat, graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0, -height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M)
\ No newline at end of file +height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M) diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp index 70424cda..71f12f5d 100644 --- a/src/Nerve_GIC/example/Nerve.cpp +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -25,7 +25,7 @@ void usage(int nbArgs, char * const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../../../data/points/human.off 2 10 0.3 --v \n"; + std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 --v \n"; exit(-1); // ----- >> } diff --git a/src/Nerve_GIC/example/VoronoiGIC.cpp b/src/Nerve_GIC/example/VoronoiGIC.cpp index 73a776a6..8ce64819 100644 --- a/src/Nerve_GIC/example/VoronoiGIC.cpp +++ b/src/Nerve_GIC/example/VoronoiGIC.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Nerve_GIC/example/km.py.COPYRIGHT b/src/Nerve_GIC/example/km.py.COPYRIGHT new file mode 100644 index 00000000..bef7b121 --- /dev/null +++ b/src/Nerve_GIC/example/km.py.COPYRIGHT @@ -0,0 +1,26 @@ +km.py is a fork of https://github.com/MLWave/kepler-mapper. +Only the visualization part has been kept (Mapper part has been removed). + +This file has te following Copyright : + +The MIT License (MIT) + +Copyright (c) 2015 Triskelion - HJ van Veen - info@mlwave.com + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 2520042d..ca8727a0 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -799,7 +799,7 @@ class Cover_complex { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) graphic << simplices[i][0] << " " << simplices[i][1] << std::endl; graphic.close(); - std::cout << "SC.txt generated. It can be visualized with e.g. python visu.py and firefox." << std::endl; + std::cout << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; } diff --git a/src/Nerve_GIC/test/CMakeLists.txt b/src/Nerve_GIC/test/CMakeLists.txt index 627778fa..77efe043 100644 --- a/src/Nerve_GIC/test/CMakeLists.txt +++ b/src/Nerve_GIC/test/CMakeLists.txt @@ -1,23 +1,14 @@ cmake_minimum_required(VERSION 2.6) project(Graph_induced_complex_tests) -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() +include(GUDHI_test_coverage) -add_executable ( graph_induced_complex_UT test_GIC.cpp ) -target_link_libraries(graph_induced_complex_UT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +add_executable ( Graph_induced_complex_test_unit test_GIC.cpp ) +target_link_libraries(Graph_induced_complex_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) if (TBB_FOUND) - target_link_libraries(graph_induced_complex_UT ${TBB_LIBRARIES}) + target_link_libraries(Graph_induced_complex_test_unit ${TBB_LIBRARIES}) endif() file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) -add_test(graph_induced_complex_UT ${CMAKE_CURRENT_BINARY_DIR}/graph_induced_complex_UT - # XML format for Jenkins xUnit plugin - --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/graph_induced_complex_UT.xml --log_level=test_suite --report_level=no) +gudhi_add_coverage_test(Graph_induced_complex_test_unit) diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp index e8c6ed4a..d7af7680 100644 --- a/src/Nerve_GIC/test/test_GIC.cpp +++ b/src/Nerve_GIC/test/test_GIC.cpp @@ -4,7 +4,7 @@ * * Author(s): Mathieu Carrière * - * Copyright (C) 2017 INRIA Saclay (France) + * Copyright (C) 2017 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index f47de4c3..e82ef04c 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -73,24 +73,18 @@ if(CGAL_FOUND) target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY}) add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY}) - add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) - target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES}) - target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) add_test(NAME Persistent_cohomology_example_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") add_test(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $<TARGET_FILE:exact_alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") - add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") install(TARGETS alpha_complex_3d_persistence DESTINATION bin) install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin) - install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) @@ -111,7 +105,7 @@ if(CGAL_FOUND) add_test(NAME Persistent_cohomology_example_alpha_complex COMMAND $<TARGET_FILE:alpha_complex_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") add_test(NAME Persistent_cohomology_example_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:periodic_alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "2" "0") + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "3" "1.0") add_test(NAME Persistent_cohomology_example_custom_persistence_sort COMMAND $<TARGET_FILE:custom_persistence_sort>) install(TARGETS alpha_complex_persistence DESTINATION bin) @@ -119,4 +113,31 @@ if(CGAL_FOUND) install(TARGETS custom_persistence_sort DESTINATION bin) endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + + if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable(weighted_periodic_alpha_complex_3d_persistence weighted_periodic_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(weighted_periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Persistent_cohomology_example_weigted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_periodic_alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "3" "1.0") + + install(TARGETS weighted_periodic_alpha_complex_3d_persistence DESTINATION bin) + + endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) + + add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin) + + add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") + endif(CGAL_FOUND) diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_helper.h b/src/Persistent_cohomology/example/alpha_complex_3d_helper.h index 7865e4ec..6b3b7d5d 100644 --- a/src/Persistent_cohomology/example/alpha_complex_3d_helper.h +++ b/src/Persistent_cohomology/example/alpha_complex_3d_helper.h @@ -23,7 +23,7 @@ #ifndef ALPHA_COMPLEX_3D_HELPER_H_ #define ALPHA_COMPLEX_3D_HELPER_H_ -template<class Vertex_list, class Cell_handle> +template <class Vertex_list, class Cell_handle> Vertex_list from_cell(const Cell_handle& ch) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -35,7 +35,7 @@ Vertex_list from_cell(const Cell_handle& ch) { return the_list; } -template<class Vertex_list, class Facet> +template <class Vertex_list, class Facet> Vertex_list from_facet(const Facet& fct) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -49,7 +49,7 @@ Vertex_list from_facet(const Facet& fct) { return the_list; } -template<class Vertex_list, class Edge_3> +template <class Vertex_list, class Edge_3> Vertex_list from_edge(const Edge_3& edg) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -63,7 +63,7 @@ Vertex_list from_edge(const Edge_3& edg) { return the_list; } -template<class Vertex_list, class Vertex_handle> +template <class Vertex_list, class Vertex_handle> Vertex_list from_vertex(const Vertex_handle& vh) { Vertex_list the_list; #ifdef DEBUG_TRACES diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp index 594efb45..26196a6f 100644 --- a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp @@ -39,6 +39,7 @@ #include <utility> #include <list> #include <vector> +#include <cstdlib> #include "alpha_complex_3d_helper.h" @@ -56,10 +57,10 @@ using Point_3 = Kernel::Point_3; // filtration with alpha values needed type definition using Alpha_value_type = Alpha_shape_3::FT; using Object = CGAL::Object; -using Dispatch = CGAL::Dispatch_output_iterator< - CGAL::cpp11::tuple<Object, Alpha_value_type>, - CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, + std::back_insert_iterator<std::vector<Alpha_value_type> > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -70,19 +71,19 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; -using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; -using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; +using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>; void usage(const std::string& progName) { - std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; + std::cerr << "Usage: " << progName + << " path_to_the_OFF_file coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); } -int main(int argc, char * const argv[]) { +int main(int argc, char* const argv[]) { // program args management if (argc != 4) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; @@ -90,13 +91,7 @@ int main(int argc, char * const argv[]) { } int coeff_field_characteristic = atoi(argv[2]); - - Filtration_value min_persistence = 0.0; - int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); - if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { - std::cerr << "Error: " << argv[3] << " is not correct\n"; - usage(argv[0]); - } + Filtration_value min_persistence = strtof(argv[3], nullptr); // Read points from file std::string offInputFile(argv[1]); @@ -143,28 +138,28 @@ int main(int argc, char * const argv[]) { Filtration_value filtration_max = 0.0; for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object - if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } - } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) { vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } @@ -210,7 +205,6 @@ int main(int argc, char * const argv[]) { std::cout << "facets \t\t" << count_facets << std::endl; std::cout << "cells \t\t" << count_cells << std::endl; - std::cout << "Information of the Simplex Tree: " << std::endl; std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; @@ -229,7 +223,7 @@ int main(int argc, char * const argv[]) { std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; // Compute the persistence diagram of the complex - PCOH pcoh(simplex_tree); + Persistent_cohomology pcoh(simplex_tree, true); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); diff --git a/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp index 08f3b334..2e2bfd2f 100644 --- a/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2014 INRIA Saclay (France) + * Copyright (C) 2014 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -39,6 +39,7 @@ #include <utility> #include <list> #include <vector> +#include <cstdlib> #include "alpha_complex_3d_helper.h" @@ -57,10 +58,10 @@ using Point_3 = Kernel::Point_3; // filtration with alpha values needed type definition using Alpha_value_type = Alpha_shape_3::FT; using Object = CGAL::Object; -using Dispatch = CGAL::Dispatch_output_iterator< - CGAL::cpp11::tuple<Object, Alpha_value_type>, - CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, + std::back_insert_iterator<std::vector<Alpha_value_type> > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -71,19 +72,19 @@ using Vertex_list = std::list<Vertex_handle>; using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; -using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; -using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; +using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>; -void usage(char * const progName) { - std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; +void usage(const std::string& progName) { + std::cerr << "Usage: " << progName + << " path_to_the_OFF_file coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); } -int main(int argc, char * const argv[]) { +int main(int argc, char* const argv[]) { // program args management if (argc != 4) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; @@ -91,13 +92,7 @@ int main(int argc, char * const argv[]) { } int coeff_field_characteristic = atoi(argv[2]); - - Filtration_value min_persistence = 0.0; - int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); - if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { - std::cerr << "Error: " << argv[3] << " is not correct\n"; - usage(argv[0]); - } + Filtration_value min_persistence = strtof(argv[3], nullptr); // Read points from file std::string offInputFile(argv[1]); @@ -144,28 +139,28 @@ int main(int argc, char * const argv[]) { Filtration_value filtration_max = 0.0; for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object - if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } - } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) { vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } @@ -212,7 +207,6 @@ int main(int argc, char * const argv[]) { std::cout << "facets \t\t" << count_facets << std::endl; std::cout << "cells \t\t" << count_cells << std::endl; - std::cout << "Information of the Simplex Tree: " << std::endl; std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; @@ -231,7 +225,7 @@ int main(int argc, char * const argv[]) { std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; // Compute the persistence diagram of the complex - PCOH pcoh(simplex_tree); + Persistent_cohomology pcoh(simplex_tree, true); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); diff --git a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp index 74940387..c6d3e236 100644 --- a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp @@ -63,10 +63,10 @@ using Point_3 = PK::Point_3; // filtration with alpha values needed type definition using Alpha_value_type = Alpha_shape_3::FT; using Object = CGAL::Object; -using Dispatch = CGAL::Dispatch_output_iterator< - CGAL::cpp11::tuple<Object, Alpha_value_type>, - CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, + std::back_insert_iterator<std::vector<Alpha_value_type> > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -77,27 +77,19 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; -using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; -using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< - ST, Gudhi::persistent_cohomology::Field_Zp >; - -void usage(char * const progName) { - std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_iso_cuboid_3_file coeff_field_characteristic[" << - "integer > 0] min_persistence[float >= -1.0]\n" << - " path_to_OFF_file is the path to your points cloud in OFF format.\n" << - " path_to_iso_cuboid_3_file is the path to the iso cuboid file with the following format :\n" << - " x_min y_min z_min x_max y_max z_max\n" << - " In this example, the periodic cube will be " << - "{ x = [x_min,x_max]; y = [y_min,y_max]; z = [z_min,z_max] }.\n" << - " For more information, please refer to\n" << - " https://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Iso__cuboid__3.html\n"; +using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>; +void usage(const std::string& progName) { + std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_iso_cuboid_3_file " + "coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); } -int main(int argc, char * const argv[]) { +int main(int argc, char* const argv[]) { // program args management if (argc != 5) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; @@ -168,29 +160,28 @@ int main(int argc, char * const argv[]) { Filtration_value filtration_max = 0.0; for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object - if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } - } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) { vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Alpha_shape_3::Vertex_handle * vertex = - CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } @@ -216,7 +207,7 @@ int main(int argc, char * const argv[]) { } } // Construction of the simplex_tree - Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); #ifdef DEBUG_TRACES std::cout << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES @@ -236,7 +227,6 @@ int main(int argc, char * const argv[]) { std::cout << "facets \t\t" << count_facets << std::endl; std::cout << "cells \t\t" << count_cells << std::endl; - std::cout << "Information of the Simplex Tree: " << std::endl; std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; diff --git a/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp index 16b44094..249a7ece 100644 --- a/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp @@ -26,12 +26,17 @@ #include <gudhi/Persistent_cohomology.h> #include <gudhi/Points_3D_off_io.h> +#include <CGAL/config.h> #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> -#include <CGAL/Regular_triangulation_euclidean_traits_3.h> #include <CGAL/Regular_triangulation_3.h> #include <CGAL/Alpha_shape_3.h> #include <CGAL/iterator.h> +// For CGAL < 4.11 +#if CGAL_VERSION_NR < 1041100000 +#include <CGAL/Regular_triangulation_euclidean_traits_3.h> +#endif // CGAL_VERSION_NR < 1041100000 + #include <fstream> #include <cmath> #include <string> @@ -44,26 +49,44 @@ #include "alpha_complex_3d_helper.h" -// Traits +// Alpha_shape_3 templates type definitions using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; + +// For CGAL < 4.11 +#if CGAL_VERSION_NR < 1041100000 using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>; using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>; using Fb = CGAL::Alpha_shape_cell_base_3<Gt>; using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; using Triangulation_3 = CGAL::Regular_triangulation_3<Gt, Tds>; -using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; // From file type definition using Point_3 = Gt::Bare_point; using Weighted_point_3 = Gt::Weighted_point; +// For CGAL >= 4.11 +#else // CGAL_VERSION_NR < 1041100000 +using Rvb = CGAL::Regular_triangulation_vertex_base_3<Kernel>; +using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel,Rvb>; +using Rcb = CGAL::Regular_triangulation_cell_base_3<Kernel>; +using Cb = CGAL::Alpha_shape_cell_base_3<Kernel,Rcb>; +using Tds = CGAL::Triangulation_data_structure_3<Vb,Cb>; +using Triangulation_3 = CGAL::Regular_triangulation_3<Kernel,Tds>; + +// From file type definition +using Point_3 = Triangulation_3::Bare_point; +using Weighted_point_3 = Triangulation_3::Weighted_point; +#endif // CGAL_VERSION_NR < 1041100000 + +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; + // filtration with alpha values needed type definition using Alpha_value_type = Alpha_shape_3::FT; using Object = CGAL::Object; -using Dispatch = CGAL::Dispatch_output_iterator< - CGAL::cpp11::tuple<Object, Alpha_value_type>, - CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, + std::back_insert_iterator<std::vector<Alpha_value_type> > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -74,24 +97,19 @@ using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; -using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; -using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< - ST, Gudhi::persistent_cohomology::Field_Zp >; - -void usage(char * const progName) { - std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; - std::cerr << " path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n"; - std::cerr << " Weights values are explained on CGAL documentation:\n"; - std::cerr << " https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n"; - std::cerr << " https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\n"; +using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>; + +void usage(const std::string& progName) { + std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_weight_file coeff_field_characteristic[integer > " + "0] min_persistence[float >= -1.0]\n"; exit(-1); } -int main(int argc, char * const argv[]) { +int main(int argc, char* const argv[]) { // program args management if (argc != 5) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; @@ -167,29 +185,28 @@ int main(int argc, char * const argv[]) { Filtration_value filtration_max = 0.0; for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object - if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } - } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) { vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Alpha_shape_3::Vertex_handle * vertex = - CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } @@ -215,7 +232,7 @@ int main(int argc, char * const argv[]) { } } // Construction of the simplex_tree - Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); #ifdef DEBUG_TRACES std::cout << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES @@ -235,7 +252,6 @@ int main(int argc, char * const argv[]) { std::cout << "facets \t\t" << count_facets << std::endl; std::cout << "cells \t\t" << count_cells << std::endl; - std::cout << "Information of the Simplex Tree: " << std::endl; std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; diff --git a/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..13634ff7 --- /dev/null +++ b/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp @@ -0,0 +1,281 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <boost/variant.hpp> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Periodic_3_regular_triangulation_traits_3.h> +#include <CGAL/Periodic_3_regular_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/iterator.h> + +#include <fstream> +#include <cmath> +#include <string> +#include <tuple> +#include <map> +#include <utility> +#include <list> +#include <vector> +#include <cstdlib> + +#include "alpha_complex_3d_helper.h" + +// Traits +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using PK = CGAL::Periodic_3_regular_triangulation_traits_3<Kernel>; + +// Vertex type +using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; +using Vb = CGAL::Regular_triangulation_vertex_base_3<PK,DsVb>; +using AsVb = CGAL::Alpha_shape_vertex_base_3<PK,Vb>; +// Cell type +using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; +using Cb = CGAL::Regular_triangulation_cell_base_3<PK,DsCb>; +using AsCb = CGAL::Alpha_shape_cell_base_3<PK,Cb>; +using Tds = CGAL::Triangulation_data_structure_3<AsVb,AsCb>; +using P3RT3 = CGAL::Periodic_3_regular_triangulation_3<PK,Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<P3RT3>; + +using Point_3 = P3RT3::Bare_point; +using Weighted_point_3 = P3RT3::Weighted_point; + +// filtration with alpha values needed type definition +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<Object> >, + std::back_insert_iterator<std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; + +// gudhi type definition +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector<Simplex_tree_vertex>; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<ST, Gudhi::persistent_cohomology::Field_Zp>; + +void usage(const std::string& progName) { + std::cerr << "Usage: " << progName << " path_to_the_OFF_file path_to_weight_file path_to_the_cuboid_file " + "coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; + exit(-1); +} + +int main(int argc, char* const argv[]) { + // program args management + if (argc != 6) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[4]); + Filtration_value min_persistence = strtof(argv[5], nullptr); + + // Read points from file + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_3> lp = off_reader.get_point_cloud(); + + // Read weights information from file + std::ifstream weights_ifstr(argv[2]); + std::vector<Weighted_point_3> wp; + if (weights_ifstr.good()) { + double weight = 0.0; + std::size_t index = 0; + wp.reserve(lp.size()); + // Attempt read the weight in a double format, return false if it fails + while ((weights_ifstr >> weight) && (index < lp.size())) { + wp.push_back(Weighted_point_3(lp[index], weight)); + index++; + } + if (index != lp.size()) { + std::cerr << "Bad number of weights in file " << argv[2] << std::endl; + usage(argv[0]); + } + } else { + std::cerr << "Unable to read file " << argv[2] << std::endl; + usage(argv[0]); + } + + // Read iso_cuboid_3 information from file + std::ifstream iso_cuboid_str(argv[3]); + double x_min, y_min, z_min, x_max, y_max, z_max; + if (iso_cuboid_str.good()) { + iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max; + } else { + std::cerr << "Unable to read file " << argv[3] << std::endl; + usage(argv[0]); + } + + // Define the periodic cube + P3RT3 prt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + prt.insert(wp.begin(), wp.end(), true); + // As prt won't be modified anymore switch to 1-sheeted cover if possible + if (prt.is_triangulation_in_1_sheet()) prt.convert_to_1_sheeted_covering(); + std::cout << "Periodic Delaunay computed." << std::endl; + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + Alpha_shape_3 as(prt, 0, Alpha_shape_3::GENERAL); + + // filtration with alpha values from alpha shape + std::vector<Object> the_objects; + std::vector<Alpha_value_type> the_alpha_values; + + Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; + + // Loop on objects vector + Vertex_list vertex_list; + ST simplex_tree; + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin(); + int dim_max = 0; + Filtration_value filtration_max = 0.0; + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { + vertex_list = from_facet<Vertex_list, Facet>(*facet); + count_facets++; + if (dim_max < 2) { + // Facet is of dim 2 + dim_max = 2; + } + } else if (const Edge_3* edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); + count_edges++; + if (dim_max < 1) { + // Edge_3 is of dim 1 + dim_max = 1; + } + } else if (const Vertex_handle* vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex_tree; + for (auto the_alpha_shape_vertex : vertex_list) { + Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Simplex_tree_vertex vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex)); + } else { + // alpha shape found + Simplex_tree_vertex vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + if (filtr > filtration_max) { + filtration_max = filtr; + } + simplex_tree.insert_simplex(the_simplex_tree, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree, true); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + pcoh.output_diagram(); + + return 0; +} diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index a5cefd6a..6f061922 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -155,7 +155,7 @@ class Tangential_complex { >::type Triangulation; typedef typename Triangulation::Geom_traits Tr_traits; typedef typename Triangulation::Weighted_point Tr_point; - typedef typename Triangulation::Bare_point Tr_bare_point; + typedef typename Tr_traits::Base::Point_d Tr_bare_point; typedef typename Triangulation::Vertex_handle Tr_vertex_handle; typedef typename Triangulation::Full_cell_handle Tr_full_cell_handle; typedef typename Tr_traits::Vector_d Tr_vector; diff --git a/src/cmake/modules/GUDHI_doxygen_target.cmake b/src/cmake/modules/GUDHI_doxygen_target.cmake index d2cb952d..f3e2d9f5 100644 --- a/src/cmake/modules/GUDHI_doxygen_target.cmake +++ b/src/cmake/modules/GUDHI_doxygen_target.cmake @@ -3,6 +3,11 @@ find_package(Doxygen) if(DOXYGEN_FOUND) # configure_file(${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile @ONLY) + #starting from cmake 3.9 the usage of DOXYGEN_EXECUTABLE is deprecated + if(TARGET Doxygen::doxygen) + get_property(DOXYGEN_EXECUTABLE TARGET Doxygen::doxygen PROPERTY IMPORTED_LOCATION) + endif() + add_custom_target(doxygen ${DOXYGEN_EXECUTABLE} ${GUDHI_USER_VERSION_DIR}/Doxyfile WORKING_DIRECTORY ${GUDHI_USER_VERSION_DIR} DEPENDS ${GUDHI_USER_VERSION_DIR}/Doxyfile ${GUDHI_DOXYGEN_DEPENDENCY} diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index f2bbafdc..8269c3bf 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -54,10 +54,12 @@ if(CGAL_FOUND) endforeach(CGAL_INCLUDE_DIR ${CGAL_INCLUDE_DIRS}) endif(NOT CGAL_VERSION VERSION_GREATER 4.9.0) - # For dev version - include_directories(BEFORE "src/common/include/gudhi_patches") - # For user version - include_directories(BEFORE "include/gudhi_patches") + if (NOT CGAL_VERSION VERSION_GREATER 4.11.0) + # For dev version + include_directories(BEFORE "src/common/include/gudhi_patches") + # For user version + include_directories(BEFORE "include/gudhi_patches") + endif (NOT CGAL_VERSION VERSION_GREATER 4.11.0) endif() endif() diff --git a/src/cmake/modules/GUDHI_user_version_target.cmake b/src/cmake/modules/GUDHI_user_version_target.cmake index cff64ad2..4abc2574 100644 --- a/src/cmake/modules/GUDHI_user_version_target.cmake +++ b/src/cmake/modules/GUDHI_user_version_target.cmake @@ -48,7 +48,11 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11) copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI) set(GUDHI_DIRECTORIES "doc;example;concept;utilities") - set(GUDHI_INCLUDE_DIRECTORIES "include/gudhi;include/gudhi_patches") + if (NOT CGAL_VERSION VERSION_GREATER 4.11.0) + set(GUDHI_INCLUDE_DIRECTORIES "include/gudhi;include/gudhi_patches") + else () + set(GUDHI_INCLUDE_DIRECTORIES "include/gudhi") + endif () foreach(GUDHI_MODULE ${GUDHI_MODULES_FULL_LIST}) foreach(GUDHI_DIRECTORY ${GUDHI_DIRECTORIES}) diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx index 45487158..32b91028 100644 --- a/src/cython/cython/simplex_tree.pyx +++ b/src/cython/cython/simplex_tree.pyx @@ -312,7 +312,7 @@ cdef class SimplexTree: """This function returns the persistence of the simplicial complex. :param homology_coeff_field: The homology coefficient field. Must be a - prime number + prime number. Default value is 11. :type homology_coeff_field: int. :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is diff --git a/src/cython/test/test_cubical_complex.py b/src/cython/test/test_cubical_complex.py index 9a365823..0e81554d 100755 --- a/src/cython/test/test_cubical_complex.py +++ b/src/cython/test/test_cubical_complex.py @@ -62,17 +62,17 @@ def test_dimension_or_perseus_file_constructor(): assert cub.__is_defined() == False assert cub.__is_persistence_defined() == False -def test_dimension_constructor(): +def test_dimension_simple_constructor(): cub = CubicalComplex(dimensions=[3, 3], top_dimensional_cells = [1,2,3,4,5,6,7,8,9]) assert cub.__is_defined() == True assert cub.__is_persistence_defined() == False - assert cub.persistence() == [(1, (0.0, 100.0)), (0, (0.0, float('inf')))] + assert cub.persistence() == [(0, (1.0, float('inf')))] assert cub.__is_persistence_defined() == True - assert cub.betti_numbers() == [1, 0] - assert cub.persistent_betti_numbers(0, 1000) == [0, 0] + assert cub.betti_numbers() == [1, 0, 0] + assert cub.persistent_betti_numbers(0, 1000) == [0, 0, 0] -def test_dimension_constructor(): +def test_dimension_file_constructor(): # Create test file test_file = open('CubicalOneSphere.txt', 'w') test_file.write('2\n3\n3\n0\n0\n0\n0\n100\n0\n0\n0\n0\n') |