1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
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;
}
|
