first version: computing coboundary indices

2 files changed, 671 insertions, 0 deletions

diff --git a/prettyprint.hpp b/prettyprint.hpp
new file mode 100644
index 0000000..6bf2543
--- /dev/null
+++ b/prettyprint.hpp
@@ -0,0 +1,445 @@
+//          Copyright Louis Delacroix 2010 - 2014.
+// Distributed under the Boost Software License, Version 1.0.
+//    (See accompanying file LICENSE_1_0.txt or copy at
+//          http://www.boost.org/LICENSE_1_0.txt)
+//
+// A pretty printing library for C++
+//
+// Usage:
+// Include this header, and operator<< will "just work".
+
+#ifndef H_PRETTY_PRINT
+#define H_PRETTY_PRINT
+
+#include <cstddef>
+#include <iterator>
+#include <memory>
+#include <ostream>
+#include <set>
+#include <tuple>
+#include <type_traits>
+#include <unordered_set>
+#include <utility>
+#include <valarray>
+
+namespace pretty_print
+{
+    namespace detail
+    {
+        // SFINAE type trait to detect whether T::const_iterator exists.
+
+        struct sfinae_base
+        {
+            using yes = char;
+            using no  = yes[2];
+        };
+
+        template <typename T>
+        struct has_const_iterator : private sfinae_base
+        {
+        private:
+            template <typename C> static yes & test(typename C::const_iterator*);
+            template <typename C> static no  & test(...);
+        public:
+            static const bool value = sizeof(test<T>(nullptr)) == sizeof(yes);
+            using type =  T;
+        };
+
+        template <typename T>
+        struct has_begin_end : private sfinae_base
+        {
+        private:
+            template <typename C>
+            static yes & f(typename std::enable_if<
+                std::is_same<decltype(static_cast<typename C::const_iterator(C::*)() const>(&C::begin)),
+                             typename C::const_iterator(C::*)() const>::value>::type *);
+
+            template <typename C> static no & f(...);
+
+            template <typename C>
+            static yes & g(typename std::enable_if<
+                std::is_same<decltype(static_cast<typename C::const_iterator(C::*)() const>(&C::end)),
+                             typename C::const_iterator(C::*)() const>::value, void>::type*);
+
+            template <typename C> static no & g(...);
+
+        public:
+            static bool const beg_value = sizeof(f<T>(nullptr)) == sizeof(yes);
+            static bool const end_value = sizeof(g<T>(nullptr)) == sizeof(yes);
+        };
+
+    }  // namespace detail
+
+
+    // Holds the delimiter values for a specific character type
+
+    template <typename TChar>
+    struct delimiters_values
+    {
+        using char_type = TChar;
+        const char_type * prefix;
+        const char_type * delimiter;
+        const char_type * postfix;
+    };
+
+
+    // Defines the delimiter values for a specific container and character type
+
+    template <typename T, typename TChar>
+    struct delimiters
+    {
+        using type = delimiters_values<TChar>;
+        static const type values; 
+    };
+
+
+    // Functor to print containers. You can use this directly if you want
+    // to specificy a non-default delimiters type. The printing logic can
+    // be customized by specializing the nested template.
+
+    template <typename T,
+              typename TChar = char,
+              typename TCharTraits = ::std::char_traits<TChar>,
+              typename TDelimiters = delimiters<T, TChar>>
+    struct print_container_helper
+    {
+        using delimiters_type = TDelimiters;
+        using ostream_type = std::basic_ostream<TChar, TCharTraits>;
+
+        template <typename U>
+        struct printer
+        {
+            static void print_body(const U & c, ostream_type & stream)
+            {
+                using std::begin;
+                using std::end;
+
+                auto it = begin(c);
+                const auto the_end = end(c);
+
+                if (it != the_end)
+                {
+                    for ( ; ; )
+                    {
+                        stream << *it;
+
+                    if (++it == the_end) break;
+
+                    if (delimiters_type::values.delimiter != NULL)
+                        stream << delimiters_type::values.delimiter;
+                    }
+                }
+            }
+        };
+
+        print_container_helper(const T & container)
+        : container_(container)
+        { }
+
+        inline void operator()(ostream_type & stream) const
+        {
+            if (delimiters_type::values.prefix != NULL)
+                stream << delimiters_type::values.prefix;
+
+            printer<T>::print_body(container_, stream);
+
+            if (delimiters_type::values.postfix != NULL)
+                stream << delimiters_type::values.postfix;
+        }
+
+    private:
+        const T & container_;
+    };
+
+    // Specialization for pairs
+
+    template <typename T, typename TChar, typename TCharTraits, typename TDelimiters>
+    template <typename T1, typename T2>
+    struct print_container_helper<T, TChar, TCharTraits, TDelimiters>::printer<std::pair<T1, T2>>
+    {
+        using ostream_type = print_container_helper<T, TChar, TCharTraits, TDelimiters>::ostream_type;
+
+        static void print_body(const std::pair<T1, T2> & c, ostream_type & stream)
+        {
+            stream << c.first;
+            if (print_container_helper<T, TChar, TCharTraits, TDelimiters>::delimiters_type::values.delimiter != NULL)
+                stream << print_container_helper<T, TChar, TCharTraits, TDelimiters>::delimiters_type::values.delimiter;
+            stream << c.second;
+        }
+    };
+
+    // Specialization for tuples
+
+    template <typename T, typename TChar, typename TCharTraits, typename TDelimiters>
+    template <typename ...Args>
+    struct print_container_helper<T, TChar, TCharTraits, TDelimiters>::printer<std::tuple<Args...>>
+    {
+        using ostream_type = print_container_helper<T, TChar, TCharTraits, TDelimiters>::ostream_type;
+        using element_type = std::tuple<Args...>;
+
+        template <std::size_t I> struct Int { };
+
+        static void print_body(const element_type & c, ostream_type & stream)
+        {
+            tuple_print(c, stream, Int<0>());
+        }
+
+        static void tuple_print(const element_type &, ostream_type &, Int<sizeof...(Args)>)
+        {
+        }
+
+        static void tuple_print(const element_type & c, ostream_type & stream,
+                                typename std::conditional<sizeof...(Args) != 0, Int<0>, std::nullptr_t>::type)
+        {
+            stream << std::get<0>(c);
+            tuple_print(c, stream, Int<1>());
+        }
+
+        template <std::size_t N>
+        static void tuple_print(const element_type & c, ostream_type & stream, Int<N>)
+        {
+            if (print_container_helper<T, TChar, TCharTraits, TDelimiters>::delimiters_type::values.delimiter != NULL)
+                stream << print_container_helper<T, TChar, TCharTraits, TDelimiters>::delimiters_type::values.delimiter;
+
+            stream << std::get<N>(c);
+
+            tuple_print(c, stream, Int<N + 1>());
+        }
+    };
+
+    // Prints a print_container_helper to the specified stream.
+
+    template<typename T, typename TChar, typename TCharTraits, typename TDelimiters>
+    inline std::basic_ostream<TChar, TCharTraits> & operator<<(
+        std::basic_ostream<TChar, TCharTraits> & stream,
+        const print_container_helper<T, TChar, TCharTraits, TDelimiters> & helper)
+    {
+        helper(stream);
+        return stream;
+    }
+
+
+    // Basic is_container template; specialize to derive from std::true_type for all desired container types
+
+    template <typename T>
+    struct is_container : public std::integral_constant<bool,
+                                                        detail::has_const_iterator<T>::value &&
+                                                        detail::has_begin_end<T>::beg_value  &&
+                                                        detail::has_begin_end<T>::end_value> { };
+
+    template <typename T, std::size_t N>
+    struct is_container<T[N]> : std::true_type { };
+
+    template <std::size_t N>
+    struct is_container<char[N]> : std::false_type { };
+
+    template <typename T>
+    struct is_container<std::valarray<T>> : std::true_type { };
+
+    template <typename T1, typename T2>
+    struct is_container<std::pair<T1, T2>> : std::true_type { };
+
+    template <typename ...Args>
+    struct is_container<std::tuple<Args...>> : std::true_type { };
+
+
+    // Default delimiters
+
+    template <typename T> struct delimiters<T, char> { static const delimiters_values<char> values; };
+    template <typename T> const delimiters_values<char> delimiters<T, char>::values = { "[", ", ", "]" };
+    template <typename T> struct delimiters<T, wchar_t> { static const delimiters_values<wchar_t> values; };
+    template <typename T> const delimiters_values<wchar_t> delimiters<T, wchar_t>::values = { L"[", L", ", L"]" };
+
+
+    // Delimiters for (multi)set and unordered_(multi)set
+
+    template <typename T, typename TComp, typename TAllocator>
+    struct delimiters< ::std::set<T, TComp, TAllocator>, char> { static const delimiters_values<char> values; };
+
+    template <typename T, typename TComp, typename TAllocator>
+    const delimiters_values<char> delimiters< ::std::set<T, TComp, TAllocator>, char>::values = { "{", ", ", "}" };
+
+    template <typename T, typename TComp, typename TAllocator>
+    struct delimiters< ::std::set<T, TComp, TAllocator>, wchar_t> { static const delimiters_values<wchar_t> values; };
+
+    template <typename T, typename TComp, typename TAllocator>
+    const delimiters_values<wchar_t> delimiters< ::std::set<T, TComp, TAllocator>, wchar_t>::values = { L"{", L", ", L"}" };
+
+    template <typename T, typename TComp, typename TAllocator>
+    struct delimiters< ::std::multiset<T, TComp, TAllocator>, char> { static const delimiters_values<char> values; };
+
+    template <typename T, typename TComp, typename TAllocator>
+    const delimiters_values<char> delimiters< ::std::multiset<T, TComp, TAllocator>, char>::values = { "{", ", ", "}" };
+
+    template <typename T, typename TComp, typename TAllocator>
+    struct delimiters< ::std::multiset<T, TComp, TAllocator>, wchar_t> { static const delimiters_values<wchar_t> values; };
+
+    template <typename T, typename TComp, typename TAllocator>
+    const delimiters_values<wchar_t> delimiters< ::std::multiset<T, TComp, TAllocator>, wchar_t>::values = { L"{", L", ", L"}" };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    struct delimiters< ::std::unordered_set<T, THash, TEqual, TAllocator>, char> { static const delimiters_values<char> values; };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    const delimiters_values<char> delimiters< ::std::unordered_set<T, THash, TEqual, TAllocator>, char>::values = { "{", ", ", "}" };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    struct delimiters< ::std::unordered_set<T, THash, TEqual, TAllocator>, wchar_t> { static const delimiters_values<wchar_t> values; };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    const delimiters_values<wchar_t> delimiters< ::std::unordered_set<T, THash, TEqual, TAllocator>, wchar_t>::values = { L"{", L", ", L"}" };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    struct delimiters< ::std::unordered_multiset<T, THash, TEqual, TAllocator>, char> { static const delimiters_values<char> values; };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    const delimiters_values<char> delimiters< ::std::unordered_multiset<T, THash, TEqual, TAllocator>, char>::values = { "{", ", ", "}" };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    struct delimiters< ::std::unordered_multiset<T, THash, TEqual, TAllocator>, wchar_t> { static const delimiters_values<wchar_t> values; };
+
+    template <typename T, typename THash, typename TEqual, typename TAllocator>
+    const delimiters_values<wchar_t> delimiters< ::std::unordered_multiset<T, THash, TEqual, TAllocator>, wchar_t>::values = { L"{", L", ", L"}" };
+
+
+    // Delimiters for pair and tuple
+
+    template <typename T1, typename T2> struct delimiters<std::pair<T1, T2>, char> { static const delimiters_values<char> values; };
+    template <typename T1, typename T2> const delimiters_values<char> delimiters<std::pair<T1, T2>, char>::values = { "(", ", ", ")" };
+    template <typename T1, typename T2> struct delimiters< ::std::pair<T1, T2>, wchar_t> { static const delimiters_values<wchar_t> values; };
+    template <typename T1, typename T2> const delimiters_values<wchar_t> delimiters< ::std::pair<T1, T2>, wchar_t>::values = { L"(", L", ", L")" };
+
+    template <typename ...Args> struct delimiters<std::tuple<Args...>, char> { static const delimiters_values<char> values; };
+    template <typename ...Args> const delimiters_values<char> delimiters<std::tuple<Args...>, char>::values = { "(", ", ", ")" };
+    template <typename ...Args> struct delimiters< ::std::tuple<Args...>, wchar_t> { static const delimiters_values<wchar_t> values; };
+    template <typename ...Args> const delimiters_values<wchar_t> delimiters< ::std::tuple<Args...>, wchar_t>::values = { L"(", L", ", L")" };
+
+
+    // Type-erasing helper class for easy use of custom delimiters.
+    // Requires TCharTraits = std::char_traits<TChar> and TChar = char or wchar_t, and MyDelims needs to be defined for TChar.
+    // Usage: "cout << pretty_print::custom_delims<MyDelims>(x)".
+
+    struct custom_delims_base
+    {
+        virtual ~custom_delims_base() { }
+        virtual std::ostream & stream(::std::ostream &) = 0;
+        virtual std::wostream & stream(::std::wostream &) = 0;
+    };
+
+    template <typename T, typename Delims>
+    struct custom_delims_wrapper : custom_delims_base
+    {
+        custom_delims_wrapper(const T & t_) : t(t_) { }
+
+        std::ostream & stream(std::ostream & s)
+        {
+            return s << print_container_helper<T, char, std::char_traits<char>, Delims>(t);
+        }
+
+        std::wostream & stream(std::wostream & s)
+        {
+            return s << print_container_helper<T, wchar_t, std::char_traits<wchar_t>, Delims>(t);
+        }
+
+    private:
+        const T & t;
+    };
+
+    template <typename Delims>
+    struct custom_delims
+    {
+        template <typename Container>
+        custom_delims(const Container & c) : base(new custom_delims_wrapper<Container, Delims>(c)) { }
+
+        std::unique_ptr<custom_delims_base> base;
+    };
+
+    template <typename TChar, typename TCharTraits, typename Delims>
+    inline std::basic_ostream<TChar, TCharTraits> & operator<<(std::basic_ostream<TChar, TCharTraits> & s, const custom_delims<Delims> & p)
+    {
+        return p.base->stream(s);
+    }
+
+
+    // A wrapper for a C-style array given as pointer-plus-size.
+    // Usage: std::cout << pretty_print_array(arr, n) << std::endl;
+
+    template<typename T>
+    struct array_wrapper_n
+    {
+        typedef const T * const_iterator;
+        typedef T value_type;
+
+        array_wrapper_n(const T * const a, size_t n) : _array(a), _n(n) { }
+        inline const_iterator begin() const { return _array; }
+        inline const_iterator end() const { return _array + _n; }
+
+    private:
+        const T * const _array;
+        size_t _n;
+    };
+
+
+    // A wrapper for hash-table based containers that offer local iterators to each bucket.
+    // Usage: std::cout << bucket_print(m, 4) << std::endl;  (Prints bucket 5 of container m.)
+
+    template <typename T>
+    struct bucket_print_wrapper
+    {
+        typedef typename T::const_local_iterator const_iterator;
+        typedef typename T::size_type size_type;
+
+        const_iterator begin() const
+        {
+            return m_map.cbegin(n);
+        }
+
+        const_iterator end() const
+        {
+            return m_map.cend(n);
+        }
+
+        bucket_print_wrapper(const T & m, size_type bucket) : m_map(m), n(bucket) { }
+
+    private:
+        const T & m_map;
+        const size_type n;
+    };
+
+}   // namespace pretty_print
+
+
+// Global accessor functions for the convenience wrappers
+
+template<typename T>
+inline pretty_print::array_wrapper_n<T> pretty_print_array(const T * const a, size_t n)
+{
+    return pretty_print::array_wrapper_n<T>(a, n);
+}
+
+template <typename T> pretty_print::bucket_print_wrapper<T>
+bucket_print(const T & m, typename T::size_type n)
+{
+    return pretty_print::bucket_print_wrapper<T>(m, n);
+}
+
+
+// Main magic entry point: An overload snuck into namespace std.
+// Can we do better?
+
+namespace std
+{
+    // Prints a container to the stream using default delimiters
+
+    template<typename T, typename TChar, typename TCharTraits>
+    inline typename enable_if< ::pretty_print::is_container<T>::value,
+                              basic_ostream<TChar, TCharTraits> &>::type
+    operator<<(basic_ostream<TChar, TCharTraits> & stream, const T & container)
+    {
+        return stream << ::pretty_print::print_container_helper<T, TChar, TCharTraits>(container);
+    }
+}
+
+
+
+#endif  // H_PRETTY_PRINT
diff --git a/ripser.cpp b/ripser.cpp
new file mode 100644
index 0000000..450e4b1
--- /dev/null
+++ b/ripser.cpp
@@ -0,0 +1,226 @@
+#include <vector>
+#include <iostream>
+#include <cassert>
+#include "prettyprint.hpp"
+//#include <boost/numeric/ublas/symmetric.hpp>
+
+
+
+class binomial_coeff_table {
+    std::vector<std::vector<int> > B;
+    int n_max, k_max;
+    
+public:
+    binomial_coeff_table(int n, int k) {
+        n_max = n;
+        k_max = k;
+    
+        B.resize(n + 1);
+        for( int i = 0; i <= n; i++ ) {
+            B[i].resize(k + 1);
+            for ( int j = 0; j <= std::min(i, k); j++ )
+            {
+                if (j == 0 || j == i)
+                    B[i][j] = 1;
+
+                else
+                    B[i][j] = B[i-1][j-1] + B[i-1][j];
+            }
+        }
+    }
+    
+    int operator()(int n, int k) const {
+//        std::cout << "B(" << n << "," << k << ")\n";
+
+        return B[n][k];
+    }
+};
+
+
+
+template<typename OutputIterator>
+OutputIterator get_simplex_vertices( int idx, int dim, int n, const binomial_coeff_table& binomial_coeff, OutputIterator out  )
+{
+    --n;
+    
+    for( int k = dim + 1; k > 0; --k ) {
+        while( binomial_coeff( n , k ) > idx ) {
+            --n;
+        }
+        *out++ = n;
+        
+        idx -= binomial_coeff( n , k );
+
+    }
+
+    return out;
+}
+
+template <typename DistanceMatrix>
+class rips_filtration_comparator {
+private:
+    std::vector<int> vertices;
+    
+    const binomial_coeff_table& binomial_coeff;
+    
+public:
+    const DistanceMatrix& dist;
+    
+    const int dim;
+
+    rips_filtration_comparator(const DistanceMatrix& _dist, const int _dim, const binomial_coeff_table& _binomial_coeff): dist(_dist), dim(_dim), vertices(_dim + 1),
+    binomial_coeff(_binomial_coeff) {};
+
+    bool operator()(const int a, const int b)
+    {
+        assert(a < binomial_coeff(dist.size(), dim + 1));
+        assert(b < binomial_coeff(dist.size(), dim + 1));
+        
+        decltype(dist(0,0)) a_diam = 0, b_diam = 0;
+
+    
+        //vertices.clear(); //std::back_inserter(vertices)
+        get_simplex_vertices(a, dim, dist.size(), binomial_coeff, vertices.begin() );
+        
+        for (int i = 0; i <= dim; ++i)
+            for (int j = i + 1; j <= dim; ++j) {
+                auto d = dist(vertices[i], vertices[j]);
+                a_diam = std::max(a_diam, dist(vertices[i], vertices[j]));
+//                std::cout << " i = " << i << " j = " << j << " d = " << d << std::endl;
+//                std::cout << " a_vertices[i] = " << vertices[i] << " a_vertices[j] = " << vertices[j] << std::endl;
+//                std::cout << " a_diam = " << a_diam << std::endl;
+            }
+        
+//        std::cout << "a_diam = " << a_diam << std::endl;
+        
+        //vertices.clear(); //std::back_inserter(vertices)
+        get_simplex_vertices(b, dim, dist.size(), binomial_coeff, vertices.begin() );
+
+        for (int i = 0; i <= dim; ++i)
+            for (int j = i + 1; j <= dim; ++j) {
+                b_diam = std::max(b_diam, dist(vertices[i], vertices[j]));
+                if (a_diam < b_diam)
+                    return true;
+            }
+        
+//        std::cout << "b_diam = " << b_diam << std::endl;
+       
+        return (a_diam == b_diam && a <= b);
+    }
+    
+};
+
+
+template<typename OutputIterator>
+void get_simplex_coboundary( int idx, int dim, int n, const binomial_coeff_table& binomial_coeff, OutputIterator coboundary )
+{
+
+    --n;
+    
+    int modified_idx = idx;
+    
+    for( int k = dim + 1; k >= 0 && n >= 0; --k ) {
+        while( binomial_coeff( n , k ) > idx ) {
+//            std::cout << "binomial_coeff(" << n << ", " << k << ") = " << binomial_coeff( n , k ) << " > " << idx << std::endl;
+            *coboundary++ = modified_idx + binomial_coeff( n , k + 1 );
+            if (n==0) break;
+            --n;
+        }
+        
+        idx -= binomial_coeff( n , k );
+            
+        modified_idx -= binomial_coeff( n , k );
+        modified_idx += binomial_coeff( n , k + 1 );
+        
+        --n;
+    }
+ 
+    return;
+}
+
+
+
+class distance_matrix  {
+public:
+    typedef double value_type;
+    std::vector<std::vector<double> > distances;
+    double operator()(const int a, const int b) const {
+        return distances[a][b];
+    }
+    
+    size_t size() const {
+        return distances.size();
+    }
+
+};
+
+
+int main() {
+
+    distance_matrix dist;
+    dist.distances = {
+        {0,1,2,3},
+        {1,0,1,2},
+        {2,1,0,1},
+        {3,2,1,0}
+    };
+    dist.distances = {
+        {0,1,1,1,1},
+        {1,0,1,1,1},
+        {1,1,0,1,1},
+        {1,1,1,0,1},
+        {1,1,1,1,0}
+    };
+    
+    binomial_coeff_table binomial_coeff(dist.size(), 4);
+    
+    
+    std::cout << dist (0,1) << std::endl;
+ 
+    rips_filtration_comparator<distance_matrix> comp1(dist, 1, binomial_coeff);
+    rips_filtration_comparator<distance_matrix> comp2(dist, 2, binomial_coeff);
+    
+    std::cout << (comp1(0,1) ? "0<1" : "0>=1") << std::endl;
+    std::cout << (comp1(1,0) ? "1<0" : "1>=0") << std::endl;
+    
+    std::cout << (comp1(0,2) ? "0<2" : "0>=2") << std::endl;
+    std::cout << (comp1(1,2) ? "1<2" : "1>=2") << std::endl;
+    
+    std::vector<int> edges = {0,1,2,3,4,5};
+    
+    std::sort(edges.begin(), edges.end(), comp1);
+    
+    std::cout << "sorted edges: " << edges << std::endl;
+   
+   
+    std::vector<int> triangles = {0,1,2,3};
+    
+    std::sort(triangles.begin(), triangles.end(), comp2);
+    
+    std::cout << "sorted triangles: " << triangles << std::endl;
+    
+    
+    int dim = 1;
+    int simplex_index = 2;
+    
+    
+    std::vector<int> vertices;
+    
+    get_simplex_vertices( simplex_index, dim, dist.size(), binomial_coeff, std::back_inserter(vertices) );
+
+    
+    std::cout << "coboundary of simplex " << vertices << ":" << std::endl;
+    
+    std::vector<int> coboundary;
+    get_simplex_coboundary( simplex_index, dim, dist.size(), binomial_coeff, std::back_inserter(coboundary) );
+    
+    
+    for (int coboundary_simplex_index: coboundary) {
+        std::vector<int> vertices;
+    
+        get_simplex_vertices( coboundary_simplex_index, dim + 1, dist.size(), binomial_coeff, std::back_inserter(vertices) );
+        std::cout << " " << vertices << std::endl;
+
+    }
+
+}
\ No newline at end of file