Skip to content

Commit

Permalink
add best_path_isclose to c source
Browse files Browse the repository at this point in the history
  • Loading branch information
@wannesm
wannesm committed Jun 14, 2023
1 parent 922027a commit 15cac70
Show file tree
Hide file tree
Showing 5 changed files with 195 additions and 22 deletions.
24 changes: 19 additions & 5 deletions dtaidistance/lib/DTAIDistanceC/DTAIDistanceC/dd_benchmark.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,17 +109,31 @@ void benchmark3() {
void benchmark4() {
// double s1[] = {0, 0, 1, 2, 1, 0, 1, 0, 0}; int l1 = 9;
// double s2[] = {0, 1, 2, 0, 0, 0, 0, 0, 0}; int l2 = 9;
double s1[] = {0., 0., 1., 2., 1., 0., 1., 0., 0., 2., 1., 0., 0.}; int l1 = 13;
double s2[] = {0., 1., 2., 3., 1., 0., 0., 0., 2., 1., 0., 0., 0.}; int l2 = 13;
// double s1[] = {0., 0., 1., 2., 1., 0., 1., 0., 0., 2., 1., 0., 0.}; int l1 = 13;
// double s2[] = {0., 1., 2., 3., 1., 0., 0., 0., 2., 1., 0., 0., 0.}; int l2 = 13;
double s1[] = {2.1, 4.1, 5.1}; int l1 = 3;
double s2[] = {1.1, 2.1, 3.1, 4.1, 5.1}; int l2 = 5;
DTWSettings settings = dtw_settings_default();
settings.use_pruning = true;
settings.inner_dist = 0;
dtw_settings_set_psi(2, &settings);
settings.psi_2b = l2;
settings.psi_2e = l2;
// dtw_settings_set_psi(2, &settings);
// double d = dtw_distance(s1, 9, s2, 9, &settings);
idx_t wps_length = dtw_settings_wps_length(l1, l2, &settings);
seq_t wps[wps_length];
double d = dtw_warping_paths(wps, s1, l1, s2, l2, true, true, false, &settings);
printf("d=%f\n", d);

idx_t i1[l1+l2];
idx_t i2[l1+l2];
for (idx_t i=0; i<(l1+l2); i++) {i1[i]=0; i2[i]=0;}
dtw_best_path_isclose(wps, i1, i2, l1, l2, 1e-05, 1e-08, &settings);
printf("[");
for (idx_t i=0; i<(l1+l2); i++) {
printf("(%zu,%zu)", i1[i], i2[i]);
}
printf("]\n");
}

void benchmark5() {
Expand Down Expand Up @@ -593,10 +607,10 @@ int main(int argc, const char * argv[]) {
time(&start_t);
clock_gettime(CLOCK_REALTIME, &start);

benchmark1();
// benchmark1();
// benchmark2();
// benchmark3();
// benchmark4();
benchmark4();
// benchmark5();
// benchmark6();
// benchmark7();
Expand Down
130 changes: 129 additions & 1 deletion dtaidistance/lib/DTAIDistanceC/DTAIDistanceC/dd_dtw.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ DTWSettings dtw_settings_default(void) {
.psi_2e = 0,
.use_pruning = false,
.only_ub = false,
.inner_dist = 0
.inner_dist = 0,
.window_type = 0
};
return s;
}
Expand Down Expand Up @@ -61,6 +62,7 @@ void dtw_settings_print(DTWSettings *settings) {
printf(" use_pruning = %d\n", settings->use_pruning);
printf(" only_ub = %d\n", settings->only_ub);
printf(" inner_dist = %d\n", settings->inner_dist);
printf(" window_type = %d\n", settings->window_type);
printf("}\n");
}

Expand Down Expand Up @@ -3206,6 +3208,7 @@ idx_t dtw_wps_max(DTWWps* p, seq_t *wps, idx_t *r, idx_t *c, idx_t l1, idx_t l2)
}



/*!
Compute best path between two series.
Expand All @@ -3219,6 +3222,7 @@ Compute best path between two series.
@param settings for Dynamic Time Warping.
@return length of path
*/

idx_t dtw_best_path(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,

DTWSettings *settings) {
Expand Down Expand Up @@ -3326,6 +3330,129 @@ idx_t dtw_best_path(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,
}


/*!
Compute best path between two series.
@param wps Array of length `(l1+1)*min(l2+1, abs(l1-l2) + 2*window-1)` with the warping paths.
@param i1 Array of length l1+l2 to store the indices for the first sequence.
Reverse ordered, last one is if i1 or i2 is zero.
@param i2 Array of length l1+l2 to store the indices for the second sequence.
Reverse ordered, last one is if i1 or i2 is zero.
@param l1 Length of first array.
@param l2 Length of second array.
@param rtol Relative tolerance for isclose, typical value is 1e-05
@param atol Absolute tolerance for isclose, typical value is 1e-08
@param settings for Dynamic Time Warping.
@return length of path
*/

idx_t dtw_best_path_isclose(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,
seq_t rtol, seq_t atol,
DTWSettings *settings) {
DTWWps p = dtw_wps_parts(l1, l2, settings);

idx_t i = 0;
idx_t rip = l1;
idx_t cip = l2;
idx_t min_ci;
idx_t wpsi_start, wpsi;
idx_t ri_widthp = p.width * (rip - 1);
idx_t ri_width = p.width * rip;

// D. ri3 <= ri < l1
min_ci = p.ri3 + 1 - p.window - p.ldiff;
wpsi_start = 2;
if (p.ri2 == p.ri3) {
wpsi_start = min_ci + 1;
} else {
min_ci = 1 + p.ri3 - p.ri2;
}
wpsi = wpsi_start + (l2 - min_ci) - 1;
while (rip > p.ri3 && cip > 0) {
if (wps[ri_width + wpsi] != -1) {
i1[i] = rip - 1;
i2[i] = cip - 1;
i++;
}
if ((wps[ri_widthp + wpsi - 1] <= wps[ri_width + wpsi - 1] || fabs(wps[ri_widthp + wpsi - 1] - wps[ri_width + wpsi - 1]) <= (atol + rtol * fabs(wps[ri_width + wpsi - 1]))) &&
(wps[ri_widthp + wpsi - 1] <= wps[ri_widthp + wpsi] || fabs(wps[ri_widthp + wpsi - 1] - wps[ri_widthp + wpsi]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi])))) {
// Go diagonal
cip--;
rip--;
wpsi = wpsi - 1;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else if ((wps[ri_width + wpsi - 1] <= wps[ri_widthp + wpsi] || fabs(wps[ri_width + wpsi - 1] - wps[ri_widthp + wpsi]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi])))) {
// Go left
cip--;
wpsi--;
} else {
// Go up
rip--;
ri_width = ri_widthp;
ri_widthp -= p.width;
}
}

