diff options
Diffstat (limited to 'src/Persistent_cohomology/include/gudhi/Persistent_cohomology')
3 files changed, 426 insertions, 0 deletions
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h new file mode 100644 index 00000000..af0d6605 --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -0,0 +1,106 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_FIELD_ZP_H +#define GUDHI_FIELD_ZP_H + +namespace Gudhi{ + +/** \brief Structure representing the coefficient field \f$\mathbb{Z}/p\mathbb{Z}\f$ + * + * \implements CoefficientField + * \ingroup persistent_cohomology + */ +class Field_Zp { +public: +typedef int Element; + +Field_Zp() +: Prime(-1) +, inverse_() {} + +void init(int charac ) { + assert(charac <= 32768); + Prime = charac; + inverse_.clear(); + inverse_.reserve(charac); + inverse_.push_back(0); + for(int i=1 ; i<Prime ; ++i) + { + int inv = 1; + while(((inv * i) % Prime) != 1) ++inv; + inverse_.push_back(inv); + } +} + +/** Set x <- x + w * y*/ +void plus_times_equal ( Element & x, Element y, Element w ) +{ x = (x + w * y) % Prime; } + +// operator= defined on Element + +/** Returns y * w */ +Element times ( Element y, int w ) { + Element res = (y * w) % Prime; + if(res < 0) return res+Prime; + else return res; +} + +void clear_coefficient(Element x) {} + +void plus_equal(Element & x, Element y) { x = ((x+y)%Prime); } + +/** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ +Element additive_identity () { return 0; } +/** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ +Element multiplicative_identity ( Element P = 0) { return 1; } +/** Returns the inverse in the field. Modifies P.*/ +std::pair<Element,Element> inverse ( Element x + , Element P ) +{ return std::pair<Element,Element>(inverse_[x],P); +} // <------ return the product of field characteristic for which x is invertible + +/** Returns -x * y.*/ +Element times_minus ( Element x, Element y ) +{ + Element out = (-x * y) % Prime; + return (out < 0) ? out + Prime : out; +} + + +bool is_one ( Element x ) { return x == 1; } +bool is_zero ( Element x ) { return x == 0; } + +//bool is_null() + +/** \brief Returns the characteristic \f$p\f$ of the field.*/ +Element characteristic() { return Prime; } + +private: + Element Prime; +/** Property map Element -> Element, which associate to an element its inverse in the field.*/ + std::vector< Element > inverse_; +}; + +} // namespace GUDHI + +#endif // GUDHI_FIELD_ZP_H diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h new file mode 100644 index 00000000..9dd0998c --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -0,0 +1,167 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_MULTI_FIELD_H +#define GUDHI_MULTI_FIELD_H + +#include <iostream> +#include <vector> +#include <gmpxx.h> + +namespace Gudhi{ + +/** \brief Structure representing coefficients in a set of finite fields simultaneously + * using the chinese remainder theorem. + * + * \implements CoefficientField + * \ingroup persistent_cohomology + + * Details on the algorithms may be found in \cite boissonnat:hal-00922572 + */ +class Multi_field { +public: + typedef mpz_class Element; + + Multi_field () {} + +/* Initialize the multi-field. The generation of prime numbers might fail with + * a very small probability.*/ + void init(int min_prime, int max_prime) + { + if(max_prime<2) + { std::cerr << "There is no prime less than " << max_prime << std::endl; } + if(min_prime > max_prime) + { std::cerr << "No prime in ["<<min_prime<<":"<<max_prime<<"]"<<std::endl; } + // fill the list of prime numbers + unsigned int curr_prime = min_prime; + mpz_t tmp_prime; mpz_init_set_ui(tmp_prime,min_prime); + //test if min_prime is prime + int is_prime = mpz_probab_prime_p(tmp_prime,25); //probabilistic primality test + + if(is_prime == 0) //min_prime is composite + { + mpz_nextprime(tmp_prime,tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + + while (curr_prime <= max_prime) + { + primes_.push_back(curr_prime); + mpz_nextprime(tmp_prime,tmp_prime); + curr_prime = mpz_get_ui(tmp_prime); + } + //set m to primorial(bound_prime) + prod_characteristics_ = 1; + for(auto p : primes_) + { mpz_mul_ui(prod_characteristics_.get_mpz_t(), + prod_characteristics_.get_mpz_t(), + p); + } + + num_primes_ = primes_.size(); + + //Uvect_ + Element Ui; Element tmp_elem; + for(auto p : primes_) + { + tmp_elem = prod_characteristics_ / p; + //Element tmp_elem_bis = 10; + mpz_powm_ui ( tmp_elem.get_mpz_t() + , tmp_elem.get_mpz_t() + , p - 1 + , prod_characteristics_.get_mpz_t() ); + Uvect_.push_back(tmp_elem); + } + mult_id_all = 0; + for(int idx = 0; idx < num_primes_; ++idx) + { mult_id_all = (mult_id_all + Uvect_[idx]) % prod_characteristics_; } + + } + + void clear_coefficient(Element & x) { mpz_clear(x.get_mpz_t()); } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + Element additive_identity () { return 0; } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + Element multiplicative_identity () { return mult_id_all; }// 1 everywhere + + Element multiplicative_identity (Element Q) + { + if(Q == prod_characteristics_) { return multiplicative_identity(); } + + Element mult_id = 0; + for(int idx = 0; idx < num_primes_; ++idx) { + if( (Q % primes_[idx]) == 0 ) + { mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; } + } + return mult_id; + } + + /** Returns y * w */ + Element times ( Element y, int w ) { + Element tmp = (y*w) % prod_characteristics_; + if(tmp < 0) return prod_characteristics_ + tmp; + return tmp; + } + + void plus_equal(Element & x, Element y) + { x += y; x %= prod_characteristics_; } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + Element characteristic() { return prod_characteristics_; } + + /** Returns the inverse in the field. Modifies P.*/ + std::pair<Element,Element> inverse ( Element x + , Element QS ) + { + Element QR; + mpz_gcd( QR.get_mpz_t(), x.get_mpz_t(), QS.get_mpz_t() ); // QR <- gcd(x,QS) + if( QR == QS ) return std::pair<Element,Element>(additive_identity() + , multiplicative_identity() ); //partial inverse is 0 + Element QT = QS / QR; + Element inv_qt; + mpz_invert(inv_qt.get_mpz_t(), x.get_mpz_t(), QT.get_mpz_t()); + + return std::pair<Element,Element>( + (inv_qt * multiplicative_identity(QT)) % prod_characteristics_ + , QT ); + } + /** Returns -x * y.*/ + Element times_minus ( Element x, Element y ) + { return prod_characteristics_ - ((x*y)%prod_characteristics_); } + + /** Set x <- x + w * y*/ + void plus_times_equal ( Element & x, Element y, Element w ) + { x = (x + w * y) % prod_characteristics_; } + + Element prod_characteristics_; //product of characteristics of the fields + //represented by the multi-field class + std::vector<int> primes_; //all the characteristics of the fields + std::vector<Element> Uvect_; + size_t num_primes_; //number of fields + Element mult_id_all; + +}; + +} // namespace GUDHI + +#endif // GUDHI_MULTI_FIELD_H diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h new file mode 100644 index 00000000..0702c58e --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h @@ -0,0 +1,153 @@ + /* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria + * + * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GUDHI_COLUMN_LIST_H +#define GUDHI_COLUMN_LIST_H + +#include "boost/tuple/tuple.hpp" +#include "boost/intrusive/set.hpp" +#include "boost/intrusive/list.hpp" + +namespace Gudhi{ + +template < typename SimplexKey + , typename ArithmeticElement + > +class Persistent_cohomology_column; + +struct cam_h_tag; // for horizontal traversal in the CAM +struct cam_v_tag; // for vertical traversal in the CAM + +typedef boost::intrusive::list_base_hook + < boost::intrusive::tag < cam_h_tag > + , boost::intrusive::link_mode < boost::intrusive::auto_unlink > //allows .unlink() + > base_hook_cam_h; + +typedef boost::intrusive::list_base_hook + < boost::intrusive::tag < cam_v_tag > + , boost::intrusive::link_mode < boost::intrusive::normal_link > //faster hook, less safe + > base_hook_cam_v; + +/** \internal + * \brief + * + */ +template < typename SimplexKey + , typename ArithmeticElement + > +class Persistent_cohomology_cell +: public base_hook_cam_h +, public base_hook_cam_v +{ + public: + template < class T1, class T2 > friend class Persistent_cohomology; + friend class Persistent_cohomology_column < SimplexKey , ArithmeticElement >; + + typedef Persistent_cohomology_column< SimplexKey, ArithmeticElement > Column; + + Persistent_cohomology_cell( SimplexKey key + , ArithmeticElement x + , Column * self_col) + : key_(key) + , coefficient_(x) + , self_col_(self_col) {} + + SimplexKey key_; + ArithmeticElement coefficient_; + Column * self_col_; +}; + + + + + +/* + * \brief Sparse column for the Compressed Annotation Matrix. + * + * The non-zero coefficients of the column are stored in a + * boost::intrusive::list. Contains a hook to be stored in a + * boost::intrusive::set. + */ +template < typename SimplexKey + , typename ArithmeticElement > +class Persistent_cohomology_column +: public boost::intrusive::set_base_hook + < boost::intrusive::link_mode< boost::intrusive::normal_link > > +{ +private: + template < class T1, class T2 > friend class Persistent_cohomology; + + typedef Persistent_cohomology_cell < SimplexKey, ArithmeticElement > Cell; + typedef boost::intrusive::list < Cell + , boost::intrusive::constant_time_size<false> + , boost::intrusive::base_hook< base_hook_cam_v > + > Col_type; + +/** \brief Creates an empty column.*/ + Persistent_cohomology_column (SimplexKey key) + { + class_key_ = key; + col_ = Col_type(); + } +public: + /** Copy constructor.*/ + Persistent_cohomology_column( Persistent_cohomology_column const &other ) + : col_() + , class_key_(other.class_key_) + { if(!other.col_.empty()) std::cerr << "Copying a non-empty column.\n"; } + +/** \brief Returns true iff the column is null.*/ + bool is_null() { return col_.empty(); } +/** \brief Returns the key of the representative simplex of + * the set of simplices having this column as annotation vector + * in the compressed annotation matrix.*/ + SimplexKey class_key () { return class_key_; } + +/** \brief Lexicographic comparison of two columns.*/ +friend bool operator< ( const Persistent_cohomology_column& c1 + , const Persistent_cohomology_column& c2) + { + typename Col_type::const_iterator it1 = c1.col_.begin(); + typename Col_type::const_iterator it2 = c2.col_.begin(); + while(it1 != c1.col_.end() && it2 != c2.col_.end()) + { + if(it1->key_ == it2->key_) + { if(it1->coefficient_ == it2->coefficient_) { ++it1; ++it2; } + else { return it1->coefficient_ < it2->coefficient_; } } + else { return it1->key_ < it2->key_; } + } + return (it2 != c2.col_.end()); + } + + // void display() + // { + // for(auto cell : col_) + // { std::cout << "(" << cell.key_ <<":"<<cell.coefficient_<<") "; } + // } + + Col_type col_; + SimplexKey class_key_; +}; + +} // namespace GUDHI + +#endif // GUDHI_COLUMN_LIST_H |