diff --git a/CMakeLists.txt b/CMakeLists.txt index 6a8169ea9d..02cdbe61c9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,6 +37,7 @@ add_gudhi_module(Hasse_complex) add_gudhi_module(Persistence_representations) add_gudhi_module(Persistent_cohomology) add_gudhi_module(Rips_complex) +add_gudhi_module(Ripser) add_gudhi_module(Simplex_tree) add_gudhi_module(Skeleton_blocker) add_gudhi_module(Spatial_searching) diff --git a/data/distance_matrix/tore3D_1307_distance_matrix.csv b/data/distance_matrix/tore3D_1307_distance_matrix.csv index afab863435..8df9e05ae8 100644 --- a/data/distance_matrix/tore3D_1307_distance_matrix.csv +++ b/data/distance_matrix/tore3D_1307_distance_matrix.csv @@ -1,4 +1,4 @@ -0. + 0.485593 0.738424 0.885478 0.188641 0.420709 0.67781 @@ -1304,4 +1304,4 @@ 1.8113 1.83101 1.40725 1.85421 1.59879 0.734918 1.35786 1.48372 1.35527 0.738936 1.87299 1.62667 1.15888 1.06351 1.01194 1.43559 1.64675 1.775 0.974258 1.48989 1.41247 0.854514 1.2493 0.271864 0.811688 1.23359 0.960687 1.67431 0.727966 0.737152 0.79253 0.683712 0.639196 0.667212 0.682323 1.14862 1.22595 1.26012 1.27812 1.37573 1.33697 1.25392 1.19951 1.39057 1.33256 1.97046 1.95145 1.80436 1.8655 1.82531 1.85551 1.88924 1.99138 2.02134 1.99613 1.96933 1.40568 1.67377 1.58185 1.7263 1.60619 1.42651 1.49596 1.49267 1.42953 1.72673 1.81775 1.93221 1.85297 1.80766 1.85765 1.96055 1.93289 2.00576 1.96405 1.86285 1.8234 1.96457 1.93582 2.02526 2.00864 2.00286 2.01824 2.02595 2.02237 2.00903 2.01265 1.7845 1.81151 1.9112 1.97035 1.98978 1.96702 1.86384 1.95143 2.01551 2.00217 2.01875 2.00244 1.98826 2.01706 1.66968 1.70363 1.68626 1.47424 1.46825 1.5299 1.63291 1.58239 1.51358 1.42563 1.4352 1.45417 1.45202 1.68629 1.75018 1.85174 1.83792 1.91947 1.9758 1.94166 1.89838 1.70397 1.74889 1.84373 1.9068 1.95718 1.96579 1.84538 1.90322 1.58901 1.52726 1.72796 1.67983 1.75682 1.76784 1.57751 1.6655 1.93873 1.87171 1.94515 1.96446 1.99388 1.99073 2.00265 2.00664 1.97896 1.96138 1.97654 1.97877 1.96627 1.95171 1.89771 1.95369 1.42061 1.43024 1.46063 1.61039 1.71344 1.55997 1.85248 1.91893 1.92685 1.86898 1.78178 1.77835 1.89548 1.84493 1.89157 1.8296 1.9435 1.93811 1.68545 1.57623 1.74553 1.75844 1.59212 1.86675 1.91556 1.90257 1.91678 1.87993 1.91166 1.90129 1.93324 1.92074 1.94322 1.9691 1.98546 1.93603 1.95789 1.72598 1.76881 1.89665 1.90364 1.86919 1.8385 1.86261 1.84555 1.80667 1.72776 1.71629 1.76025 1.98268 1.92722 1.93867 1.9793 1.89189 1.88906 1.78215 1.85157 1.36816 1.44855 1.43386 1.64106 1.54476 1.87489 1.86016 1.75828 1.67845 1.7556 1.82112 1.80523 1.70747 1.532 1.45531 1.65536 1.62925 1.77515 1.76142 1.54888 1.74794 1.75403 1.44406 1.52462 1.36109 1.36598 1.31429 1.25768 1.33653 1.34458 1.42112 1.61191 1.58118 1.54378 1.61754 1.32259 1.26893 1.39248 1.3482 1.35749 1.41012 1.43227 1.62882 1.5242 1.5245 1.6006 1.83633 1.70594 1.63731 1.59898 1.77485 1.71038 1.0668 1.05746 1.16435 1.10789 1.15585 1.09257 0.969334 1.01022 0.991799 0.943252 0.887349 0.964135 1.79384 1.76752 1.77415 1.80313 1.64305 1.44631 1.45004 1.25433 1.25515 1.35517 1.76848 1.76948 1.68333 1.71417 1.81893 1.71534 1.79983 1.35689 1.09527 1.1546 1.43486 1.34593 1.27261 1.71953 1.64948 1.67351 1.59911 1.51536 1.57753 1.62504 1.6471 1.67441 1.45793 1.54317 1.47098 1.66935 1.64055 1.50134 1.70298 1.73519 1.72391 1.659 1.70996 1.58477 1.56706 1.37889 1.39795 1.27385 1.18131 1.2227 1.27722 1.5255 1.57151 1.6708 1.6555 1.60871 1.57448 1.53082 1.62625 1.6071 1.57204 1.51411 1.55341 0.873697 1.06431 0.997937 0.971373 0.960128 1.0965 1.02087 1.52002 1.50205 1.54809 1.09773 1.0766 1.24163 1.09525 1.24611 1.18575 1.03101 0.93392 0.912656 0.976079 1.0423 1.07896 1.41822 1.29191 1.3674 1.42345 1.49788 1.38438 1.3147 1.41888 1.44542 1.16303 1.12012 1.16752 1.18438 1.51145 1.57805 1.60372 1.57632 1.50968 1.5669 1.23548 1.23372 1.22455 1.12944 1.05352 1.0912 1.17433 1.36255 1.2636 1.38944 1.33043 1.51052 1.47754 1.4605 1.37448 1.35176 1.32133 1.17785 1.22136 1.21835 1.17859 1.16144 1.47635 1.41482 1.49736 1.52457 1.36895 1.31967 1.2967 1.35086 1.04093 0.916437 0.929274 0.79795 0.859095 0.695163 1.14828 1.05845 0.949859 1.1118 1.13279 1.26424 1.26561 1.19127 1.132 1.09552 1.20813 1.177 0.930393 0.91014 0.989762 0.980488 1.43423 1.4442 1.29022 1.37393 1.43576 1.40385 1.30426 1.37744 1.21372 1.22102 1.14048 1.03928 1.19242 1.30868 1.33738 1.2659 1.12811 1.02895 1.11244 1.1618 1.1638 1.10346 1.15013 1.06763 1.01931 1.04773 1.133 1.18901 1.00868 0.984821 1.07414 1.08318 1.2976 1.21636 1.15483 1.15539 1.08717 1.00725 1.07828 1.01964 1.04081 0.913429 0.979598 0.988745 0.938433 0.892117 0.963176 0.85066 0.83968 0.796097 0.890168 0.934945 0.889839 0.788952 0.715328 0.748381 0.80791 0.86644 0.809344 0.850741 0.663995 0.742685 1.2793 1.21112 1.30508 1.24767 1.21934 1.22921 1.29952 1.34351 0.808688 0.740942 0.92246 0.897505 0.650656 0.642831 0.569524 0.465291 0.604859 0.48917 0.40529 0.612826 0.588742 1.07992 1.01122 0.872725 0.776768 0.913865 0.931939 0.863055 0.888094 0.861099 0.811209 0.799629 0.850497 0.72409 0.68026 0.552046 0.455959 0.579881 0.601431 0.867201 0.837127 0.692761 0.798585 1.20295 1.18306 1.11343 1.16146 0.476312 0.574951 0.703217 0.76394 0.597126 0.621667 0.764241 0.765179 0.688236 0.624801 0.678343 0.712825 1.06596 0.986578 1.09304 1.11838 0.796594 0.773031 0.667449 0.74529 0.66354 0.715221 0.423705 0.502116 0.396115 0.482944 0.625745 0.523279 0.771323 0.816987 0.738525 0.683228 0.785948 1.02157 1.03283 0.97919 0.126503 0.211205 0.291484 0.272291 0.290479 0.396017 0.326337 0.37907 0.36418 0.121886 0.2308 0.300104 0.314049 0.30922 0.121067 0.235705 0.609519 0.632652 0.590344 0.564364 0.496134 0.423377 0.453326 0.557537 0.659466 0.70029 0.719593 0.685578 0.642449 0.648685 0.627906 0.597771 0.614412 0.84144 0.783122 0.749231 0.613206 0.606897 0.598438 0.564964 0.602098 0.587168 0.427922 0.380698 0.52745 0.497667 0.503691 0.537161 0.516879 0.442823 0.384398 0.431576 0.680797 0.636572 0.511195 0.417376 0.53458 0.130915 0.240805 0.389947 0.45873 0.336373 0.311546 0.703108 0.671967 0.754937 0.763575 0.69352 0.589324 0.570873 0.510442 0.437927 0.441057 0.527846 0.605883 0.656054 0.596806 0.578645 0.311169 0.133036 0.626541 0.430062 0.472189 0.563565 0.49529 0.568428 0.536673 0.566394 0.594521 0.623209 0.654513 0.609384 0.610891 0.227386 0.11966 0.118102 0.226183 0.277791 0.283032 0.723218 0.728508 0.661142 0.62411 0.641573 0.676011 0.527855 0.438778 0.427401 0.5025 0.586256 0.601419 0.843361 0.893317 0.868497 0.779933 0.706338 0.753861 1.00819 1.03844 0.954903 1.66379 1.59253 1.5083 1.49317 1.57343 1.6681 1.53212 1.5951 1.55187 1.46302 1.41187 1.44564 1.45099 1.4179 1.43501 1.50856 1.56921 1.53176 1.57713 1.5509 1.64333 1.71746 1.68028 1.47928 1.55795 1.55689 1.47512 0.662094 0.742351 0.760192 0.681638 1.18916 1.22139 1.29108 0.843934 0.832544 0.896884 0.915947 0.802527 0.83659 0.825165 0.609357 0.604068 0.550072 0.544923 0.919838 0.914328 0.861435 0.913015 0.882821 0.792672 0.771442 0.847163 1.53605 1.60794 1.64368 1.57263 1.7035 1.61916 1.62335 1.70878 1.39883 1.41615 1.38875 1.74505 1.72135 1.78414 1.45472 1.44102 1.49336 1.30846 1.32934 1.35699 1.37155 1.35146 1.04037 1.06041 0.99604 1.38613 1.3326 1.38906 1.43806 1.83752 1.83194 1.8553 1.86498 1.73835 1.78416 1.73976 1.41622 1.39727 1.3805 1.40672 0.937361 0.837279 0.794736 0.853857 0.934672 1.36452 1.4044 1.49135 1.51055 1.43757 1.48617 1.44294 1.38089 1.41762 1.32601 1.36839 1.44329 1.39769 1.55685 1.62659 1.63069 1.56643 1.70919 1.62524 1.68616 0.394551 0.387226 0.296244 0.310481 0.705857 0.658655 0.703345 0.746256 1.27347 1.27913 1.22296 1.1966 1.24175 1.37872 1.42792 1.37406 1.289 1.29399 0.984593 1.03172 1.69945 1.613 1.68874 1.45192 1.42745 1.42022 1.84972 1.80964 1.81454 1.87124 1.78626 1.83223 1.12101 1.14246 1.21939 1.193 1.39589 1.34172 1.46308 1.21451 1.26527 1.22129 1.16622 0.238796 0.305216 0.35907 0.305265 0.495571 0.435718 0.489897 0.543145 0.837877 0.805183 0.850178 0.885425 0.575302 0.552404 0.469227 0.499093 1.0018 0.972536 1.0347 1.05513 0.534851 0.530583 0.567363 0.586379 0.578315 0.883779 0.938304 0.965812 0.914616 0.873708 0.9183 0.557698 1.01313 0.522482 0.846213 0.488525 0.302506 1.21519 1.17006 1.3516 1.23849 0.702858 0.343768 1.59976 1.38259 1.43091 1.43892 0.867951 1.39702 1.8504 1.38613 1.33755 1.75607 1.66274 1.5882 0.838831 0.578259 0.844001 0.875418 0.709374 1.51594 1.62831 1.476 1.49588 1.58597 0.961611 0.80177 0.510185 0.680305 0.189329 0.638388 0.582976 0.509809 0.573873 0.210348 0.620164 0.512624 0.676809 0.640611 0.393932 0.227141 0.41504 0.583688 0.467249 0.461949 0.554513 0.592052 0.593682 0.774987 0.774444 0.777393 0.569256 0.613528 0.677729 0.506252 0.525899 0.209868 0.20758 0.35218 0.43679 0.462101 0.346819 0.348812 0.200512 0.201095 0.966618 0.712409 0.747531 0.846645 0.579578 0.474287 0.401887 0.658253 0.729157 1.04843 0.656196 0.722866 0.682722 0.528892 1.14933 0.766365 0.520436 0.650643 0.810672 0.831851 0.840755 1.00022 0.512294 0.396505 0.487879 0.556706 0.831526 1.28064 1.24413 0.649285 0.758711 0.723748 0.844258 0.778549 0.92804 0.911689 0.968153 0.994372 1.08416 1.22427 0.989852 0.925608 1.10676 1.09703 1.04067 1.22705 1.12678 1.37002 1.29923 1.36018 0.919901 0.958454 0.851302 1.15144 1.21746 1.10454 1.04888 0.844568 0.694397 0.767827 0.86556 0.980617 1.08262 1.31725 1.46937 1.09044 1.11985 1.26921 1.40823 1.42798 1.28868 1.15187 1.15685 1.5536 1.51233 1.10716 1.37959 1.46217 1.31709 1.30689 0.986266 1.13962 1.30143 0.999937 1.46598 1.03615 1.05236 0.938564 0.900284 1.47898 1.57527 1.63715 1.56484 1.53493 1.19577 1.28398 1.48085 1.65486 1.69909 1.59229 1.48292 1.57024 1.5572 1.6255 1.69598 1.33858 1.15377 1.13921 1.29261 1.2931 1.74361 1.73675 1.26317 1.35356 1.54501 1.74707 1.80651 0.915041 0.922697 1.06789 1.16358 1.11405 1.00361 1.69392 1.76369 1.61341 1.53719 1.33235 1.29067 1.32992 1.49801 1.65659 1.51305 1.43062 1.29705 1.44477 1.29123 1.8146 1.66462 1.52032 1.65892 1.54019 1.7309 1.78485 1.53059 1.37613 1.37891 1.83538 1.93643 1.68674 1.76378 1.82047 1.8784 1.81127 1.66094 1.9693 1.92927 1.92276 1.88183 1.88465 1.84905 1.51395 1.65576 1.89773 1.83057 1.85675 1.67984 1.52724 1.39222 1.93219 1.9389 1.96286 1.99783 1.91415 1.97858 1.67713 1.62578 1.91097 1.83746 1.63779 1.95608 1.90004 1.77654 1.64913 1.41657 1.51088 1.5599 1.60043 1.99721 2.0142 1.92637 1.88933 1.98024 2.02369 2.00853 1.9033 1.98674 1.98617 1.89529 1.76604 1.77442 1.43972 1.4305 1.64259 1.69519 1.57292 1.46244 1.41279 1.39296 2.01166 1.94185 1.78671 1.75622 1.90229 2.0042 1.30602 1.39861 1.32163 1.21035 1.16667 0.629723 0.731106 0.792859 0.732444 1.28782 1.92318 1.49696 1.89986 2.01196 2.0223 1.74832 1.45826 1.97937 1.81801 1.947 1.647 1.76425 1.98477 1.97858 1.49632 1.79124 1.93008 1.97107 1.78628 1.93483 1.95442 1.8012 1.99941 1.89196 1.47112 1.46997 1.60294 1.83798 1.75768 1.8013 1.3839 1.6879 1.37723 1.46375 1.04013 1.83361 1.88389 1.46488 1.73401 1.86193 1.4093 1.2111 1.7433 1.67131 1.61855 1.69166 1.74829 1.56965 1.36641 1.43024 1.6412 1.60272 1.13546 1.45595 1.2152 1.61296 1.03545 1.44008 1.2017 1.27208 1.52851 1.39262 1.46009 0.93886 1.44622 1.21179 1.22769 1.08537 1.14407 1.03808 0.851871 1.06503 0.96875 0.833911 1.35211 1.41836 1.22869 1.23445 0.893315 0.980055 0.686534 0.585535 0.949308 0.894314 0.629484 0.921146 1.20268 0.650772 0.750682 1.16368 1.31811 1.13942 0.820043 0.581968 0.895501 0.359771 0.665476 0.651668 0.889529 0.628317 0.476551 0.593194 0.326165 0.616672 0.584721 0.347747 1.06882 0.748825 0.559261 0.401563 0.604459 0.570964 0.341348 0.818873 0.697751 0.519713 1.89976 1.92976 1.61101 1.97065 1.84013 0.585472 1.02538 1.74986 1.02428 1.00864 1.74242 1.74854 1.47774 1.38265 1.33031 1.37387 1.81119 1.53898 0.663134 1.75274 1.32086 0.566972 1.40414 0.49327 1.06903 1.50145 1.19234 1.79895 0.427452 0.447907 0.490767 0.386815 0.349789 0.344431 0.372258 1.43033 1.50653 1.56281 1.56863 1.66472 1.62173 1.52247 1.45446 1.65399 1.60459 1.9136 1.88231 1.74451 1.794 1.74042 1.73967 1.79181 1.91282 1.9491 1.89131 1.87567 1.28006 1.62927 1.52955 1.64067 1.53192 1.35271 1.45772 1.47867 1.39523 1.71658 1.80559 1.88417 1.80446 1.77105 1.83228 1.92312 1.9068 1.98958 1.94939 1.86319 1.83898 1.9749 1.93437 1.96846 1.9574 1.96384 1.98913 1.99435 1.97647 1.90977 1.93084 1.72726 1.73498 1.82562 1.89546 1.94379 1.93521 1.8144 1.89398 2.01066 2.01051 2.01219 2.0075 1.96757 1.99725 1.64951 1.66962 1.62765 1.40822 1.42284 1.49334 1.58895 1.57039 1.4958 1.41332 1.43548 1.47286 1.44836 1.52041 1.59905 1.67564 1.68123 1.7927 1.85458 1.78863 1.7544 1.57145 1.63265 1.81128 1.88052 1.94765 1.96917 1.84283 1.89071 1.59374 1.5492 1.74671 1.68447 1.73631 1.76299 1.5699 1.64128 1.96316 1.89973 1.99404 1.99982 2.02657 2.03396 2.02084 2.0263 2.03484 2.02034 2.00661 2.02117 2.03029 2.03006 1.95532 2.00884 1.44603 1.4941 1.5006 1.62267 1.72357 1.59278 1.86392 1.93375 1.95595 1.91141 1.81986 1.7992 1.9496 1.88935 1.96095 1.90881 2.02676 2.00984 1.72238 1.61689 1.81567 1.80916 1.43101 1.94174 2.00203 2.00682 2.01477 2.00609 2.02172 1.72144 1.76551 1.74086 1.77433 1.82264 1.85363 1.78284 1.80164 1.52368 1.55045 1.70179 1.71093 1.66884 1.62372 1.65228 1.62291 1.56787 1.49323 1.45483 1.51395 1.87231 1.82343 1.80081 1.85236 1.71386 1.72473 1.59391 1.66179 1.24183 1.51694 1.52362 1.71799 1.60652 1.74288 1.74592 1.65295 1.56358 1.58448 1.6613 1.66272 1.57335 1.39835 1.31497 1.4916 1.4853 1.8558 1.86182 1.68288 1.90048 1.89009 1.27983 1.35248 1.16943 1.19204 1.14586 1.07344 1.19054 1.19582 1.60647 1.81001 1.40464 1.33065 1.41987 1.15531 1.08315 1.55114 1.52279 1.49976 1.51783 1.55971 1.72462 1.62569 1.64361 1.72718 1.96312 1.81521 1.77463 1.75046 1.91326 1.84496 0.756405 0.754721 0.844978 0.799579 0.852077 0.792356 0.689739 0.723022 0.705148 0.670649 0.62237 0.683976 1.95464 1.92143 1.94997 1.96384 1.37781 1.17582 1.19759 1.04831 1.00593 1.11503 1.94433 1.95226 1.87995 1.89678 1.6089 1.52938 1.60855 1.45742 1.28581 1.32617 1.56321 1.46672 1.40823 1.45107 1.36517 1.39217 1.31991 1.23366 1.27783 1.33695 1.44752 1.45367 1.25766 1.35258 1.25222 1.40934 1.36598 1.24795 1.90763 1.93145 1.90813 1.86796 1.90756 1.79067 1.75457 1.15254 1.15487 1.05488 0.951015 1.02226 1.07438 1.78251 1.82044 1.88888 1.87632 1.84145 1.82128 1.78172 1.85365 1.84504 1.79768 1.75347 1.76425 1.11898 1.26733 1.21648 1.2098 1.20882 1.32049 1.2499 1.73983 1.75037 1.78251 1.35274 1.34568 1.48312 1.29417 1.4564 1.41164 1.27404 1.18252 1.15252 1.20017 1.25862 1.30644 1.66804 1.55355 1.63046 1.70595 1.76519 1.66033 1.60115 1.70417 1.72119 0.946521 0.890032 0.927185 0.963718 1.24217 1.29937 1.30916 1.27289 1.2119 1.27535 1.51732 1.48814 1.46366 1.36528 1.30282 1.35346 1.43581 1.08518 0.994765 1.09715 1.02855 1.20686 1.18347 1.1466 1.0603 1.02164 0.998503 0.851113 0.901407 0.921543 0.899154 0.854026 1.15703 1.08719 1.18751 1.21219 1.66355 1.61673 1.60297 1.64888 1.37442 1.24641 1.25272 1.11145 1.17401 0.993017 1.43531 1.33418 1.22901 1.43297 1.43266 1.57708 1.57268 1.50056 1.45793 1.42724 1.52877 1.49952 1.23545 1.1734 1.26133 1.26474 1.15959 1.15487 1.00704 1.08057 1.12059 1.07793 0.996749 1.071 0.929315 0.9256 0.846718 0.757281 0.948462 1.0588 1.07252 0.995482 0.854151 0.75486 0.857745 0.896512 1.48772 1.43443 1.46701 1.3844 1.32558 1.3432 1.42807 1.49526 0.750264 0.698615 0.801338 0.829396 0.962377 0.882478 0.815843 0.819338 0.762241 0.674601 0.746523 1.34417 1.37288 1.24667 1.31563 1.32621 1.27422 1.23113 1.30215 1.1847 1.17636 1.10568 1.20283 1.25681 1.21668 1.11827 1.04357 1.06502 1.13228 1.1251 1.06766 1.12248 0.944634 1.01238 0.941713 0.871932 0.971484 0.916614 0.896463 0.913364 0.981336 1.01631 0.492834 0.415509 0.593734 0.578687 0.968102 0.966598 0.882431 0.785052 0.906137 0.786625 0.702988 0.906156 0.880837 0.741953 0.675398 0.54136 0.44693 0.585741 0.59868 1.20057 1.22422 1.18852 1.13046 1.12511 1.18136 1.02837 0.978844 0.84232 0.742937 0.864328 0.892447 0.596583 0.574617 0.42249 0.527373 0.870338 0.846219 0.791284 0.835812 0.798432 0.898876 1.0332 1.09755 0.916463 0.947185 1.07983 1.07861 0.986908 0.911964 0.969386 1.01487 0.730831 0.658742 0.766762 0.784179 1.12945 1.09878 0.988422 1.07609 0.973469 1.02862 0.740115 0.811251 0.663768 0.770546 0.88984 0.74808 0.53223 0.49758 0.413304 0.401242 0.48617 0.716332 0.716269 0.66304 0.378291 0.426847 0.549959 0.448938 0.404577 0.585684 0.519325 0.303656 0.123078 0.437357 0.527613 0.608858 0.627273 0.593238 0.453747 0.563523 0.313843 0.311934 0.262998 0.258096 0.205533 0.109298 0.119418 0.226016 0.936666 1.00032 1.02289 0.975581 0.9024 0.919348 0.903643 0.872793 0.874827 0.593511 0.57138 0.597714 0.839721 0.811968 0.753321 0.767503 0.829072 0.718385 0.570167 0.471562 0.594751 0.626645 0.73919 0.798959 0.789079 0.715111 0.637813 0.656496 0.464339 0.394907 0.245398 0.134382 0.315482 0.403578 0.491496 0.612007 0.654785 0.440084 0.505489 0.584998 0.421342 0.506934 0.530626 0.492162 0.496829 0.535255 0.528074 0.458478 0.39342 0.42726 0.586811 0.599234 0.667859 0.621129 0.125434 0.296582 0.471014 0.500414 0.600203 0.551609 0.484497 0.611954 0.64387 0.713888 0.717654 0.657736 0.624754 0.634343 0.672646 0.302234 0.289002 0.223569 0.115406 0.127413 0.242525 0.595509 0.569054 0.534066 0.56121 0.604826 0.609968 0.378872 0.29279 0.222576 0.259429 0.357068 0.409877 1.11695 1.17826 1.16421 1.07905 0.999101 1.03134 1.31766 1.33901 1.24647 1.83887 1.77612 1.68576 1.65499 1.72549 1.82814 1.41985 1.46476 1.40651 1.31551 1.27563 1.33033 1.37211 1.32428 1.32541 1.39357 1.46428 1.44677 1.48052 1.47345 1.57108 1.63165 1.58015 1.55633 1.63568 1.6519 1.57021 0.949919 1.02962 1.03947 0.96345 1.35655 1.39481 1.45312 0.550942 0.559435 0.612787 0.614794 0.6056 0.615993 0.602859 0.901921 0.879346 0.811607 0.826357 0.636692 0.640365 0.583055 1.21763 1.17872 1.08659 1.07157 1.15427 1.28998 1.35595 1.40561 1.34168 1.77017 1.68489 1.67125 1.75809 1.30316 1.31085 1.28033 1.48678 1.46657 1.53947 1.72774 1.70776 1.76142 1.41644 1.42312 1.43071 1.44829 1.45031 0.723714 0.753728 0.696957 1.63666 1.57472 1.62139 1.67859 1.98934 1.9841 1.9939 2.00304 1.50949 1.56136 1.4975 1.45146 1.43462 1.43733 1.46153 0.666788 0.563562 0.514905 0.578304 0.667117 1.52618 1.58507 1.66777 1.66694 1.58435 1.315 1.28521 1.21195 1.23605 1.11656 1.14652 1.23257 1.1992 1.3262 1.39431 1.38255 1.32038 1.87201 1.80365 1.8625 0.710437 0.70826 0.622863 0.631167 0.561919 0.524092 0.525357 0.562161 1.40093 1.42168 1.38303 1.3541 1.37932 1.61245 1.64604 1.57449 1.49621 1.52197 1.3079 1.35827 1.82523 1.71883 1.79832 1.40542 1.38091 1.35788 1.97823 1.9373 1.92693 1.97235 1.89224 1.945 1.32175 1.35262 1.41986 1.38254 1.6148 1.55166 1.67248 1.39028 1.44576 1.41669 1.35727 0.285157 0.342956 0.293416 0.224778 0.338607 0.333127 0.415699 0.416083 0.612896 0.590918 0.60905 0.634768 0.868003 0.851834 0.775783 0.799003 0.746745 0.710876 0.782773 0.811097 0.752243 0.70099 0.735328 0.79719 0.811916 0.634516 0.665414 0.677781 0.638161 0.622599 0.649066 0.758977 0.760111 0.822159 0.612584 0.373601 0.283994 1.40085 1.37005 1.56899 1.38317 0.542007 0.66555 1.36018 1.17226 1.26134 1.60281 0.595102 1.44299 1.99583 1.62775 1.42762 1.50805 1.72 1.3465 1.1394 0.857072 0.616846 0.58641 0.993357 1.60223 1.54211 1.37895 1.36552 1.75429 1.26185 1.08855 0.316263 0.581165 0.19967 0.636431 0.669517 0.564404 0.491533 0.196831 0.627311 0.471639 0.45119 0.55519 0.549196 0.388766 0.23404 0.397092 0.720979 0.542088 0.666232 0.771291 0.777359 0.593581 0.591085 0.559983 0.823001 0.858332 0.962947 0.178561 0.213594 0.515162 0.535007 0.209786 0.174983 0.356915 0.499957 0.572928 0.463258 0.369152 0.666897 0.40318 0.480283 0.592262 0.825462 0.726339 0.700272 0.951583 1.05311 0.720672 0.945463 1.02684 1.01067 0.849326 0.818267 0.50429 0.803434 0.945275 1.13392 1.16477 0.511371 0.665732 0.80533 0.705093 0.799179 0.877294 0.505089 0.954394 0.907533 0.936592 1.02919 1.04729 1.16305 1.11142 1.26696 1.24914 1.29822 0.670114 0.747587 0.886176 0.721367 0.658375 1.41373 1.4225 0.775288 0.96983 0.843286 1.05191 1.0037 1.08074 1.19503 1.25385 1.15562 1.47857 1.53358 1.41591 1.33816 1.1297 0.98756 1.07385 1.18785 1.31213 1.41 1.62122 1.14955 0.798081 0.802813 0.939656 1.08525 1.1235 1.00211 1.40197 1.43004 1.25074 1.22815 0.871748 1.66644 1.73828 1.56493 1.58057 1.21905 1.34901 1.52774 1.25763 1.70096 1.27325 1.26882 1.16995 1.14654 1.70397 1.81363 1.86715 1.80902 1.79315 0.980084 1.04945 1.68007 1.84729 1.90788 1.32968 1.24647 1.35749 1.26499 1.3324 1.42561 1.48111 1.33194 1.32417 1.47562 1.40938 1.53476 1.92722 1.03609 1.09381 1.27993 1.91647 1.97319 0.649024 0.647993 0.776775 0.851091 0.798371 0.710234 1.8426 1.88617 1.7264 1.67456 1.4543 1.45187 1.14107 1.30601 1.84406 1.70793 1.60196 1.13093 1.26311 1.1051 1.95808 1.80583 1.63389 1.7466 1.38547 1.57795 1.66055 1.61309 1.45308 1.25329 1.6596 1.81507 1.43634 1.50948 1.59293 1.67144 1.60804 1.44159 1.82705 1.76215 1.74159 1.99478 2.00765 1.94493 1.35627 1.71245 1.98327 1.89238 1.88311 1.70976 1.58421 1.44452 2.00205 1.98187 2.03214 2.02873 1.95331 2.02441 1.67145 1.64719 1.91187 1.82066 1.52298 1.81853 1.73685 1.60549 1.50148 1.43049 1.51269 1.50309 1.57241 1.98916 2.02082 1.88522 1.82203 1.8873 1.98123 1.97555 1.91603 1.98363 1.95911 1.8581 1.74233 1.77707 1.4132 1.37615 1.54754 1.62991 1.54323 1.33795 1.33388 1.28405 1.92316 1.8321 1.68686 1.68422 1.843 1.93983 1.56276 1.67533 1.61806 1.50402 1.43475 0.318625 0.417016 0.505004 0.451636 1.56732 1.84056 1.4344 1.88641 1.94936 2.00427 1.70104 1.43322 1.93273 1.68303 1.92663 1.6397 1.72596 2.00654 2.03088 1.50915 1.79831 1.99088 1.98788 1.81869 2.02862 1.79588 1.57079 1.88342 1.70328 1.38003 1.51326 1.47853 1.74234 1.87648 1.85995 1.22808 1.50441 1.23017 1.60936 0.745087 1.97653 1.73506 1.23415 1.90213 1.66767 1.49578 1.36404 1.83403 1.39593 1.43949 1.44992 1.93624 1.73842 1.15728 1.53987 1.85358 1.824 1.37634 1.71942 1.01323 1.32754 1.24348 1.15631 1.42168 0.960762 1.21836 1.68314 1.71425 1.25327 1.14325 0.936667 1.54218 0.841101 0.812869 1.37439 1.18183 1.34791 1.2321 1.0925 1.0188 1.15925 0.922342 1.53114 0.592153 1.29513 0.993443 0.875952 0.616202 1.22721 0.924001 0.647294 0.875602 0.977511 1.06027 0.903469 1.51169 0.800863 1.15651 0.886894 0.580977 0.683731 0.335251 0.92988 0.639245 0.874124 0.666531 0.330129 0.572623 0.465032 0.315736 0.631683 0.751818 1.06444 0.573822 0.393734 0.701527 0.847469 0.362623 0.567206 0.604514 0.690462 0.339427 2.42085 2.4499 2.09679 2.49161 2.30944 0.398883 1.09951 2.19391 1.07661 1.44905 2.17696 2.26402 1.82802 1.71608 1.65548 1.86769 2.31644 1.85537 0.897208 2.20131 1.79876 0.897271 1.91374 0.822265 1.52322 1.92892 1.6637 2.31338 0.762653 0.781448 0.810082 0.728869 0.692905 0.65486 0.708024 1.84621 1.92616 1.95576 1.9787 2.08173 2.04153 1.95613 1.90055 2.09754 2.0376 2.39868 2.36031 2.22551 2.26907 2.20719 2.18548 2.25118 2.38665 2.42806 2.3488 2.33945 1.73843 2.12268 2.02179 2.11271 2.0129 1.8396 1.95863 1.98763 1.8996 2.21879 2.3073 2.37266 2.29144 2.26315 2.32948 2.41722 2.40505 2.49388 2.45286 2.36954 2.34907 2.48662 2.44166 2.45582 2.44661 2.45857 2.48882 2.49366 2.46954 2.37232 2.40468 2.21378 2.21142 2.29686 2.37327 2.43663 2.43395 2.30455 2.38092 2.51948 2.52298 2.52093 2.51973 2.47098 2.50132 2.14883 2.16317 2.11011 1.89532 1.91992 1.99062 2.07951 2.07437 2.00052 1.92488 1.95023 1.98915 1.96022 1.92171 2.01456 2.06777 2.09225 2.23135 2.29948 2.2062 2.17829 2.00395 2.07781 2.30862 2.38074 2.45457 2.4803 2.35055 2.39587 2.10224 2.0631 2.25695 2.19142 2.23803 2.2697 2.07803 2.14256 2.47784 2.41358 2.51389 2.51779 2.54513 2.55472 2.53663 2.54226 2.5572 2.54232 2.5247 2.54168 2.5535 2.55406 2.47621 2.53074 1.96583 2.01608 2.0206 2.13493 2.23433 2.10945 2.37557 2.44718 2.47283 2.42965 2.33626 2.31241 2.46871 2.40579 2.48172 2.42857 2.55051 2.53241 2.23617 2.13231 2.33335 2.32519 1.8394 2.46342 2.52547 2.52988 2.53792 2.5267 2.54467 2.11175 2.16943 2.13448 2.17859 2.24759 2.29131 2.20334 2.21855 1.89541 1.90085 2.07662 2.09002 2.04012 1.97718 2.00753 1.96402 1.88555 1.82163 1.73796 1.82249 2.32745 2.28311 2.23479 2.29493 2.1103 2.13514 1.9808 2.04593 1.70134 2.03715 2.04313 2.23522 2.12419 2.18184 2.19835 2.11227 2.01782 1.98976 2.0766 2.09355 2.01254 1.84382 1.75922 1.90638 1.91857 2.37587 2.38142 2.19557 2.41289 2.40583 1.69326 1.75185 1.55743 1.60208 1.56961 1.48863 1.63676 1.63635 2.10154 2.30337 1.80949 1.69513 1.80019 1.58508 1.49887 2.05711 2.0228 2.01173 2.03675 2.07437 2.24139 2.14228 2.15777 2.24027 2.48184 2.3319 2.28625 2.25817 2.42828 2.35865 0.628148 0.672942 0.807961 0.769041 0.916378 0.823412 0.834948 0.785506 0.629618 0.524458 0.666409 0.693211 2.46671 2.43334 2.45688 2.47583 1.65691 1.45384 1.50892 1.44366 1.34039 1.45394 2.45141 2.45709 2.37651 2.39969 1.96165 1.90853 1.98276 1.97826 1.78422 1.83252 2.07636 1.98389 1.9241 1.71792 1.59361 1.63251 1.57219 1.4842 1.47849 1.56642 1.81021 1.78785 1.62553 1.72903 1.59189 1.68161 1.61005 1.52584 2.40158 2.43003 2.41011 2.35733 2.40418 2.27808 2.24951 1.48488 1.45465 1.41146 1.30175 1.42022 1.45734 2.23646 2.2832 2.37511 2.36004 2.31615 2.28602 2.24163 2.33302 2.31714 2.27579 2.22073 2.24968 1.5828 1.75819 1.69875 1.67236 1.65926 1.79166 1.71971 2.21742 2.20891 2.25275 1.79873 1.77613 1.94163 1.78746 1.93804 1.88285 1.7337 1.64056 1.6208 1.67975 1.7416 1.77877 2.12095 1.99007 2.0695 2.13256 2.20959 2.08843 2.01299 2.12677 2.15429 1.33432 1.25951 1.28676 1.34798 1.48767 1.52951 1.51168 1.45711 1.39149 1.47731 1.92872 1.92971 1.9196 1.82446 1.74894 1.78354 1.86828 1.35486 1.29211 1.32967 1.23919 1.40098 1.4062 1.30976 1.22795 1.11663 1.13449 0.976811 1.05717 1.15799 1.18496 1.06503 1.29382 1.18485 1.33625 1.36605 2.07205 2.01709 1.99256 2.05306 1.67288 1.50653 1.54111 1.37607 1.46814 1.31886 1.83419 1.74531 1.63255 1.77115 1.81081 1.95286 1.95496 1.87048 1.79188 1.74451 1.8862 1.84752 1.61005 1.61523 1.68958 1.67539 1.38921 1.34766 1.20465 1.25377 1.24047 1.16168 1.10644 1.20289 1.11636 1.06575 0.971171 0.928664 1.26544 1.34882 1.32487 1.23314 1.0819 0.984495 1.15957 1.15561 1.83639 1.7596 1.82951 1.74133 1.70034 1.73727 1.82547 1.87865 1.07064 1.00486 1.11877 1.1795 1.02266 0.958791 0.821315 0.874752 0.90723 0.764858 0.844179 1.67972 1.68923 1.54158 1.60911 1.601 1.5293 1.47891 1.57406 1.46274 1.42815 1.46948 1.55979 1.59246 1.53315 1.4104 1.32642 1.40816 1.4499 1.56789 1.51875 1.5427 1.37316 1.45517 0.961522 0.859776 0.992026 0.903773 0.892224 0.952101 1.04578 1.06933 0.763599 0.666302 0.737845 0.801072 1.30966 1.27198 1.25106 1.13876 1.30579 1.17703 1.07682 1.25384 1.2608 0.722126 0.699549 0.646352 0.601173 0.524282 0.617731 1.4237 1.4534 1.42353 1.36364 1.30899 1.39534 1.29346 1.27392 1.17251 1.08414 1.13451 1.18999 0.782609 0.708623 0.714718 0.786426 0.822143 0.793198 0.681415 0.766574 1.10403 1.18384 1.28346 1.32966 1.11817 1.19102 1.25363 1.28034 1.22111 1.14975 1.14222 1.21673 0.649797 0.537568 0.636974 0.691589 1.3242 1.26801 1.16107 1.27422 1.10492 1.16266 1.00681 1.0132 0.822079 0.920105 0.937478 0.769093 0.125517 0.445846 0.500331 0.301496 0.337182 0.53178 0.536444 0.469929 0.757387 0.780372 0.889111 0.769575 0.7161 0.827248 0.810389 0.616894 0.572715 0.815022 0.899045 0.988793 0.912122 0.82217 0.807378 0.904815 0.378175 0.466745 0.601391 0.631371 0.601945 0.55833 0.509163 0.525009 1.0124 1.11438 1.1637 1.11656 0.977977 1.03485 1.07875 1.08292 1.0131 0.596665 0.125757 0.320407 0.88238 0.80473 0.700344 0.741578 0.829495 0.689357 0.764093 0.7076 0.655075 0.739381 0.941097 1.02495 1.05829 1.02233 0.939454 0.899611 0.119091 0.226388 0.371523 0.448588 0.29239 0.711634 0.724041 0.717246 0.682676 0.535245 0.651771 0.409574 0.667887 0.664429 0.606692 0.572535 0.548691 0.562907 0.623703 0.64894 0.617521 0.583604 0.532394 0.450896 0.597162 0.595339 0.480692 0.619081 0.125671 0.473789 0.557772 0.328009 0.352669 0.42427 0.531819 0.626568 0.609887 0.504587 0.420592 0.414891 0.50621 0.671836 0.684226 0.619618 0.556637 0.581017 0.624618 0.244717 0.124104 0.125432 0.24305 0.314468 0.317892 0.591064 0.601254 0.595208 0.610875 0.626046 0.595641 1.5482 1.58951 1.55726 1.46825 1.40717 1.4629 1.66309 1.7068 1.62583 2.34256 2.27478 2.18504 2.15977 2.23403 2.33642 1.87435 1.9019 1.83211 1.74458 1.71965 1.7885 1.85692 1.80255 1.79308 1.85382 1.92862 1.92461 1.9439 1.94994 2.04717 2.09822 2.03792 2.07673 2.15485 2.17024 2.08946 1.34343 1.41655 1.44979 1.37881 1.86306 1.89674 1.95815 0.306262 0.236141 0.312723 0.368711 0.300691 0.362474 0.230918 1.00462 0.947753 0.875501 0.931107 0.963231 0.979866 0.91828 1.5613 1.54301 1.45044 1.40694 1.47715 1.60712 1.6602 1.73105 1.68041 2.28893 2.20354 2.18874 2.27552 1.78058 1.7791 1.74977 1.76919 1.75139 1.84522 2.16551 2.15062 2.20549 1.93842 1.94635 1.95482 1.97184 1.9723 0.918867 0.984319 0.966593 2.09089 2.03408 2.08927 2.14215 2.50475 2.49939 2.51238 2.52181 1.83391 1.89615 1.8026 1.97183 1.95639 1.96068 1.98341 0.939336 0.848575 0.834229 0.902683 0.973756 2.02986 2.08039 2.1656 2.1721 2.09248 1.72995 1.71543 1.63603 1.64488 1.50102 1.51102 1.60599 1.59083 1.64468 1.71052 1.67301 1.61278 2.37967 2.30378 2.36542 1.09859 1.07934 0.989402 1.01351 0.45701 0.510953 0.522539 0.477575 1.91973 1.93567 1.8923 1.86542 1.89618 2.07597 2.12152 2.05954 1.97688 1.98795 1.61387 1.66848 2.34124 2.23619 2.31604 1.90499 1.88133 1.85168 2.49788 2.45582 2.44705 2.49374 2.41089 2.46466 1.81095 1.8338 1.90582 1.87631 2.09161 2.0333 2.15615 1.89252 1.94427 1.908 1.85275 0.531152 0.493263 0.419397 0.463708 0.290096 0.346309 0.303473 0.231337 0.430251 0.471513 0.52903 0.494761 1.28486 1.26132 1.17665 1.20371 1.09137 1.05846 1.1399 1.17137 0.896742 0.7977 0.776821 0.863551 0.932482 0.359976 0.446967 0.447118 0.341776 0.295 0.375365 0.853748 1.11393 1.23131 0.479127 0.292731 0.472145 1.89788 1.85773 2.04329 1.89787 0.491107 1.04377 1.66422 1.55072 1.6813 2.10516 0.897373 1.96514 2.51297 2.08885 1.95122 1.80599 2.23809 1.66832 1.48481 0.941533 0.303228 0.304571 1.39557 2.12157 2.01034 1.85233 1.80778 2.2582 1.62253 1.50081 0.603648 0.211077 0.620443 0.408331 0.5166 0.437369 0.227503 0.540241 0.513615 0.596979 0.622485 0.494174 0.598681 0.608157 0.3829 0.20234 0.979522 0.690338 0.690223 0.719009 0.785011 0.21221 0.398819 0.537536 0.998148 0.951074 1.07107 0.569488 0.446569 0.812727 0.909783 0.598587 0.587288 0.601848 0.76635 0.865524 0.82675 0.729842 0.425188 0.398255 0.217283 0.222495 0.852563 0.814462 0.914905 1.04494 1.21687 0.582439 1.15036 1.19767 1.22501 1.09924 0.731753 0.722463 1.10765 1.21167 1.33971 1.36551 0.550215 0.621463 1.16878 1.09723 1.18784 1.20662 0.70053 0.977453 0.890308 1.36026 1.46254 1.3619 1.5019 1.37512 1.51649 1.52122 1.61267 0.83442 0.787328 0.917788 1.06123 0.992426 1.79133 1.76253 1.05194 1.24769 1.02483 1.14296 1.15914 1.29432 1.62033 1.64669 1.53129 1.81621 1.89737 1.77125 1.72739 1.5204 1.34929 1.37386 1.4507 1.59542 1.73064 2.01405 1.27085 1.06669 0.982728 1.04379 1.2145 1.31989 1.25747 1.84656 1.84803 1.42719 1.44032 1.24391 2.08364 2.17311 2.01524 2.00566 1.69052 1.83348 1.99846 1.70259 2.16686 1.73374 1.749 1.64426 1.60507 2.17912 2.28423 2.34529 2.27426 2.24748 1.35725 1.37497 2.16745 2.34337 2.39953 1.59268 1.55611 1.70334 1.48839 1.5452 1.68363 1.99229 1.83519 1.82327 1.97057 1.92998 1.88645 2.4282 1.40364 1.39962 1.5633 2.42391 2.48416 0.592115 0.74738 0.877537 0.860233 0.720028 0.556952 2.35319 2.40296 2.24162 2.18511 1.97249 1.95861 1.549 1.69838 2.34318 2.19997 2.10306 1.56194 1.65631 1.51014 2.47444 2.3185 2.15032 2.2654 1.81298 2.00085 2.10611 2.13103 1.97588 1.71041 2.05911 2.26218 1.7407 1.80217 1.92414 2.03327 1.97706 1.79335 2.2568 2.16923 2.13247 2.5163 2.529 2.46668 1.77178 2.22868 2.50535 2.41084 2.39822 2.22461 2.10364 1.96761 2.52499 2.50161 2.55558 2.54735 2.47047 2.54498 2.17837 2.1587 2.42138 2.32376 1.9712 2.24985 2.14337 2.00146 1.92149 1.94829 2.02303 1.9902 2.07038 2.49721 2.53352 2.37952 2.30346 2.35427 2.4754 2.47432 2.42682 2.49192 2.45846 2.35088 2.23961 2.28321 1.92044 1.8723 2.0157 2.11344 2.04462 1.79389 1.82009 1.75417 2.39202 2.28412 2.14374 2.15866 2.32524 2.42213 2.00926 2.10676 2.02304 1.90673 1.86695 0.669275 0.726899 0.838679 0.767602 1.99062 2.31006 1.92515 2.3888 2.43412 2.50878 2.18835 1.93855 2.42316 2.11317 2.42959 2.14366 2.22096 2.52206 2.55291 2.0243 2.30834 2.51194 2.50282 2.3321 2.55263 2.20959 1.90315 2.33382 2.08788 1.85307 2.03188 1.92809 2.20734 2.3946 2.37917 1.6563 1.89837 1.67384 2.11908 0.729667 2.49315 2.15985 1.5789 2.41026 2.04079 2.01736 1.87634 2.35221 1.65274 1.82751 1.7533 2.4382 2.24041 1.51924 2.05688 2.34149 2.30453 1.83562 2.16355 1.41538 1.54798 1.73289 1.40618 1.89379 1.14583 1.38829 2.0989 2.16851 1.57345 1.29898 1.15901 1.91429 1.18295 0.911735 1.67227 1.48115 1.76071 1.66181 1.54467 1.07505 1.42472 1.01091 1.92385 0.883755 1.65025 1.3693 1.28361 0.692144 1.46512 1.24326 0.873691 0.839257 1.25063 1.29336 1.23065 2.00274 0.756375 1.37981 1.03491 0.417383 1.01604 0.569144 1.07799 0.427723 0.909686 0.839438 0.309316 0.746737 0.557928 0.65004 0.981991 0.60237 1.21484 0.593454 0.409931 0.599778 1.09035 0.658104 0.649401 0.354569 0.649851 0.712571 0.524131 -2.40859 2.42069 1.99303 2.45586 2.16306 0.466146 1.40423 2.02177 1.37759 1.29686 2.32016 2.22623 1.59334 1.47611 1.41387 1.95853 2.23253 2.07792 1.14877 2.03789 1.9058 1.10339 1.8427 0.676623 1.38191 1.75157 1.53567 2.26161 0.956311 0.966318 1.01918 0.909112 0.856988 0.857398 0.900106 1.67252 1.75152 1.7529 1.79088 1.89184 1.85869 1.79385 1.75354 1.93707 1.86946 2.49051 2.46053 2.31417 2.36798 2.31475 2.31693 2.37013 2.49597 2.53526 2.47881 2.45904 1.87341 2.20969 2.11226 2.2306 2.11977 1.93475 2.03631 2.04817 1.96867 2.26966 2.36125 2.45637 2.37304 2.33487 2.39411 2.49409 2.47243 2.55713 2.51288 2.41525 2.38257 2.52804 2.49078 2.55158 2.53629 2.53868 2.56266 2.57115 2.55841 2.50065 2.52058 2.31186 2.32399 2.4169 2.48524 2.52681 2.5138 2.39764 2.48009 2.57525 2.56793 2.57967 2.57047 2.54208 2.57071 2.20628 2.23106 2.19544 1.98421 1.9949 2.05896 2.15363 2.12532 2.05571 1.979 1.99809 2.01907 2.00641 2.08417 2.16842 2.24508 2.25381 2.37373 2.44083 2.37051 2.33289 2.14266 2.20605 2.38898 2.45671 2.51762 2.53333 2.40858 2.46147 2.14113 2.0896 2.28657 2.23092 2.30878 2.32897 2.13647 2.21471 2.50776 2.43884 2.52585 2.54065 2.57189 2.57443 2.57683 2.57969 2.56663 2.54772 2.55659 2.56349 2.55781 2.54686 2.48776 2.54314 1.9958 2.02021 2.04234 2.17967 2.28199 2.13847 2.42305 2.49192 2.50644 2.45335 2.36339 2.35256 2.47501 2.41813 2.47642 2.41425 2.53744 2.52716 2.25183 2.1448 2.32456 2.33129 1.99522 2.46151 2.51322 2.50295 2.51245 2.47894 2.51096 2.29513 2.34527 2.32223 2.35741 2.40919 2.44255 2.3725 2.38875 2.09904 2.11709 2.27275 2.28665 2.24331 2.19081 2.21423 2.18192 2.11514 2.04941 1.98537 2.0589 2.46512 2.41657 2.39298 2.44442 2.29761 2.31266 2.17496 2.24182 1.83326 2.02998 2.01973 2.2204 2.12095 2.33578 2.33914 2.24512 2.15591 2.17084 2.25039 2.25441 2.16538 1.99095 1.9078 2.07963 2.07652 2.37014 2.36009 2.14839 2.34855 2.35541 1.84708 1.91335 1.7272 1.76054 1.72434 1.6534 1.7805 1.77951 2.01338 2.20135 1.98987 1.90169 1.99757 1.74498 1.6683 1.99087 1.94407 1.95811 2.00919 2.03228 2.21199 2.10861 2.11019 2.18658 2.43135 2.2925 2.22415 2.18448 2.36715 2.2994 0.921016 0.943901 1.09546 1.03911 1.17128 1.07429 1.04511 1.0134 0.876638 0.767961 0.863664 0.916852 2.39376 2.35994 2.3672 2.39864 1.90683 1.7022 1.74251 1.62829 1.56003 1.67271 2.36538 2.36328 2.27406 2.30938 2.16464 2.08718 2.16942 1.9485 1.69115 1.75296 2.02227 1.93796 1.86885 1.96905 1.85376 1.89308 1.83237 1.74451 1.75432 1.83346 1.99633 1.99173 1.80452 1.9043 1.78479 1.9187 1.85908 1.74748 2.2902 2.3243 2.31181 2.23751 2.29515 2.15917 2.14591 1.67918 1.66245 1.59607 1.4898 1.59297 1.63189 2.07354 2.12954 2.25118 2.23118 2.17587 2.13621 2.08986 2.20208 2.17537 2.14874 2.08169 2.13715 1.44875 1.65314 1.583 1.53536 1.51111 1.66272 1.58971 2.08589 2.05227 2.10906 1.65424 1.61819 1.80663 1.69011 1.83085 1.76312 1.60197 1.5061 1.49281 1.56286 1.63139 1.65948 1.9617 1.81894 1.89735 1.94828 2.03783 1.90435 1.81724 1.93696 1.97262 1.51643 1.44627 1.49163 1.53784 1.72372 1.77838 1.77939 1.7342 1.65573 1.7386 1.73664 1.76618 1.77155 1.68214 1.59522 1.61425 1.69817 1.60519 1.52902 1.59464 1.51047 1.67987 1.67601 1.59695 1.51335 1.41873 1.4274 1.26842 1.34222 1.41759 1.4249 1.3315 1.58535 1.4839 1.61391 1.64984 1.87165 1.81139 1.78224 1.85143 1.41377 1.23572 1.28436 1.11923 1.22059 1.1015 1.63605 1.56124 1.44706 1.53121 1.5974 1.7319 1.73759 1.64764 1.54837 1.4935 1.65526 1.61082 1.41521 1.46847 1.5319 1.50437 1.62565 1.59893 1.44043 1.5051 1.52111 1.45616 1.37077 1.46936 1.34843 1.31064 1.20984 1.14583 1.46275 1.55703 1.54689 1.4555 1.29772 1.19223 1.35931 1.36659 1.60386 1.5149 1.6101 1.52475 1.50057 1.55061 1.63655 1.6735 1.26593 1.23053 1.34254 1.38705 1.33082 1.26316 1.1371 1.18163 1.19009 1.05688 1.13846 1.4515 1.44519 1.29613 1.35584 1.33296 1.25427 1.20696 1.30827 1.2133 1.16483 1.27431 1.35545 1.37122 1.30322 1.1747 1.09224 1.20363 1.22706 1.41469 1.37369 1.37265 1.21425 1.30594 1.27381 1.17723 1.29725 1.20674 1.18131 1.22749 1.32525 1.36381 0.99314 0.895302 1.01117 1.05181 1.10974 1.05492 1.06791 0.954517 1.13588 1.00348 0.902975 1.05319 1.07678 1.03979 1.00584 0.924274 0.857178 0.838316 0.919027 1.14097 1.17011 1.14528 1.09003 1.01359 1.10807 1.04498 1.0437 0.969828 0.895438 0.906761 0.969393 0.97212 0.890419 0.875994 0.965695 1.13218 1.1118 0.987432 1.06901 0.900416 0.963048 1.03413 1.06441 0.859519 0.947445 0.958639 0.997954 0.963872 0.903585 0.862093 0.942282 0.974858 0.865043 0.955997 1.01493 1.03529 0.969639 0.879974 0.994124 0.806821 0.857301 0.795112 0.76853 0.588516 0.66225 0.617074 0.46742 0.463771 0.745808 0.762521 0.571364 0.646343 0.834444 0.856978 0.800609 0.641365 0.653902 0.729197 0.640529 0.611183 0.651506 0.657849 0.606901 0.635317 0.67961 0.746359 0.826262 0.721338 0.623145 0.662814 0.737936 0.602237 0.692315 0.783984 0.78705 0.726866 0.665266 0.644363 0.709107 0.696187 0.804021 0.861776 0.822207 0.670157 0.740309 0.810351 0.832509 0.736739 0.776815 0.416347 0.403869 0.577968 0.485719 0.389377 0.416903 0.502974 0.416868 0.595791 0.587723 0.488449 0.544219 0.714945 0.792584 0.840766 0.827026 0.756213 0.702387 0.404749 0.471008 0.535968 0.561959 0.453718 0.575201 0.552031 0.482822 0.411567 0.382131 0.463977 0.424512 0.811639 0.821861 0.756482 0.683141 0.538913 0.488643 0.513511 0.560952 0.58337 0.586828 0.41882 0.367017 0.362679 0.417801 0.520671 0.526956 0.288852 0.281494 0.289711 0.0942476 0.207717 0.0973594 0.215389 0.294184 0.283215 0.263503 0.266901 0.0922615 0.200502 0.602362 0.601147 0.559238 0.556836 0.597546 0.60068 0.292273 0.297969 0.215011 0.0974517 0.0948727 0.207196 0.635848 0.63623 0.660427 0.713264 0.738437 0.683064 1.39401 1.42125 1.37743 1.28879 1.23998 1.30855 1.43454 1.49138 1.42302 2.26013 2.18604 2.10263 2.09048 2.17216 2.26692 1.99481 2.03609 1.97572 1.88764 1.85469 1.91015 1.96103 1.91495 1.91829 1.98653 2.05636 2.03675 2.05428 2.04674 2.1415 2.20365 2.15276 2.07244 2.1509 2.15328 2.07189 1.16289 1.22872 1.27569 1.21043 1.78402 1.81036 1.8792 0.638157 0.572696 0.641422 0.702138 0.475627 0.559772 0.484245 0.705091 0.637532 0.5764 0.645757 1.17618 1.18222 1.1157 1.33794 1.33426 1.24654 1.18857 1.24633 1.83904 1.89903 1.95996 1.90173 2.29408 2.20902 2.20754 2.29389 1.89181 1.8945 1.8684 2.01069 1.98722 2.07688 1.99171 1.98491 2.03513 1.90775 1.92414 1.94815 1.96501 1.94982 1.1855 1.24079 1.20707 1.9439 1.89656 1.95992 2.00388 2.43651 2.43245 2.45693 2.46533 2.04787 2.10773 2.02825 1.99085 1.97563 1.96533 1.98784 1.1353 1.03524 1.02071 1.09492 1.16813 1.94941 1.98308 2.07126 2.09439 2.0233 1.90262 1.87612 1.80151 1.82191 1.69379 1.71607 1.8064 1.77954 1.85069 1.92131 1.89726 1.83177 2.29793 2.20822 2.27242 0.930677 0.906179 0.821882 0.849448 0.509934 0.548438 0.613512 0.58519 1.87495 1.88107 1.8243 1.79784 1.84237 1.93149 1.99114 1.94532 1.85787 1.84954 1.36274 1.41643 2.2999 2.21104 2.28779 1.9863 1.96043 1.94078 2.45187 2.41154 2.41545 2.46386 2.37517 2.42473 1.70392 1.71651 1.7948 1.77639 1.97454 1.92507 2.04656 1.81311 1.86124 1.81338 1.76171 0.462452 0.406968 0.403334 0.459287 0.379784 0.363519 0.276845 0.296857 0.614411 0.626711 0.714204 0.70578 1.12179 1.09629 1.00991 1.03663 1.28599 1.25728 1.33893 1.36408 0.653365 0.56334 0.514732 0.583156 0.665286 0.610661 0.712619 0.740626 0.648766 0.574461 0.657646 0.595651 1.30992 1.0662 0.664572 0.329527 0.428537 1.81096 1.74899 1.91401 1.84057 0.563963 0.875866 1.87959 1.74764 1.85011 2.02174 1.08788 1.97705 2.45099 1.95067 1.93357 2.03814 2.25021 1.89852 1.26739 0.642227 0.53429 0.639755 1.21813 2.11089 2.11427 1.97073 1.93993 2.18268 1.40687 1.33487 0.678511 0.184162 0.583467 0.167915 0.179722 0.178038 0.187225 0.508262 0.341121 0.54484 0.74929 0.466582 0.377821 0.476812 0.459126 0.374403 0.770702 0.546569 0.466401 0.388917 0.490698 0.391096 0.521396 0.681341 0.745317 0.661602 0.764545 0.718939 0.623743 0.645095 0.755011 0.62248 0.680102 0.610338 0.623412 0.696206 0.687785 0.625295 0.744738 0.6678 0.532793 0.539812 0.535555 0.543322 0.691639 0.732928 0.922597 0.909668 0.886712 0.908162 0.960647 0.868711 1.04681 0.890709 0.899559 0.97269 1.053 1.0745 0.837778 0.939272 0.982485 0.930167 1.01621 1.00659 0.953345 1.27029 1.20318 1.20099 1.30424 1.14345 1.28908 1.12478 1.2406 1.26179 1.37459 1.10736 1.09451 1.23036 1.27025 1.18798 1.58782 1.53034 1.25373 1.45668 1.24997 1.42444 1.40855 1.53083 1.46153 1.46335 1.34096 1.57503 1.6691 1.54312 1.5268 1.33039 1.14996 1.13741 1.18608 1.33171 1.4816 1.80321 1.56032 1.31092 1.25725 1.34331 1.50961 1.59761 1.51189 1.68887 1.66549 1.70046 1.6905 1.43731 1.88969 1.9938 1.85763 1.82253 1.56803 1.72826 1.87631 1.55828 2.02148 1.5941 1.6304 1.52026 1.46497 2.0549 2.14404 2.20967 2.12319 2.08626 1.5375 1.57057 2.05365 2.23678 2.28216 1.82742 1.76307 1.89616 1.75842 1.81459 1.93296 1.92855 1.75312 1.7307 1.87325 1.89076 2.08431 2.32944 1.6058 1.63434 1.81077 2.3379 2.40412 0.806985 0.952027 1.11007 1.13192 1.00907 0.831896 2.28255 2.35405 2.19875 2.12287 1.93347 1.89043 1.72615 1.88708 2.24988 2.10249 2.02642 1.71714 1.82191 1.67498 2.41639 2.26481 2.11909 2.25432 1.97578 2.16801 2.25329 2.11207 1.96435 1.83932 2.24471 2.40847 1.97985 2.04304 2.14272 2.23998 2.18252 2.00834 2.41712 2.34817 2.31939 2.47805 2.48613 2.4481 1.92292 2.22931 2.48812 2.40975 2.43435 2.25728 2.11461 1.97793 2.52602 2.52415 2.55413 2.57719 2.48855 2.56103 2.23722 2.18549 2.4768 2.39244 2.0955 2.40275 2.31311 2.17136 2.06945 1.9838 2.06428 2.07279 2.13286 2.55835 2.58159 2.46662 2.40993 2.4794 2.55997 2.55421 2.46438 2.5444 2.52875 2.42529 2.30104 2.32541 1.98693 1.95463 2.13852 2.21592 2.11871 1.93134 1.92075 1.87553 2.51072 2.41458 2.26113 2.25501 2.41627 2.52149 1.8579 1.93116 1.82598 1.71718 1.70748 0.848009 0.935968 1.03321 0.944095 1.81543 2.41923 2.02 2.44573 2.53607 2.57528 2.26661 1.99966 2.50777 2.2575 2.50095 2.19157 2.30487 2.556 2.5664 2.06613 2.35891 2.5144 2.54528 2.35235 2.53268 2.37988 2.12806 2.47392 2.28261 1.95878 2.0433 2.07101 2.33426 2.35843 2.39052 1.80358 2.08713 1.82271 2.06295 0.978232 2.43053 2.32613 1.79583 2.33252 2.2334 2.00508 1.80971 2.32761 1.91035 1.99718 1.97552 2.33834 2.15245 1.70089 2.01789 2.22526 2.1731 1.70349 1.99253 1.59201 1.80607 1.62738 1.66567 1.76492 1.42381 1.67232 1.90509 2.01583 1.33412 1.56516 1.38195 1.69517 1.37503 1.20865 1.41488 1.24411 1.58827 1.49578 1.40368 1.38017 1.64656 1.27087 1.71601 1.11123 1.43896 1.1853 1.1141 0.983238 1.18515 1.02856 1.07072 1.13856 1.01757 1.027 1.45215 1.90873 1.07791 1.10304 0.757211 0.744648 0.832295 0.790778 0.794591 0.657492 0.590052 0.629502 0.524636 0.540425 0.61166 0.787705 0.809366 0.911711 0.911612 0.463617 0.324233 0.310889 0.855359 0.593239 0.822955 0.319989 0.341667 0.602168 0.593494 0.339723 \ No newline at end of file +2.40859 2.42069 1.99303 2.45586 2.16306 0.466146 1.40423 2.02177 1.37759 1.29686 2.32016 2.22623 1.59334 1.47611 1.41387 1.95853 2.23253 2.07792 1.14877 2.03789 1.9058 1.10339 1.8427 0.676623 1.38191 1.75157 1.53567 2.26161 0.956311 0.966318 1.01918 0.909112 0.856988 0.857398 0.900106 1.67252 1.75152 1.7529 1.79088 1.89184 1.85869 1.79385 1.75354 1.93707 1.86946 2.49051 2.46053 2.31417 2.36798 2.31475 2.31693 2.37013 2.49597 2.53526 2.47881 2.45904 1.87341 2.20969 2.11226 2.2306 2.11977 1.93475 2.03631 2.04817 1.96867 2.26966 2.36125 2.45637 2.37304 2.33487 2.39411 2.49409 2.47243 2.55713 2.51288 2.41525 2.38257 2.52804 2.49078 2.55158 2.53629 2.53868 2.56266 2.57115 2.55841 2.50065 2.52058 2.31186 2.32399 2.4169 2.48524 2.52681 2.5138 2.39764 2.48009 2.57525 2.56793 2.57967 2.57047 2.54208 2.57071 2.20628 2.23106 2.19544 1.98421 1.9949 2.05896 2.15363 2.12532 2.05571 1.979 1.99809 2.01907 2.00641 2.08417 2.16842 2.24508 2.25381 2.37373 2.44083 2.37051 2.33289 2.14266 2.20605 2.38898 2.45671 2.51762 2.53333 2.40858 2.46147 2.14113 2.0896 2.28657 2.23092 2.30878 2.32897 2.13647 2.21471 2.50776 2.43884 2.52585 2.54065 2.57189 2.57443 2.57683 2.57969 2.56663 2.54772 2.55659 2.56349 2.55781 2.54686 2.48776 2.54314 1.9958 2.02021 2.04234 2.17967 2.28199 2.13847 2.42305 2.49192 2.50644 2.45335 2.36339 2.35256 2.47501 2.41813 2.47642 2.41425 2.53744 2.52716 2.25183 2.1448 2.32456 2.33129 1.99522 2.46151 2.51322 2.50295 2.51245 2.47894 2.51096 2.29513 2.34527 2.32223 2.35741 2.40919 2.44255 2.3725 2.38875 2.09904 2.11709 2.27275 2.28665 2.24331 2.19081 2.21423 2.18192 2.11514 2.04941 1.98537 2.0589 2.46512 2.41657 2.39298 2.44442 2.29761 2.31266 2.17496 2.24182 1.83326 2.02998 2.01973 2.2204 2.12095 2.33578 2.33914 2.24512 2.15591 2.17084 2.25039 2.25441 2.16538 1.99095 1.9078 2.07963 2.07652 2.37014 2.36009 2.14839 2.34855 2.35541 1.84708 1.91335 1.7272 1.76054 1.72434 1.6534 1.7805 1.77951 2.01338 2.20135 1.98987 1.90169 1.99757 1.74498 1.6683 1.99087 1.94407 1.95811 2.00919 2.03228 2.21199 2.10861 2.11019 2.18658 2.43135 2.2925 2.22415 2.18448 2.36715 2.2994 0.921016 0.943901 1.09546 1.03911 1.17128 1.07429 1.04511 1.0134 0.876638 0.767961 0.863664 0.916852 2.39376 2.35994 2.3672 2.39864 1.90683 1.7022 1.74251 1.62829 1.56003 1.67271 2.36538 2.36328 2.27406 2.30938 2.16464 2.08718 2.16942 1.9485 1.69115 1.75296 2.02227 1.93796 1.86885 1.96905 1.85376 1.89308 1.83237 1.74451 1.75432 1.83346 1.99633 1.99173 1.80452 1.9043 1.78479 1.9187 1.85908 1.74748 2.2902 2.3243 2.31181 2.23751 2.29515 2.15917 2.14591 1.67918 1.66245 1.59607 1.4898 1.59297 1.63189 2.07354 2.12954 2.25118 2.23118 2.17587 2.13621 2.08986 2.20208 2.17537 2.14874 2.08169 2.13715 1.44875 1.65314 1.583 1.53536 1.51111 1.66272 1.58971 2.08589 2.05227 2.10906 1.65424 1.61819 1.80663 1.69011 1.83085 1.76312 1.60197 1.5061 1.49281 1.56286 1.63139 1.65948 1.9617 1.81894 1.89735 1.94828 2.03783 1.90435 1.81724 1.93696 1.97262 1.51643 1.44627 1.49163 1.53784 1.72372 1.77838 1.77939 1.7342 1.65573 1.7386 1.73664 1.76618 1.77155 1.68214 1.59522 1.61425 1.69817 1.60519 1.52902 1.59464 1.51047 1.67987 1.67601 1.59695 1.51335 1.41873 1.4274 1.26842 1.34222 1.41759 1.4249 1.3315 1.58535 1.4839 1.61391 1.64984 1.87165 1.81139 1.78224 1.85143 1.41377 1.23572 1.28436 1.11923 1.22059 1.1015 1.63605 1.56124 1.44706 1.53121 1.5974 1.7319 1.73759 1.64764 1.54837 1.4935 1.65526 1.61082 1.41521 1.46847 1.5319 1.50437 1.62565 1.59893 1.44043 1.5051 1.52111 1.45616 1.37077 1.46936 1.34843 1.31064 1.20984 1.14583 1.46275 1.55703 1.54689 1.4555 1.29772 1.19223 1.35931 1.36659 1.60386 1.5149 1.6101 1.52475 1.50057 1.55061 1.63655 1.6735 1.26593 1.23053 1.34254 1.38705 1.33082 1.26316 1.1371 1.18163 1.19009 1.05688 1.13846 1.4515 1.44519 1.29613 1.35584 1.33296 1.25427 1.20696 1.30827 1.2133 1.16483 1.27431 1.35545 1.37122 1.30322 1.1747 1.09224 1.20363 1.22706 1.41469 1.37369 1.37265 1.21425 1.30594 1.27381 1.17723 1.29725 1.20674 1.18131 1.22749 1.32525 1.36381 0.99314 0.895302 1.01117 1.05181 1.10974 1.05492 1.06791 0.954517 1.13588 1.00348 0.902975 1.05319 1.07678 1.03979 1.00584 0.924274 0.857178 0.838316 0.919027 1.14097 1.17011 1.14528 1.09003 1.01359 1.10807 1.04498 1.0437 0.969828 0.895438 0.906761 0.969393 0.97212 0.890419 0.875994 0.965695 1.13218 1.1118 0.987432 1.06901 0.900416 0.963048 1.03413 1.06441 0.859519 0.947445 0.958639 0.997954 0.963872 0.903585 0.862093 0.942282 0.974858 0.865043 0.955997 1.01493 1.03529 0.969639 0.879974 0.994124 0.806821 0.857301 0.795112 0.76853 0.588516 0.66225 0.617074 0.46742 0.463771 0.745808 0.762521 0.571364 0.646343 0.834444 0.856978 0.800609 0.641365 0.653902 0.729197 0.640529 0.611183 0.651506 0.657849 0.606901 0.635317 0.67961 0.746359 0.826262 0.721338 0.623145 0.662814 0.737936 0.602237 0.692315 0.783984 0.78705 0.726866 0.665266 0.644363 0.709107 0.696187 0.804021 0.861776 0.822207 0.670157 0.740309 0.810351 0.832509 0.736739 0.776815 0.416347 0.403869 0.577968 0.485719 0.389377 0.416903 0.502974 0.416868 0.595791 0.587723 0.488449 0.544219 0.714945 0.792584 0.840766 0.827026 0.756213 0.702387 0.404749 0.471008 0.535968 0.561959 0.453718 0.575201 0.552031 0.482822 0.411567 0.382131 0.463977 0.424512 0.811639 0.821861 0.756482 0.683141 0.538913 0.488643 0.513511 0.560952 0.58337 0.586828 0.41882 0.367017 0.362679 0.417801 0.520671 0.526956 0.288852 0.281494 0.289711 0.0942476 0.207717 0.0973594 0.215389 0.294184 0.283215 0.263503 0.266901 0.0922615 0.200502 0.602362 0.601147 0.559238 0.556836 0.597546 0.60068 0.292273 0.297969 0.215011 0.0974517 0.0948727 0.207196 0.635848 0.63623 0.660427 0.713264 0.738437 0.683064 1.39401 1.42125 1.37743 1.28879 1.23998 1.30855 1.43454 1.49138 1.42302 2.26013 2.18604 2.10263 2.09048 2.17216 2.26692 1.99481 2.03609 1.97572 1.88764 1.85469 1.91015 1.96103 1.91495 1.91829 1.98653 2.05636 2.03675 2.05428 2.04674 2.1415 2.20365 2.15276 2.07244 2.1509 2.15328 2.07189 1.16289 1.22872 1.27569 1.21043 1.78402 1.81036 1.8792 0.638157 0.572696 0.641422 0.702138 0.475627 0.559772 0.484245 0.705091 0.637532 0.5764 0.645757 1.17618 1.18222 1.1157 1.33794 1.33426 1.24654 1.18857 1.24633 1.83904 1.89903 1.95996 1.90173 2.29408 2.20902 2.20754 2.29389 1.89181 1.8945 1.8684 2.01069 1.98722 2.07688 1.99171 1.98491 2.03513 1.90775 1.92414 1.94815 1.96501 1.94982 1.1855 1.24079 1.20707 1.9439 1.89656 1.95992 2.00388 2.43651 2.43245 2.45693 2.46533 2.04787 2.10773 2.02825 1.99085 1.97563 1.96533 1.98784 1.1353 1.03524 1.02071 1.09492 1.16813 1.94941 1.98308 2.07126 2.09439 2.0233 1.90262 1.87612 1.80151 1.82191 1.69379 1.71607 1.8064 1.77954 1.85069 1.92131 1.89726 1.83177 2.29793 2.20822 2.27242 0.930677 0.906179 0.821882 0.849448 0.509934 0.548438 0.613512 0.58519 1.87495 1.88107 1.8243 1.79784 1.84237 1.93149 1.99114 1.94532 1.85787 1.84954 1.36274 1.41643 2.2999 2.21104 2.28779 1.9863 1.96043 1.94078 2.45187 2.41154 2.41545 2.46386 2.37517 2.42473 1.70392 1.71651 1.7948 1.77639 1.97454 1.92507 2.04656 1.81311 1.86124 1.81338 1.76171 0.462452 0.406968 0.403334 0.459287 0.379784 0.363519 0.276845 0.296857 0.614411 0.626711 0.714204 0.70578 1.12179 1.09629 1.00991 1.03663 1.28599 1.25728 1.33893 1.36408 0.653365 0.56334 0.514732 0.583156 0.665286 0.610661 0.712619 0.740626 0.648766 0.574461 0.657646 0.595651 1.30992 1.0662 0.664572 0.329527 0.428537 1.81096 1.74899 1.91401 1.84057 0.563963 0.875866 1.87959 1.74764 1.85011 2.02174 1.08788 1.97705 2.45099 1.95067 1.93357 2.03814 2.25021 1.89852 1.26739 0.642227 0.53429 0.639755 1.21813 2.11089 2.11427 1.97073 1.93993 2.18268 1.40687 1.33487 0.678511 0.184162 0.583467 0.167915 0.179722 0.178038 0.187225 0.508262 0.341121 0.54484 0.74929 0.466582 0.377821 0.476812 0.459126 0.374403 0.770702 0.546569 0.466401 0.388917 0.490698 0.391096 0.521396 0.681341 0.745317 0.661602 0.764545 0.718939 0.623743 0.645095 0.755011 0.62248 0.680102 0.610338 0.623412 0.696206 0.687785 0.625295 0.744738 0.6678 0.532793 0.539812 0.535555 0.543322 0.691639 0.732928 0.922597 0.909668 0.886712 0.908162 0.960647 0.868711 1.04681 0.890709 0.899559 0.97269 1.053 1.0745 0.837778 0.939272 0.982485 0.930167 1.01621 1.00659 0.953345 1.27029 1.20318 1.20099 1.30424 1.14345 1.28908 1.12478 1.2406 1.26179 1.37459 1.10736 1.09451 1.23036 1.27025 1.18798 1.58782 1.53034 1.25373 1.45668 1.24997 1.42444 1.40855 1.53083 1.46153 1.46335 1.34096 1.57503 1.6691 1.54312 1.5268 1.33039 1.14996 1.13741 1.18608 1.33171 1.4816 1.80321 1.56032 1.31092 1.25725 1.34331 1.50961 1.59761 1.51189 1.68887 1.66549 1.70046 1.6905 1.43731 1.88969 1.9938 1.85763 1.82253 1.56803 1.72826 1.87631 1.55828 2.02148 1.5941 1.6304 1.52026 1.46497 2.0549 2.14404 2.20967 2.12319 2.08626 1.5375 1.57057 2.05365 2.23678 2.28216 1.82742 1.76307 1.89616 1.75842 1.81459 1.93296 1.92855 1.75312 1.7307 1.87325 1.89076 2.08431 2.32944 1.6058 1.63434 1.81077 2.3379 2.40412 0.806985 0.952027 1.11007 1.13192 1.00907 0.831896 2.28255 2.35405 2.19875 2.12287 1.93347 1.89043 1.72615 1.88708 2.24988 2.10249 2.02642 1.71714 1.82191 1.67498 2.41639 2.26481 2.11909 2.25432 1.97578 2.16801 2.25329 2.11207 1.96435 1.83932 2.24471 2.40847 1.97985 2.04304 2.14272 2.23998 2.18252 2.00834 2.41712 2.34817 2.31939 2.47805 2.48613 2.4481 1.92292 2.22931 2.48812 2.40975 2.43435 2.25728 2.11461 1.97793 2.52602 2.52415 2.55413 2.57719 2.48855 2.56103 2.23722 2.18549 2.4768 2.39244 2.0955 2.40275 2.31311 2.17136 2.06945 1.9838 2.06428 2.07279 2.13286 2.55835 2.58159 2.46662 2.40993 2.4794 2.55997 2.55421 2.46438 2.5444 2.52875 2.42529 2.30104 2.32541 1.98693 1.95463 2.13852 2.21592 2.11871 1.93134 1.92075 1.87553 2.51072 2.41458 2.26113 2.25501 2.41627 2.52149 1.8579 1.93116 1.82598 1.71718 1.70748 0.848009 0.935968 1.03321 0.944095 1.81543 2.41923 2.02 2.44573 2.53607 2.57528 2.26661 1.99966 2.50777 2.2575 2.50095 2.19157 2.30487 2.556 2.5664 2.06613 2.35891 2.5144 2.54528 2.35235 2.53268 2.37988 2.12806 2.47392 2.28261 1.95878 2.0433 2.07101 2.33426 2.35843 2.39052 1.80358 2.08713 1.82271 2.06295 0.978232 2.43053 2.32613 1.79583 2.33252 2.2334 2.00508 1.80971 2.32761 1.91035 1.99718 1.97552 2.33834 2.15245 1.70089 2.01789 2.22526 2.1731 1.70349 1.99253 1.59201 1.80607 1.62738 1.66567 1.76492 1.42381 1.67232 1.90509 2.01583 1.33412 1.56516 1.38195 1.69517 1.37503 1.20865 1.41488 1.24411 1.58827 1.49578 1.40368 1.38017 1.64656 1.27087 1.71601 1.11123 1.43896 1.1853 1.1141 0.983238 1.18515 1.02856 1.07072 1.13856 1.01757 1.027 1.45215 1.90873 1.07791 1.10304 0.757211 0.744648 0.832295 0.790778 0.794591 0.657492 0.590052 0.629502 0.524636 0.540425 0.61166 0.787705 0.809366 0.911711 0.911612 0.463617 0.324233 0.310889 0.855359 0.593239 0.822955 0.319989 0.341667 0.602168 0.593494 0.339723 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 53a22f8fd5..bae5a0ee7f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -36,6 +36,7 @@ add_gudhi_module(Hasse_complex) add_gudhi_module(Persistence_representations) add_gudhi_module(Persistent_cohomology) add_gudhi_module(Rips_complex) +add_gudhi_module(Ripser) add_gudhi_module(Simplex_tree) add_gudhi_module(Skeleton_blocker) add_gudhi_module(Spatial_searching) diff --git a/src/Ripser/include/gudhi/ripser.h b/src/Ripser/include/gudhi/ripser.h new file mode 100644 index 0000000000..1e712b397e --- /dev/null +++ b/src/Ripser/include/gudhi/ripser.h @@ -0,0 +1,1419 @@ +/* Based on Ripser commit 140670f2c76997404601e43d8054151f46be9fd7 + * Modification(s): + * - YYYY/MM Author: Description of the modification + * - 2024 Marc Glisse: Heavy refactoring +*/ + +/* + + Ripser: a lean C++ code for computation of Vietoris-Rips persistence barcodes + + MIT License + + Copyright (c) 2015–2021 Ulrich Bauer + + 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. + + You are under no obligation whatsoever to provide any bug fixes, patches, or + upgrades to the features, functionality or performance of the source code + ("Enhancements") to anyone; however, if you choose to make your Enhancements + available either publicly, or directly to the author of this software, without + imposing a separate written license agreement for such Enhancements, then you + hereby grant the following license: a non-exclusive, royalty-free perpetual + license to install, use, modify, prepare derivative works, incorporate into + other computer software, distribute, and sublicense such enhancements or + derivative works thereof, in binary and source code form. + +*/ + +//#define GUDHI_INDICATE_PROGRESS + +// #define GUDHI_RIPSER_USE_BOOST_HEAP // only useful when heap::push dominates? +// #define GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT // only useful in cases where edge-collapse is more important? + +/* TODO: + * from branch representative-cocycles + - for vertices: brute force check all vertices that are in the same connected component (could do better, but there is already quadraticness elsewhere in most cases, even the output for dim 0 may be quadratic) + - for others: print_chain just calls iteratively get_pivot on working_reduction_column + * from branch representative-cycles + - dim 0: trivial + - parametrize a number of functions (like add_coboundary) by the (co)boundary iterator so they can be also used for homology + - once cohomology is computed for some dim, assemble the relevant simplices (don't forget the essential ones) and reduce them homology-style. I think we have 2 choices: reduce the birth and check which columns we add (only possibility for essential classes, but do we really need to do cohomology first if we are going to do that?), or reduce the death and look at the column after reduction (Ripser(er) chose this). + * check out the persistence image branch + + * allow non-0 filtration value on vertices, so we can handle all flag-type filtrations, not just plain Rips. cf scikit-tda +*/ + +#ifndef GUDHI_RIPSER_H +#define GUDHI_RIPSER_H + +#include +#include // sqrt +#include +#include +#include +#include +#ifdef GUDHI_INDICATE_PROGRESS +#include +#endif + +#include +#include +#if BOOST_VERSION >= 108100 +#include +#else +#include +#endif + +#ifdef GUDHI_RIPSER_USE_BOOST_HEAP +#include +#endif + +#include +#include + + +namespace Gudhi::ripser { + +#define GUDHI_assert(X) GUDHI_CHECK(X, std::logic_error("")) + +#ifdef GUDHI_RIPSER_USE_BOOST_HEAP +template +using Heap = boost::heap::d_ary_heap, boost::heap::compare>; +#else +template +struct Heap : std::priority_queue { + typedef std::priority_queue Base; + using Base::Base; + void clear() { this->c.clear(); } +}; +#endif + +template +class Union_find { + public: + typedef vertex_t_ vertex_t; + private: + std::vector parent; + std::vector rank; + public: + Union_find(const vertex_t n) : parent(n), rank(n, 0) { + for (vertex_t i = 0; i < n; ++i) parent[i] = i; + } + + // TODO: check which one is faster (probably doesn't matter) + vertex_t find(vertex_t x) { +#if 0 + vertex_t y = x, z; + while ((z = parent[y]) != y) y = z; + while ((z = parent[x]) != y) { + parent[x] = y; + x = z; + } + return z; +#else + // Path halving + vertex_t y = parent[x]; + vertex_t z = parent[y]; + while (y != z) + { + parent[x] = z; + x = z; + y = parent[x]; + z = parent[y]; + } + return y; +#endif + } + + void link(vertex_t x, vertex_t y) { + // this line is redundant, the caller already does it + if ((x = find(x)) == (y = find(y))) return; + if (rank[x] > rank[y]) + parent[y] = x; + else { + parent[x] = y; + if (rank[x] == rank[y]) ++rank[y]; + } + } +}; + +template +bool is_prime(const coefficient_t n) { + if (!(n & 1) || n < 2) return n == 2; + for (coefficient_t p = 3; p * p <= n; p += 2) + if (!(n % p)) return false; + return true; +} + +// For pretty printing, modulo 11, we prefer -1 to 10. +template +coefficient_t normalize(const coefficient_t n, const coefficient_t modulus) { + return n > modulus / 2 ? n - modulus : n; +} + +template +std::vector multiplicative_inverse_vector(const coefficient_t m) { + std::vector inverse(m); + if (!is_prime(m)) + throw std::domain_error("Modulus must be a prime number"); + if ((m - 1) != (coefficient_storage_t)(m - 1)) + throw std::overflow_error("Modulus is too large"); + inverse[1] = 1; + // m = a * (m / a) + m % a + // Multipying with inverse(a) * inverse(m % a): + // 0 = inverse(m % a) * (m / a) + inverse(a) (mod m) + for (coefficient_t a = 2; a < m; ++a) inverse[a] = m - (inverse[m % a] * (m / a)) % m; + return inverse; +} + +#ifdef GUDHI_INDICATE_PROGRESS +constexpr std::chrono::milliseconds time_step(40); +constexpr const char* clear_line="\r\033[K"; +#endif + +/* concept DistanceMatrix +struct DistanceMatrix { + typedef Category; + typedef vertex_t; + typedef value_t; + value_t operator()(vertex_t i, vertex_t j) const; + vertex_t size() const; + // Optional: + static const bool is_sparse; + vector> neighbors; +}; +*/ + +class Tag_dense {}; // Use directly, iterate on vertices +class Tag_sparse {}; // Use directly, iterate on edges +class Tag_other {}; // Could use directly as dense, but prefer converting it. + +template +struct Full_distance_matrix { + public: + typedef Tag_dense Category; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::value_t value_t; + std::vector distances; // TODO: private + private: + vertex_t n; + public: + //Full_distance_matrix(std::vector&& _distances) + // : distances(std::move(_distances)), n(std::sqrt(_distances.size())) {} + + template + Full_distance_matrix(const DistanceMatrix& mat) + : distances(static_cast(mat.size()) * mat.size()), n(mat.size()) { // vertex_t is meant for vertices. Using it for edges could be unsafe, so we cast to size_t. + for (vertex_t i = 0; i < size(); ++i) + for (vertex_t j = 0; j < size(); ++j) + distances[i * n + j] = mat(i, j); // If mat::operator() involves computations, should we try to take advantage of the symmetry? + } + + // Confusing: why is the order j,i significantly faster than i,j? + value_t operator()(const vertex_t j, const vertex_t i) const { + return distances[i * n + j]; + } + vertex_t size() const { return n; } +}; + +enum Compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; + +template +struct Compressed_distance_matrix { + public: + typedef Tag_dense Category; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::value_t value_t; + std::vector distances; // TODO: private + private: + std::vector rows; // Surprisingly, this is more efficient than computing i*(i-1)/2 + + public: + Compressed_distance_matrix(std::vector&& _distances) + : distances(std::move(_distances)), rows((1 + std::sqrt(1 + 8 * distances.size())) / 2) { + GUDHI_assert(distances.size() == (std::size_t)size() * (size() - 1) / 2); + init_rows(); + } + + template + Compressed_distance_matrix(const DistanceMatrix& mat) + : distances(static_cast(mat.size()) * (mat.size() - 1) / 2), rows(mat.size()) { // vertex_t is meant for vertices. Using it for edges could be unsafe, so we cast to size_t. + init_rows(); + + for (vertex_t i = 1; i < size(); ++i) + for (vertex_t j = 0; j < i; ++j) rows[i][j] = mat(i, j); + } + + value_t operator()(const vertex_t i, const vertex_t j) const { + if (i == j) return 0; + if ((Layout == LOWER_TRIANGULAR) ? (i < j) : (i > j)) + return rows[j][i]; + else + return rows[i][j]; + } + vertex_t size() const { return rows.size(); } + private: + void init_rows() { + if constexpr (Layout == LOWER_TRIANGULAR) { + value_t* pointer = &distances[0]; + for (vertex_t i = 1; i < size(); ++i) { + rows[i] = pointer; + pointer += i; + } + } else { // UPPER_TRIANGULAR + value_t* pointer = &distances[0] - 1; + for (vertex_t i = 0; i < size() - 1; ++i) { + rows[i] = pointer; + pointer += size() - i - 2; + } + } + } +}; + +template +struct Sparse_distance_matrix { + public: + typedef Tag_sparse Category; + static constexpr bool is_sparse = true; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::value_t value_t; + struct vertex_diameter_t { + vertex_diameter_t() =default; + vertex_diameter_t(vertex_t i_, value_t d_) : i(i_), d(d_) {} +#if 1 + // gcc is faster with those 2 lines (or with a std::pair) :-( + vertex_diameter_t(vertex_diameter_t const&o)noexcept :i(o.i),d(o.d){} + vertex_diameter_t&operator=(vertex_diameter_t const&o) noexcept{i=o.i;d=o.d;return*this;} +#endif + vertex_t i; value_t d; + friend vertex_t get_vertex(const vertex_diameter_t& i) { return i.i; } + friend value_t get_diameter(const vertex_diameter_t& i) { return i.d; } + friend bool operator<(vertex_diameter_t const& a, vertex_diameter_t const& b) { + if (a.i < b.i) return true; + if (a.i > b.i) return false; + return a.d < b.d; + } + }; + + std::vector> neighbors; + std::size_t num_edges; // TODO: useless, remove + + private: +#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT +#if BOOST_VERSION >= 108100 + boost::unordered_flat_map,value_t> m; +#else + boost::unordered_map,value_t> m; +#endif + // Would a vector (one map per vertex) have any advantage? +#endif + + public: + Sparse_distance_matrix(std::vector>&& _neighbors, + std::size_t _num_edges = 0) + : neighbors(std::move(_neighbors)), num_edges(_num_edges) {init();} + + template + Sparse_distance_matrix(const DistanceMatrix& mat, const value_t threshold) + : neighbors(mat.size()), num_edges(0) { + + for (vertex_t i = 0; i < size(); ++i) + for (vertex_t j = 0; j < size(); ++j) + if (i != j) { + auto d = mat(i, j); + if (d <= threshold) { + ++num_edges; + neighbors[i].emplace_back(j, d); + } + } + init(); + } + + value_t operator()(const vertex_t i, const vertex_t j) const { +#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT + // We could insert in both orders to save the minmax for each query. + return m.at(std::minmax(i,j)); + //// We never hit the infinity case? + // auto it = m.find(std::minmax(i,j)); + // return (it != m.end()) ? *it : std::numeric_limits::infinity(); +#else + auto neighbor = + std::lower_bound(neighbors[i].begin(), neighbors[i].end(), vertex_diameter_t{j, 0}); + return (neighbor != neighbors[i].end() && get_vertex(*neighbor) == j) + ? get_diameter(*neighbor) + : std::numeric_limits::infinity(); +#endif + } + vertex_t size() const { return neighbors.size(); } + private: + void init() { +#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT + for(vertex_t i=0; i +struct Euclidean_distance_matrix { + public: + typedef Tag_other Category; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::value_t value_t; + + Euclidean_distance_matrix(std::vector>&& _points) + : points(std::move(_points)) { + for (auto p : points) { GUDHI_assert(p.size() == points.front().size()); } + } + + value_t operator()(const vertex_t i, const vertex_t j) const { + GUDHI_assert((std::size_t)i < points.size()); + GUDHI_assert((std::size_t)j < points.size()); + return std::sqrt(std::inner_product( + points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus(), + [](value_t u, value_t v) { return (u - v) * (u - v); })); + } + + vertex_t size() const { return points.size(); } + private: + std::vector> points; +}; + +// The gratuitous restrictions on what can be specialized in C++ are annoying. +template > struct Is_sparse_impl : std::bool_constant {}; +template struct Is_sparse_impl> : std::bool_constant {}; +template constexpr bool is_sparse () { return Is_sparse_impl::value; } + +template class Compressed_sparse_matrix_ { + std::vector bounds; + std::vector entries; + + typedef typename std::vector::const_iterator iterator; + typedef boost::iterator_range iterator_pair; + + public: + iterator_pair subrange(const size_t index) const { + return {entries.begin() + (index == 0 ? 0 : bounds[index - 1]), + entries.begin() + bounds[index]}; + } + + // Close the current column, the next push_back will be for a new column + void append_column() { bounds.push_back(entries.size()); } + + void push_back(const ValueType e) { + GUDHI_assert(0 < bounds.size()); + entries.push_back(e); + ++bounds.back(); + } +}; + +/* concept SimplexEncoding +struct SimplexEncoding { + typedef dimension_t; + typedef vertex_t; + typedef simplex_t; // has to be an integer + SimplexEncoding(vertex_t n_vertices, dimension_t max_dim_plus_one); // max_vertices_per_simplex + simplex_t operator()(vertex_t v, dimension_t position) const; // position starts at 1? + vertex_t get_max_vertex(simplex_t s, dimension_t position, vertex_t n_vertices) const; + int num_extra_bits() const; +}; +*/ + +// number of bits necessary to store x with 0 <= x < n +template +constexpr int log2up(vertex_t n) { + --n; + int k = 0; + while(n>0) { n>>=1; ++k; } + return k; +} + +template +class Cns_encoding { + public: + typedef typename Params::dimension_t dimension_t; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::simplex_t simplex_t; + + private: + std::vector> B; // table of binomial coefficients + int extra_bits; + + public: + Cns_encoding(vertex_t n, dimension_t k) : B(k + 1, std::vector(n + 1, 0)) { + static_assert(std::numeric_limits::radix == 2); + const int available_bits = std::numeric_limits::digits; + simplex_t max_simplex_index = 0; + for (vertex_t i = 0; i <= n; ++i) { + B[0][i] = 1; + for (dimension_t j = 1; (vertex_t)j < std::min(i, k + 1); ++j) + B[j][i] = B[j - 1][i - 1] + B[j][i - 1]; + if (i <= k) B[i][i] = 1; + vertex_t mi = std::min(i >> 1, (vertex_t)k); // max + max_simplex_index = B[mi][i]; + if (i > 1 && max_simplex_index < B[mi][i-1]) { // overflow + throw std::overflow_error("cannot encode all simplices of dimension " + std::to_string(k) + " with " + std::to_string(n) + " vertices using only " + std::to_string(available_bits) + " bits"); + } + } + extra_bits = available_bits - log2up(max_simplex_index + 1); + } + + simplex_t operator()(vertex_t n, dimension_t k) const { + GUDHI_assert(n >= k - 1); + return B[k][n]; + } + + // We could get `n` from B and avoid passing it as argument + vertex_t get_max_vertex(const simplex_t idx, const dimension_t k, const vertex_t n) const { + return get_max(n, k - 1, [&](vertex_t w) -> bool { return (*this)(w, k) <= idx; }); + } + + int num_extra_bits() const { return extra_bits; } + + private: + template + static vertex_t get_max(vertex_t top, const vertex_t bottom, const Predicate pred) { + if (!pred(top)) { + vertex_t count = top - bottom; + while (count > 0) { + vertex_t step = count >> 1, mid = top - step; + if (!pred(mid)) { + top = mid - 1; + count -= step + 1; + } else + count = step; + } + } + return top; + } +}; + +template +class Bitfield_encoding { + public: + typedef typename Params::dimension_t dimension_t; + typedef typename Params::vertex_t vertex_t; + typedef typename Params::simplex_t simplex_t; + + private: + int bits_per_vertex; + int extra_bits; + + public: + Bitfield_encoding(vertex_t n, dimension_t k) : bits_per_vertex(log2up(n)) { + static_assert(std::numeric_limits::radix == 2); + const int available_bits = std::numeric_limits::digits; + extra_bits = available_bits - bits_per_vertex * k; + if (extra_bits < 0) + throw std::overflow_error("cannot encode all simplices of dimension " + std::to_string(k - 1) + " with " + std::to_string(n) + " vertices using only " + std::to_string(available_bits) + " bits"); + // The message is a bit misleading, it is tuples that we cannot encode, and just with this representation. + } + + simplex_t operator()(vertex_t n, dimension_t k) const { + if(k==0) return 1; // because of odd use in (co)boundary... + --k; + // USE_N_MINUS_K only useful if it somehow helps remove the test k==0 above +#ifdef USE_N_MINUS_K + return (simplex_t)(n - k) << (bits_per_vertex * k); +#else + return (simplex_t)n << (bits_per_vertex * k); +#endif + } + + vertex_t get_max_vertex(const simplex_t idx, dimension_t k, const vertex_t) const { + GUDHI_assert(k > 0); + --k; +#ifdef USE_N_MINUS_K + return static_cast(idx >> (bits_per_vertex * k)) + k; +#else + return static_cast(idx >> (bits_per_vertex * k)); +#endif + } + + int num_extra_bits() const { return extra_bits; } +}; + +template struct Rips_filtration { + using size_t = typename Params::size_t; + using vertex_t = typename SimplexEncoding::vertex_t; + // static_assert(std::is_same_v); // too strict + using simplex_t = typename SimplexEncoding::simplex_t; + using dimension_t = typename SimplexEncoding::dimension_t; + using value_t = typename DistanceMatrix::value_t; + using coefficient_storage_t = typename Params::coefficient_storage_t; + using coefficient_t = typename Params::coefficient_t; + static constexpr bool use_coefficients = Params::use_coefficients; + + // The definition of entry_t could be added in some intermediate layer between SimplexEncoding and here + struct entry_with_coeff_t { + simplex_t content; + entry_with_coeff_t(simplex_t _index, coefficient_t _coefficient, int bits_for_coeff) + : content((_index << bits_for_coeff) | (_coefficient - 1)) { GUDHI_assert(_coefficient != 0); } + entry_with_coeff_t() {} + // We never store a coefficient of 0, so we can store coef-1 so %2 requires 0 bit and %3 only 1 bit + friend const entry_with_coeff_t& get_entry(const entry_with_coeff_t& e) { return e; } + }; + simplex_t get_index(const entry_with_coeff_t& e) const { return e.content >> num_bits_for_coeff(); } + coefficient_t get_coefficient(const entry_with_coeff_t& e) const { return static_cast(e.content & (((simplex_t)1 << num_bits_for_coeff()) - 1)) + 1; } + void set_coefficient(entry_with_coeff_t& e, const coefficient_t c) const { GUDHI_assert(c!=0); e.content = (e.content & ((simplex_t)(-1) << num_bits_for_coeff())) | (c - 1); } + // Should we cache the masks derived from num_bits_for_coeff? + + struct entry_plain_t { + simplex_t index; + entry_plain_t(simplex_t _index, coefficient_t, int) : index(_index) {} + entry_plain_t() {} + friend const entry_plain_t& get_entry(const entry_plain_t& e) { return e; } + }; + static simplex_t get_index(const entry_plain_t& i) { return i.index; } + static coefficient_t get_coefficient(const entry_plain_t& i) { return 1; } + static void set_coefficient(entry_plain_t& e, const coefficient_t c) { GUDHI_assert(c==1); } + + typedef std::conditional_t entry_t; + entry_t make_entry(simplex_t i, coefficient_t c) const { return entry_t(i, c, num_bits_for_coeff()); } + + static_assert(sizeof(entry_t) == sizeof(simplex_t), "size of entry_t is not the same as simplex_t"); + + // TODO: avoid storing filtp when !use_coefficients (same for Equal_index and greater_*) + struct Entry_hash { + Entry_hash(Rips_filtration const& filt) : filtp(&filt) {} + Rips_filtration const* filtp; + std::size_t operator()(const entry_t& e) const { return boost::hash()(filtp->get_index(e)); } + }; + + struct Equal_index { + Equal_index(Rips_filtration const& filt) : filtp(&filt) {} + Rips_filtration const* filtp; + bool operator()(const entry_t& e, const entry_t& f) const { + return filtp->get_index(e) == filtp->get_index(f); + } + }; + + struct diameter_simplex_t { + value_t diameter; + simplex_t index; + friend value_t get_diameter(const diameter_simplex_t& i) { return i.diameter; } + }; + static simplex_t get_index(const diameter_simplex_t& i) { return i.index; } + + struct diameter_entry_t : std::pair { + using std::pair::pair; + friend const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } + friend entry_t& get_entry(diameter_entry_t& p) { return p.second; } + friend value_t get_diameter(const diameter_entry_t& p) { return p.first; } + }; + simplex_t get_index(const diameter_entry_t& p) const { return get_index(get_entry(p)); } + coefficient_t get_coefficient(const diameter_entry_t& p) const { + return get_coefficient(get_entry(p)); + } + void set_coefficient(diameter_entry_t& p, const coefficient_t c) const { + set_coefficient(get_entry(p), c); + } + diameter_entry_t make_diameter_entry(value_t diameter, simplex_t index, coefficient_t coefficient) const { + return diameter_entry_t(diameter, make_entry(index, coefficient)); + } + diameter_entry_t make_diameter_entry(const diameter_simplex_t& diameter_index, coefficient_t coefficient) const { + return diameter_entry_t(get_diameter(diameter_index), make_entry(get_index(diameter_index), coefficient)); + } + + template struct Greater_diameter_or_smaller_index { + Greater_diameter_or_smaller_index(Rips_filtration const& filt) : filtp(&filt) {} + Rips_filtration const* filtp; + bool operator()(const Entry& a, const Entry& b) const { + return (get_diameter(a) > get_diameter(b)) || + ((get_diameter(a) == get_diameter(b)) && (filtp->get_index(a) < filtp->get_index(b))); + } + }; + + const DistanceMatrix dist; // only store a reference instead? + const vertex_t n; // redundant with dist? + const dimension_t dim_max; + const value_t threshold; // It would be nice if this was only in DistanceMatrix, but inconvenient. Only used to list the edges of dense distance matrices. + const coefficient_t modulus; + const SimplexEncoding simplex_encoding; // only store a reference instead? + mutable std::vector vertices; // we must not have several threads looking at the same complex + int bits_for_coeff; + + Rips_filtration(DistanceMatrix&& _dist, dimension_t _dim_max, value_t _threshold, coefficient_t _modulus) + : dist(std::move(_dist)), n(dist.size()), + dim_max(std::min(_dim_max, dist.size() - 2)), threshold(_threshold), + modulus(_modulus), simplex_encoding(n, dim_max + 2), bits_for_coeff(log2up(modulus-1)) { + // See entry_with_coeff_t for the logic for log2up(modulus-1) (storing coeff-1) + if (use_coefficients && simplex_encoding.num_extra_bits() < num_bits_for_coeff()) + // TODO: include relevant numbers in the message + throw std::overflow_error("Not enough spare bits in the simplex encoding to store a coefficient"); + } + + vertex_t num_vertices() const { return n; } + int num_bits_for_coeff() const { return bits_for_coeff; } + + simplex_t get_edge_index(const vertex_t i, const vertex_t j) const { + return simplex_encoding(i, 2) + j; + } + + template + OutputIterator get_simplex_vertices(simplex_t idx, const dimension_t dim, vertex_t n, + OutputIterator out) const { + --n; + for (dimension_t k = dim + 1; k > 1; --k) { + n = simplex_encoding.get_max_vertex(idx, k, n); + *out++ = n; + idx -= simplex_encoding(n, k); + } + *out = static_cast(idx); + return out; + } + + value_t compute_diameter(const simplex_t index, const dimension_t dim) const { + value_t diam = -std::numeric_limits::infinity(); + + vertices.resize(dim + 1); + get_simplex_vertices(index, dim, dist.size(), vertices.rbegin()); + + for (dimension_t i = 0; i <= dim; ++i) + for (dimension_t j = 0; j < i; ++j) { + diam = std::max(diam, dist(vertices[i], vertices[j])); + } + return diam; + } + + std::vector get_edges() { + if constexpr (!std::is_same_v) { // Compressed_lower_distance_matrix + // TODO: it would be convenient to have DistanceMatrix provide a range of neighbors at dist<=threshold even in the dense case (as a filtered_range) + std::vector edges; + for (vertex_t i = 0; i < n; ++i) { + for (vertex_t j = 0; j < i; ++j) { + value_t length = dist(i, j); + if (length <= threshold) edges.push_back({length, get_edge_index(i, j)}); + } + } + return edges; + } else { // Sparse_distance_matrix + std::vector edges; + for (vertex_t i = 0; i < n; ++i) + for (auto n : dist.neighbors[i]) { + vertex_t j = get_vertex(n); + if (i > j) edges.push_back({get_diameter(n), get_edge_index(i, j)}); + } + return edges; + } + } + + // TODO: document in what way (if any) the order matters + template class Simplex_coboundary_enumerator_ { // Compressed_lower_distance_matrix + simplex_t idx_below, idx_above; + vertex_t j; + dimension_t k; + std::vector vertices; + diameter_entry_t simplex; + const DistanceMatrix2& dist; + const SimplexEncoding& simplex_encoding; + const Rips_filtration& parent; // for n and get_simplex_vertices + // at least dist and simplex_encoding are redundant with parent, but using parent.dist and parent.simplex_encoding seems to have a bad impact on performance. + + public: + Simplex_coboundary_enumerator_(const Rips_filtration& _parent) : dist(_parent.dist), + simplex_encoding(_parent.simplex_encoding), parent(_parent) {} + + void set_simplex(const diameter_entry_t _simplex, const dimension_t _dim) { + idx_below = parent.get_index(_simplex); + idx_above = 0; + j = dist.size() - 1; + k = _dim + 1; + simplex = _simplex; + vertices.resize(_dim + 1); + parent.get_simplex_vertices(parent.get_index(_simplex), _dim, dist.size(), vertices.rbegin()); + } + + bool has_next(bool all_cofacets = true) { + return (j >= k && (all_cofacets || simplex_encoding(j, k) > idx_below)); + } + + std::optional next_raw(bool all_cofacets = true) { + if (!has_next(all_cofacets)) return std::nullopt; + // this requires simplex_encoding(x,0)>0 + while (simplex_encoding(j, k) <= idx_below) { + idx_below -= simplex_encoding(j, k); + idx_above += simplex_encoding(j, k + 1); + --j; + --k; + GUDHI_assert(k != -1); + } + value_t cofacet_diameter = get_diameter(simplex); + // The order of j and i matters for performance + for (vertex_t i : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(j, i)); + simplex_t cofacet_index = idx_above + simplex_encoding(j--, k + 1) + idx_below; + coefficient_t cofacet_coefficient = parent.get_coefficient(simplex); + if (k & 1) cofacet_coefficient = parent.modulus - cofacet_coefficient; + return parent.make_diameter_entry(cofacet_diameter, cofacet_index, cofacet_coefficient); + } + std::optional next(bool all_cofacets = true) { + while(true) { + std::optional res = next_raw(all_cofacets); + if(!res || get_diameter(*res) <= parent.threshold) return res; + } + } + }; + + template class Simplex_coboundary_enumerator_ { + typedef typename DistanceMatrix2::vertex_diameter_t vertex_diameter_t; + simplex_t idx_below, idx_above; + dimension_t k; + std::vector vertices; + diameter_entry_t simplex; + const DistanceMatrix2& dist; + const SimplexEncoding& simplex_encoding; + std::vector::const_reverse_iterator> neighbor_it; + std::vector::const_reverse_iterator> neighbor_end; + vertex_diameter_t neighbor; + const Rips_filtration& parent; // for n and get_simplex_vertices + + public: + Simplex_coboundary_enumerator_(const Rips_filtration& _parent) + : dist(_parent.dist), + simplex_encoding(_parent.simplex_encoding), parent(_parent) {} + + void set_simplex(const diameter_entry_t _simplex, const dimension_t _dim) { + idx_below = parent.get_index(_simplex); + idx_above = 0; + k = _dim + 1; + simplex = _simplex; + vertices.resize(_dim + 1); + parent.get_simplex_vertices(idx_below, _dim, dist.size(), vertices.rbegin()); + + neighbor_it.resize(_dim + 1); + neighbor_end.resize(_dim + 1); + for (dimension_t i = 0; i <= _dim; ++i) { + auto v = vertices[i]; + neighbor_it[i] = dist.neighbors[v].rbegin(); + neighbor_end[i] = dist.neighbors[v].rend(); + } + } + + bool has_next(bool all_cofacets = true) { + for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { + neighbor = *it0; + for (size_t idx = 1; idx < neighbor_it.size(); ++idx) { + auto &it = neighbor_it[idx], end = neighbor_end[idx]; + while (get_vertex(*it) > get_vertex(neighbor)) + if (++it == end) return false; + if (get_vertex(*it) != get_vertex(neighbor)) + goto continue_outer; + else + neighbor = std::max(neighbor, *it); + } + while (k > 0 && vertices[k - 1] > get_vertex(neighbor)) { + if (!all_cofacets) return false; + idx_below -= simplex_encoding(vertices[k - 1], k); + idx_above += simplex_encoding(vertices[k - 1], k + 1); + --k; + } + return true; +continue_outer:; + } + return false; + } + + std::optional next_raw(bool all_cofacets = true) { + if (!has_next(all_cofacets)) return std::nullopt; + ++neighbor_it[0]; + value_t cofacet_diameter = std::max(get_diameter(simplex), get_diameter(neighbor)); + simplex_t cofacet_index = idx_above + simplex_encoding(get_vertex(neighbor), k + 1) + idx_below; + const coefficient_t modulus = parent.modulus; + coefficient_t cofacet_coefficient = + (k & 1 ? modulus - 1 : 1) * parent.get_coefficient(simplex) % modulus; + return parent.make_diameter_entry(cofacet_diameter, cofacet_index, cofacet_coefficient); + } + std::optional next(bool all_cofacets = true) { + return next_raw(all_cofacets); + } + }; + + typedef Simplex_coboundary_enumerator_ Simplex_coboundary_enumerator; + + class Simplex_boundary_enumerator { + private: + simplex_t idx_below, idx_above; + vertex_t j; + dimension_t k; + diameter_entry_t simplex; + dimension_t dim; + const SimplexEncoding& simplex_encoding; + const Rips_filtration& parent; // for n, get_max_vertex, compute_diameter + + public: + Simplex_boundary_enumerator(const Rips_filtration& _parent) + : simplex_encoding(_parent.simplex_encoding), parent(_parent) {} + + void set_simplex(const diameter_entry_t _simplex, const dimension_t _dim) { + idx_below = parent.get_index(_simplex); + idx_above = 0; + j = parent.n - 1; + k = _dim; + simplex = _simplex; + dim = _dim; + } + + bool has_next() { return (k >= 0); } + + std::optional next() { + if (!has_next()) return std::nullopt; + j = simplex_encoding.get_max_vertex(idx_below, k + 1, j); + + simplex_t face_index = idx_above - simplex_encoding(j, k + 1) + idx_below; + + // It would make sense to extract the vertices once in set_simplex + // and pass the proper subset to compute_diameter, but even in cases + // where this dominates it does not seem to help (probably because we + // stop at the first coface). + value_t face_diameter = parent.compute_diameter(face_index, dim - 1); + + const coefficient_t modulus = parent.modulus; + coefficient_t face_coefficient = + (k & 1 ? -1 + modulus : 1) * parent.get_coefficient(simplex) % modulus; + + idx_below -= simplex_encoding(j, k + 1); + idx_above += simplex_encoding(j, k); + + --k; + + return parent.make_diameter_entry(face_diameter, face_index, face_coefficient); + } + }; +}; + +#if BOOST_VERSION >= 108100 +template using Hash_map = boost::unordered_flat_map; +#else +template using Hash_map = boost::unordered_map; +#endif + +template class Persistent_cohomology { + using size_t = typename Filtration::size_t; + using coefficient_t = typename Filtration::coefficient_t; + using coefficient_storage_t = typename Filtration::coefficient_storage_t; + using dimension_t = typename Filtration::dimension_t; + using value_t = typename Filtration::value_t; + using vertex_t = typename Filtration::vertex_t; + using simplex_t = typename Filtration::simplex_t; + using diameter_simplex_t = typename Filtration::diameter_simplex_t; + using entry_t = typename Filtration::entry_t; + using diameter_entry_t = typename Filtration::diameter_entry_t; + using Entry_hash = typename Filtration::Entry_hash; + using Equal_index = typename Filtration::Equal_index; + using Simplex_boundary_enumerator = typename Filtration::Simplex_boundary_enumerator; + using Simplex_coboundary_enumerator = typename Filtration::Simplex_coboundary_enumerator; + templateusing Greater_diameter_or_smaller_index = typename Filtration::template Greater_diameter_or_smaller_index; + + typedef Compressed_sparse_matrix_ Compressed_sparse_matrix; + typedef Hash_map entry_hash_map; + + Filtration filt; + const vertex_t n; + const dimension_t dim_max; + const coefficient_t modulus; + const std::vector multiplicative_inverse_; + mutable std::vector cofacet_entries; + mutable std::vector vertices; + Simplex_boundary_enumerator facets; + // Creating a new one in each function that needs it wastes a bit of time, but we may need 2 at the same time. + Simplex_coboundary_enumerator cofacets1, cofacets2; + + coefficient_t multiplicative_inverse(coefficient_t c) const { + return multiplicative_inverse_[c]; + } + public: + Persistent_cohomology(Filtration&& _filt, dimension_t _dim_max, coefficient_t _modulus) + : filt(std::move(_filt)), n(filt.num_vertices()), + dim_max(std::min(_dim_max, n - 2)), + modulus(_modulus), + multiplicative_inverse_(multiplicative_inverse_vector(_modulus)), + facets(filt), cofacets1(filt), cofacets2(filt) + { + } + + + std::optional get_zero_pivot_facet(const diameter_entry_t simplex, const dimension_t dim) { + facets.set_simplex(simplex, dim); + while(true) { + std::optional facet = facets.next(); + if (!facet) break; + if (get_diameter(*facet) == get_diameter(simplex)) return *facet; + } + return std::nullopt; + } + + std::optional get_zero_pivot_cofacet(const diameter_entry_t simplex, const dimension_t dim) { + cofacets1.set_simplex(simplex, dim); + while(true) { + std::optional cofacet = cofacets1.next_raw(); + if (!cofacet) break; + if (get_diameter(*cofacet) == get_diameter(simplex)) return *cofacet; + } + return std::nullopt; + } + + // Apparent pairs are implicit in Ripser. + // pro: we don't need to store them + // con: we may have to recompute them many times, and each test is more expensive than emergent pairs + std::optional get_zero_apparent_facet(const diameter_entry_t simplex, const dimension_t dim) { + std::optional facet = get_zero_pivot_facet(simplex, dim); + if (!facet) return std::nullopt; + std::optional cofacet = get_zero_pivot_cofacet(*facet, dim - 1); + if (!cofacet || filt.get_index(*cofacet) != filt.get_index(simplex)) return std::nullopt; + return *facet; + } + + std::optional get_zero_apparent_cofacet(const diameter_entry_t simplex, const dimension_t dim) { + std::optional cofacet = get_zero_pivot_cofacet(simplex, dim); + if (!cofacet) return std::nullopt; + std::optional facet = get_zero_pivot_facet(*cofacet, dim + 1); + if (!facet || filt.get_index(*facet) != filt.get_index(simplex)) return std::nullopt; + return *cofacet; + } + + bool is_in_zero_apparent_pair(const diameter_entry_t simplex, const dimension_t dim) { + return get_zero_apparent_cofacet(simplex, dim) || get_zero_apparent_facet(simplex, dim); + } + + void assemble_columns_to_reduce(std::vector& simplices, + std::vector& columns_to_reduce, + entry_hash_map& pivot_column_index, dimension_t dim) { + +#ifdef GUDHI_INDICATE_PROGRESS + std::cerr << clear_line << "assembling columns" << std::flush; + std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; +#endif + + columns_to_reduce.clear(); + std::vector next_simplices; + + for (diameter_simplex_t& simplex : simplices) { + cofacets2.set_simplex(filt.make_diameter_entry(simplex, 1), dim - 1); + + while(true) { + std::optional cofacet = cofacets2.next(false); + if (!cofacet) break; +#ifdef GUDHI_INDICATE_PROGRESS + if (std::chrono::steady_clock::now() > next) { + std::cerr << clear_line << "assembling " << next_simplices.size() + << " columns (processing " << std::distance(&simplices[0], &simplex) + << "/" << simplices.size() << " simplices)" << std::flush; + next = std::chrono::steady_clock::now() + time_step; + } +#endif + if (dim < dim_max) next_simplices.push_back({get_diameter(*cofacet), filt.get_index(*cofacet)}); + // Wouldn't it be cheaper in the reverse order? Seems negligible + if (!is_in_zero_apparent_pair(*cofacet, dim) && + (pivot_column_index.find(get_entry(*cofacet)) == pivot_column_index.end())) + columns_to_reduce.push_back({get_diameter(*cofacet), filt.get_index(*cofacet)}); + } + } + + if (dim < dim_max) simplices.swap(next_simplices); + +#ifdef GUDHI_INDICATE_PROGRESS + std::cerr << clear_line << "sorting " << columns_to_reduce.size() << " columns" + << std::flush; +#endif + + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + Greater_diameter_or_smaller_index(filt)); +#ifdef GUDHI_INDICATE_PROGRESS + std::cerr << clear_line << std::flush; +#endif + } + + template + void compute_dim_0_pairs(std::vector& edges, + std::vector& columns_to_reduce, OutPair& output_pair) { + Union_find dset(n); + + edges = filt.get_edges(); + std::sort(edges.rbegin(), edges.rend(), + Greater_diameter_or_smaller_index(filt)); + std::vector vertices_of_edge(2); + for (auto e : edges) { + // Should we work with pairs of vertices instead of edges, to skip get_simplex_vertices? + filt.get_simplex_vertices(filt.get_index(e), 1, n, vertices_of_edge.rbegin()); + vertex_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); + + if (u != v) { + if (get_diameter(e) != 0) + output_pair(0, get_diameter(e)); + dset.link(u, v); + } else if ((dim_max > 0) && !get_zero_apparent_cofacet(filt.make_diameter_entry(e, 1), 1)) + columns_to_reduce.push_back(e); + } + if (dim_max > 0) std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); + + for (vertex_t i = 0; i < n; ++i) + if (dset.find(i) == i) output_pair(0, std::numeric_limits::infinity()); + } + + template std::optional pop_pivot(Column& column) { + while(!column.empty()) { // At this stage the partial sum is 0 + diameter_entry_t pivot = column.top(); + column.pop(); + while(true) { // At this stage the partial sum is led by pivot + if (column.empty() || filt.get_index(column.top()) != filt.get_index(pivot)) return pivot; + coefficient_t sum = (filt.get_coefficient(pivot) + filt.get_coefficient(column.top())) % modulus; + column.pop(); + if (sum == 0) { + break; + } + filt.set_coefficient(pivot, sum); + } + } + return std::nullopt; + } + + template std::optional get_pivot(Column& column) { + // We could look for the pivot without needing to push it back, but it does not seem to help in benchmarks. + std::optional result = pop_pivot(column); + if (result) column.push(*result); + return result; + } + + // Either return the pivot in an emergent pair, or fill working_coboundary with the full coboundary (and return its pivot) + template + std::optional init_coboundary_and_get_pivot(const diameter_entry_t simplex, + Column& working_coboundary, const dimension_t dim, + entry_hash_map& pivot_column_index) { + bool check_for_emergent_pair = true; + cofacet_entries.clear(); + cofacets2.set_simplex(simplex, dim); + while(true) { + std::optional cofacet = cofacets2.next(); + if (!cofacet) break; + cofacet_entries.push_back(*cofacet); + if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(*cofacet))) { + if ((pivot_column_index.find(get_entry(*cofacet)) == pivot_column_index.end()) && + !get_zero_apparent_facet(*cofacet, dim + 1)) + return *cofacet; + check_for_emergent_pair = false; + } + } + for (auto cofacet : cofacet_entries) working_coboundary.push(cofacet); + return get_pivot(working_coboundary); + } + + // Beyond apparent/emergent pairs, Ripser does the column reduction eagerly. + // To keep with the lazy paradigm, it would be possible to represent a column as a heap of coboundary_iterator, using the order of the current simplices pointed to by the iterators. The equivalent of `pop` would be ++ on the top iterator, which decreases its priority (or removes it if it reached the end). If we do it right, it could also help us notice if we have twice the same iterator and avoid computing its full coboundary twice (although we may still end up computing the first simplex of the coboundary for both, since it may be the easiest way to detect duplicates without maintaining yet another structure on the side). Another advantage is that its size would be bounded by the number of simplices of dimension d instead of d+1 currently. (Dory seems to do something related but too complicated for me.) + template + void add_simplex_coboundary(const diameter_entry_t simplex, const dimension_t dim, + Column& working_reduction_column, Column& working_coboundary) { + working_reduction_column.push(simplex); + cofacets1.set_simplex(simplex, dim); + while(true) { // stupid C++ + std::optional cofacet = cofacets1.next(); + if (!cofacet) break; + working_coboundary.push(*cofacet); + } + } + + // add an already reduced column, i.e. add all the simplex coboundaries that were involved in that reduction + template + void add_coboundary(Compressed_sparse_matrix& reduction_matrix, + const std::vector& columns_to_reduce, + const size_t index_column_to_add, const coefficient_t factor, + const dimension_t dim, Column& working_reduction_column, + Column& working_coboundary) { + diameter_entry_t column_to_add = filt.make_diameter_entry(columns_to_reduce[index_column_to_add], factor); + add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary); + + for (diameter_entry_t simplex : reduction_matrix.subrange(index_column_to_add)) { + filt.set_coefficient(simplex, filt.get_coefficient(simplex) * factor % modulus); + add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary); + } + } + + template + void compute_pairs(const std::vector& columns_to_reduce, + entry_hash_map& pivot_column_index, const dimension_t dim, OutPair& output_pair) { + Compressed_sparse_matrix reduction_matrix; + Greater_diameter_or_smaller_index cmp(filt); + Heap, + Greater_diameter_or_smaller_index> + working_reduction_column(cmp), working_coboundary(cmp); + +#ifdef GUDHI_INDICATE_PROGRESS + std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; +#endif + for (size_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); + ++index_column_to_reduce) { + + diameter_entry_t column_to_reduce = filt.make_diameter_entry(columns_to_reduce[index_column_to_reduce], 1); + value_t diameter = get_diameter(column_to_reduce); + + reduction_matrix.append_column(); + + working_reduction_column.clear(); working_coboundary.clear(); + + std::optional pivot = init_coboundary_and_get_pivot( + column_to_reduce, working_coboundary, dim, pivot_column_index); + // When we found an emergent pair, we could avoid checking again below, but it does not seem to gain anything in practice. + + while (true) { +#ifdef GUDHI_INDICATE_PROGRESS + if (std::chrono::steady_clock::now() > next) { + std::cerr << clear_line << "reducing column " << index_column_to_reduce + 1 + << "/" << columns_to_reduce.size() << " (diameter " << diameter << ")" + << std::flush; + next = std::chrono::steady_clock::now() + time_step; + } +#endif + if (pivot) { + auto pair = pivot_column_index.find(get_entry(*pivot)); + if (pair != pivot_column_index.end()) { + entry_t other_pivot = pair->first; + size_t index_column_to_add = pair->second; + coefficient_t factor = + modulus - filt.get_coefficient(*pivot) * + multiplicative_inverse(filt.get_coefficient(other_pivot)) % + modulus; + + // It saves a little bit (3% on an example, 0% on another) if we pass pivot to add_coboundary and avoid pushing entries smaller than pivot in working_coboundary + add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, + factor, dim, working_reduction_column, working_coboundary); + + pivot = get_pivot(working_coboundary); + } else if (std::optional e = get_zero_apparent_facet(*pivot, dim + 1); e) { + filt.set_coefficient(*e, modulus - filt.get_coefficient(*e)); + + add_simplex_coboundary(*e, dim, working_reduction_column, working_coboundary); + + pivot = get_pivot(working_coboundary); + } else { + value_t death = get_diameter(*pivot); + output_pair(diameter, death); + pivot_column_index.insert({get_entry(*pivot), index_column_to_reduce}); + // CubicalRipser suggests caching the column here, at least if it took many operations to reduce it. + + while (true) { + std::optional e = pop_pivot(working_reduction_column); + if (!e) break; + GUDHI_assert(filt.get_coefficient(*e) > 0); + reduction_matrix.push_back(*e); + } + break; + } + } else { + output_pair(diameter, std::numeric_limits::infinity()); + break; + } + } + } +#ifdef GUDHI_INDICATE_PROGRESS + std::cerr << clear_line << std::flush; +#endif + } + + // Add a separate output_essential? + // TODO: Should output_pair also take a simplex_t argument? + template + void compute_barcodes(OutDim&& output_dim, OutPair&& output_pair) { + std::vector simplices, columns_to_reduce; + + output_dim(0); + compute_dim_0_pairs(simplices, columns_to_reduce, output_pair); + + for (dimension_t dim = 1; dim <= dim_max; ++dim) { + entry_hash_map pivot_column_index(0, filt, filt); + pivot_column_index.reserve(columns_to_reduce.size()); + + output_dim(dim); + compute_pairs(columns_to_reduce, pivot_column_index, dim, output_pair); + + if (dim < dim_max) + assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, + dim + 1); + // If for some odd reason one wanted all the infinite intervals in the last dimension, assemble_columns_to_reduce should give us that. + } + } +}; + +// An example of what Params can be +#if 0 +struct Params1 { + // size_t (not always from Params...) is used for counting (ok) and storage for the index of columns in a Hash_map. + // To gain on a pair by reducing size_t, simplex_t has to be smaller than size_t I guess, which is very small, not worth it. + typedef std::size_t size_t; + typedef float value_t; + typedef int8_t dimension_t; // Does it need to be signed? Experimentally no. + typedef int vertex_t; // Currently needs to be signed for Simplex_coboundary_enumerator_::has_next. Reducing to int16_t helps perf a bit. + typedef Gudhi::numbers::uint128_t simplex_t; + // We could introduce a smaller edge_t, but it is not trivial to separate and probably not worth it + typedef uint16_t coefficient_storage_t; // For the table of multiplicative inverses + typedef uint_least32_t coefficient_t; // We need x * y % z to work, but promotion from uint16_t to int is not enough + static const bool use_coefficients = false; + + // Assumptions used in the code: + // * dimension_t smaller than vertex_t + // Check which ones need to be signed, and which ones could be unsigned instead +}; +#endif + +// Trying to write a magic function +template +void help2(DistanceMatrix&& dist, int dim_max, typename DistanceMatrix::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + typedef Rips_filtration Filt; + typedef Persistent_cohomology Pcoh; + Filt filt(std::move(dist), dim_max, threshold, modulus); + Pcoh(std::move(filt), dim_max, modulus).compute_barcodes(output_dim, output_pair); +} +template struct TParams { + // hardcode most options + typedef std::size_t size_t; + typedef value_t_ value_t; + typedef int8_t dimension_t; + typedef int vertex_t; + typedef simplex_t_ simplex_t; + typedef uint16_t coefficient_storage_t; + typedef uint_least32_t coefficient_t; + static const bool use_coefficients = use_coefficients_; +}; +template struct TParams2 { + typedef int vertex_t; + typedef value_t_ value_t; +}; +// Choose simplex encoding +template +void help1(DistanceMatrix&& dist, int dim_max, typename DistanceMatrix::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + auto n = dist.size(); + // TODO: if n>=3 we could even go to n-3, because the Rips cannot have homology in dimension n-2 + // (e.g. for 3 points, there is no 1-homology) + if (dim_max > n - 2) dim_max = n - 2; // FIXME: duplicated. problem if n unsigned and 0 or 1. + int bits_per_vertex = log2up(n); + int bits_for_coeff = log2up(modulus - 1); // duplicating logic :-( Also, if modulus is something absurd, the diagnostic is too late + int bitfield_size = bits_per_vertex * (dim_max + 2) + bits_for_coeff; + if (bitfield_size <= 64) { // use bitfield-64 + typedef TParams P; + help2>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else if (bitfield_size <= 128) { // use bitfield-128 + typedef TParams P; + help2>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else { // use cns-128 + typedef TParams P; + help2>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } + // Does cns-64 have its place? +} +// Select hardcoded Z/2Z or runtime Z/pZ +template +void ripser(DistanceMatrix dist, int dim_max, typename DistanceMatrix::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + if (modulus == 2) + help1(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + else + help1(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); +} +#if 0 +template +void ripser_auto(Sparse_distance_matrix dist, int dim_max, typename DMParams::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); +} +template +void ripser_auto(Compressed_distance_matrix dist, int dim_max, typename DMParams::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + typedef typename DMParams::value_t value_t; + typedef typename DMParams::vertex_t vertex_t; + if (threshold < std::numeric_limits::max()) { // or infinity() + Sparse_distance_matrix new_dist(dist, threshold); + ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair); + } else { + for (vertex_t i = 0; i < dist.size(); ++i) { + value_t r_i = -std::numeric_limits::infinity(); + for (vertex_t j = 0; j < dist.size(); ++j) + r_i = std::max(r_i, dist(i, j)); + threshold = std::min(threshold, r_i); + // Should we also compute this when an explicit threshold is passed, in case the user passed an unnecessarily large threshold? + } + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } +} +// Other = Euclidean. TODO: more robust dispatching... (how?) +template +void ripser_auto(DistanceMatrix dist, int dim_max, typename DistanceMatrix::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + typedef typename DistanceMatrix::value_t value_t; + typedef TParams2 P; + if (threshold < std::numeric_limits::max()) { // or infinity() + Sparse_distance_matrix