// C. ri2 <= ri < ri3
while (rip > p.ri2 && cip > 0) {
if (wps[ri_width + wpsi] != -1) {
i1[i] = rip - 1;
i2[i] = cip - 1;
i++;
}
if ((wps[ri_widthp + wpsi] <= wps[ri_width + wpsi - 1] || fabs(wps[ri_widthp + wpsi] - wps[ri_width + wpsi - 1]) <= (atol + rtol * fabs(wps[ri_width + wpsi - 1]))) &&
(wps[ri_widthp + wpsi] <= wps[ri_widthp + wpsi + 1] || fabs(wps[ri_widthp + wpsi] - wps[ri_widthp + wpsi + 1]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi + 1])))) {
// Go diagonal
cip--;
rip--;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else if ((wps[ri_width + wpsi - 1] <= wps[ri_widthp + wpsi + 1] || fabs(wps[ri_width + wpsi - 1] - wps[ri_widthp + wpsi + 1]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi + 1])))) {
// Go left
cip--;
wpsi--;
} else {
// Go up
rip--;
wpsi++;
ri_width = ri_widthp;
ri_widthp -= p.width;
}
}

// A-B. 0 <= ri < ri2
while (rip > 0 && cip > 0) {
if (wps[ri_width + wpsi] != -1) {
i1[i] = rip - 1;
i2[i] = cip - 1;
i++;
}
if ((wps[ri_widthp + wpsi - 1] <= wps[ri_width + wpsi - 1] || fabs(wps[ri_widthp + wpsi - 1] - wps[ri_width + wpsi - 1]) <= (atol + rtol * fabs(wps[ri_width + wpsi - 1]))) &&
(wps[ri_widthp + wpsi - 1] <= wps[ri_widthp + wpsi] || fabs(wps[ri_widthp + wpsi - 1] - wps[ri_widthp + wpsi]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi])))) {
// Go diagonal
cip--;
rip--;
wpsi--;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else {
if ((wps[ri_width + wpsi - 1] <= wps[ri_widthp + wpsi] || fabs(wps[ri_width + wpsi - 1] - wps[ri_widthp + wpsi]) <= (atol + rtol * fabs(wps[ri_widthp + wpsi])))) {
// Go left
cip--;
wpsi--;
} else {
// Go up
rip--;
ri_width = ri_widthp;
ri_widthp -= p.width;
}
}
}
return i;
}


