Initial commit.

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;
+}