new_dist(dist, threshold); + ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair); + } else { + Compressed_distance_matrix new_dist(dist); + ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair); + } +} +#endif +template +void ripser_auto(DistanceMatrix dist, int dim_max, typename DistanceMatrix::value_t threshold, unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) { + typedef typename DistanceMatrix::vertex_t vertex_t; + typedef typename DistanceMatrix::value_t value_t; + typedef TParams2 P; + if constexpr (std::is_same_v) { + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else if (threshold < std::numeric_limits::max()) { // or infinity() + Sparse_distance_matrix

new_dist(dist, threshold); + ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair); + } else if constexpr (std::is_same_v) { + for (vertex_t i = 0; i < dist.size(); ++i) { + value_t r_i = -std::numeric_limits::infinity(); + for (vertex_t j = 0; j < dist.size(); ++j) + r_i = std::max(r_i, dist(i, j)); + threshold = std::min(threshold, r_i); + // Should we also compute this when an explicit threshold is passed, in case the user passed an unnecessarily large threshold? + } + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else { + Compressed_distance_matrix new_dist(dist); + ripser_auto(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair); + } +} +// - sparse input -> sparse matrix +// - euclidean input & threshold -> sparse matrix (don't build dense matrix) +// - euclidean input & !threshold -> dense matrix +// - dense matrix & threshold -> sparse matrix +// - dense matrix & !threshold -> compute minmax, keep dense +} +#undef GUDHI_assert +#endif // GUDHI_RIPSER_H + +/* Relevant benchmarks where different functions dominate the profile +# push +ripser --format point-cloud ripser/examples/o3_1024.txt +# pop +ripser --format point-cloud ripser/examples/o3_4096.txt --threshold 1.6 +# apparent_pair (get_simplex_vertices + dist) - edge_collapse would help +ripser --format point-cloud --dim 2 --threshold .7 tore +# apparent_pair (get_simplex_vertices) - no edge collapse possible +ripser --format point-cloud circle24 --dim 25 + +where equivalent datasets can be obtained through + tore -> tore3D_1307.off + circle24: + t=np.linspace(0, 2*np.pi, 24, endpoint=False) + np.stack((np.cos(t),np.sin(t))).T + OR + gudhi.datasets.generators.points.torus(24,1,'grid') + o3 (?): + scipy.stats.ortho_group.rvs(3,4096).reshape(-1,9) +*/ diff --git a/src/Ripser/utilities/CMakeLists.txt b/src/Ripser/utilities/CMakeLists.txt new file mode 100644 index 0000000000..fac8513523 --- /dev/null +++ b/src/Ripser/utilities/CMakeLists.txt @@ -0,0 +1,3 @@ +add_executable(ripser ripser.cc) +install(TARGETS ripser DESTINATION bin) +add_test(NAME ripser_lower COMMAND $ --format lower --dim 2 --threshold .7 "${CMAKE_SOURCE_DIR}/data/distance_matrix/tore3D_1307_distance_matrix.csv") diff --git a/src/Ripser/utilities/ripser.cc b/src/Ripser/utilities/ripser.cc new file mode 100644 index 0000000000..6dabba826b --- /dev/null +++ b/src/Ripser/utilities/ripser.cc @@ -0,0 +1,391 @@ +/* Based on Ripser commit 140670f2c76997404601e43d8054151f46be9fd7 + * Modification(s): + * - YYYY/MM Author: Description of the modification + * - 2024 Marc Glisse: Heavy refactoring +*/ + +/* + + Ripser: a lean C++ code for computation of Vietoris-Rips persistence barcodes + + MIT License + + Copyright (c) 2015–2021 Ulrich Bauer + + 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. + + You are under no obligation whatsoever to provide any bug fixes, patches, or + upgrades to the features, functionality or performance of the source code + ("Enhancements") to anyone; however, if you choose to make your Enhancements + available either publicly, or directly to the author of this software, without + imposing a separate written license agreement for such Enhancements, then you + hereby grant the following license: a non-exclusive, royalty-free perpetual + license to install, use, modify, prepare derivative works, incorporate into + other computer software, distribute, and sublicense such enhancements or + derivative works thereof, in binary and source code form. + +*/ + +//#define GUDHI_INDICATE_PROGRESS + +#include +#include +#include +#include + +using Gudhi::ripser::ripser; + +struct Params { + typedef float value_t; + typedef int8_t dimension_t; + typedef int vertex_t; + typedef uint_least32_t coefficient_t; +}; + +typedef Params::value_t value_t; +typedef Params::dimension_t dimension_t; +typedef Params::vertex_t vertex_t; +typedef Params::coefficient_t coefficient_t; + +typedef Gudhi::ripser::Full_distance_matrix Full_distance_matrix; +typedef Gudhi::ripser::Compressed_distance_matrix Compressed_lower_distance_matrix; +typedef Gudhi::ripser::Compressed_distance_matrix Compressed_upper_distance_matrix; +typedef Gudhi::ripser::Sparse_distance_matrix Sparse_distance_matrix; +typedef Gudhi::ripser::Euclidean_distance_matrix Euclidean_distance_matrix; + +enum file_format { + LOWER_DISTANCE_MATRIX, + UPPER_DISTANCE_MATRIX, + DISTANCE_MATRIX, + POINT_CLOUD, + DIPHA, + SPARSE, + BINARY +}; + +static constexpr uint16_t endian_check=0xff00; +static bool is_big_endian() { return *reinterpret_cast(&endian_check); } + +template T read(std::istream& input_stream) { + T result; + char* p = reinterpret_cast(&result); + if (input_stream.read(p, sizeof(T)).gcount() != sizeof(T)) return T(); + if (is_big_endian()) std::reverse(p, p + sizeof(T)); + return result; +} + +Euclidean_distance_matrix read_point_cloud(std::istream& input_stream) { + std::vector> points; + + std::string line; + value_t value; + while (std::getline(input_stream, line)) { + std::vector point; + std::istringstream s(line); + while (s >> value) { + point.push_back(value); + s.ignore(); + } + if (!point.empty()) points.push_back(point); + assert(point.size() == points.front().size()); + } + + std::size_t d = points.front().size(); + Euclidean_distance_matrix eucl_dist(std::move(points)); + vertex_t n = eucl_dist.size(); + std::cout << "point cloud with " << n << " points in dimension " + << d << std::endl; + + return eucl_dist; +} + +Sparse_distance_matrix read_sparse_distance_matrix(std::istream& input_stream) { + typedef typename Sparse_distance_matrix::vertex_diameter_t vertex_diameter_t; + std::vector> neighbors; + std::size_t num_edges = 0; + + std::string line; + while (std::getline(input_stream, line)) { + std::istringstream s(line); + std::size_t i, j; + value_t value; + s >> i; + s.ignore(); + s >> j; + s.ignore(); + s >> value; + s.ignore(); + if (i != j) { + neighbors.resize(std::max({neighbors.size(), i + 1, j + 1})); + neighbors[i].emplace_back(j, value); + neighbors[j].emplace_back(i, value); + ++num_edges; + } + } + + for (std::size_t i = 0; i < neighbors.size(); ++i) + std::sort(neighbors[i].begin(), neighbors[i].end()); + + return Sparse_distance_matrix(std::move(neighbors), num_edges); +} + +Compressed_lower_distance_matrix read_lower_distance_matrix(std::istream& input_stream) { + std::vector distances; + value_t value; + while (input_stream >> value) { + distances.push_back(value); + input_stream.ignore(); + } + + return Compressed_lower_distance_matrix(std::move(distances)); +} + +Compressed_lower_distance_matrix read_upper_distance_matrix(std::istream& input_stream) { + std::vector distances; + value_t value; + while (input_stream >> value) { + distances.push_back(value); + input_stream.ignore(); + } + + return Compressed_lower_distance_matrix(Compressed_upper_distance_matrix(std::move(distances))); +} + +Compressed_lower_distance_matrix read_distance_matrix(std::istream& input_stream) { + std::vector distances; + + std::string line; + value_t value; + for (int i = 0; std::getline(input_stream, line); ++i) { + std::istringstream s(line); + for (int j = 0; j < i && s >> value; ++j) { + distances.push_back(value); + s.ignore(); + } + } + + return Compressed_lower_distance_matrix(std::move(distances)); +} + +Compressed_lower_distance_matrix read_dipha(std::istream& input_stream) { + if (read(input_stream) != 8067171840) { + std::cerr << "input is not a Dipha file (magic number: 8067171840)" << std::endl; + exit(-1); + } + + if (read(input_stream) != 7) { + std::cerr << "input is not a Dipha distance matrix (file type: 7)" << std::endl; + exit(-1); + } + + vertex_t n = read(input_stream); + + std::vector distances; + + for (vertex_t i = 0; i < n; ++i) + for (vertex_t j = 0; j < n; ++j) + if (i > j) + distances.push_back(read(input_stream)); + else + read(input_stream); + + return Compressed_lower_distance_matrix(std::move(distances)); +} + +Compressed_lower_distance_matrix read_binary(std::istream& input_stream) { + std::vector distances; + while (!input_stream.eof()) distances.push_back(read(input_stream)); + return Compressed_lower_distance_matrix(std::move(distances)); +} + +Compressed_lower_distance_matrix read_file(std::istream& input_stream, const file_format format) { + switch (format) { + case LOWER_DISTANCE_MATRIX: + return read_lower_distance_matrix(input_stream); + case UPPER_DISTANCE_MATRIX: + return read_upper_distance_matrix(input_stream); + case DISTANCE_MATRIX: + return read_distance_matrix(input_stream); + case POINT_CLOUD: + return read_point_cloud(input_stream); + case DIPHA: + return read_dipha(input_stream); + default: + return read_binary(input_stream); + } +} + +void print_usage_and_exit(int exit_code) { + std::cerr + << "Usage: " + << "ripser " + << "[options] [filename]" << std::endl + << std::endl + << "Options:" << std::endl + << std::endl + << " --help print this screen" << std::endl + << " --format use the specified file format for the input. Options are:" + << std::endl + << " lower-distance (lower triangular distance matrix)" + << std::endl + << " upper-distance (upper triangular distance matrix)" << std::endl + << " (default:) distance (distance matrix; only lower triangular part is read)" << std::endl + << " point-cloud (point cloud in Euclidean space)" << std::endl + << " dipha (distance matrix in DIPHA file format)" << std::endl + << " sparse (sparse distance matrix in sparse triplet format)" + << std::endl + << " binary (lower triangular distance matrix in binary format)" + << std::endl + << " --dim compute persistent homology up to dimension k" << std::endl + << " --threshold compute Rips complexes up to diameter t" << std::endl + << " --modulus

compute homology with coefficients in the prime field Z/pZ" + << std::endl + << " --ratio only show persistence pairs with death/birth ratio > r" << std::endl + << std::endl; + exit(exit_code); +} + +int main(int argc, char** argv) { + const char* filename = nullptr; + + file_format format = DISTANCE_MATRIX; + + dimension_t dim_max = 1; + value_t threshold = std::numeric_limits::max(); + float ratio = 1; + coefficient_t modulus = 2; + + for (int i = 1; i < argc; ++i) { + const std::string arg(argv[i]); + if (arg == "--help") { + print_usage_and_exit(0); + } else if (arg == "--dim") { + std::string parameter = std::string(argv[++i]); + std::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]); + std::size_t next_pos; + threshold = std::stof(parameter, &next_pos); + if (next_pos != parameter.size()) print_usage_and_exit(-1); + } else if (arg == "--ratio") { + std::string parameter = std::string(argv[++i]); + std::size_t next_pos; + ratio = std::stof(parameter, &next_pos); + if (next_pos != parameter.size()) print_usage_and_exit(-1); + } else if (arg == "--format") { + std::string parameter = std::string(argv[++i]); + if (parameter.rfind("lower", 0) == 0) + format = LOWER_DISTANCE_MATRIX; + else if (parameter.rfind("upper", 0) == 0) + format = UPPER_DISTANCE_MATRIX; + else if (parameter.rfind("dist", 0) == 0) + format = DISTANCE_MATRIX; + else if (parameter.rfind("point", 0) == 0) + format = POINT_CLOUD; + else if (parameter == "dipha") + format = DIPHA; + else if (parameter == "sparse") + format = SPARSE; + else if (parameter == "binary") + format = BINARY; + else + print_usage_and_exit(-1); + } else if (arg == "--modulus") { + std::string parameter = std::string(argv[++i]); + std::size_t next_pos; + modulus = std::stol(parameter, &next_pos); + if (next_pos != parameter.size() || !Gudhi::ripser::is_prime(modulus)) print_usage_and_exit(-1); + } else { + if (filename) { print_usage_and_exit(-1); } + filename = argv[i]; + } + } + + std::ifstream file_stream(filename); + if (filename && file_stream.fail()) { + std::cerr << "couldn't open file " << filename << std::endl; + exit(-1); + } + + auto output_dim = [](dimension_t dim) { + std::cout << "persistence intervals in dim " << (int)dim << ":" << std::endl; + }; + auto output_pair = [ratio](value_t birth, value_t death) { +#ifdef GUDHI_INDICATE_PROGRESS + // Not necessary if we redirect stdout + std::cerr << Gudhi::ripser::clear_line << std::flush; +#endif + if (death == std::numeric_limits::infinity()) + std::cout << " [" << birth << ", )" << std::endl; + else if (death > birth * ratio) + std::cout << " [" << birth << "," << death << ")" << std::endl; + }; + if (format == SPARSE) { + Sparse_distance_matrix dist = + read_sparse_distance_matrix(filename ? file_stream : std::cin); + std::cout << "sparse distance matrix with " << dist.size() << " points and " + << dist.num_edges << "/" << (dist.size() * (dist.size() - 1)) / 2 << " entries" + << std::endl; + + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else if (format == POINT_CLOUD && threshold < std::numeric_limits::max()) { + Sparse_distance_matrix dist(read_point_cloud(filename ? file_stream : std::cin), threshold); + ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair); + } else { + Compressed_lower_distance_matrix dist = + read_file(filename ? file_stream : std::cin, format); + + value_t min = std::numeric_limits::infinity(), + max = -std::numeric_limits::infinity(), max_finite = max; + std::size_t num_edges = 0; + + value_t enclosing_radius = std::numeric_limits::infinity(); + if (threshold >= std::numeric_limits::max()) { + for (vertex_t i = 0; i < dist.size(); ++i) { + value_t r_i = -std::numeric_limits::infinity(); + for (vertex_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j)); + enclosing_radius = std::min(enclosing_radius, r_i); + } + } + + for (auto d : dist.distances) { + min = std::min(min, d); + max = std::max(max, d); + if (d != std::numeric_limits::infinity()) max_finite = std::max(max_finite, d); + if (d <= threshold) ++num_edges; + } + std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; + + if (threshold >= std::numeric_limits::max()) { + std::cout << "distance matrix with " << dist.size() + << " points, using threshold at enclosing radius " << enclosing_radius + << std::endl; + ripser(std::move(dist), dim_max, enclosing_radius, modulus, output_dim, output_pair); + } else { + std::cout << "sparse distance matrix with " << dist.size() << " points and " + << num_edges << "/" << (dist.size() * (dist.size() - 1)) / 2 << " entries" + << std::endl; + + ripser(Sparse_distance_matrix(std::move(dist), threshold), dim_max, threshold, modulus, output_dim, output_pair); + } + } + return 0; +} diff --git a/src/common/include/gudhi/uint128.h b/src/common/include/gudhi/uint128.h new file mode 100644 index 0000000000..1907314b96 --- /dev/null +++ b/src/common/include/gudhi/uint128.h @@ -0,0 +1,157 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GUDHI_UINT128_H_ +#define GUDHI_UINT128_H_ +#include +#include +#include +#include + +/* What about boost::multiprecision::uint128_t? + * It works on linux, mac, windows, etc. + * On linux x86_64, it has low overhead since it delegates to unsigned __int128. + * It is a bit slower than Fake_uint128, probably because of checks in << and >>, but I only noticed a 10% + * overhead compared to unsigned __int128 on 1 testcase. + * On windows, it is implemented as { size_t; uint64_t[2]; }, that's 50% overhead on storage. On the same + * testcase, forcing linux-gcc to use the windows code led to 4x slowdown. Other testcases were not affected. + * The pathological testcase is computing homology in all dimensions for the vertices of a regular 24-gon. + */ + +// GUDHI_FORCE_FAKE_UINT128 is only used for tests +#if !defined __SIZEOF_INT128__ || defined GUDHI_FORCE_FAKE_UINT128 +namespace Gudhi::numbers { +class Fake_uint128 { + // Debug + #ifdef __SIZEOF_INT128__ + unsigned __int128 native() const { return ((unsigned __int128)high << 64) + low; } + #define GUDHI_VERIF(X) GUDHI_CHECK(res.native() == (X), "") + #else + #define GUDHI_VERIF(X) + #endif + public: + constexpr Fake_uint128(): high(0), low(0) {} + constexpr Fake_uint128(std::uint64_t a): high(0), low(a) {} + // Arithmetic + // (multiplication and division are not needed for now) + friend Fake_uint128 operator+(Fake_uint128 a, Fake_uint128 b){ + Fake_uint128 res; + res.low = a.low + b.low; + res.high = a.high + b.high + (res.low < a.low); + GUDHI_VERIF (a.native() + b.native()); + return res; + } + friend Fake_uint128 operator-(Fake_uint128 a, Fake_uint128 b){ + Fake_uint128 res; + res.low = a.low - b.low; + res.high = a.high - b.high - (res.low > a.low); + GUDHI_VERIF (a.native() - b.native()); + return res; + } + friend Fake_uint128 operator<<(Fake_uint128 a, uint8_t b){ + Fake_uint128 res; + GUDHI_CHECK(b < 128, ""); + if (b >= 64) { res.low = 0; res.high = a.low << (b-64); } + else if (b == 0) { res = a; } + else { res.low = a.low << b; res.high = a.high << b | a.low >> (64-b); } + GUDHI_VERIF (a.native() << b); + return res; + } + friend Fake_uint128 operator>>(Fake_uint128 a, uint8_t b){ + Fake_uint128 res; + GUDHI_CHECK(b < 128, ""); + if (b >= 64) { res.high = 0; res.low = a.high >> (b-64); } + else if (b == 0) { res = a; } + else { res.high = a.high >> b; res.low = a.low >> b | a.high << (64-b); } + GUDHI_VERIF (a.native() >> b); + return res; + } + friend Fake_uint128 operator&(Fake_uint128 a, Fake_uint128 b){ + Fake_uint128 res; + res.low = a.low & b.low; + res.high = a.high & b.high; + GUDHI_VERIF (a.native() & b.native()); + return res; + } + friend Fake_uint128 operator|(Fake_uint128 a, Fake_uint128 b){ + Fake_uint128 res; + res.low = a.low | b.low; + res.high = a.high | b.high; + GUDHI_VERIF (a.native() | b.native()); + return res; + } + friend Fake_uint128 operator~(Fake_uint128 a){ + Fake_uint128 res; + res.low = ~a.low; + res.high = ~a.high; + return res; + } + // In-place arithmetic + Fake_uint128& operator+=(Fake_uint128 a) { *this = *this + a; return *this; } + Fake_uint128& operator-=(Fake_uint128 a) { *this = *this - a; return *this; } + Fake_uint128& operator++() { if (++low == 0) ++high; return *this; } + Fake_uint128& operator--() { if (low-- == 0) --high; return *this; } + Fake_uint128& operator<<=(uint8_t a) { *this = *this << a; return *this; } + Fake_uint128& operator>>=(uint8_t a) { *this = *this >> a; return *this; } + Fake_uint128& operator&=(Fake_uint128 a) { *this = *this & a; return *this; } + Fake_uint128& operator|=(Fake_uint128 a) { *this = *this | a; return *this; } + // Comparisons + friend bool operator==(Fake_uint128 a, Fake_uint128 b){ + return a.low == b.low && a.high == b.high; + } + friend bool operator!=(Fake_uint128 a, Fake_uint128 b){ + return a.low != b.low || a.high != b.high; + } + friend bool operator<(Fake_uint128 a, Fake_uint128 b){ + return a.high < b.high || (a.high == b.high && a.low < b.low); + } + friend bool operator>(Fake_uint128 a, Fake_uint128 b){ + return a.high > b.high || (a.high == b.high && a.low > b.low); + } + friend bool operator<=(Fake_uint128 a, Fake_uint128 b){ + return a.high < b.high || (a.high == b.high && a.low <= b.low); + } + friend bool operator>=(Fake_uint128 a, Fake_uint128 b){ + return a.high > b.high || (a.high == b.high && a.low >= b.low); + } + // Misc + friend std::size_t hash_value(Fake_uint128 a) { + typedef std::pair P; + return boost::hash_value(P(a.high, a.low)); + } + template >> + explicit operator T() const { + GUDHI_CHECK(high == 0 && low <= std::numeric_limits::max(), ""); + return static_cast(low); + } + private: + std::uint64_t high, low; // does the order matter? + #undef GUDHI_VERIF +}; +typedef Fake_uint128 uint128_t; +} // namespace Gudhi::numbers +template<> class std::numeric_limits { + public: + static constexpr bool is_specialized = true; + static constexpr bool is_signed = false; + static constexpr bool is_integer = true; + static constexpr bool is_exact = true; + static constexpr bool has_infinity = false; + static constexpr bool is_modulo = true; + static constexpr int digits = 128; + static constexpr int radix = 2; + // etc +}; +#else +namespace Gudhi::numbers { +typedef unsigned __int128 uint128_t; +} // namespace Gudhi::numbers +#endif +#endif // GUDHI_UINT128_H_ diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index cce1ea77de..88aec7593f 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -245,6 +245,7 @@ if (TARGET CGAL::CGAL AND TARGET Eigen3::Eigen) endif () set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_pers_cub_low_dim', ") set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_edge_collapse', ") +set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'_ripser', ") # from windows vcpkg eigen 3.4.0#2 : build fails with # error C2440: '': cannot convert from 'Eigen::EigenBase::Index' to '__gmp_expr' diff --git a/src/python/gudhi/_ripser.cc b/src/python/gudhi/_ripser.cc new file mode 100644 index 0000000000..36e84ae9ec --- /dev/null +++ b/src/python/gudhi/_ripser.cc @@ -0,0 +1,200 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Marc Glisse + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include + +#include +#include +#include +#include + +#include + +using namespace Gudhi::ripser; + +namespace py = pybind11; +typedef std::vector< int> Vi; +typedef std::vector Vd; +typedef std::vector> V2f; +typedef std::vector> V2d; +PYBIND11_MAKE_OPAQUE(Vi); +PYBIND11_MAKE_OPAQUE(Vd); +PYBIND11_MAKE_OPAQUE(V2f); +PYBIND11_MAKE_OPAQUE(V2d); + +templatestruct Full { + typedef Tag_dense Category; + typedef int vertex_t; + typedef T value_t; + decltype(std::declval&>().template unchecked<2>()) data; + int size() const { return data.shape(0); } + T operator()(int i, int j) const { + return data(i, j); + } +}; + +templatestruct DParams { + typedef vertex_t_ vertex_t; + typedef value_t_ value_t; +}; + +template +py::list doit(DistanceMatrix&& dist, int max_dimension, typename DistanceMatrix::value_t max_edge_length, unsigned homology_coeff_field) { + //static_assert(!std::is_lvalue_reference_v); + typedef typename DistanceMatrix::value_t T; + // We could put everything in a single vector, and return slices of it, but I don't think there is much advantage. + std::vector>> dgms; + { + py::gil_scoped_release release; + auto output = [&](T birth, T death){ + // Skip empty intervals + if (birth < death) + dgms.back().push_back({birth, death}); + }; + auto switch_dim = [&](int new_dim){ + dgms.emplace_back(); + }; + ripser_auto(std::move(dist), max_dimension, max_edge_length, homology_coeff_field, switch_dim, output); + } + py::list ret; + for (auto&& dgm : dgms) + ret.append(py::array(py::cast(std::move(dgm)))); + return ret; +} + +template +py::list full(py::array_t matrix, int max_dimension, T max_edge_length, unsigned homology_coeff_field) { + Full dist{matrix.template unchecked<2>()}; + if(dist.data.ndim() != 2 || dist.data.shape(0) != dist.data.shape(1)) + throw std::runtime_error("Distance matrix must be a square 2-dimensional array"); + return doit(std::move(dist), max_dimension, max_edge_length, homology_coeff_field); +} + +py::list lower(py::object low_mat, int max_dimension, double max_edge_length, unsigned homology_coeff_field) { + using Dist = Compressed_distance_matrix, LOWER_TRIANGULAR>; + std::vector distances; + int rowi = 0; + for (auto&& row : low_mat) { + if (rowi == 0) { ++rowi; continue; } + int coli = 0; + for (auto&& elem : row) { + distances.push_back(elem.cast()); // need a cast? + if (++coli == rowi) break; + } + if (coli < rowi) throw std::invalid_argument("Not enough elements for a lower triangular matrix"); + ++rowi; + }; + + // optional as a trick to allow destruction where I want it + std::optional release_local(std::in_place); + Dist dist(std::move(distances)); + release_local.reset(); + + return doit(std::move(dist), max_dimension, max_edge_length, homology_coeff_field); +} + +template +py::list sparse(py::array_t is_, py::array_t js_, py::array_t fs_, int num_vertices, int max_dimension, T max_edge_length, unsigned homology_coeff_field) { + // Duplicate entries and self loops are forbidden + auto is = is_.unchecked(); + auto js = js_.unchecked(); + auto fs = fs_.unchecked(); + if (is.ndim() != 1 || js.ndim() != 1 || fs.ndim() != 1) + throw std::runtime_error("vertices and filtrations must be 1-dimensional arrays"); + if (is.shape(0) != js.shape(0) || is.shape(0) != js.shape(0)) + throw std::runtime_error("vertices and filtrations must have the same shape"); + + typedef DParams P; + typedef Sparse_distance_matrix

Dist; + typedef typename Dist::vertex_diameter_t vertex_diameter_t; + + std::optional release_local(std::in_place); + std::vector> neighbors(num_vertices); + for (py::ssize_t e = 0; e < is.shape(0); ++e) { + neighbors[is(e)].emplace_back(js(e), fs(e)); + neighbors[js(e)].emplace_back(is(e), fs(e)); + } + // We could easily parallelize this loop, but it is unlikely to be worth it. + for (size_t i = 0; i < neighbors.size(); ++i) + std::sort(neighbors[i].begin(), neighbors[i].end()); + Dist dist(std::move(neighbors)); + release_local.reset(); + + return doit(std::move(dist), max_dimension, max_edge_length, homology_coeff_field); +} + +py::list lower_to_coo(py::object low_mat, double max_edge_length) { + // Cannot release the GIL since we keep accessing Python objects. + // TODO: full_to_coo for numpy arrays? + // Should we compute the cone radius at the same time? + std::vector is, js; + std::vector fs; + int rowi = 0; + for (auto&& row : low_mat) { + if (rowi == 0) { ++rowi; continue; } + int coli = 0; + for (auto&& elem : row) { + double d = elem.cast(); // need a cast? + if (d <= max_edge_length) { + is.push_back(rowi); + js.push_back(coli); + fs.push_back(d); + } + if (++coli == rowi) break; + } + if (coli < rowi) throw std::invalid_argument("Not enough elements for a lower triangular matrix"); + ++rowi; + }; + return py::make_tuple( + py::array(py::cast(std::move(is))), + py::array(py::cast(std::move(js))), + py::array(py::cast(std::move(fs)))); +} + +double lower_cone_radius(py::object low_mat) { + // It would be more efficient to read the matrix only once + auto n = py::len(low_mat); + std::vector maxs(n, -std::numeric_limits::infinity()); + int rowi = 0; + for (auto&& row : low_mat) { + if (rowi == 0) { ++rowi; continue; } + int coli = 0; + for (auto&& elem : row) { + double d = elem.cast(); + maxs[rowi] = std::max(maxs[rowi], d); + maxs[coli] = std::max(maxs[coli], d); + if (++coli == rowi) break; + } + if (coli < rowi) throw std::invalid_argument("Not enough elements for a lower triangular matrix"); + ++rowi; + }; + return *std::min_element(maxs.begin(), maxs.end()); +} + +PYBIND11_MODULE(_ripser, m) { + py::bind_vector(m, "VectorInt" , py::buffer_protocol()); + py::bind_vector(m, "VectorDouble" , py::buffer_protocol()); + py::bind_vector(m, "VectorPairFloat" , py::buffer_protocol()); + py::bind_vector(m, "VectorPairDouble", py::buffer_protocol()); + // Remove the default for max_dimension? + m.def("_full", full, py::arg("matrix").noconvert(), py::arg("max_dimension") = std::numeric_limits::max(), py::arg("max_edge_length") = std::numeric_limits::infinity(), py::arg("homology_coeff_field") = 2); + m.def("_full", full, py::arg("matrix"), py::arg("max_dimension") = std::numeric_limits::max(), py::arg("max_edge_length") = std::numeric_limits::infinity(), py::arg("homology_coeff_field") = 2); + m.def("_lower", lower, py::arg("matrix"), py::arg("max_dimension") = std::numeric_limits::max(), py::arg("max_edge_length") = std::numeric_limits::infinity(), py::arg("homology_coeff_field") = 2); + // We could do a version with long, but copying the arrays of integers shouldn't be too costly + m.def("_sparse", sparse, py::arg("row"), py::arg("col"), py::arg("data").noconvert(), py::arg("num_vertices"), py::arg("max_dimension") = std::numeric_limits::max(), py::arg("max_edge_length") = std::numeric_limits::infinity(), py::arg("homology_coeff_field") = 2); + m.def("_sparse", sparse, py::arg("row"), py::arg("col"), py::arg("data"), py::arg("num_vertices"), py::arg("max_dimension") = std::numeric_limits::max(), py::arg("max_edge_length") = std::numeric_limits::infinity(), py::arg("homology_coeff_field") = 2); + // Not directly an interface to Ripser... + m.def("_lower_to_coo", lower_to_coo, py::arg("matrix"), py::arg("max_edge_length")); + m.def("_lower_cone_radius", lower_cone_radius, py::arg("matrix")); +} + +// We could also create a RipsComplex class, that allows looking at a simplex, querying its (co)boundary, etc. But I am not convinced it is worth the effort. diff --git a/src/python/gudhi/sklearn/rips_persistence.py b/src/python/gudhi/sklearn/rips_persistence.py index f34af5a464..68f0646e14 100644 --- a/src/python/gudhi/sklearn/rips_persistence.py +++ b/src/python/gudhi/sklearn/rips_persistence.py @@ -7,8 +7,16 @@ # Modification(s): # - YYYY/MM Author: Description of the modification -from .. import RipsComplex +from .._ripser import _lower, _full, _sparse, _lower_to_coo, _lower_cone_radius +from ..flag_filtration.edge_collapse import reduce_graph +from .. import SimplexTree +import math +import numpy as np +from typing import Union, Iterable, Literal, Optional from sklearn.base import BaseEstimator, TransformerMixin +from scipy.sparse import coo_matrix +from scipy.spatial import cKDTree +from scipy.spatial.distance import pdist, squareform # joblib is required by scikit-learn from joblib import Parallel, delayed @@ -35,12 +43,14 @@ class RipsPersistence(BaseEstimator, TransformerMixin): def __init__( self, - homology_dimensions, - threshold=float('inf'), - input_type='point cloud', - num_collapses=1, - homology_coeff_field=11, - n_jobs=None, + homology_dimensions: Union[int, Iterable[int]], + threshold: float = float('inf'), + input_type: Literal[ + "point cloud", "full distance matrix", "lower distance matrix", "distance coo_matrix" + ] = 'point cloud', + num_collapses: Union[int, Literal["auto"]] = 'auto', + homology_coeff_field: int = 11, + n_jobs: Optional[int] = None, ): """ Constructor for the RipsPersistence class. @@ -49,14 +59,16 @@ def __init__( homology_dimensions (int or list of int): The returned persistence diagrams dimension(s). Short circuit the use of :class:`~gudhi.representations.preprocessing.DimensionSelector` when only one dimension matters (in other words, when `homology_dimensions` is an int). - threshold (float): Rips maximal edge length value. Default is +Inf. - input_type (str): Can be 'point cloud' when inputs are point clouds, or 'lower distance matrix', when - inputs are lower triangular distance matrix (can be full square, but the upper part of the distance - matrix will not be considered). Default is 'point cloud'. - num_collapses (int): Specify the number of :func:`~gudhi.SimplexTree.collapse_edges` iterations to perform - on the SimplexTree. Default value is 1 (a relatively good enough number of iterations). + threshold (float): Rips maximal edge length value. Default is +Inf. Ignored if input_type is 'distance coo_matrix'. + input_type (str): Can be 'point cloud' when inputs are point clouds, 'full distance matrix', + 'lower distance matrix' when inputs are lower triangular distance matrix (can be full square, + but the upper part will not be considered), or 'distance coo_matrix' for a distance matrix in SciPy's + sparse format, which should contain each edge at most once (avoid the symmetric) and no diagonal entry. + Default is 'point cloud'. + num_collapses (int|str): Specify the number of iterations of :func:`~gudhi.flag_filtration.edge_collapse.reduce_graph` + (edge collapse) to perform on the graph. Default value is 'auto'. homology_coeff_field (int): The homology coefficient field. Must be a prime number. Default value is 11. - n_jobs (int): Number of jobs to run in parallel. `None` (default value) means `n_jobs = 1` unless in a + n_jobs (Optional[int]): Number of jobs to run in parallel. `None` (default value) means `n_jobs = 1` unless in a joblib.parallel_backend context. `-1` means using all processors. cf. https://joblib.readthedocs.io/en/latest/generated/joblib.Parallel.html for more details. """ @@ -73,36 +85,89 @@ def fit(self, X, Y=None): """ return self - def __transform(self, inputs): - max_dimension = max(self.dim_list_) + 1 - - if self.input_type == 'point cloud': - rips = RipsComplex(points=inputs, max_edge_length = self.threshold) - elif self.input_type == 'lower distance matrix': - rips = RipsComplex(distance_matrix=inputs, max_edge_length = self.threshold) + def __transform(self, inp): + # TODO: give the user more control over the strategy + # Should we use threshold in the sparse case? + num_collapses = self.num_collapses + input_type = self.input_type + threshold = self.threshold + n = inp.shape[0] if input_type == 'distance coo_matrix' else len(inp) + max_dimension = min(max(self.dim_list_), max(0, n-3)) + # Ripser needs to encode simplices and coefficients in 128 bits, which may not always fit + # Instead of a 256 bit version which may not always suffice either, fall back to SimplexTree + use_simplex_tree = math.comb(n, min(n // 2, max_dimension + 2)) >= (1 << (128 - (self.homology_coeff_field - 2).bit_length())) + if num_collapses == 'auto': + num_collapses = 1 if max_dimension > (not use_simplex_tree) else 0 + # or num_collapses=max_dimension-1 maybe? + elif max_dimension == 0: + num_collapses = 0 + + # Points -> distance matrix + if input_type == 'point cloud': + if threshold < float('inf'): + # Hope that the user gave a useful threshold + tree = cKDTree(inp) + + ## V1: Returns self-loops and every edge twice (symmetry) + # inp = tree.sparse_distance_matrix(tree, max_distance=threshold, output_type="coo_matrix") + # mask = inp.row < inp.col + # inp = coo_matrix((inp.data[mask], (inp.row[mask], inp.col[mask])), shape=inp.shape) + + # V2: Gets the right edges, but forgets the distances + pairs = tree.query_pairs(r=threshold, output_type='ndarray') + data = np.ravel(np.linalg.norm(np.diff(inp[pairs], axis=1), axis=-1)) + inp = coo_matrix((data, (pairs[:,0], pairs[:,1])), shape=(n,)*2) + + input_type = 'distance coo_matrix' + else: + inp = squareform(pdist(inp)) + input_type = 'full distance matrix' + + # Dense -> sparse + if input_type in ('full distance matrix', 'lower distance matrix'): + # After this filtration value, all complexes are cones, nothing happens + if input_type == 'full distance matrix': + inp = np.asarray(inp) + cone_radius = inp.max(-1).min() + else: + cone_radius = _lower_cone_radius(inp) + sparsify = use_simplex_tree or num_collapses > 0 or threshold < cone_radius # heuristic + threshold = min(threshold, cone_radius) + if sparsify: + # For 'full' we could use i, j = np.triu_indices_from(inp, k=1), etc + i, j, f = _lower_to_coo(inp, threshold) + inp = coo_matrix((f, (i, j)), shape=(n,) * 2) + input_type = 'distance coo_matrix' + + if num_collapses > 0: + assert input_type == 'distance coo_matrix' + inp = reduce_graph(inp, num_collapses) + + if use_simplex_tree: + st = SimplexTree() + # Use create_from_array in case of full matrix? + # (not important since this fallback mostly matters in high dimension, where we use edge-collapse anyway) + st.insert_batch(np.arange(n).reshape(1,-1), np.zeros(n)) + st.insert_edges_from_coo_matrix(inp) + st.expansion(max_dimension + 1) + st.compute_persistence(homology_coeff_field=self.homology_coeff_field, persistence_dim_max=max_dimension>=st.dimension()) + return [ st.persistence_intervals_in_dimension(dim) for dim in self.dim_list_ ] + + if input_type == 'full distance matrix': + ## Possibly transpose for performance? + # if inp.strides[0] > inp.strides[1]: # or the reverse? + # inp = inp.T + dgm = _full(inp, max_dimension=max_dimension, max_edge_length=threshold, homology_coeff_field=self.homology_coeff_field) + elif input_type == 'lower distance matrix': + dgm = _lower(inp, max_dimension=max_dimension, max_edge_length=threshold, homology_coeff_field=self.homology_coeff_field) + elif input_type == 'distance coo_matrix': + # Switch to coo_array (danger: row/col seem deprecated)? + dgm = _sparse(inp.row, inp.col, inp.data, inp.shape[0], max_dimension=max_dimension, max_edge_length=threshold, homology_coeff_field=self.homology_coeff_field) else: - raise ValueError("Only 'point cloud' and 'lower distance matrix' are valid input_type") + raise ValueError("Only 'point cloud', 'lower distance matrix', 'full distance matrix' and 'distance coo_matrix' are valid input_type") # move to __init__? - if max_dimension > 1: - stree = rips.create_simplex_tree(max_dimension=1) - stree.collapse_edges(nb_iterations = self.num_collapses) - stree.expansion(max_dimension) - else: - stree = rips.create_simplex_tree(max_dimension=max_dimension) - - persistence_dim_max = False - # Specific case where, despite expansion(max_dimension), stree has a lower dimension - if max_dimension > stree.dimension(): - persistence_dim_max = True - - stree.compute_persistence( - homology_coeff_field=self.homology_coeff_field, - persistence_dim_max=persistence_dim_max - ) - - return [ - stree.persistence_intervals_in_dimension(dim) for dim in self.dim_list_ - ] + # dgm stops at n-2 + return [dgm[dim] if dim < len(dgm) else np.empty((0,2)) for dim in self.dim_list_] def transform(self, X, Y=None): """Compute all the Vietoris-Rips complexes and their associated persistence diagrams. @@ -117,7 +182,7 @@ def transform(self, X, Y=None): `[[array( Hi(X[0]) ), array( Hj(X[0]) )], [array( Hi(X[1]) ), array( Hj(X[1]) )], ...]` :rtype: list of numpy ndarray of shape (,2) or list of list of numpy ndarray of shape (,2) """ - # Depends on homology_dimensions is an integer or a list of integer (else case) + # Depends if homology_dimensions is an integer or a list of integers (else case) if isinstance(self.homology_dimensions, int): unwrap = True self.dim_list_ = [ self.homology_dimensions ] diff --git a/src/python/test/test_sklearn_rips_persistence.py b/src/python/test/test_sklearn_rips_persistence.py index 2ca0c313e0..928c82c313 100644 --- a/src/python/test/test_sklearn_rips_persistence.py +++ b/src/python/test/test_sklearn_rips_persistence.py @@ -10,6 +10,12 @@ from gudhi.datasets.generators import points from gudhi.sklearn.rips_persistence import RipsPersistence +from gudhi import RipsComplex, SimplexTree +from gudhi._ripser import _lower, _full, _sparse, _lower_to_coo, _lower_cone_radius +from gudhi import bottleneck_distance +from scipy.sparse import coo_matrix +from scipy.spatial.distance import cdist +import numpy as np import random import pytest @@ -36,12 +42,12 @@ def test_h1_only_rips_persistence_of_points_on_a_circle(): rips = RipsPersistence(homology_dimensions=1, n_jobs=-2) diags = rips.fit_transform([points.sphere(n_samples=150, ambient_dim=2)])[0] assert len(diags) == 1 - assert 0. < diags[0][0] < 0.6 - assert 1. < diags[0][1] < 2. + assert 0.0 < diags[0][0] < 0.6 + assert 1.0 < diags[0][1] < 2.0 def test_invalid_input_type(): - rips = RipsPersistence(homology_dimensions=0, input_type='whatsoever') + rips = RipsPersistence(homology_dimensions=0, input_type="whatsoever") with pytest.raises(ValueError): rips.fit_transform([points.sphere(n_samples=10, ambient_dim=2)]) @@ -49,13 +55,14 @@ def test_invalid_input_type(): def test_distance_matrix_rips_persistence_of_points_on_a_circle(): try: from scipy.spatial.distance import cdist + pts = points.sphere(n_samples=150, ambient_dim=2) distance_matrix = cdist(pts, pts) - rips = RipsPersistence(homology_dimensions=1, input_type='lower distance matrix') + rips = RipsPersistence(homology_dimensions=1, input_type="lower distance matrix") diags = rips.fit_transform([distance_matrix])[0] assert len(diags) == 1 - assert 0. < diags[0][0] < 0.6 - assert 1. < diags[0][1] < 2. + assert 0.0 < diags[0][0] < 0.6 + assert 1.0 < diags[0][1] < 2.0 except ValueError: pass @@ -63,13 +70,95 @@ def test_distance_matrix_rips_persistence_of_points_on_a_circle(): def test_set_output(): try: import pandas + NB_PC = 5 point_clouds = [points.sphere(n_samples=random.randint(100, 150), ambient_dim=2) for _ in range(NB_PC)] rips = RipsPersistence(homology_dimensions=[0, 2], n_jobs=-2) diags_pandas = rips.set_output(transform="pandas").fit_transform(point_clouds) - assert 'H0' == diags_pandas.columns[0] - assert 'H2' == diags_pandas.columns[1] + assert "H0" == diags_pandas.columns[0] + assert "H2" == diags_pandas.columns[1] assert len(diags_pandas.index) == NB_PC except ImportError: print("Missing pandas, skipping set_output test") + + +def test_big(): + # A circle + many isolated points + n = 1000000 + X = np.zeros((n, 2)) + X[:, 0] = np.arange(n) * 100 + X[:24] = points.torus(24, 1, "grid") + # Ripser cannot handle it, have to fall back to SimplexTree + # Computing the full distance matrix would require too much memory -> kd-tree + RipsPersistence(range(25), threshold=10).fit_transform([X]) + + +def cmp_rips(point_cloud): + primes = [2, 3, 11, 17, 29] # Small list so 2 is often selected + field = random.choice(primes) + print(f"random prime = {field}") + + print(f"nb points = {len(point_cloud)}, dim = {point_cloud.shape[1]}") + dists = cdist(point_cloud, point_cloud) + + # Check cone radius + cr = dists.max(-1).min() + assert cr == _lower_cone_radius(dists) + assert cr < 2.0 + + ## Compute with the SimplexTree + stree = RipsComplex(distance_matrix=dists).create_simplex_tree(max_dimension=2) + stree.compute_persistence(homology_coeff_field=field, persistence_dim_max=True) + dgm0 = stree.persistence_intervals_in_dimension(0) + dgm1 = stree.persistence_intervals_in_dimension(1) + + # Compute with Ripser + dgm = _full(dists, max_dimension=1, max_edge_length=float("inf"), homology_coeff_field=field) + # The order of the intervals may differ, so we cannot compare the arrays with np.testing.assert_almost_equal + assert bottleneck_distance(dgm0, dgm[0]) < 1e-8 + assert bottleneck_distance(dgm1, dgm[1]) < 1e-8 + + dgm = _lower(dists, max_dimension=1, max_edge_length=float("inf"), homology_coeff_field=field) + assert bottleneck_distance(dgm0, dgm[0]) < 1e-8 + assert bottleneck_distance(dgm1, dgm[1]) < 1e-8 + + # Convert to coo matrix + n = len(point_cloud) + dists_copy = np.array(dists, copy=True) + dists_copy[np.triu_indices_from(dists_copy)] = 0 # Keep only the lower entries + dists_sparse = coo_matrix(dists_copy) + s1 = _lower_to_coo(dists, float("inf")) + s2 = _lower_to_coo(dists_copy, float("inf")) + s1 = coo_matrix((s1[2], (s1[0], s1[1])), shape=(n, n)) + s2 = coo_matrix((s2[2], (s2[0], s2[1])), shape=(n, n)) + assert np.array_equal(s1.toarray(), s2.toarray()) + assert np.array_equal(dists_sparse.toarray(), s1.toarray()) + ## Compute with the SimplexTree + stree = SimplexTree() + stree.insert_batch(np.arange(n).reshape(1, -1), np.zeros(n)) + stree.insert_edges_from_coo_matrix(dists_sparse) + stree.expansion(2) + stree.compute_persistence(homology_coeff_field=field, persistence_dim_max=True) + sp_dgm0 = stree.persistence_intervals_in_dimension(0) + sp_dgm1 = stree.persistence_intervals_in_dimension(1) + + # Compute with Ripser + sp_dgm = _sparse( + dists_sparse.row, + dists_sparse.col, + dists_sparse.data, + dists_sparse.shape[0], + max_dimension=1, + max_edge_length=float("inf"), + homology_coeff_field=field, + ) + assert bottleneck_distance(sp_dgm0, sp_dgm[0]) < 1e-8 + assert bottleneck_distance(sp_dgm1, sp_dgm[1]) < 1e-8 + assert bottleneck_distance(sp_dgm0, dgm0) < 1e-8 + assert bottleneck_distance(sp_dgm1, dgm1) < 1e-8 + + +def test_ripser_interfaces(): + cmp_rips(points.sphere(n_samples=random.randint(100, 150), ambient_dim=2)) + cmp_rips(np.random.rand(random.randint(100, 150), random.randint(2, 3)))