diff options
| author | Ulrich Bauer <ulrich.bauer@tum.de> | 2016-05-24 16:23:28 +0200 |
|---|---|---|
| committer | Ulrich Bauer <ulrich.bauer@tum.de> | 2016-05-24 16:23:51 +0200 |
| commit | 05747807349f86b372f746edc384d7a26770d6f1 (patch) | |
| tree | f72505575e77b38c97ea4c9690a47d3d43d94c78 | |
| parent | 8b7b6c13ee106480c15134e3d3397021ae750ac5 (diff) | |
| parent | 7a6a20ef738a1abe182e775dc8c9a1267ee11e63 (diff) | |
updated Google hash map support
| -rw-r--r-- | .clang-format | 11 | ||||
| -rw-r--r-- | benchmarks/benchmarks.txt | 10 | ||||
| -rw-r--r-- | benchmarks/examples.txt | 4 | ||||
| -rw-r--r-- | benchmarks/sphere_3_192/dionysus.txt | 28 | ||||
| -rw-r--r-- | benchmarks/sphere_3_192/dipha.txt | 58 | ||||
| -rw-r--r-- | benchmarks/sphere_3_192/gudhi.txt | 18 | ||||
| -rw-r--r-- | benchmarks/sphere_3_192/phat.txt | 17 | ||||
| -rw-r--r-- | benchmarks/sphere_3_192/ripser.txt | 268 | ||||
| -rw-r--r-- | benchmarks/sphere_3_64/dionysus.txt | 28 | ||||
| -rw-r--r-- | benchmarks/sphere_3_64/dipha.txt | 58 | ||||
| -rw-r--r-- | benchmarks/sphere_3_64/gudhi.txt | 18 | ||||
| -rw-r--r-- | benchmarks/sphere_3_64/phat.txt | 17 | ||||
| -rw-r--r-- | benchmarks/sphere_3_64/ripser.txt | 103 | ||||
| -rw-r--r-- | random16.compressed_lower_distance_matrix | 1 | ||||
| -rw-r--r-- | random20.compressed_lower_distance_matrix | 1 | ||||
| -rw-r--r-- | ripser.cpp | 1877 | ||||
| -rw-r--r-- | ripser.xcodeproj/project.pbxproj | 2 |
17 files changed, 1352 insertions, 1167 deletions
diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..1975b19 --- /dev/null +++ b/.clang-format @@ -0,0 +1,11 @@ +BasedOnStyle: LLVM +IndentWidth: 4 +TabWidth: 4 +ColumnLimit: 100 +AccessModifierOffset: -4 +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AllowShortBlocksOnASingleLine: true +AllowShortFunctionsOnASingleLine: true +DerivePointerAlignment: true +UseTab: ForIndentation diff --git a/benchmarks/benchmarks.txt b/benchmarks/benchmarks.txt new file mode 100644 index 0000000..9fbc688 --- /dev/null +++ b/benchmarks/benchmarks.txt @@ -0,0 +1,10 @@ +/usr/bin/time -l ~/Source/Dionysus/examples/cohomology/rips-pairwise-cohomology -s3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_192_points.dat + +/usr/bin/time -l ~/Bitbucket/phat/benchmark --primal --bit_tree_pivot_column --twist ~/Bitbucket/phat-paper/benchmark/sphere_3_192_coboundary.phat + +/usr/bin/time -l ~/Bitbucket/dipha/dipha --benchmark --upper_dim 3 --dual ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex /dev/null + +/usr/bin/time -l ~/Source/Gudhi_library_1.3.0/example/Persistent_cohomology/rips_persistence -d3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_192_points.dat -o/dev/null + +c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D FILE_FORMAT_DIPHA -D PRINT_PERSISTENCE_PAIRS && /usr/bin/time -l ./ripser --top_dim 2 ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex + diff --git a/benchmarks/examples.txt b/benchmarks/examples.txt new file mode 100644 index 0000000..39afeb6 --- /dev/null +++ b/benchmarks/examples.txt @@ -0,0 +1,4 @@ +c++ -std=c++14 ripser.cpp -o ripser -Ofast -D NDEBUG -D STORE_DIAMETERS -D FILE_FORMAT_UPPER_TRIANGULAR_CSV -D PRINT_PERSISTENCE_PAIRS && /usr/bin/time -l ./ripser --top_dim 1 ~/Downloads/avian_all_nt_concat_jukes_cantor.csv + +c++ -std=c++14 ripser.cpp -o ripser -Ofast -D NDEBUG -D STORE_DIAMETERS -D FILE_FORMAT_DIPHA -D PRINT_PERSISTENCE_PAIRS && /usr/bin/time -l ./ripser --top_dim 2 ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex + diff --git a/benchmarks/sphere_3_192/dionysus.txt b/benchmarks/sphere_3_192/dionysus.txt new file mode 100644 index 0000000..56e23d1 --- /dev/null +++ b/benchmarks/sphere_3_192/dionysus.txt @@ -0,0 +1,28 @@ +uli:Dionysus uli$ /usr/bin/time -l ~/Source/Dionysus/examples/cohomology/rips-pairwise-cohomology -s3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_192_points.dat +Boundary matrix: +Cocycles: *.ccl +Vertices: +Diagram: +Simplex vector generated, size: 56050288 + +0% 10 20 30 40 50 60 70 80 90 100% +|----|----|----|----|----|----|----|----|----|----| +*************************************************** +Rips timer : Elapsed time [393.20] seconds +Persistence timer : Elapsed time [322.41] seconds +Total timer : Elapsed time [807.98] seconds + 843.40 real 787.80 user 26.13 sys +3240779776 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 4808144 page reclaims + 0 page faults + 0 swaps + 11 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 11585 voluntary context switches + 691823 involuntary context switches diff --git a/benchmarks/sphere_3_192/dipha.txt b/benchmarks/sphere_3_192/dipha.txt new file mode 100644 index 0000000..38cc112 --- /dev/null +++ b/benchmarks/sphere_3_192/dipha.txt @@ -0,0 +1,58 @@ +uli:dipha uli$ /usr/bin/time -l ~/Bitbucket/dipha/dipha --benchmark --upper_dim 3 --dual ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex /dev/null + +Input filename: +/Users/uli/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex + +upper_dim: 3 + +Number of processes used: +1 + +Detailed information for rank 0: + time prior mem peak mem bytes recv + 0.0s 3 MB 4 MB 0 MB complex.load_binary( input_filename, upper_dim ); + +Number of cells in input: +56050288 + 20.1s 4 MB 1714 MB 0 MB get_filtration_to_cell_map( complex, dualize, filtration_to_cell_map ); + 11.5s 431 MB 2756 MB 855 MB get_cell_to_filtration_map( complex.get_num_cells(), filtration_to_cell_map, cell_to_filtration_map ); + 1.9s 1713 MB 2756 MB 0 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 0.0s 2144 MB 2756 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 1.7s 2145 MB 2756 MB 53 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 0.1s 2226 MB 2756 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 35.9s 2244 MB 2809 MB 3349 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 17.1s 1359 MB 2809 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 14.9s 943 MB 3638 MB 106 MB dipha::outputs::save_persistence_diagram( output_filename, complex, filtration_to_cell_map, reduced_columns, dualize, upper_dim ); + +Overall running time in seconds: +104.5 + +Reduction kernel running time in seconds: +17.2 + +Overall peak mem in GB of all ranks: +3.6 + +Individual peak mem in GB of per rank: +3.6 + +Maximal communication traffic (without sorting) in GB between any pair of nodes: +4.3 + +Total communication traffic (without sorting) in GB between all pairs of nodes: +4.3 + 104.51 real 71.06 user 25.97 sys +3815686144 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 8788553 page reclaims + 0 page faults + 0 swaps + 0 block input operations + 4 block output operations + 0 messages sent + 0 messages received + 0 signals received + 2337 voluntary context switches + 238131 involuntary context switches diff --git a/benchmarks/sphere_3_192/gudhi.txt b/benchmarks/sphere_3_192/gudhi.txt new file mode 100644 index 0000000..87763c8 --- /dev/null +++ b/benchmarks/sphere_3_192/gudhi.txt @@ -0,0 +1,18 @@ +uli:Gudhi_library_1.3.0 uli$ /usr/bin/time -l ~/Source/Gudhi_library_1.3.0/example/Persistent_cohomology/rips_persistence -d3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_192_points.dat -o/dev/null +The complex contains 56050288 simplices + and has dimension 3 + 71.98 real 120.27 user 5.86 sys +2905890816 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 1535169 page reclaims + 129 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 1957 voluntary context switches + 168798 involuntary context switches diff --git a/benchmarks/sphere_3_192/phat.txt b/benchmarks/sphere_3_192/phat.txt new file mode 100644 index 0000000..3ebcfff --- /dev/null +++ b/benchmarks/sphere_3_192/phat.txt @@ -0,0 +1,17 @@ +uli:phat uli$ /usr/bin/time -l ~/Bitbucket/phat/benchmark --primal --bit_tree_pivot_column --twist ~/Bitbucket/phat-paper/benchmark/sphere_3_192_coboundary.phat +/Users/uli/Bitbucket/phat-paper/benchmark/sphere_3_192_coboundary.phat, bit_tree_pivot_column, twist, primal, Reduction time: 2.649s + 16.91 real 12.22 user 4.38 sys +3742277632 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 1568913 page reclaims + 8 page faults + 0 swaps + 10 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 20 voluntary context switches + 13471 involuntary context switches diff --git a/benchmarks/sphere_3_192/ripser.txt b/benchmarks/sphere_3_192/ripser.txt new file mode 100644 index 0000000..60b7939 --- /dev/null +++ b/benchmarks/sphere_3_192/ripser.txt @@ -0,0 +1,268 @@ +uli:ripser uli$ c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D FILE_FORMAT_DIPHA -D PRINT_PERSISTENCE_PAIRS && /usr/bin/time -l ./ripser --top_dim 2 ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_192.complex +distance matrix with 192 points +distance matrix transformed to lower triangular form +value range: [0.00367531,2] +persistence intervals in dim 0: + [0,0.145018) + [0,0.0831521) + [0,0.115939) + [0,0.202564) + [0,0.0217315) + [0,0.127383) + [0,0.0387451) + [0,0.0623736) + [0,0.282169) + [0,0.126595) + [0,0.195722) + [0,0.0792212) + [0,0.149226) + [0,0.2286) + [0,0.197377) + [0,0.0651112) + [0,0.138674) + [0,0.147427) + [0,0.124774) + [0,0.293466) + [0,0.17493) + [0,0.131438) + [0,0.151643) + [0,0.208063) + [0,0.129741) + [0,0.127864) + [0,0.131638) + [0,0.0959033) + [0,0.150747) + [0,0.119187) + [0,0.0870489) + [0,0.165592) + [0,0.121955) + [0,0.1569) + [0,0.0433426) + [0,0.280629) + [0,0.0487061) + [0,0.0746054) + [0,0.0728647) + [0,0.0471186) + [0,0.0988752) + [0,0.111131) + [0,0.141823) + [0,0.249345) + [0,0.283881) + [0,0.0522132) + [0,0.157293) + [0,0.0496187) + [0,0.143151) + [0,0.135402) + [0,0.192733) + [0,0.057488) + [0,0.0311057) + [0,0.159282) + [0,0.109526) + [0,0.234112) + [0,0.084341) + [0,0.0593071) + [0,0.133033) + [0,0.240632) + [0,0.248879) + [0,0.134886) + [0,0.228917) + [0,0.249339) + [0,0.204544) + [0,0.048826) + [0,0.143085) + [0,0.231675) + [0,0.191661) + [0,0.0637053) + [0,0.0560443) + [0,0.142537) + [0,0.00367531) + [0,0.257997) + [0,0.037954) + [0,0.204208) + [0,0.201961) + [0,0.215759) + [0,0.162464) + [0,0.0692368) + [0,0.145725) + [0,0.205389) + [0,0.203812) + [0,0.15587) + [0,0.123527) + [0,0.108768) + [0,0.165314) + [0,0.0782455) + [0,0.0173134) + [0,0.176306) + [0,0.0741239) + [0,0.0831358) + [0,0.266592) + [0,0.0515157) + [0,0.223232) + [0,0.0120098) + [0,0.122089) + [0,0.147203) + [0,0.283562) + [0,0.304062) + [0,0.253227) + [0,0.164809) + [0,0.149017) + [0,0.133525) + [0,0.237571) + [0,0.0743236) + [0,0.0672609) + [0,0.218526) + [0,0.0559886) + [0,0.165907) + [0,0.0800348) + [0,0.192356) + [0,0.164815) + [0,0.142834) + [0,0.174172) + [0,0.0714516) + [0,0.227933) + [0,0.0702075) + [0,0.18428) + [0,0.0970981) + [0,0.133722) + [0,0.0753955) + [0,0.0884133) + [0,0.15244) + [0,0.157373) + [0,0.234969) + [0,0.173412) + [0,0.193791) + [0,0.147289) + [0,0.111798) + [0,0.277191) + [0,0.167455) + [0,0.189106) + [0,0.179609) + [0,0.297053) + [0,0.158413) + [0,0.165494) + [0,0.199673) + [0,0.227524) + [0,0.076909) + [0,0.332695) + [0,0.0802451) + [0,0.0994522) + [0,0.154095) + [0,0.274611) + [0,0.24348) + [0,0.207016) + [0,0.218561) + [0,0.270996) + [0,0.214018) + [0,0.323964) + [0,0.184757) + [0,0.193804) + [0,0.276431) + [0,0.175016) + [0,0.142152) + [0,0.0857401) + [0,0.207714) + [0,0.00391946) + [0,0.168587) + [0,0.192406) + [0,0.246491) + [0,0.148506) + [0,0.25725) + [0,0.165202) + [0,0.0814685) + [0,0.247042) + [0,0.0633425) + [0,0.253163) + [0,0.157251) + [0,0.285416) + [0,0.173586) + [0,0.210372) + [0,0.266635) + [0,0.316023) + [0,0.287316) + [0,0.264334) + [0,0.208851) + [0,0.277582) + [0,0.288591) + [0,0.231179) + [0,0.190642) + [0,0.24509) + [0,0.319493) + [0,0.207261) + [0,0.173684) + [0,0.264483) + [0,0.289189) + [0,0.23879) + [0,0.297933) + [0,0.268806) + [0, ) +persistence intervals in dim 1: + [0.542696,0.558863) + [0.531636,0.578093) + [0.530723,0.576869) + [0.463389,0.505345) + [0.445398,0.448892) + [0.443911,0.54761) + [0.431628,0.477277) + [0.413789,0.487379) + [0.412572,0.46308) + [0.411549,0.471715) + [0.409968,0.478461) + [0.386278,0.4922) + [0.381084,0.421374) + [0.377729,0.514046) + [0.377147,0.414788) + [0.377019,0.434385) + [0.374531,0.477153) + [0.370051,0.483155) + [0.361715,0.403181) + [0.354747,0.454956) + [0.352356,0.541947) + [0.350913,0.369543) + [0.35058,0.580726) + [0.347806,0.638039) + [0.347388,0.559978) + [0.344005,0.387193) + [0.34298,0.65758) + [0.339394,0.463666) + [0.33836,0.350411) + [0.33392,0.35111) + [0.33364,0.3797) + [0.331259,0.356381) + [0.329066,0.386148) + [0.324571,0.343083) + [0.324159,0.399094) + [0.322683,0.482958) + [0.318032,0.417462) + [0.317521,0.349622) + [0.316534,0.682559) + [0.309419,0.37151) + [0.309336,0.310216) + [0.308779,0.346297) + [0.304761,0.647803) + [0.301986,0.478334) + [0.301045,0.303943) + [0.29997,0.311993) + [0.298341,0.360737) + [0.279016,0.301331) + [0.2576,0.373749) + [0.253337,0.292286) + [0.24546,0.248903) + [0.240869,0.288747) + [0.189652,0.197515) +persistence intervals in dim 2: + [0.720484,1.65562) + 2.22 real 2.09 user 0.11 sys + 182439936 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 44546 page reclaims + 9 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 1 voluntary context switches + 1307 involuntary context switches diff --git a/benchmarks/sphere_3_64/dionysus.txt b/benchmarks/sphere_3_64/dionysus.txt new file mode 100644 index 0000000..52b879d --- /dev/null +++ b/benchmarks/sphere_3_64/dionysus.txt @@ -0,0 +1,28 @@ +uli:Dionysus uli$ /usr/bin/time -l ~/Source/Dionysus/examples/cohomology/rips-pairwise-cohomology -s3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_64_points.dat +Boundary matrix: +Cocycles: *.ccl +Vertices: +Diagram: +Simplex vector generated, size: 679120 + +0% 10 20 30 40 50 60 70 80 90 100% +|----|----|----|----|----|----|----|----|----|----| +*************************************************** +Rips timer : Elapsed time [3.01] seconds +Persistence timer : Elapsed time [1.98] seconds +Total timer : Elapsed time [5.85] seconds + 5.94 real 5.85 user 0.05 sys + 76812288 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 18766 page reclaims + 0 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 0 voluntary context switches + 2160 involuntary context switches diff --git a/benchmarks/sphere_3_64/dipha.txt b/benchmarks/sphere_3_64/dipha.txt new file mode 100644 index 0000000..67534d6 --- /dev/null +++ b/benchmarks/sphere_3_64/dipha.txt @@ -0,0 +1,58 @@ +uli:dipha uli$ /usr/bin/time -l ~/Bitbucket/dipha/dipha --benchmark --upper_dim 3 --dual ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_64.complex /dev/null + +Input filename: +/Users/uli/Bitbucket/phat-paper/benchmark/dipha/sphere_3_64.complex + +upper_dim: 3 + +Number of processes used: +1 + +Detailed information for rank 0: + time prior mem peak mem bytes recv + 0.0s 3 MB 3 MB 0 MB complex.load_binary( input_filename, upper_dim ); + +Number of cells in input: +679120 + 0.2s 3 MB 24 MB 0 MB get_filtration_to_cell_map( complex, dualize, filtration_to_cell_map ); + 0.1s 24 MB 69 MB 10 MB get_cell_to_filtration_map( complex.get_num_cells(), filtration_to_cell_map, cell_to_filtration_map ); + 0.0s 84 MB 85 MB 0 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 0.0s 85 MB 85 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 0.0s 85 MB 88 MB 1 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 0.0s 88 MB 89 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 0.3s 89 MB 122 MB 38 MB generate_unreduced_columns( complex, filtration_to_cell_map, cell_to_filtration_map, cur_dim, dualize, unreduced_columns ); + 0.1s 102 MB 142 MB 0 MB reduction_kernel( complex.get_num_cells(), unreduced_columns, reduced_columns ); + 0.1s 141 MB 142 MB 3 MB dipha::outputs::save_persistence_diagram( output_filename, complex, filtration_to_cell_map, reduced_columns, dualize, upper_dim ); + +Overall running time in seconds: +0.8 + +Reduction kernel running time in seconds: +0.1 + +Overall peak mem in GB of all ranks: +0.1 + +Individual peak mem in GB of per rank: +0.1 + +Maximal communication traffic (without sorting) in GB between any pair of nodes: +0.1 + +Total communication traffic (without sorting) in GB between all pairs of nodes: +0.1 + 0.81 real 0.65 user 0.13 sys + 149577728 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 51467 page reclaims + 0 page faults + 0 swaps + 0 block input operations + 1 block output operations + 0 messages sent + 0 messages received + 0 signals received + 7 voluntary context switches + 1069 involuntary context switches diff --git a/benchmarks/sphere_3_64/gudhi.txt b/benchmarks/sphere_3_64/gudhi.txt new file mode 100644 index 0000000..ee69c4a --- /dev/null +++ b/benchmarks/sphere_3_64/gudhi.txt @@ -0,0 +1,18 @@ +uli:Gudhi_library_1.3.0 uli$ /usr/bin/time -l ~/Source/Gudhi_library_1.3.0/example/Persistent_cohomology/rips_persistence -d3 -p2 ~/Bitbucket/phat-paper/benchmark/point\ cloud/sphere_3_64_points.dat -o/dev/null +The complex contains 679120 simplices + and has dimension 3 + 0.45 real 0.75 user 0.02 sys + 37740544 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 9124 page reclaims + 122 page faults + 0 swaps + 1 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 14 voluntary context switches + 2324 involuntary context switches diff --git a/benchmarks/sphere_3_64/phat.txt b/benchmarks/sphere_3_64/phat.txt new file mode 100644 index 0000000..0dff324 --- /dev/null +++ b/benchmarks/sphere_3_64/phat.txt @@ -0,0 +1,17 @@ +uli:phat uli$ /usr/bin/time -l ~/Bitbucket/phat/benchmark --primal --bit_tree_pivot_column --twist ~/Bitbucket/phat-paper/benchmark/sphere_3_64_coboundary.phat +/Users/uli/Bitbucket/phat-paper/benchmark/sphere_3_64_coboundary.phat, bit_tree_pivot_column, twist, primal, Reduction time: 0.019s + 0.22 real 0.15 user 0.04 sys + 47742976 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 11669 page reclaims + 2 page faults + 0 swaps + 2 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 4 voluntary context switches + 106 involuntary context switches diff --git a/benchmarks/sphere_3_64/ripser.txt b/benchmarks/sphere_3_64/ripser.txt new file mode 100644 index 0000000..016143d --- /dev/null +++ b/benchmarks/sphere_3_64/ripser.txt @@ -0,0 +1,103 @@ +uli:ripser uli$ c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D FILE_FORMAT_DIPHA -D PRINT_PERSISTENCE_PAIRS && /usr/bin/time -l ./ripser --top_dim 2 ~/Bitbucket/phat-paper/benchmark/dipha/sphere_3_64.complex +distance matrix with 64 points +distance matrix transformed to lower triangular form +value range: [0.00391946,1.99927] +persistence intervals in dim 0: + [0,0.147289) + [0,0.111798) + [0,0.279354) + [0,0.379849) + [0,0.217485) + [0,0.179609) + [0,0.297053) + [0,0.227524) + [0,0.334523) + [0,0.205129) + [0,0.266643) + [0,0.076909) + [0,0.424333) + [0,0.0802451) + [0,0.141103) + [0,0.154095) + [0,0.280747) + [0,0.561426) + [0,0.255364) + [0,0.218561) + [0,0.344005) + [0,0.457212) + [0,0.339394) + [0,0.184757) + [0,0.193804) + [0,0.276431) + [0,0.276383) + [0,0.273482) + [0,0.162738) + [0,0.208851) + [0,0.00391946) + [0,0.168587) + [0,0.231288) + [0,0.370337) + [0,0.151879) + [0,0.336643) + [0,0.176636) + [0,0.0814685) + [0,0.287316) + [0,0.0633425) + [0,0.46675) + [0,0.169115) + [0,0.387245) + [0,0.354112) + [0,0.236525) + [0,0.377019) + [0,0.396994) + [0,0.58991) + [0,0.374531) + [0,0.374138) + [0,0.439977) + [0,0.487279) + [0,0.231179) + [0,0.200333) + [0,0.401576) + [0,0.574945) + [0,0.252863) + [0,0.173684) + [0,0.284319) + [0,0.355419) + [0,0.23879) + [0,0.381084) + [0,0.532104) + [0, ) +persistence intervals in dim 1: + [0.84172,1.04196) + [0.763446,0.863406) + [0.727244,0.74155) + [0.708398,0.757279) + [0.655113,0.849319) + [0.62495,0.753672) + [0.60904,0.682559) + [0.607973,0.672219) + [0.606907,0.829538) + [0.588045,1.04527) + [0.57698,0.769491) + [0.559498,0.586113) + [0.525224,0.645181) + [0.445398,0.448892) + [0.412572,0.478334) + [0.376457,0.63549) +persistence intervals in dim 2: + [1.10328,1.67029) + 0.05 real 0.05 user 0.00 sys + 4300800 maximum resident set size + 0 average shared memory size + 0 average unshared data size + 0 average unshared stack size + 1054 page reclaims + 9 page faults + 0 swaps + 0 block input operations + 0 block output operations + 0 messages sent + 0 messages received + 0 signals received + 1 voluntary context switches + 10 involuntary context switches diff --git a/random16.compressed_lower_distance_matrix b/random16.compressed_lower_distance_matrix new file mode 100644 index 0000000..74e5c29 --- /dev/null +++ b/random16.compressed_lower_distance_matrix @@ -0,0 +1 @@ +85,120,48,10,73,44,37,47,76,42,14,104,116,97,58,57,111,78,93,22,46,91,55,45,56,107,7,23,62,105,81,72,94,34,52,30,77,117,17,103,19,28,99,29,89,108,101,43,25,53,65,40,98,75,13,51,92,112,36,100,115,2,82,70,32,38,12,90,24,68,86,64,8,59,49,69,114,3,95,35,61,15,71,31,74,41,5,4,16,84,50,27,102,87,39,118,109,21,11,9,83,88,80,60,66,18,54,67,79,26,96,6,20,63,1,33,110,113,106,119
\ No newline at end of file diff --git a/random20.compressed_lower_distance_matrix b/random20.compressed_lower_distance_matrix new file mode 100644 index 0000000..241408f --- /dev/null +++ b/random20.compressed_lower_distance_matrix @@ -0,0 +1 @@ +68,6,155,10,173,171,168,52,32,63,136,16,78,163,8,175,28,70,107,165,18,97,99,118,49,76,48,133,58,92,44,57,190,101,25,94,148,37,2,146,56,95,125,121,142,31,21,152,154,124,75,120,112,45,39,115,170,179,157,183,15,14,38,182,151,164,185,96,127,41,53,180,166,122,134,100,20,169,36,40,87,46,140,82,187,71,12,178,159,184,189,3,86,61,137,93,116,167,150,186,126,172,29,34,135,27,144,177,73,19,147,17,26,1,139,22,69,9,138,160,181,105,24,129,91,11,64,103,132,130,4,47,66,85,65,149,109,161,81,128,7,113,67,77,54,35,30,74,98,114,5,176,158,62,89,79,156,102,119,141,51,143,108,131,106,33,123,43,174,55,83,80,145,104,153,13,188,90,117,111,60,84,23,59,88,110,42,50,72,162
\ No newline at end of file @@ -1,25 +1,28 @@ -#include <vector> #include <iostream> -#include <iomanip> #include <fstream> -#include <iterator> -#include <string> #include <cassert> #include <queue> +#include <cmath> + +#ifdef USE_GOOGLE_HASHMAP + #include <sparsehash/sparse_hash_map> -#include "prettyprint.hpp" +template <class Key, class T> class hash_map : public google::sparse_hash_map<Key, T> { + public: inline void reserve(size_t hint) { this->resize(hint); } }; + +#else + +#include <unordered_map> +template <class Key, class T> class hash_map : public std::unordered_map<Key, T> {}; + +#endif typedef float value_t; -//typedef uint16_t value_t; +// typedef uint16_t value_t; typedef long index_t; typedef long coefficient_t; -//#define PRECOMPUTE_DIAMETERS -//#define PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION - -//#define STORE_DIAMETERS - #define USE_BINARY_SEARCH //#define USE_EXPONENTIAL_SEARCH @@ -34,1279 +37,821 @@ typedef long coefficient_t; //#define FILE_FORMAT_LOWER_TRIANGULAR_CSV class binomial_coeff_table { - std::vector<std::vector<index_t> > B; - index_t n_max, k_max; - + std::vector<std::vector<index_t>> B; + index_t n_max, k_max; + public: - binomial_coeff_table(index_t n, index_t k) { - n_max = n; - k_max = k; - - B.resize(n + 1); - for( index_t i = 0; i <= n; i++ ) { - B[i].resize(k + 1); - for ( index_t j = 0; j <= std::min(i, k); j++ ) - { - if (j == 0 || j == i) - B[i][j] = 1; - else - B[i][j] = B[i-1][j-1] + B[i-1][j]; - } - } - } - - inline index_t operator()(index_t n, index_t k) const { - assert(n <= n_max); - assert(k <= k_max); - return B[n][k]; - } + binomial_coeff_table(index_t n, index_t k) { + n_max = n; + k_max = k; + + B.resize(n + 1); + for (index_t i = 0; i <= n; i++) { + B[i].resize(k + 1); + for (index_t j = 0; j <= std::min(i, k); j++) { + if (j == 0 || j == i) + B[i][j] = 1; + else + B[i][j] = B[i - 1][j - 1] + B[i - 1][j]; + } + } + } + + inline index_t operator()(index_t n, index_t k) const { + assert(n <= n_max); + assert(k <= k_max); + return B[n][k]; + } }; -// -// https://comeoncodeon.wordpress.com/2011/10/09/modular-multiplicative-inverse/ -// - -// a * (m / a) + m % a = m -// m % a = -a * (m / a) (mod m) -//Dividing by (a * (m % a)): -// inverse(a) = - (m / a) * inverse(m % a) (mod m) - -std::vector<coefficient_t> multiplicative_inverse_vector (const coefficient_t m) { - std::vector<coefficient_t> mod_inverse(m); - mod_inverse[1] = 1; - for (coefficient_t i = 2; i < m; ++i) { - mod_inverse[i] = (-(m/i) * mod_inverse[m % i]) % m + m; - } - return mod_inverse; +std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) { + std::vector<coefficient_t> inverse(m); + inverse[1] = 1; + for (coefficient_t a = 2; a < m; ++a) { + // m = a * (m / a) + m % a + // Multipying with inverse(a) * inverse(m % a): + // 0 = (m / a) * inverse(m % a) + inverse(a) = 0 (mod m) + inverse[a] = m - ((m / a) * inverse[m % a]) % m; + } + return inverse; } -template<typename OutputIterator> -inline OutputIterator get_simplex_vertices( index_t idx, const index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out ) -{ - --n; - - for( index_t k = dim + 1; k > 0; --k ) { - - #ifdef USE_BINARY_SEARCH - if ( binomial_coeff( n , k ) > idx ) { - index_t count; - - #ifdef USE_EXPONENTIAL_SEARCH - for (count = 1; (binomial_coeff( n - count , k ) > idx); count = std::min(count << 1, n)); - #else - count = n; - #endif - - while (count > 0) { - index_t i = n; - index_t step = count >> 1; - i -= step; - if (binomial_coeff( i , k ) > idx) { - n = --i; - count -= step + 1; - } else count = step; - } - } - #else - while( binomial_coeff( n , k ) > idx ) - --n; - #endif - - assert( binomial_coeff( n , k ) <= idx ); - assert( binomial_coeff( n + 1, k ) > idx ); - - *out++ = n; - - idx -= binomial_coeff( n , k ); +template <typename OutputIterator> +inline OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n, + const binomial_coeff_table& binomial_coeff, + OutputIterator out) { + --n; - } + for (index_t k = dim + 1; k > 0; --k) { - return out; -} +#ifdef USE_BINARY_SEARCH + if (binomial_coeff(n, k) > idx) { + index_t count; -std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim, const index_t n, const binomial_coeff_table& binomial_coeff) { - std::vector<index_t> vertices; - get_simplex_vertices( simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices) ); - return vertices; -} +#ifdef USE_EXPONENTIAL_SEARCH + for (count = 1; (binomial_coeff(n - count, k) > idx); count = std::min(count << 1, n)) + ; +#else + count = n; +#endif + while (count > 0) { + index_t i = n; + index_t step = count >> 1; + i -= step; + if (binomial_coeff(i, k) > idx) { + n = --i; + count -= step + 1; + } else + count = step; + } + } +#else + while (binomial_coeff(n, k) > idx) --n; +#endif + assert(binomial_coeff(n, k) <= idx); + assert(binomial_coeff(n + 1, k) > idx); -#ifdef USE_COEFFICIENTS -struct entry_t { - index_t index; - coefficient_t coefficient; - - entry_t(index_t _index, coefficient_t _coefficient) : - index(_index), coefficient(_coefficient) {} - - entry_t(index_t _index) : - index(_index), coefficient(1) {} - - entry_t() : - index(0), coefficient(1) {} -}; + *out++ = n; -inline entry_t make_entry(index_t _index, coefficient_t _coefficient) { - return entry_t(_index, _coefficient); -} + idx -= binomial_coeff(n, k); + } -inline index_t get_index(entry_t e) { - return e.index; + return out; } -inline index_t get_coefficient(entry_t e) { - return e.coefficient; +std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim, + const index_t n, + const binomial_coeff_table& binomial_coeff) { + std::vector<index_t> vertices; + get_simplex_vertices(simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices)); + return vertices; } -inline void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; } +#ifdef USE_COEFFICIENTS +struct entry_t { + index_t index; + coefficient_t coefficient; + entry_t(index_t _index, coefficient_t _coefficient) + : index(_index), coefficient(_coefficient) {} + entry_t(index_t _index) : index(_index), coefficient(1) {} + entry_t() : index(0), coefficient(1) {} +}; +inline entry_t make_entry(index_t _index, coefficient_t _coefficient) { + return entry_t(_index, _coefficient); +} +inline index_t get_index(entry_t e) { return e.index; } +inline index_t get_coefficient(entry_t e) { return e.coefficient; } +inline void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; } -bool operator== (const entry_t& e1, const entry_t& e2) { - return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2); +bool operator==(const entry_t& e1, const entry_t& e2) { + return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2); } -std::ostream& operator<< (std::ostream& stream, const entry_t& e) { - stream << get_index(e) << ":" << get_coefficient(e); - return stream; +std::ostream& operator<<(std::ostream& stream, const entry_t& e) { + stream << get_index(e) << ":" << get_coefficient(e); + return stream; } #else typedef index_t entry_t; - -inline const index_t get_index(entry_t i) { - return i; -} - -inline index_t get_coefficient(entry_t i) { - return 1; -} - -inline entry_t make_entry(index_t _index, coefficient_t _value) { - return entry_t(_index); -} - +inline const index_t get_index(entry_t i) { return i; } +inline index_t get_coefficient(entry_t i) { return 1; } +inline entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); } inline void set_coefficient(index_t& e, const coefficient_t c) { e = c; } #endif inline const entry_t& get_entry(const entry_t& e) { return e; } -template <typename Entry> -struct greater_index { - bool operator() (const Entry& a, const Entry& b) { - return get_index(a) > get_index(b); - } +template <typename Entry> struct smaller_index { + bool operator()(const Entry& a, const Entry& b) { return get_index(a) < get_index(b); } }; - -#ifdef STORE_DIAMETERS typedef std::pair<value_t, index_t> diameter_index_t; -inline value_t get_diameter(diameter_index_t i) { - return i.first; -} -inline index_t get_index(diameter_index_t i) { - return i.second; -} -#else -typedef index_t diameter_index_t; -#endif +inline value_t get_diameter(diameter_index_t i) { return i.first; } +inline index_t get_index(diameter_index_t i) { return i.second; } -#ifdef STORE_DIAMETERS -class diameter_entry_t: public std::pair<value_t, entry_t> { +class diameter_entry_t : public std::pair<value_t, entry_t> { public: - diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {} -#ifdef USE_COEFFICIENTS - diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {} -#endif - diameter_entry_t(index_t i) : std::pair<value_t, entry_t>(0, i) {} + diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {} + diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {} + diameter_entry_t() : diameter_entry_t(0) {} }; inline const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } inline entry_t& get_entry(diameter_entry_t& p) { return p.second; } inline const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); } -inline const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); } +inline const coefficient_t get_coefficient(const diameter_entry_t& p) { + return get_coefficient(get_entry(p)); +} inline const value_t& get_diameter(const diameter_entry_t& p) { return p.first; } -inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } -inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) { - return std::make_pair(_diameter, make_entry(_index, _coefficient)); +inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) { + set_coefficient(get_entry(p), c); } -inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) { - return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient)); +inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, + coefficient_t _coefficient) { + return std::make_pair(_diameter, make_entry(_index, _coefficient)); } - - -struct greater_diameter_or_index { - inline bool operator() (const diameter_entry_t& a, const diameter_entry_t& b) { - return (get_diameter(a) > get_diameter(b)) || - ((get_diameter(a) == get_diameter(b)) && (get_index(a) > get_index(b))); - } -}; -#else -typedef entry_t diameter_entry_t; -inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) { - return make_entry(_index, _coefficient); +inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, + coefficient_t _coefficient) { + return std::make_pair(get_diameter(_diameter_index), + make_entry(get_index(_diameter_index), _coefficient)); } -inline diameter_entry_t make_diameter_entry(index_t _index, coefficient_t _coefficient) { - return make_entry(_index, _coefficient); -} -#endif - +template <typename Entry> struct greater_diameter_or_smaller_index { + inline bool operator()(const Entry& a, const Entry& b) { + return (get_diameter(a) > get_diameter(b)) || + ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b))); + } +}; -template <typename DistanceMatrix> -class rips_filtration_comparator { +template <typename DistanceMatrix> class rips_filtration_comparator { public: - const DistanceMatrix& dist; - - const index_t dim; + const DistanceMatrix& dist; + const index_t dim; private: - mutable std::vector<index_t> vertices; - - const binomial_coeff_table& binomial_coeff; - -public: - rips_filtration_comparator( - const DistanceMatrix& _dist, - const index_t _dim, - const binomial_coeff_table& _binomial_coeff - ): dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff) {}; - - inline value_t diameter(const index_t index) const { - value_t diam = 0; - get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin() ); - - for (index_t i = 0; i <= dim; ++i) - for (index_t j = 0; j < i; ++j) { - diam = std::max(diam, dist(vertices[i], vertices[j])); - } - return diam; - } - - inline bool operator()(const index_t a, const index_t b) const - { - assert(a < binomial_coeff(dist.size(), dim + 1)); - assert(b < binomial_coeff(dist.size(), dim + 1)); - - value_t a_diam = 0, b_diam = 0; - - b_diam = diameter(b); - - get_simplex_vertices(a, dim, dist.size(), binomial_coeff, vertices.begin() ); - for (index_t i = 0; i <= dim; ++i) - for (index_t j = i + 1; j <= dim; ++j) { - a_diam = std::max(a_diam, dist(vertices[i], vertices[j])); - if (a_diam > b_diam) - return true; - } - - return (a_diam == b_diam) && (a > b); - } - - template <typename Entry> - inline bool operator()(const Entry& a, const Entry& b) const - { - return operator()(get_index(a), get_index(b)); - } + mutable std::vector<index_t> vertices; + const binomial_coeff_table& binomial_coeff; +public: + rips_filtration_comparator(const DistanceMatrix& _dist, const index_t _dim, + const binomial_coeff_table& _binomial_coeff) + : dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff){}; + + inline value_t diameter(const index_t index) const { + value_t diam = 0; + get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin()); + + for (index_t i = 0; i <= dim; ++i) + for (index_t j = 0; j < i; ++j) { + diam = std::max(diam, dist(vertices[i], vertices[j])); + } + return diam; + } + + inline bool operator()(const index_t a, const index_t b) const { + assert(a < binomial_coeff(dist.size(), dim + 1)); + assert(b < binomial_coeff(dist.size(), dim + 1)); + + value_t a_diam = 0, b_diam = diameter(b); + + get_simplex_vertices(a, dim, dist.size(), binomial_coeff, vertices.begin()); + for (index_t i = 0; i <= dim; ++i) + for (index_t j = i + 1; j <= dim; ++j) { + a_diam = std::max(a_diam, dist(vertices[i], vertices[j])); + if (a_diam > b_diam) return true; + } + + return (a_diam == b_diam) && (a < b); + } + + template <typename Entry> inline bool operator()(const Entry& a, const Entry& b) const { + return operator()(get_index(a), get_index(b)); + } }; - class simplex_coboundary_enumerator { private: - index_t idx, modified_idx, dim, n, k; - - const binomial_coeff_table& binomial_coeff; - + index_t idx, modified_idx, dim, v, k; + const binomial_coeff_table& binomial_coeff; + public: - inline simplex_coboundary_enumerator ( - index_t _idx, index_t _dim, index_t _n, - const binomial_coeff_table& _binomial_coeff): - idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff) {} - - inline bool has_next() - { - while ( (k != -1 && n != -1) && (binomial_coeff( n , k ) <= idx) ) { - idx -= binomial_coeff( n , k ); - - modified_idx -= binomial_coeff( n , k ); - modified_idx += binomial_coeff( n , k + 1 ); - - --n; + inline simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n, + const binomial_coeff_table& _binomial_coeff) + : idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), v(_n - 1), + binomial_coeff(_binomial_coeff) {} + + inline bool has_next() { + while ((v != -1) && (binomial_coeff(v, k) <= idx)) { + idx -= binomial_coeff(v, k); + modified_idx += binomial_coeff(v, k + 1) - binomial_coeff(v, k); + --v; --k; - } - return k != -1 && n != -1; - } - - inline std::pair<entry_t, index_t> next() - { - while ( binomial_coeff( n , k ) <= idx ) { - idx -= binomial_coeff( n , k ); - - modified_idx -= binomial_coeff( n , k ); - modified_idx += binomial_coeff( n , k + 1 ); - - --n; - } - auto result = std::make_pair(make_entry( - modified_idx + binomial_coeff( n , k + 1 ), - k & 1 ? -1 : 1), - n); - - --n; - return result; - } + assert(k != -1); + } + return v != -1; + } + + inline std::pair<entry_t, index_t> next() { + auto result = + std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v); + --v; + return result; + } }; -template <typename DistanceMatrix> -std::vector<value_t> get_diameters ( - const DistanceMatrix& dist, const index_t dim, - const std::vector<value_t>& previous_diameters, - const binomial_coeff_table& binomial_coeff - ) -{ - index_t n = dist.size(); - - std::vector<value_t> diameters(binomial_coeff(n, dim + 1), 0); - - std::vector<index_t> coboundary; - - for (index_t simplex = 0; simplex < previous_diameters.size(); ++simplex) { - coboundary.clear(); - - #ifdef INDICATE_PROGRESS - std::cout << "\033[Kpropagating diameter of simplex " << simplex + 1 << "/" << previous_diameters.size() << std::flush << "\r"; - #endif - - simplex_coboundary_enumerator coboundaries(simplex, dim - 1, n, binomial_coeff); - while (coboundaries.has_next()) { - index_t coface = get_index(coboundaries.next().first); - diameters[coface] = std::max( diameters[coface], previous_diameters[simplex]); - } - } - - #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; - #endif - - return diameters; -} - - -class rips_filtration_diameter_comparator { -private: - const std::vector<value_t>& diameters; - - const index_t dim; - -public: - std::vector<index_t> vertices; - - typedef value_t dist_t; - - const binomial_coeff_table& binomial_coeff; - +class distance_matrix { public: - rips_filtration_diameter_comparator( - const std::vector<value_t>& _diameters, - const index_t _dim, - const binomial_coeff_table& _binomial_coeff - ): - diameters(_diameters), dim(_dim), - binomial_coeff(_binomial_coeff) {} - - inline value_t diameter(const index_t a) const { - assert(a < diameters.size()); - return diameters[a]; - } - - inline bool operator()(const index_t a, const index_t b) const - { - assert(a < diameters.size()); - assert(b < diameters.size()); - - dist_t a_diam = diameters[a], b_diam = diameters[b]; - - return ((a_diam > b_diam) || ((a_diam == b_diam) && (a > b))); - } - - template <typename Entry> - inline bool operator()(const Entry& a, const Entry& b) const - { - return operator()(get_index(a), get_index(b)); - } + typedef value_t value_type; + std::vector<std::vector<value_t>> distances; + inline value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; } + + inline size_t size() const { return distances.size(); } }; +enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; -class distance_matrix { +template <compressed_matrix_layout Layout> class compressed_distance_matrix_adapter { public: - typedef value_t value_type; - std::vector<std::vector<value_t> > distances; - inline value_t operator()(const index_t a, const index_t b) const { - return distances[a][b]; - } - - inline size_t size() const { - return distances.size(); - } -}; + typedef value_t value_type; + std::vector<value_t>& distances; + std::vector<value_t*> rows; + void init_distances() { distances.resize(size() * (size() - 1) / 2); } -class compressed_upper_distance_matrix_adapter { -public: - typedef value_t value_type; - std::vector<value_t>& distances; - std::vector<value_t*> rows; - - index_t n; - - void init_distances () { - distances.resize(n * (n-1) / 2); - } - - void init_rows () { - rows.resize(n); - value_t* pointer = &distances[0] - 1; - for (index_t i = 0; i < n - 1; ++i) { - rows[i] = pointer; - pointer += n - i - 2; - } - } - - compressed_upper_distance_matrix_adapter(std::vector<value_t>& _distances) : - distances(_distances) - { - n = (1 + std::sqrt(1+ 8 * _distances.size())) / 2; - - assert( distances.size() == n * (n-1) / 2 ); - - init_rows(); - } - - inline value_t operator()(const index_t a, const index_t b) const { - if (a < b) - return rows[a][b]; - else if (a > b) - return rows[b][a]; - else - return 0; - } - - inline size_t size() const { - return n; - } -}; + void init_rows(); + compressed_distance_matrix_adapter(std::vector<value_t>& _distances) + : distances(_distances), rows((1 + std::sqrt(1 + 8 * _distances.size())) / 2) { + assert(distances.size() == size() * (size() - 1) / 2); + init_rows(); + } -class compressed_lower_distance_matrix_adapter { -public: - typedef value_t value_type; - std::vector<value_t>& distances; - std::vector<value_t*> rows; - - index_t n; - - void init_distances () { - distances.resize(n * (n-1) / 2); - } - - void init_rows () { - rows.resize(n); - value_t* pointer = &distances[0]; - for (index_t i = 1; i < n; ++i) { - rows[i] = pointer; - pointer += i; - } - } - - compressed_lower_distance_matrix_adapter( - std::vector<value_t>& _distances, const index_t _n) : - distances(_distances), n(_n) - { - init_distances(); - init_rows(); - } - - compressed_lower_distance_matrix_adapter(std::vector<value_t>& _distances) : - distances(_distances) - { - n = (1 + std::sqrt(1+ 8 * distances.size())) / 2; - assert( distances.size() == n * (n-1) / 2 ); - - init_rows(); - } - - template <typename DistanceMatrix> - compressed_lower_distance_matrix_adapter( - std::vector<value_t>& _distances, const DistanceMatrix& mat) : - distances(_distances), n(mat.size()) { - init_distances(); - init_rows(); - - for (index_t i = 1; i < n; ++i) { - for (index_t j = 0; j < i; ++j) { - rows[i][j] = mat(i, j); - } - } - } - - inline value_t operator()(const index_t i, const index_t j) const { - if (i > j) - return rows[i][j]; - else if (i < j) - return rows[j][i]; - else - return 0; - } - - inline size_t size() const { - return n; - } -}; + template <typename DistanceMatrix> + compressed_distance_matrix_adapter(std::vector<value_t>& _distances, const DistanceMatrix& mat) + : distances(_distances), rows(mat.size()) { + init_distances(); + init_rows(); + for (index_t i = 1; i < size(); ++i) { + for (index_t j = 0; j < i; ++j) { rows[i][j] = mat(i, j); } + } + } + inline value_t operator()(const index_t i, const index_t j) const; + inline size_t size() const { return rows.size(); } +}; +template <> void compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::init_rows() { + value_t* pointer = &distances[0]; + for (index_t i = 1; i < size(); ++i) { + rows[i] = pointer; + pointer += i; + } +} -#ifdef USE_COEFFICIENTS -template <typename Heap> -inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) -{ - - if( column.empty() ) - return -1; - else { - auto pivot = column.top(); - - coefficient_t coefficient = 0; - do { - coefficient = (coefficient + get_coefficient(column.top())) % modulus; - column.pop(); - - if( coefficient == 0 ) { - if (column.empty()) { - return -1; - } else { - pivot = column.top(); - } - } - } while ( !column.empty() && get_index(column.top()) == get_index(pivot) ); - if( get_index(pivot) != -1 ) { - set_coefficient(pivot, coefficient); - } - return pivot; - } +template <> void compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::init_rows() { + value_t* pointer = &distances[0] - 1; + for (index_t i = 0; i < size() - 1; ++i) { + rows[i] = pointer; + pointer += size() - i - 2; + } } -#else +template <> +inline value_t compressed_distance_matrix_adapter<UPPER_TRIANGULAR>:: +operator()(const index_t i, const index_t j) const { + if (i < j) + return rows[i][j]; + else if (i > j) + return rows[j][i]; + else + return 0; +} -template <typename Heap> -inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) -{ - - if( column.empty() ) - return -1; - else { - auto pivot = column.top(); - column.pop(); - while( !column.empty() && get_index(column.top()) == get_index(pivot) ) { - column.pop(); - if( column.empty() ) - return -1; - else { - pivot = column.top(); - column.pop(); - } - } - return pivot; - } +template <> +inline value_t compressed_distance_matrix_adapter<LOWER_TRIANGULAR>:: +operator()(const index_t i, const index_t j) const { + if (i > j) + return rows[i][j]; + else if (i < j) + return rows[j][i]; + else + return 0; } -#endif +typedef compressed_distance_matrix_adapter<LOWER_TRIANGULAR> + compressed_lower_distance_matrix_adapter; +typedef compressed_distance_matrix_adapter<UPPER_TRIANGULAR> + compressed_upper_distance_matrix_adapter; -template <typename Heap> -inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) -{ - diameter_entry_t result = pop_pivot(column, modulus); - if( get_index(result) != -1 ) { - column.push(result); - } - return result; + +template <typename Heap> inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { + if (column.empty()) + return diameter_entry_t(-1); + else { + auto pivot = column.top(); + +#ifdef USE_COEFFICIENTS + coefficient_t coefficient = 0; + do { + coefficient = (coefficient + get_coefficient(column.top())) % modulus; + column.pop(); + + if (coefficient == 0) { + if (column.empty()) { + return diameter_entry_t(-1); + } else { + pivot = column.top(); + } + } + } while (!column.empty() && get_index(column.top()) == get_index(pivot)); + if (get_index(pivot) != -1) { set_coefficient(pivot, coefficient); } +#else + column.pop(); + while (!column.empty() && get_index(column.top()) == get_index(pivot)) { + column.pop(); + if (column.empty()) + return diameter_entry_t(-1); + else { + pivot = column.top(); + column.pop(); + } + } +#endif + return pivot; + } } +template <typename Heap> inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) { + diameter_entry_t result = pop_pivot(column, modulus); + if (get_index(result) != -1) { column.push(result); } + return result; +} template <typename Comparator> -void assemble_columns_to_reduce ( - std::vector<diameter_index_t>& columns_to_reduce, - google::sparse_hash_map<index_t, index_t>& pivot_column_index, - const Comparator& comp, - index_t dim, index_t n, - value_t threshold, - const binomial_coeff_table& binomial_coeff -) { - index_t num_simplices = binomial_coeff(n, dim + 2); - - columns_to_reduce.clear(); - - for (index_t index = 0; index < num_simplices; ++index) { - if (pivot_column_index.find(index) == pivot_column_index.end()) { - value_t diameter = comp.diameter(index); - if (diameter <= threshold) - #ifdef STORE_DIAMETERS - columns_to_reduce.push_back(std::make_pair(diameter, index)); - #else - columns_to_reduce.push_back(index); - #endif - - } - } - - #ifdef STORE_DIAMETERS - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), std::greater<diameter_index_t>()); - #else - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp); - #endif +void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce, + hash_map<index_t, index_t>& pivot_column_index, + const Comparator& comp, index_t dim, index_t n, value_t threshold, + const binomial_coeff_table& binomial_coeff) { + index_t num_simplices = binomial_coeff(n, dim + 2); + + columns_to_reduce.clear(); + + for (index_t index = 0; index < num_simplices; ++index) { + if (pivot_column_index.find(index) == pivot_column_index.end()) { + value_t diameter = comp.diameter(index); + if (diameter <= threshold) + columns_to_reduce.push_back(std::make_pair(diameter, index)); + } + } + + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index<diameter_index_t>()); } - -template <typename ValueType> -class compressed_sparse_matrix { +template <typename ValueType> class compressed_sparse_matrix { public: - std::vector<size_t> bounds; - std::vector<ValueType> entries; - - - inline size_t size() const { - return bounds.size(); - } - - inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const { - assert(index < size()); - return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; - } - - inline typename std::vector<ValueType>::const_iterator cend(size_t index) const { - assert(index < size()); - return entries.cbegin() + bounds[index]; - } - - template <typename Iterator> - inline void append(Iterator begin, Iterator end) { - for (Iterator it = begin; it != end; ++it) { - entries.push_back(*it); - } - bounds.push_back(entries.size()); - } - - inline void append() { - bounds.push_back(entries.size()); - } - - inline void push_back(ValueType e) { - assert(0 < size()); - entries.push_back(e); - ++bounds.back(); - } - - inline void pop_back() { - assert(0 < size()); - entries.pop_back(); - --bounds.back(); - } - - template <typename Collection> - inline void append(const Collection collection) { - append(collection.cbegin(), collection.cend()); - } - + std::vector<size_t> bounds; + std::vector<ValueType> entries; + + inline size_t size() const { return bounds.size(); } + + inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const { + assert(index < size()); + return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; + } + + inline typename std::vector<ValueType>::const_iterator cend(size_t index) const { + assert(index < size()); + return entries.cbegin() + bounds[index]; + } + + template <typename Iterator> inline void append(Iterator begin, Iterator end) { + for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); } + bounds.push_back(entries.size()); + } + + inline void append() { bounds.push_back(entries.size()); } + + inline void push_back(ValueType e) { + assert(0 < size()); + entries.push_back(e); + ++bounds.back(); + } + + inline void pop_back() { + assert(0 < size()); + entries.pop_back(); + --bounds.back(); + } + + template <typename Collection> inline void append(const Collection collection) { + append(collection.cbegin(), collection.cend()); + } }; - template <typename Heap> -inline std::vector<entry_t> get_column_vector(Heap column, coefficient_t modulus) -{ - std::vector<entry_t> temp_col; - entry_t pivot = pop_pivot( column, modulus ); - while( get_index(pivot) != -1 ) { - temp_col.push_back( pivot ); - pivot = pop_pivot( column, modulus ); - } - return temp_col; +inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) { + entry_t e = make_entry(i, c); + column.push(std::make_pair(diameter, e)); } +template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator> +void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, + hash_map<index_t, index_t>& pivot_column_index, + const DistanceMatrix& dist, const ComparatorCofaces& comp, + const Comparator& comp_prev, index_t dim, index_t n, value_t threshold, + coefficient_t modulus, const std::vector<coefficient_t>& multiplicative_inverse, + const binomial_coeff_table& binomial_coeff) { + +#ifdef PRINT_PERSISTENCE_PAIRS + std::cout << "persistence intervals in dim " << dim << ":" << std::endl; +#endif -template <typename Heap> -inline std::vector<entry_t> get_heap_vector(Heap heap) -{ - std::vector<entry_t> temp_col; - while( !heap.empty() ) { - temp_col.push_back( heap.top() ); - heap.pop(); - } - return temp_col; -} +#ifdef ASSEMBLE_REDUCTION_MATRIX + compressed_sparse_matrix<diameter_entry_t> reduction_matrix; +#else +#ifdef USE_COEFFICIENTS + std::vector<diameter_entry_t> reduction_coefficients; +#endif +#endif + std::vector<diameter_entry_t> coface_entries; + std::vector<index_t> vertices; -template <typename Heap> -inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) { - entry_t e = make_entry(i, c); -#ifdef STORE_DIAMETERS - column.push(std::make_pair(diameter, e)); -#else - column.push(e); + for (index_t i = 0; i < columns_to_reduce.size(); ++i) { + auto column_to_reduce = columns_to_reduce[i]; + +#ifdef ASSEMBLE_REDUCTION_MATRIX + std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, + smaller_index<diameter_entry_t>> reduction_column; #endif -} + std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, + greater_diameter_or_smaller_index<diameter_entry_t>> working_coboundary; -template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator> -void compute_pairs( - std::vector<diameter_index_t>& columns_to_reduce, - google::sparse_hash_map<index_t, index_t>& pivot_column_index, - const DistanceMatrix& dist, - const ComparatorCofaces& comp, const Comparator& comp_prev, - index_t dim, index_t n, - value_t threshold, - coefficient_t modulus, const std::vector<coefficient_t>& multiplicative_inverse, - const binomial_coeff_table& binomial_coeff -) { - - #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; - #endif - - #ifdef ASSEMBLE_REDUCTION_MATRIX - compressed_sparse_matrix <diameter_entry_t> reduction_matrix; - #else - #ifdef USE_COEFFICIENTS - std::vector <diameter_entry_t> reduction_coefficients; - #endif - #endif - -// size_t boundary_additions = 0; - - for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - auto column_to_reduce = columns_to_reduce[i]; - - #ifdef ASSEMBLE_REDUCTION_MATRIX - - std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_index<diameter_entry_t>> reduction_column; - - #endif - - - #ifdef STORE_DIAMETERS - std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_diameter_or_index > working_coboundary; - #else - std::priority_queue<entry_t, std::vector<entry_t>, decltype(comp) > working_coboundary(comp); - #endif - - #ifdef STORE_DIAMETERS - value_t diameter = get_diameter(columns_to_reduce[i]); - #else - value_t diameter = comp_prev.diameter(get_index(column_to_reduce)); - #endif + value_t diameter = get_diameter(columns_to_reduce[i]); - - #ifdef INDICATE_PROGRESS - std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size() - << " (diameter " << diameter << ")" - << std::flush << "\r"; - #endif - - index_t j = i; - - auto column_to_add = column_to_reduce; - - - // start with a dummy pivot entry with coefficient -1 in order to initialize working_coboundary - // with the coboundary of the simplex with index column_to_reduce - - diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1); - -// std::cout << "reducing " << column_to_reduce << ":" << std::endl; - - #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_matrix.append(); - #endif - +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " + << diameter << ")" << std::flush << "\r"; +#endif - // initialize reduction_matrix as identity matrix - #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1)); - #else - #ifdef USE_COEFFICIENTS - reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1)); - #endif - #endif - -// --boundary_additions; - - do { - - const coefficient_t factor = modulus - get_coefficient(pivot); + index_t j = i; + auto column_to_add = column_to_reduce; + // start with a dummy pivot entry with coefficient -1 in order to initialize + // working_coboundary with the coboundary of the simplex with index column_to_reduce + diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1); -// std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, decltype(comp) > eliminating_coboundary(comp); +#ifdef ASSEMBLE_REDUCTION_MATRIX + reduction_matrix.append(); +#endif -// std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl; +#ifdef ASSEMBLE_REDUCTION_MATRIX + // initialize reduction_matrix as identity matrix + reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1)); +#else +#ifdef USE_COEFFICIENTS + reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1)); +#endif +#endif - #ifdef ASSEMBLE_REDUCTION_MATRIX - for (auto it = reduction_matrix.cbegin(j); it != reduction_matrix.cend(j); ++it) - #endif - { - -// ++boundary_additions; - - - #ifdef ASSEMBLE_REDUCTION_MATRIX - auto simplex = *it; - coefficient_t simplex_coefficient = get_coefficient(simplex); - simplex_coefficient *= factor; - simplex_coefficient %= modulus; - #else - #ifdef USE_COEFFICIENTS - const auto& simplex = reduction_coefficients[j]; - #else - const auto simplex = column_to_add; - #endif - #endif - - #ifdef STORE_DIAMETERS - value_t simplex_diameter = get_diameter(simplex); - #else - value_t simplex_diameter = comp_prev.diameter(get_index(simplex)); - #endif - assert(simplex_diameter == comp_prev.diameter(get_index(simplex))); - - #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_column.push( make_diameter_entry( simplex_diameter, get_index(simplex), simplex_coefficient ) ); - #endif - - std::vector<index_t> vertices = vertices_of_simplex(get_index(simplex), dim, n, binomial_coeff); + bool might_be_apparent_pair = true; + + do { + const coefficient_t factor = modulus - get_coefficient(pivot); + +#ifdef ASSEMBLE_REDUCTION_MATRIX + for (auto it = reduction_matrix.cbegin(j); it != reduction_matrix.cend(j); ++it) +#endif + { +#ifdef ASSEMBLE_REDUCTION_MATRIX + auto simplex = *it; + coefficient_t simplex_coefficient = get_coefficient(simplex); + simplex_coefficient *= factor; + simplex_coefficient %= modulus; +#else +#ifdef USE_COEFFICIENTS + const auto& simplex = reduction_coefficients[j]; +#else + const auto simplex = column_to_add; +#endif +#endif + value_t simplex_diameter = get_diameter(simplex); + assert(simplex_diameter == comp_prev.diameter(get_index(simplex))); + +#ifdef ASSEMBLE_REDUCTION_MATRIX + reduction_column.push( + make_diameter_entry(simplex_diameter, get_index(simplex), simplex_coefficient)); +#endif + + vertices.clear(); + get_simplex_vertices(get_index(simplex), dim, n, binomial_coeff, std::back_inserter(vertices)); - simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff); - while (cofaces.has_next()) { - auto coface_descriptor = cofaces.next(); - - entry_t coface = coface_descriptor.first; - index_t covertex = coface_descriptor.second; - - index_t coface_index = get_index(coface); - - value_t coface_diameter = simplex_diameter; - for (index_t v : vertices) { - coface_diameter = std::max(coface_diameter, dist(v, covertex)); - } - - assert(comp.diameter(coface_index) == coface_diameter); - - if (coface_diameter <= threshold) { - coefficient_t coface_coefficient = get_coefficient(coface) * get_coefficient(simplex) * factor; - coface_coefficient %= modulus; - if (coface_coefficient < 0) coface_coefficient += modulus; + coface_entries.clear(); + simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff); + while (cofaces.has_next()) { + auto coface_descriptor = cofaces.next(); + entry_t coface = coface_descriptor.first; + index_t covertex = coface_descriptor.second; + index_t coface_index = get_index(coface); + value_t coface_diameter = simplex_diameter; + for (index_t v : vertices) { + coface_diameter = std::max(coface_diameter, dist(v, covertex)); + } + assert(comp.diameter(coface_index) == coface_diameter); + + if (coface_diameter <= threshold) { + coefficient_t coface_coefficient = + get_coefficient(coface) * get_coefficient(simplex) * factor; + coface_coefficient %= modulus; + if (coface_coefficient < 0) coface_coefficient += modulus; + assert(coface_coefficient >= 0); - assert(coface_coefficient >= 0); + diameter_entry_t coface_entry = make_diameter_entry(coface_diameter, coface_index, coface_coefficient); + coface_entries.push_back(coface_entry); - push_entry(working_coboundary, coface_index, coface_coefficient, coface_diameter); -// push_entry(eliminating_coboundary, coface_index, coface_coefficient, coface_diameter); - } - } - } + if (might_be_apparent_pair && (simplex_diameter == coface_diameter)) { + if (pivot_column_index.find(coface_index) == pivot_column_index.end()) + { + pivot = coface_entry; + goto found_pivot; + } + might_be_apparent_pair = false; + } + } + } + for (auto e: coface_entries) working_coboundary.push(e); + } + pivot = get_pivot(working_coboundary, modulus); + found_pivot: + if (get_index(pivot) != -1) { + auto pair = pivot_column_index.find(get_index(pivot)); - -// std::cout << get_heap_vector(working_coboundary) << std::endl; + if (pair == pivot_column_index.end()) { + pivot_column_index.insert(std::make_pair(get_index(pivot), i)); -// std::cout << "e:" << get_column_vector(eliminating_coboundary, modulus) << std::endl; -// std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl << std::endl; - - pivot = get_pivot(working_coboundary, modulus); - -// std::cout << "pivot " << get_index(pivot) << std::endl; - - if (get_index(pivot) != -1) { - auto pair = pivot_column_index.find(get_index(pivot)); - - if (pair == pivot_column_index.end()) { -// std::cout << std::endl; - - pivot_column_index.insert(std::make_pair(get_index(pivot), i)); - - #ifdef USE_COEFFICIENTS - const coefficient_t inverse = multiplicative_inverse[ get_coefficient( pivot ) ]; - #else - #ifdef ASSEMBLE_REDUCTION_MATRIX - const coefficient_t inverse = 1; - #endif - #endif - - // replace current column of reduction_matrix (with a single diagonal 1 entry) - // by reduction_column (possibly with a different entry on the diagonal) - #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_matrix.pop_back(); - while (true) { - diameter_entry_t e = pop_pivot(reduction_column, modulus); - index_t index = get_index(e); - if (index == -1) - break; - const coefficient_t coefficient = inverse * get_coefficient(e) % modulus; - assert(coefficient > 0); - reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient)); - } - #else - #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); - reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse)); - #endif - #endif - - #ifdef PRINT_PERSISTENCE_PAIRS - #ifdef STORE_DIAMETERS - value_t death = get_diameter(pivot); - #else - value_t death = comp.diameter(get_index(pivot)); - #endif - if (diameter != death) - { - #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; - #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; - -// std::cout << " (" << vertices_of_simplex(column_to_reduce, dim, n, binomial_coeff) << ", " << vertices_of_simplex(get_index(pivot), dim + 1, n, binomial_coeff) << ")" << std::endl; - } - #endif +#ifdef USE_COEFFICIENTS + const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; +#endif - break; - } - - j = pair->second; - column_to_add = columns_to_reduce[j]; - } +#ifdef ASSEMBLE_REDUCTION_MATRIX + // replace current column of reduction_matrix (with a single diagonal 1 entry) + // by reduction_column (possibly with a different entry on the diagonal) + reduction_matrix.pop_back(); + while (true) { + diameter_entry_t e = pop_pivot(reduction_column, modulus); + index_t index = get_index(e); + if (index == -1) break; +#ifdef USE_COEFFICIENTS + const coefficient_t coefficient = inverse * get_coefficient(e) % modulus; + assert(coefficient > 0); +#else + const coefficient_t coefficient = 1; +#endif + reduction_matrix.push_back( + make_diameter_entry(get_diameter(e), index, coefficient)); + } +#else +#ifdef USE_COEFFICIENTS + reduction_coefficients.pop_back(); + reduction_coefficients.push_back( + make_diameter_entry(column_to_reduce, inverse)); +#endif +#endif - } while ( get_index(pivot) != -1 ); - - #ifdef PRINT_PERSISTENCE_PAIRS - if ( get_index(pivot) == -1 ) { -// std::cout << std::endl; - - value_t birth = diameter; - #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; - #endif - std::cout << " [" << birth << ", )" << std::endl << std::flush; - } - #endif - - -// #ifdef ASSEMBLE_REDUCTION_MATRIX -// std::cout << "reduction matrix fill-in: " << i + 1 << " + " << reduction_matrix.entries.size() - (i + 1) << std::endl; -// #endif - - } - - #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; - #endif - -// std::cout << boundary_additions << " boundary additions" << std::endl; +#ifdef PRINT_PERSISTENCE_PAIRS + value_t death = get_diameter(pivot); + if (diameter != death) { +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif + std::cout << " [" << diameter << "," << death << ")" << std::endl + << std::flush; + } +#endif + break; + } + j = pair->second; + column_to_add = columns_to_reduce[j]; + } + } while (get_index(pivot) != -1); + +#ifdef PRINT_PERSISTENCE_PAIRS + if (get_index(pivot) == -1) { + value_t birth = diameter; +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif + std::cout << " [" << birth << ", )" << std::endl << std::flush; + } +#endif + } +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif } -bool is_prime(const long n) { - bool is_prime = true; - for (int i = 2; i <= n/2; ++i) - if (n%i == 0) { - is_prime = false; - break; - } - return is_prime; +bool is_prime(const coefficient_t n) { + if (!(n & 1) || n < 2) return n == 2; + for (coefficient_t p = 3, q = n / p, r = n % p; p <= q; p += 2, q = n / p, r = n % p) + if (!r) return false; + return true; } -void print_usage_and_exit(int exit_code) -{ - std::cerr << "Usage: " << "ripser " << "[options] filename" << std::endl; - std::cerr << std::endl; - std::cerr << "Options:" << std::endl; - std::cerr << std::endl; - std::cerr << " --help print this screen" << std::endl; - std::cerr << " --top_dim <k> compute persistent homology up to dimension <k>" << std::endl; - std::cerr << " --threshold <t> compute Rips complexes up to diameter <t>" << std::endl; - #ifdef USE_COEFFICIENTS - std::cerr << " --modulus <p> compute homology with coefficients in the prime field Z/<p>Z" << std::endl; - #endif - - exit(exit_code); -} +void print_usage_and_exit(int exit_code) { + std::cerr << "Usage: " + << "ripser " + << "[options] filename" << std::endl; + std::cerr << std::endl; + std::cerr << "Options:" << std::endl; + std::cerr << std::endl; + std::cerr << " --help print this screen" << std::endl; + std::cerr << " --top_dim <k> compute persistent homology up to dimension <k>" << std::endl; + std::cerr << " --threshold <t> compute Rips complexes up to diameter <t>" << std::endl; +#ifdef USE_COEFFICIENTS + std::cerr << " --modulus <p> compute homology with coefficients in the prime field Z/<p>Z" + << std::endl; +#endif -int main( int argc, char** argv ) { - - if( argc < 2 ) print_usage_and_exit(-1); - - const char *filename = nullptr; - - index_t dim_max = 1; - value_t threshold = std::numeric_limits<value_t>::max(); - - #ifdef USE_COEFFICIENTS - coefficient_t modulus = 2; - #else - const coefficient_t modulus = 2; - #endif - - for( index_t i = 1; i < argc; ++i ) { - const std::string arg(argv[ i ]); - if( arg == "--help" ) { - print_usage_and_exit(0); - } else if( arg == "--top_dim" ) { - std::string parameter = std::string( argv[ ++i ] ); - size_t next_pos; - dim_max = std::stol( parameter, &next_pos ); - if( next_pos != parameter.size() ) - print_usage_and_exit( -1 ); - } else if( arg == "--threshold" ) { - std::string parameter = std::string( argv[ ++i ] ); - size_t next_pos; - threshold = std::stof( parameter, &next_pos ); - if( next_pos != parameter.size() ) - print_usage_and_exit( -1 ); - #ifdef USE_COEFFICIENTS - } else if( arg == "--modulus" ) { - std::string parameter = std::string( argv[ ++i ] ); - size_t next_pos; - modulus = std::stol( parameter, &next_pos ); - if( next_pos != parameter.size() || !is_prime(modulus) ) - print_usage_and_exit( -1 ); - #endif - } else { - if (filename) { - print_usage_and_exit( -1 ); - } - filename = argv[i]; - } - } - - - std::ifstream input_stream( filename, std::ios_base::binary | std::ios_base::in ); - if( input_stream.fail( ) ) { - std::cerr << "couldn't open file " << filename << std::endl; - exit(-1); - } - - std::vector<std::vector<value_t>> diameters(dim_max + 2); - - #ifdef FILE_FORMAT_DIPHA - - int64_t magic_number; - input_stream.read( reinterpret_cast<char*>(&magic_number), sizeof( int64_t ) ); - if( magic_number != 8067171840 ) { - std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; - exit(-1); - } - - int64_t file_type; - input_stream.read( reinterpret_cast<char*>(&file_type), sizeof( int64_t ) ); - if( file_type != 7 ) { - std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; - exit(-1); - } - - int64_t n; - input_stream.read( reinterpret_cast<char*>(&n), sizeof( int64_t ) ); - - distance_matrix dist_full; - dist_full.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n)); - - for( int i = 0; i < n; ++i ) { - for( int j = 0; j < n; ++j ) { - double val; - input_stream.read( reinterpret_cast<char*>(&val), sizeof(double) ); - dist_full.distances[i][j] = val; - } - } - - std::cout << "distance matrix with " << n << " points" << std::endl; - - compressed_lower_distance_matrix_adapter dist(diameters[1], dist_full); - - std::cout << "distance matrix transformed to lower triangular form" << std::endl; - - #endif - - - #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV - - std::vector<value_t> distances; - std::string value_string; - while(std::getline(input_stream, value_string, ',')) - distances.push_back(std::stod(value_string)); - - compressed_upper_distance_matrix_adapter dist_upper(distances); - - index_t n = dist_upper.size(); - - std::cout << "distance matrix with " << n << " points" << std::endl; - - compressed_lower_distance_matrix_adapter dist(diameters[1], dist_upper); - - std::cout << "distance matrix transformed to lower triangular form" << std::endl; - - #endif - - - #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV - - std::vector<value_t>& distances = diameters[1]; - std::string value_string; - while(std::getline(input_stream, value_string, ',')) - distances.push_back(std::stod(value_string)); - - compressed_lower_distance_matrix_adapter dist(distances); - - index_t n = dist.size(); - - std::cout << "distance matrix with " << n << " points" << std::endl; - - #endif - - - - auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); - std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; - - - assert(dim_max <= n - 2); - - binomial_coeff_table binomial_coeff(n, dim_max + 2); - - std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus)); - - std::vector<diameter_index_t> columns_to_reduce; - - - for (index_t index = n; index-- > 0; ) { - #ifdef STORE_DIAMETERS - columns_to_reduce.push_back(std::make_pair(0, index)); - #else - columns_to_reduce.push_back(index); - #endif - - } - - - - index_t dim = 0; - - { - rips_filtration_diameter_comparator comp(diameters[1], dim + 1, binomial_coeff); - rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff); - - google::sparse_hash_map<index_t, index_t> pivot_column_index; - - compute_pairs( - columns_to_reduce, - pivot_column_index, - dist, - comp, comp_prev, - dim, n, - threshold, - modulus, multiplicative_inverse, - binomial_coeff - ); - - assemble_columns_to_reduce( - columns_to_reduce, - pivot_column_index, - comp, - dim, n, - threshold, - binomial_coeff - ); - } - - #ifdef PRECOMPUTE_DIAMETERS - #ifdef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION - for (dim = 1; dim <= dim_max; ++dim) { - #else - for (dim = 1; dim < dim_max; ++dim) { - #endif - #else - for (dim = 1; dim <= dim_max; ++dim) { - #endif - - #ifdef PRECOMPUTE_DIAMETERS - - #ifdef INDICATE_PROGRESS - std::cout << "precomputing " << dim + 1 << "-simplex diameters" << std::endl; - #endif - diameters[dim + 1] = get_diameters( dist, dim + 1, diameters[dim], binomial_coeff ); - - rips_filtration_diameter_comparator comp(diameters[dim + 1], dim + 1, binomial_coeff); - rips_filtration_diameter_comparator comp_prev(diameters[dim], dim, binomial_coeff); - - #else - - rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff); - rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff); - - #endif - - google::sparse_hash_map<index_t, index_t> pivot_column_index; - pivot_column_index.resize(columns_to_reduce.size()); - - compute_pairs( - columns_to_reduce, - pivot_column_index, - dist, - comp, comp_prev, - dim, n, - threshold, - modulus, - multiplicative_inverse, - binomial_coeff - ); - - if (dim < dim_max) - assemble_columns_to_reduce( - columns_to_reduce, - pivot_column_index, - comp, - dim, n, - threshold, - binomial_coeff - ); - -// if ( dim > 1 ) -// diameters[dim] = std::vector<value_t>(); - } - - #ifdef PRECOMPUTE_DIAMETERS - #ifndef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION - { - dim = dim_max; - - rips_filtration_diameter_comparator comp_prev(diameters[dim], dim, binomial_coeff); - rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff); - - std::unordered_map<index_t, index_t> pivot_column_index; - pivot_column_index.resize(columns_to_reduce.size()); - - compute_pairs( - columns_to_reduce, - pivot_column_index, - comp, comp_prev, - dim, n, - threshold, - modulus, - multiplicative_inverse, - binomial_coeff - ); - } - #endif - #endif - + exit(exit_code); } + +int main(int argc, char** argv) { + + if (argc < 2) print_usage_and_exit(-1); + + const char* filename = nullptr; + + index_t dim_max = 1; + value_t threshold = std::numeric_limits<value_t>::max(); + +#ifdef USE_COEFFICIENTS + coefficient_t modulus = 2; +#else + const coefficient_t modulus = 2; +#endif + + for (index_t i = 1; i < argc; ++i) { + const std::string arg(argv[i]); + if (arg == "--help") { + print_usage_and_exit(0); + } else if (arg == "--top_dim") { + std::string parameter = std::string(argv[++i]); + size_t next_pos; + dim_max = std::stol(parameter, &next_pos); + if (next_pos != parameter.size()) print_usage_and_exit(-1); + } else if (arg == "--threshold") { + std::string parameter = std::string(argv[++i]); + size_t next_pos; + threshold = std::stof(parameter, &next_pos); + if (next_pos != parameter.size()) print_usage_and_exit(-1); +#ifdef USE_COEFFICIENTS + } else if (arg == "--modulus") { + std::string parameter = std::string(argv[++i]); + size_t next_pos; + modulus = std::stol(parameter, &next_pos); + if (next_pos != parameter.size() || !is_prime(modulus)) print_usage_and_exit(-1); +#endif + } else { + if (filename) { print_usage_and_exit(-1); } + filename = argv[i]; + } + } + + std::ifstream input_stream(filename, std::ios_base::binary | std::ios_base::in); + if (input_stream.fail()) { + std::cerr << "couldn't open file " << filename << std::endl; + exit(-1); + } + + std::vector<value_t> diameters; + +#ifdef FILE_FORMAT_DIPHA + + int64_t magic_number; + input_stream.read(reinterpret_cast<char*>(&magic_number), sizeof(int64_t)); + if (magic_number != 8067171840) { + std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl; + exit(-1); + } + + int64_t file_type; + input_stream.read(reinterpret_cast<char*>(&file_type), sizeof(int64_t)); + if (file_type != 7) { + std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl; + exit(-1); + } + + int64_t n; + input_stream.read(reinterpret_cast<char*>(&n), sizeof(int64_t)); + + distance_matrix dist_full; + dist_full.distances = std::vector<std::vector<value_t>>(n, std::vector<value_t>(n)); + + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + double val; + input_stream.read(reinterpret_cast<char*>(&val), sizeof(double)); + dist_full.distances[i][j] = val; + } + } + std::cout << "distance matrix with " << n << " points" << std::endl; + + compressed_lower_distance_matrix_adapter dist(diameters, dist_full); + std::cout << "distance matrix transformed to lower triangular form" << std::endl; + +#endif + +#ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV + + std::vector<value_t> distances; + std::string value_string; + while (std::getline(input_stream, value_string, ',')) + distances.push_back(std::stod(value_string)); + + compressed_upper_distance_matrix_adapter dist_upper(distances); + + index_t n = dist_upper.size(); + std::cout << "distance matrix with " << n << " points" << std::endl; + + compressed_lower_distance_matrix_adapter dist(diameters, dist_upper); + std::cout << "distance matrix transformed to lower triangular form" << std::endl; + +#endif + +#ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV + + std::vector<value_t>& distances = diameters; + std::string value_string; + while (std::getline(input_stream, value_string, ',')) + distances.push_back(std::stod(value_string)); + + compressed_lower_distance_matrix_adapter dist(distances); + + index_t n = dist.size(); + + std::cout << "distance matrix with " << n << " points" << std::endl; + +#endif + + auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); + std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; + + assert(dim_max <= n - 2); + + binomial_coeff_table binomial_coeff(n, dim_max + 2); + std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus)); + + std::vector<diameter_index_t> columns_to_reduce; + for (index_t index = n; index-- > 0;) columns_to_reduce.push_back(diameter_index_t(0, index)); + + for (index_t dim = 0; dim <= dim_max; ++dim) { + + rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff); + rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff); + + hash_map<index_t, index_t> pivot_column_index; + pivot_column_index.reserve(columns_to_reduce.size()); + + compute_pairs(columns_to_reduce, pivot_column_index, dist, comp, comp_prev, dim, n, + threshold, modulus, multiplicative_inverse, binomial_coeff); + + if (dim < dim_max) + assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, + threshold, binomial_coeff); + } +}
\ No newline at end of file diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index f3e4722..d924596 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -209,7 +209,7 @@ buildSettings = { GCC_PREPROCESSOR_DEFINITIONS = ( "$(inherited)", - FILE_FORMAT_LOWER_TRIANGULAR_CSV, + FILE_FORMAT_UPPER_TRIANGULAR_CSV, USE_COEFFICIENTS, STORE_DIAMETERS, PRINT_PERSISTENCE_PAIRS, |