diff options
Diffstat (limited to 'include/pd.hpp')
-rw-r--r-- | include/pd.hpp | 126 |
1 files changed, 126 insertions, 0 deletions
diff --git a/include/pd.hpp b/include/pd.hpp new file mode 100644 index 0000000..04b8609 --- /dev/null +++ b/include/pd.hpp @@ -0,0 +1,126 @@ +#pragma once + +#include <cstdint> +#include <fstream> +#include <limits> +#include <string> +#include <vector> + +#include "misc.hpp" + +template <typename T> class PD +{ +public: + class Interval + { + public: + T birth; + T death; + }; + + + PD() : intervals() + { + }; + + void add(T b, T d, unsigned int m) + { + if (m > 0 && b < d) + { + Interval interval; + interval.birth = b; + interval.death = d; + intervals.push_back(std::make_pair(interval, m)); + } + } + + void clear() { intervals.clear(); } + + void finitize(T scale) + { + for (auto it = intervals.begin(); it != intervals.end(); ++it) + { + if (std::isinf(it->first.death)) + it->first.death = scale; + } + } + + typename std::vector<std::pair<PD<T>::Interval, unsigned int> >::size_type size() const { return intervals.size(); } + + using Iterator = typename std::vector<std::pair<PD<T>::Interval, unsigned int> >::const_iterator; + + inline Iterator cbegin() const { return intervals.cbegin(); } + inline Iterator cend() const { return intervals.cend(); } + +private: + std::vector<std::pair<PD<T>::Interval, unsigned int> > intervals; +}; + +inline double sqdist(PD<double>::Interval x, PD<double>::Interval y) +{ + return (x.birth - y.birth)*(x.birth - y.birth) + (x.death - y.death)*(x.death - y.death); +} + +template <typename T> T heat_kernel(T sigma, const PD<T> & a, const PD<T> & b) +{ + T ret = 0; + + for (auto it = a.cbegin(); it != a.cend(); ++it) + { + auto x = *it; + typename PD<T>::Interval xbar; + xbar.birth = x.first.death; + xbar.death = x.first.birth; + + for (auto jt = b.cbegin(); jt != b.cend(); ++jt) + { + auto y = *jt; + ret += x.second*y.second*(std::exp(-sqdist(x.first, y.first)/(8*sigma)) - std::exp(-sqdist(xbar, y.first)/(8*sigma))); + } + } + + ret /= 8*sigma*PI; + + return ret; +} + +int read_dipha_degree(const std::string & fname, unsigned int degree, PD<double> & pd) +{ + std::ifstream f; + f.open(fname, std::ios::in | std::ios::binary); + if (!f.is_open()) + { + std::cerr << "Failed to open " << fname << ". Bailing." << std::endl; + return 1; + } + + if (read_le<int64_t>(f) != DIPHA_MAGIC) + { + std::cerr << "File " << fname << " is not DIPHA file. Bailing." << std::endl; + f.close(); + return 1; + } + + if (read_le<int64_t>(f) != DIPHA_PERSISTENCE_DIAGRAM) + { + std::cerr << "File " << fname << " is not persistence diagram. Bailing." << std::endl; + f.close(); + return 1; + } + + int64_t n = read_le<int64_t>(f); + + for (int64_t i = 0; i < n; ++i) + { + int64_t d = read_le<int64_t>(f); + double birth = read_le<double>(f); + double death = read_le<double>(f); + if (d == degree) + pd.add(birth, death, 1); + else if (-d - 1 == degree) + pd.add(birth, std::numeric_limits<double>::infinity(), 1); + } + + f.close(); + return 0; +} |