diff options
Diffstat (limited to 'src/Persistent_cohomology/include/gudhi/Persistent_cohomology')
3 files changed, 405 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..0673625c --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Field_Zp.h @@ -0,0 +1,104 @@ +/* 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): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ +#define PERSISTENT_COHOMOLOGY_FIELD_ZP_H_ + +#include <utility> +#include <vector> + +namespace Gudhi { + +namespace persistent_cohomology { + +/** \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(0), + inverse_() { + } + + void init(int charac) { + assert(charac > 0); // division by zero + non negative values + 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*/ + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(Prime > 0); // division by zero + non negative values + Element result = (x + w * y) % Prime; + if (result < 0) + result += Prime; + return result; + } + +// operator= defined on Element + + /** Returns y * w */ + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, (Element)w); + } + + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + Element additive_identity() const { + return 0; + } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + Element multiplicative_identity(Element = 0) const { + 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) { + assert(Prime > 0); // division by zero + non negative values + Element out = (-x * y) % Prime; + return (out < 0) ? out + Prime : out; + } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + int characteristic() const { + return Prime; + } + + private: + int Prime; + /** Property map Element -> Element, which associate to an element its inverse in the field.*/ + std::vector<Element> inverse_; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_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..1754a2ec --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Multi_field.h @@ -0,0 +1,173 @@ +/* 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): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ +#define PERSISTENT_COHOMOLOGY_MULTI_FIELD_H_ + +#include <gmpxx.h> + +#include <vector> +#include <utility> + +namespace Gudhi { + +namespace persistent_cohomology { + +/** \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() + : prod_characteristics_(0), + mult_id_all(0), + add_id_all(0) { + } + + /* 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 + 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); + } + mpz_clear(tmp_prime); + // set m to primorial(bound_prime) + prod_characteristics_ = 1; + for (auto p : primes_) { + prod_characteristics_ *= p; + } + + // Uvect_ + Element Ui; + Element tmp_elem; + for (auto p : primes_) { + assert(p > 0); // division by zero + non negative values + 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 (auto uvect : Uvect_) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + mult_id_all = (mult_id_all + uvect) % prod_characteristics_; + } + } + + /** \brief Returns the additive idendity \f$0_{\Bbbk}\f$ of the field.*/ + const Element& additive_identity() const { + return add_id_all; + } + /** \brief Returns the multiplicative identity \f$1_{\Bbbk}\f$ of the field.*/ + const Element& multiplicative_identity() const { + return mult_id_all; + } // 1 everywhere + + Element multiplicative_identity(Element Q) { + if (Q == prod_characteristics_) { + return multiplicative_identity(); + } + + assert(prod_characteristics_ > 0); // division by zero + non negative values + Element mult_id = 0; + for (unsigned int idx = 0; idx < primes_.size(); ++idx) { + assert(primes_[idx] > 0); // division by zero + non negative values + if ((Q % primes_[idx]) == 0) { + mult_id = (mult_id + Uvect_[idx]) % prod_characteristics_; + } + } + return mult_id; + } + + /** Returns y * w */ + Element times(const Element& y, const Element& w) { + return plus_times_equal(0, y, w); + } + + Element plus_equal(const Element& x, const Element& y) { + return plus_times_equal(x, y, (Element)1); + } + + /** \brief Returns the characteristic \f$p\f$ of the field.*/ + const Element& characteristic() const { + 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()); + + assert(prod_characteristics_ > 0); // division by zero + non negative values + return { (inv_qt * multiplicative_identity(QT)) % prod_characteristics_, QT }; + } + /** Returns -x * y.*/ + Element times_minus(const Element& x, const Element& y) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + /* This assumes that (x*y)%pc cannot be zero, but Field_Zp has specific code for the 0 case ??? */ + return prod_characteristics_ - ((x * y) % prod_characteristics_); + } + + /** Set x <- x + w * y*/ + Element plus_times_equal(const Element& x, const Element& y, const Element& w) { + assert(prod_characteristics_ > 0); // division by zero + non negative values + Element result = (x + w * y) % prod_characteristics_; + if (result < 0) + result += prod_characteristics_; + return result; + } + + 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_; + Element mult_id_all; + const Element add_id_all; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_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..480be389 --- /dev/null +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h @@ -0,0 +1,128 @@ +/* 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): Clément Maria + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ +#define PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ + +#include <boost/intrusive/set.hpp> +#include <boost/intrusive/list.hpp> + +#include <list> + +namespace Gudhi { + +namespace persistent_cohomology { + +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. + * + * Movable but not Copyable. + */ +template<typename SimplexKey, typename ArithmeticElement> +class Persistent_cohomology_column : public boost::intrusive::set_base_hook< + boost::intrusive::link_mode<boost::intrusive::normal_link> > { + template<class T1, class T2> friend class Persistent_cohomology; + + public: + 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.*/ + explicit Persistent_cohomology_column(SimplexKey key) + : col_(), + class_key_(key) {} + + /** \brief Returns true iff the column is null.*/ + bool is_null() const { + 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() const { + 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()); + } + + Col_type col_; + SimplexKey class_key_; +}; + +} // namespace persistent_cohomology + +} // namespace Gudhi + +#endif // PERSISTENT_COHOMOLOGY_PERSISTENT_COHOMOLOGY_COLUMN_H_ |