forked from jMotif/jmotif-R
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathdiscord.cpp
More file actions
149 lines (114 loc) · 4.61 KB
/
Copy pathdiscord.cpp
File metadata and controls
149 lines (114 loc) · 4.61 KB
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include <RcppArmadillo.h>
using namespace Rcpp ;
//
#include <jmotif.h>
//
// NumericVector subseries(const NumericVector& ts, int start, int end) {
// if(start < 0 || end > ts.length()){
// stop("provided start and stop indexes are invalid.");
// }
// NumericVector res(end-start);
// for (int i=start; i<end; i++) {
// res[i-start] = ts[i];
// }
// return res;
// }
discord_record find_best_discord_brute_force(const NumericVector& series,
int w_size, VisitRegistry* globalRegistry) {
// Rcout << "looking for the best discord, series length " << series.size() << "\n";
double best_so_far_distance = -1.0;
int best_so_far_index = -1;
VisitRegistry outerRegistry(series.size() - w_size);
int outer_idx = outerRegistry.getNextUnvisited();
while(!(-1==outer_idx)){
outerRegistry.markVisited(outer_idx);
if(globalRegistry->isVisited(outer_idx)){
// Rcout << " skipping " << outer_idx << ", marked as visited in global\n";
outer_idx = outerRegistry.getNextUnvisited();
continue;
}
// Rcout << " outer unvisited candidate at " << outer_idx << "\n";
NumericVector candidate_seq = subseries(series, outer_idx, outer_idx + w_size);
double nnDistance = std::numeric_limits<double>::max();
VisitRegistry innerRegistry(series.size() - w_size);
int inner_idx = innerRegistry.getNextUnvisited();
while(!(-1==inner_idx)){
innerRegistry.markVisited(inner_idx);
// Rcout << "examining the subsequences starting at outer " << outer_idx << " and inner " << inner_idx << "\n";
if(std::abs(inner_idx - outer_idx) > w_size) {
NumericVector curr_seq = subseries(series, inner_idx, inner_idx + w_size);
double dist = early_abandoned_dist(candidate_seq, curr_seq, nnDistance);
// Rcout << " .. dist: " << dist << " best dist " << nnDistance << "\n";
if ( (!std::isnan(dist)) && dist < nnDistance) {
nnDistance = dist;
}
}
inner_idx = innerRegistry.getNextUnvisited();
}
if (!(std::numeric_limits<double>::max() == nnDistance)
&& nnDistance > best_so_far_distance) {
// Rcout << "** updating discord " << nnDistance << " at " << outer_idx << "\n";
best_so_far_distance = nnDistance;
best_so_far_index = outer_idx;
}
outer_idx = outerRegistry.getNextUnvisited();
}
struct discord_record res;
res.index = best_so_far_index;
res.nn_distance = best_so_far_distance;
return res;
}
//' Finds a discord using brute force algorithm.
//'
//' @param ts the input timeseries.
//' @param w_size the sliding window size.
//' @param discords_num the number of discords to report.
//' @useDynLib jmotif
//' @export
//' @references Keogh, E., Lin, J., Fu, A.,
//' HOT SAX: Efficiently finding the most unusual time series subsequence.
//' Proceeding ICDM '05 Proceedings of the Fifth IEEE International Conference on Data Mining
//' @examples
//' discords = find_discords_brute_force(ecg0606[1:600], 100, 1)
//' plot(ecg0606[1:600], type = "l", col = "cornflowerblue", main = "ECG 0606")
//' lines(x=c(discords[1,2]:(discords[1,2]+100)),
//' y=ecg0606[discords[1,2]:(discords[1,2]+100)], col="red")
// [[Rcpp::export]]
Rcpp::DataFrame find_discords_brute_force(
NumericVector ts, int w_size, int discords_num) {
std::map<int, double> res;
VisitRegistry registry(ts.length());
registry.markVisited(ts.length() - w_size, ts.length());
// Rcout << "starting search of " << discords_num << " discords..." << "\n";
int discord_counter = 0;
while(discord_counter < discords_num){
discord_record rec = find_best_discord_brute_force(ts, w_size, ®istry);
// Rcout << "found a discord " << discord_counter << " at " << rec.index;
// Rcout << ", NN distance: " << rec.nn_distance << "\n";
if(rec.nn_distance == 0 || rec.index == -1){ break; }
res.insert(std::make_pair(rec.index, rec.nn_distance));
int start = rec.index - w_size;
if(start<0){
start = 0;
}
int end = rec.index + w_size;
// it can't be greater
// if(end>=ts.length()){
// end = ts.length();
//}
// Rcout << "marking as visited from " << start << " to " << end << "\n";
registry.markVisited(start, end);
discord_counter = discord_counter + 1;
}
std::vector<int> positions;
std::vector<double > distances;
for(std::map<int, double>::iterator it = res.begin(); it != res.end(); it++) {
positions.push_back(it->first);
distances.push_back(it->second);
}
// make results
return Rcpp::DataFrame::create(
Named("nn_distance") = distances,
Named("position") = positions
);
}