/*!
Compute best path in affinity matrix for two series.
Expand All @@ -3341,6 +3468,7 @@ Compute best path in affinity matrix for two series.
@param settings for Dynamic Time Warping.
@return length of path
*/

idx_t dtw_best_path_affinity(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,
idx_t rs, idx_t cs,
DTWSettings *settings) {
Expand Down
10 changes: 8 additions & 2 deletions dtaidistance/lib/DTAIDistanceC/DTAIDistanceC/dd_dtw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,12 @@ static char printFormat[5];


// Inner distance options
//const int kSquaredEuclidean = 0;
//const int kEuclidean = 1;
//const int kSquaredEuclideanInnerDist = 0;
//const int kEuclideanInnerDist = 1;

// Band type
//const int kSakoeChibaBand = 0;
//const int kSlantedBand = 1;


/**
Expand Down Expand Up @@ -76,6 +80,7 @@ struct DTWSettings_s {
bool use_pruning;
bool only_ub;
int inner_dist; // 0=squared euclidean, 1=euclidean
int window_type; // 0=band around two diagonals, 1=band around slanted diagonal

};
typedef struct DTWSettings_s DTWSettings;
Expand Down Expand Up @@ -153,6 +158,7 @@ idx_t dtw_wps_loc(DTWWps* p, idx_t r, idx_t c, idx_t l1, idx_t l2);
idx_t dtw_wps_loc_columns(DTWWps* p, idx_t r, idx_t *cb, idx_t *ce, idx_t l1, idx_t l2);
idx_t dtw_wps_max(DTWWps* p, seq_t *wps, idx_t *r, idx_t *c, idx_t l1, idx_t l2);
idx_t dtw_best_path(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2, DTWSettings *settings);
idx_t dtw_best_path_isclose(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2, seq_t rtol, seq_t atol, DTWSettings *settings);
idx_t dtw_best_path_affinity(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2, idx_t s1, idx_t s2, DTWSettings *settings);
idx_t dtw_best_path_prob(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2, seq_t avg, DTWSettings *settings);
seq_t dtw_warping_path(seq_t *from_s, idx_t from_l, seq_t* to_s, idx_t to_l, idx_t *from_i, idx_t *to_i, idx_t * length_i, DTWSettings * settings);
Expand Down
9 changes: 8 additions & 1 deletion dtaidistance/lib/DTAIDistanceC/DTAIDistanceC/jinja/dd_dtw.jinja.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ DTWSettings dtw_settings_default(void) {
.psi_2e = 0,
.use_pruning = false,
.only_ub = false,
.inner_dist = 0
.inner_dist = 0,
.window_type = 0
};
return s;
}
Expand Down Expand Up @@ -61,6 +62,7 @@ void dtw_settings_print(DTWSettings *settings) {
printf(" use_pruning = %d\n", settings->use_pruning);
printf(" only_ub = %d\n", settings->only_ub);
printf(" inner_dist = %d\n", settings->inner_dist);
printf(" window_type = %d\n", settings->window_type);
printf("}\n");
}

Expand Down Expand Up @@ -605,10 +607,15 @@ idx_t dtw_wps_max(DTWWps* p, seq_t *wps, idx_t *r, idx_t *c, idx_t l1, idx_t l2)


{% set suffix = '' %}
{% set use_isclose = 0 %}
{%- include 'dtw_bestpath.jinja.c' %}

{% set suffix = '' %}
{% set use_isclose = 1 %}
{%- include 'dtw_bestpath.jinja.c' %}

{% set suffix = '_affinity' %}
{% set use_isclose = 0 %}
{%- include 'dtw_bestpath.jinja.c' %}


Expand Down
44 changes: 31 additions & 13 deletions dtaidistance/lib/DTAIDistanceC/DTAIDistanceC/jinja/dtw_bestpath.jinja.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,17 +16,35 @@ Compute best path between two series.
@param rs Start position row.
@param cs Start position column.
{%- endif %}
{%- if use_isclose == 1 %}
@param rtol Relative tolerance for isclose, typical value is 1e-05
@param atol Absolute tolerance for isclose, typical value is 1e-08
{%- endif %}
@param settings for Dynamic Time Warping.
@return length of path
*/
{%- if "affinity" in suffix %}
{%- set cmp=">=" %}
{% macro cmpfn(a, b) -%}
{%- if use_isclose == 1 -%}
({{a}} <= {{b}} || fabs({{a}} - {{b}}) <= (atol + rtol * fabs({{b}})))
{%- else %}
{%- if "affinity" in suffix -%}
{{a}} >= {{b}}
{%- else -%}
{{a}} <= {{b}}
{%- endif %}
{%- endif %}
{%- endmacro -%}
{%- if use_isclose == 1 %}
{%- set suffix2 = "_isclose" %}
{%- else %}
{%- set cmp="<=" %}
{%- set suffix2 = "" %}
{%- endif %}
idx_t dtw_best_path{{suffix}}(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,
idx_t dtw_best_path{{suffix}}{{suffix2}}(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t l2,
{% if "affinity" in suffix -%}
idx_t rs, idx_t cs,
{%- endif -%}
{%- if use_isclose == 1 -%}
seq_t rtol, seq_t atol,
{%- endif %}
DTWSettings *settings) {
DTWWps p = dtw_wps_parts(l1, l2, settings);
Expand Down Expand Up @@ -75,15 +93,15 @@ idx_t dtw_best_path{{suffix}}(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t
i++;
}
{%- endif %}
if (wps[ri_widthp + wpsi - 1] {{cmp}} wps[ri_width + wpsi - 1] &&
wps[ri_widthp + wpsi - 1] {{cmp}} wps[ri_widthp + wpsi]) {
if ({{cmpfn("wps[ri_widthp + wpsi - 1]", "wps[ri_width + wpsi - 1]")}} &&
{{cmpfn("wps[ri_widthp + wpsi - 1]", "wps[ri_widthp + wpsi]")}}) {
// Go diagonal
cip--;
rip--;
wpsi = wpsi - 1;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else if (wps[ri_width + wpsi - 1] {{cmp}} wps[ri_widthp + wpsi]) {
} else if ({{cmpfn("wps[ri_width + wpsi - 1]","wps[ri_widthp + wpsi]")}}) {
// Go left
cip--;
wpsi--;
Expand Down Expand Up @@ -113,14 +131,14 @@ idx_t dtw_best_path{{suffix}}(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t
i++;
}
{%- endif %}
if (wps[ri_widthp + wpsi] {{cmp}} wps[ri_width + wpsi - 1] &&
wps[ri_widthp + wpsi] {{cmp}} wps[ri_widthp + wpsi + 1]) {
if ({{cmpfn("wps[ri_widthp + wpsi]","wps[ri_width + wpsi - 1]")}} &&
{{cmpfn("wps[ri_widthp + wpsi]","wps[ri_widthp + wpsi + 1]")}}) {
// Go diagonal
cip--;
rip--;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else if (wps[ri_width + wpsi - 1] {{cmp}} wps[ri_widthp + wpsi + 1]) {
} else if ({{cmpfn("wps[ri_width + wpsi - 1]","wps[ri_widthp + wpsi + 1]")}}) {
// Go left
cip--;
wpsi--;
Expand Down Expand Up @@ -151,16 +169,16 @@ idx_t dtw_best_path{{suffix}}(seq_t *wps, idx_t *i1, idx_t *i2, idx_t l1, idx_t
i++;
}
{%- endif %}
if (wps[ri_widthp + wpsi - 1] {{cmp}} wps[ri_width + wpsi - 1] &&
wps[ri_widthp + wpsi - 1] {{cmp}} wps[ri_widthp + wpsi]) {
if ({{cmpfn("wps[ri_widthp + wpsi - 1]","wps[ri_width + wpsi - 1]")}} &&
{{cmpfn("wps[ri_widthp + wpsi - 1]","wps[ri_widthp + wpsi]")}}) {
// Go diagonal
cip--;
rip--;
wpsi--;
ri_width = ri_widthp;
ri_widthp -= p.width;
} else {
if (wps[ri_width + wpsi - 1] {{cmp}} wps[ri_widthp + wpsi]) {
if ({{cmpfn("wps[ri_width + wpsi - 1]","wps[ri_widthp + wpsi]")}}) {
// Go left
cip--;
wpsi--;
Expand Down

0 comments on commit 15cac70

Please sign in to comment.