diff options
Diffstat (limited to 'include')
-rw-r--r-- | include/misc.hpp | 66 | ||||
-rw-r--r-- | include/pd.hpp | 126 |
2 files changed, 192 insertions, 0 deletions
diff --git a/include/misc.hpp b/include/misc.hpp new file mode 100644 index 0000000..f8afe77 --- /dev/null +++ b/include/misc.hpp @@ -0,0 +1,66 @@ +#pragma once + +#include <algorithm> +#include <cmath> +#include <cstdint> +#include <sstream> + +constexpr int64_t DIPHA_MAGIC = 8067171840; +constexpr int64_t DIPHA_DISTANCE_MATRIX = 7; +constexpr int64_t DIPHA_PERSISTENCE_DIAGRAM = 2; +constexpr int64_t DIPHA_WEIGHTED_BOUNDARY_MATRIX = 0; +inline constexpr double PI = 4*std::atan(1); + +template <typename T> inline void reverse_endianness(T & x) +{ + uint8_t * p = reinterpret_cast<uint8_t *>(&x); + std::reverse(p, p + sizeof(T)); +} + +template <typename T> inline T read_le(std::istream & s) +{ + T result; + s.read(reinterpret_cast<char *>(&result), sizeof(T)); + #ifdef BIGENDIAN + reverse_endianness(result); + #endif + return result; +} + +template <typename T> inline void write_le(std::ostream & s, T value) +{ + #ifdef BIGENDIAN + reverse_endianness(value); + #endif + s.write(reinterpret_cast<char *>(&value), sizeof(T)); +} + +template <typename T> inline T read_be(std::istream & s) +{ + T result; + s.read(reinterpret_cast<char *>(&result), sizeof(T)); + #ifndef BIGENDIAN + reverse_endianness(result); + #endif + return result; +} + +template <typename T> inline void write_be(std::ostream & s, T value) +{ + #ifndef BIGENDIAN + reverse_endianness(value); + #endif + s.write(reinterpret_cast<char *>(&value), sizeof(T)); +} + +enum class Message_tag : int { result, work }; + +template <typename T> inline T unroll_upper_tri(T n, T i, T j) +{ + return i*(n-1) - ((i-1)*i)/2 + j; +} + +template <typename T> inline T unroll_full(T n, T i, T j) +{ + return i*n + j; +} 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; +} |