Newer
Older
#include "dtw.h"
#include "calcul.h"
#include "print.h"
#include "help.h"
///sumary: entry point function for dtw method
///parameters:
///input_method - input for the dtw method
///input_info - info about input data
///parameter - parameter for dtw method
result_dtw dtw::main(input_method const &input, input_info const &info, parameter const ¶ms)
result_dtw result;
if (params.segmented)
result = main_segment(input, info, params);
return main_pair(input, info, params);
///sumary: entry point function for dtw method
///parameters:
///input_method - input for the dtw method
///input_info - info about input data
///parameter - parameter for dtw method
result_dtw dtw::main_pair(input_method const &input, input_info const &info, parameter const ¶ms)
if ((int)((input.A.size() * input.B.size()) / 131072) > params.ram) //131072 to convert bytes to MB
{
cout << "size A: " << input.A.size() << ", size B: " << input.B.size() << endl;
//throw runtime_error("DTW aborted. Input too large: " + to_string(input.A.size() * input.B.size()) + "B");
cout << "DTW aborted. Input too large: " << (input.A.size() * input.B.size() * 8) / 1024 / 1024 << "MB" << endl;
cout << "For overriding RAM limit use -ram [GBs] switch." << endl;
exit(0);
}
auto result = configure(input, info, params);
if (params.drawOut.size() > 0)
draw::plot_pair(result, input, info, params);
if (params.isRatioReversed())
{
for (size_t i = 0; i < result.score.size(); i++)
result.score[params.scoreType] = 1 - result.score[params.scoreType];
///sumary: entry point function for dtw method
///parameters:
///input_method - input for the dtw method
///input_info - info about input data
///parameter - parameter for dtw method
result_dtw dtw::main_segment(input_method const &input, input_info const &info, parameter const ¶ms)
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
{
//segmentig
int win = params.segmented + 1;
tree trieA1, trieA2, trieB1, trieB2;
trieA1 = tree(input.A, win);
trieB1 = tree(input.B, win);
vtr2<double> revA(input.A);
reverse(revA.begin(), revA.end());
trieA2 = tree(revA, win);
vtr2<double> revB(input.B);
reverse(revB.begin(), revB.end());
trieB2 = tree(revB, win);
vtr3<double> setA, setB;
//setA = trieA.process(input.A, params.segmented);
//setB = trieA.process(input.B, params.segmented);
setA = trieA1.getSegments(input.A, trieA1.getVotes(input.A, win - 1), trieA2.getVotes(revA, win- 1), win);
setB = trieA1.getSegments(input.B, trieB1.getVotes(input.B, win - 1), trieB2.getVotes(revB, win- 1), win);
vtr2<result_dtw> result(setA.size());
for (size_t i = 0; i < setA.size(); i++)
{
result[i] = vtr<result_dtw>(setB.size());
for (size_t j = 0; j < setB.size(); j++)
{
input_method subInput(setA[i], setB[j]);
result[i][j] = configure(subInput, info, params);
}
}
if (params.drawOut.size() > 0) {
draw::plot_segmets(input, result, info, params);
}
if (params.isRatioReversed())
{
for (size_t i = 0; i < result.size(); i++)
for (size_t j = 0; j < result[i].size(); j++)
result[i][j].score[j] -= 1;
return result[0][0];
/*
sumary: entry point function for dtw method
parameters:
input_method - input for the dtw method
input_info - info about input data
parameter - parameter for dtw method
*/
result_dtw dtw::configure(input_method const &input, input_info const &info, parameter const ¶ms)
DISTANCE d;
d.classic = calcul::distance_dtw_euklid;
d.classic = calcul::distance_dtw_manhattan;
d.csi_chroma = calcul::distance_dtw_csiChroma;
else if (params.distance == 4)
d.csi_chord = calcul::distance_dtw_csiChord;
else if (params.distance == 5)
d.classic = calcul::distance_dtw_euklid_mean;
if (params.isMemoization())
back = matrix_memoized(input.A, input.B, params);
else if(params.block > 0)
back = dtw::matrix_tiled(input.A, input.B, params);
else if(params.isSimd())
back = dtw::matrix_simd(input.A, input.B, params);
else
{
if (params.localAlignment)
result = dtw::alignment_local(input, info, d, params);
result = dtw::alignment(input, info, d, params);
}
///sumary: entry point function for dtw method
///parameters:
///input_method - input for the dtw method
///input_info - info about input data
///DISTANCE - pointer to function used for distance calculations
///parameter - parameter for dtw method
result_dtw dtw::alignment(input_method const &input, input_info const &info, DISTANCE distance, parameter const ¶ms)
auto m = dtw::matrix(input, distance, params);
auto end = get_relaxedEnds(m, params);
//auto warping = get_warping(m, end, params);
vtr<result_path> warping(1);
warping[0] = get_warping(m, end, params);
warping[0].scoreRaw = m[end.row][end.col].value;
warping[0].end = end;
vtr2<range> ranges(1);
ranges[0].push_back(range(0, (int)warping[0].path.size()));
result_dtw result;
if (params.drawOut.size() > 0) {
result.matrix_acc = draw::draw_combine(m, warping, params);
result.matrix_noacc = draw::draw_combine(dtw::matrix_noaccumulation(input, distance, params), warping, params);
cout << endl << print::distanceMatrix(m);
cout << endl << warping[0].path;
//cout << endl << print::printPathShape(back.path, end, (int)A.size() + 1, (int)B.size() + 1);
//print::write(print::printHtmlDistanceMatrix<T>(m), "c:\\code\\data\\sc\\dm.html", false);
}
//vtr<double> resultScore;
//resultScore.push_back(warping[0].scoreRaw);
//resultScore.push_back(calcul::score_dtw_s2(warping[0].scoreRaw, warping[0].path.size())); // back.path.size());
//resultScore.push_back(calcul::score_dtw_s3(warping[0].end.row - warping[0].start.row, warping[0].end.col - warping[0].start.col, warping[0].path.size()));
//resultScore.push_back(calcul::score_dtw_s4(warping[0].scoreRaw, calcul::score_dtw_max(input.A, input.B, warping[0].start, warping[0].end)));
//resultScore.push_back(calcul::score_dtw_s5(warping[0].scoreRaw,
// calcul::score_dtw_max(input.A, input.B, warping[0].start, warping[0].end),
// calcul::lenRatio(warping[0].end.row - warping[0].start.row, warping[0].end.col - warping[0].start.row)));
result.score = score(input, warping);
result.warpings = warping;
return result;
}
vtr<double> dtw::score(input_method const &input, vtr<result_path> const &warpings)
{
vtr<double> score(5);
for (auto &i : warpings) {
score[0] += i.scoreRaw;
score[1] += calcul::score_dtw_s2(i.scoreRaw, i.path.size()); // back.path.size());
score[2] += calcul::score_dtw_s3(i.end.row - i.start.row, i.end.col - i.start.col, i.path.size());
score[3] += calcul::score_dtw_s4(i.scoreRaw, calcul::score_dtw_max(input.A, input.B, i.start, i.end));
score[4] += calcul::score_dtw_s5(i.scoreRaw,
calcul::score_dtw_max(input.A, input.B, i.start, i.end),
calcul::lenRatio(i.end.row - i.start.row, i.end.col - i.start.row));
}
for (auto &i : score)
i /= warpings.size();
return score;
}
result_dtw dtw::alignment_local(input_method const &input, input_info const &info, DISTANCE distance, parameter const ¶ms)
{
auto m = dtw::matrix_noaccumulation(input, distance, params);
auto mMin = dtw::get_minimums(m, params);
for (auto &i : mMin)
m[i.row][i.col] = 0;
dtw::accumulate_mod(m, mMin, params);
result_dtw result;
result.warpings = get_warpings(m, mMin, params);
filterPaths(m, result.warpings, params);
if (params.drawOut.size() > 0) {
result.matrix_acc = draw::draw_combine(m, result.warpings/*, ranges*/, params);
result.matrix_noacc = draw::draw_combine(dtw::matrix_noaccumulation(input, distance, params), result.warpings, params);
}
if (params.isDebugInfo()) {
cout << endl << print::distanceMatrix(m);
//cout << endl << warpings[0].path;
//cout << endl << print::printPathShape(back.path, end, (int)A.size() + 1, (int)B.size() + 1);
//print::write(print::printHtmlDistanceMatrix<T>(m), "c:\\code\\data\\sc\\dm.html", false);
}
//result.score = vtr<double>(5);
//for (auto &i : result.warpings) {
// result.score[0] += i.scoreRaw;
// result.score[1] += calcul::score_dtw_s2(i.scoreRaw, i.path.size()); // back.path.size());
// result.score[2] += calcul::score_dtw_s3(i.end.row - i.start.row, i.end.col - i.start.col, i.path.size());
// result.score[3] += calcul::score_dtw_s4(i.scoreRaw, calcul::score_dtw_max(input.A, input.B, i.start, i.end));
// result.score[4] += calcul::score_dtw_s5(i.scoreRaw, calcul::score_dtw_max(input.A, input.B, i.start, i.end),
// calcul::lenRatio(i.end.row - i.start.row, i.end.col - i.start.row));
//}
//for (auto &i : result.score)
// i /= result.warpings.size();
result.score = score(input, result.warpings);
///sumary: entry point function for dtw method
///parameters:
///input_method - input for the dtw method
///DISTANCE - pointer to function used for distance calculations
///parameter - parameter for dtw method
vtr2<node> dtw::matrix(input_method const &input, DISTANCE distance, parameter const ¶ms)
int lenA = (int)input.A.size();
int lenB = (int)input.B.size();
vtr2<node> m(lenA + 1);
for (int i = 0; i < lenA + 1; i++)
m[i] = vtr<node>(lenB + 1);
for (int i = 0; i < min(lenA, params.relax + 1); i++) //include 0,0 = 0 !!
for (int i = 0; i < min(lenB, params.relax + 1); i++)
if (params.isSubsequence() && calcul::lenRatio(lenA, lenB) < params.subsequence)
if(lenA < lenB)
for (int i = 0; i < lenB + 1; i++)
if (lenB < lenA)
for (int i = 0; i < lenA + 1; i++)
const double coef = lenB / static_cast<double>(lenA);
const int w = (int)(lenB * params.w);
for (int i = 1; i < lenA + 1; i++) //row - y
const size_t start = calcul::dtw_wpassStart(i, w, coef);
const size_t end = calcul::dtw_wpassEnd(i, w, coef, lenB);
for (size_t j = start; j < end; j++) //col - x
{
const double u = m[i - 1][j].value;
const double l = m[i][j - 1].value;
const double d = m[i - 1][j - 1].value;
m[i][j].value = distance.csi_chroma(input.A[i - 1], input.B[j - 1], 0.07) + min;
else if (params.distance == 4)
m[i][j].value = distance.csi_chord(input.A[i - 1], input.B[j - 1], input.A2[i - 1], input.B2[j - 1]) + min;
m[i][j].value = distance.classic(input.A[i - 1], input.B[j - 1]) + min;
}
}
return m;
}
///sumary:
///parameters:
vtr2<node> dtw::matrix_noaccumulation(input_method const &input, DISTANCE distance, parameter const ¶ms)
{
int lenA = (int)input.A.size();
int lenB = (int)input.B.size();
vtr2<node> m(lenA + 1);
for (int i = 0; i < lenA + 1; i++)
m[i] = vtr<node>(lenB + 1);
for (int i = 0; i < min(lenA, params.relax + 1); i++)
for (int i = 0; i < min(lenB, params.relax + 1); i++)
m[0][i].value = 0;
if (params.isSubsequence() && calcul::lenRatio(lenA, lenB) < params.subsequence)
{
if (lenA < lenB)
for (int i = 0; i < lenB + 1; i++)
m[0][i].value = 0;
if (lenB < lenA)
for (int i = 0; i < lenA + 1; i++)
m[i][0].value = 0;
}
const double coef = lenB / (double)lenA;
const int w = (int)(lenB * params.w);
for (int i = 1; i < lenA + 1; i++) //row = y
{
const size_t start = calcul::dtw_wpassStart(i, w, coef);
const size_t end = calcul::dtw_wpassEnd(i, w, coef, lenB);
for (size_t j = start; j < end; j++) //col = x
{
if (params.distance == 3)
m[i][j].value = distance.csi_chroma(input.A[i - 1], input.B[j - 1], 0.07);
else if (params.distance == 4)
m[i][j].value = distance.csi_chord(input.A[i - 1], input.B[j - 1], input.A2[i - 1], input.B2[j - 1]);
m[i][j].value = distance.classic(input.A[i - 1], input.B[j - 1]);
void dtw::accumulate(vtr2<node> &m, parameter const ¶ms)
{
int lenA = (int)m.size();
int lenB = (int)m[0].size();
const double coef = lenB / (double)lenA;
const int w = (int)(lenB * params.w);
for (int i = 1; i < lenA; i++) //row = y
{
const size_t start = calcul::dtw_wpassStart(i, w, coef);
const size_t end = calcul::dtw_wpassEnd(i, w, coef, lenB);
for (size_t j = start; j < end; j++) //col = x
{
m[i][j].value += std::min({ m[i - 1][j - 1].value, m[i - 1][j].value, m[i][j - 1].value });
}
}
}
void dtw::accumulate_mod(vtr2<node> &m, vtr<coord> const &minims, parameter const ¶ms)
{
int lenA = (int)m.size();
int lenB = (int)m[0].size();
const double coef = lenB / (double)lenA;
const int w = (int)(lenB * params.w);
for (int i = 1; i < lenA; i++) //row = y
{
const size_t start = calcul::dtw_wpassStart(i, w, coef);
const size_t end = calcul::dtw_wpassEnd(i, w, coef, lenB);
if(m[i][j].value != 0)
m[i][j].value += std::min({ m[i - 1][j - 1].value, m[i - 1][j].value, m[i][j - 1].value });
vtr<coord> dtw::get_minimums(vtr2<node> const &m, parameter const ¶ms)
vtr<coord> minims;
int lenA = (int)m.size();
int lenB = (int)m[0].size();
const double coef = lenB / (double)lenA;
const int w = (int)(lenB * params.w);
for (int i = 1; i < lenA; i++) //row = y
{
const int start = calcul::dtw_wpassStart(i, w, coef);
const int end = calcul::dtw_wpassEnd(i, w, coef, lenB);
for (int j = start; j < end; j++) //col = x
{
double current = m[i][j].value;
if (i + 1 < static_cast<int>(m.size()) && j + 1 < static_cast<int>(m[i].size()) &&
m[i - 1][j].value > current && m[i - 1][j + 1].value > current && m[i][j - 1].value > current &&
m[i][j + 1].value > current && m[i + 1][j - 1].value > current && m[i + 1][j].value > current)
{
minims.push_back(coord(i, j));
vtr<coord> filtered;
for (size_t i = 0; i < minims.size(); i++)
{
for (size_t j = i + 1; j < minims.size(); j++)
{
if (minims[j].row - minims[i].row < 2) // if dist in rows is larger than 1 ...cant find neigh
{
if (minims[j].col - minims[i].col == 1)
break;
else
continue;
}
filtered.push_back(minims[i]);
break;
}
}
return filtered;
///sumary:
///parameters:
///vtr2<node<T>> -
///i
///j
///params
result_path dtw::get_warping(vtr2<node> const &m, coord coords, parameter const ¶ms)
result_path warping;
double lenA = coords.row;
double lenB = coords.col;
const bool flexible = params.isPassFlexible();
while (coords.row > 0 && coords.col > 0)
warping.pathCoords.push_back(coord(coords.row, coords.col));
double d = 0, u = 0, l = 0;
if (flexible) {
coord past = (int)warping.pathCoords.size() - params.d >= 0 ? warping.pathCoords[warping.pathCoords.size() - params.d] : coord(coords.row, coords.col);
d = calcul::dtw_isFlexiblePass(coords.row - 1, coords.col - 1, past.row, past.col, params.wd, params.d) ? m[coords.row - 1][coords.col - 1].value : constant::MAX_double;
u = calcul::dtw_isFlexiblePass(coords.row - 1, coords.col, past.row, past.col, params.wd, params.d) ? m[coords.row - 1][coords.col].value : constant::MAX_double;
l = calcul::dtw_isFlexiblePass(coords.row, coords.col - 1, past.row, past.col, params.wd, params.d) ? m[coords.row][coords.col - 1].value : constant::MAX_double;
}
else{
d = m[coords.row - 1][coords.col - 1].value;
u = m[coords.row - 1][coords.col].value;
l = m[coords.row][coords.col - 1].value;
}
if (min({ d, u, l }) == d){
warping.path = "M" + warping.path;
coords.row--;
coords.col--;
if (l < u){
warping.path = "L" + warping.path;
coords.col--;
else if (u < l){
warping.path = "U" + warping.path;
coords.row--;
if(lenA / coords.row > lenB / coords.col)
warping.path = "L" + warping.path;
coords.col--;
else{
warping.path = "U" + warping.path;
coords.row--;
//warping.pathCoords.push_back(coord(coords.row, coords.col));
if (!params.isSubsequence() || (params.isSubsequence() && m.size() < m[0].size()))
while (coords.row > params.relax)
warping.path = "U" + warping.path;
coords.row--;
warping.pathCoords.push_back(coord(coords.row, coords.col));
if (!params.isSubsequence() || (params.isSubsequence() && m[0].size() < m.size()))
while (coords.col > params.relax)
warping.path = "L" + warping.path;
coords.col--;
warping.pathCoords.push_back(coord(coords.row, coords.col));
warping.start.row = coords.row; // start is upper left...so its where path is finished building(sort of end)
warping.start.col = coords.col;
std::reverse(warping.pathCoords.begin(), warping.pathCoords.end());
return warping;
vtr<result_path> dtw::get_warpings(vtr2<node> const &m, vtr<coord> const &minims, parameter const ¶ms)
{
vtr<result_path> outPaths;
for (auto i : minims) {//no ref!
outPaths.push_back(get_warping(m, i, params));
outPaths.back().end = i;
outPaths.back().scoreRaw = m[i.row][i.col].value;
}
return outPaths;
}
void dtw::filterPaths(vtr2<node> const &m, vtr<result_path> &warpings, parameter const ¶ms)
vtr2<range> ranges(warpings.size());
for (int k = 0; k < (int)warpings.size(); k++)
for (int i = 0; i < (int)warpings[k].path.size(); i++)
{
double cTotal = 0;
double max = constant::MIN_double;
range rangeTmp = range(-1, 0);
bool insert = true;
bool accept = false;
if ((int)warpings[k].path.size() >= params.treshold_l) {
for (int j = i; j < (int)warpings[k].path.size(); j++)
int len = j - i + 1;
double current = m[warpings[k].pathCoords[j].row][warpings[k].pathCoords[j].col].value - cTotal;
cTotal += current;
if (max < current)
max = current;
if (max <= params.treshold_e /*&& cTotal <= params.treshold_t*/ && cTotal / len <= params.treshold_a && len >= params.treshold_l)
{
accept = true;
if (!insert)
if (i - rangeTmp.end < 2)
rangeTmp.end = j;
if (insert) {
rangeTmp = range(i, j);
insert = false;
}
}
}
if (accept) //filtr for excluding subsequences // true subsequences
{
for (auto &i : ranges[k])
if (i.start <= rangeTmp.start && rangeTmp.end <= i.end)
accept = false;
else if (rangeTmp.start - i.end < 2) {
i.end = rangeTmp.end;
accept = false;
break;
}
if (accept)
ranges[k].push_back(rangeTmp);
}
vtr<result_path> result;
for (size_t i = 0; i < warpings.size(); i++)
result_path warp;
if (!ranges[i].empty()) {
warp = warpings[i];
warp.path = warpings[i].path.substr(ranges[i][0].start, ranges[i][0].end - ranges[i][0].start + 1);
warp.pathCoords = vtr<coord>(warpings[i].pathCoords.begin() + ranges[i][0].start, warpings[i].pathCoords.begin() + ranges[i][0].end + 1);
result.push_back(warp);
}
warpings = result;
///sumary:
///parameters:
coord dtw::get_relaxedEnds(vtr2<node> const &m, parameter const ¶ms)
double min = constant::MAX_double;
coord coordMin;
int lenA = (int)m.size() - 1;
int lenB = (int)m[0].size() - 1;
int startA = (lenA - params.relax) < 0 ? 0 : lenA - params.relax;
int startB = (lenB - params.relax) < 0 ? 0 : lenB - params.relax;
if (lenA < lenB)
startB = 0;
else if (lenB < lenA)
startA = 0;
}
for (size_t i = startA; i < m.size(); i++)
{
if (m[i][lenB].value <= min)
min = m[i][lenB].value;
coordMin.row = (int)i;
coordMin.col = lenB;
for (size_t i = startB; i < m[0].size(); i++)
min = m[lenA][i].value;
coordMin.row = lenA;
coordMin.col = (int)i;
///sumary:
///parameters:
result_path dtw::matrix_tiled(vtr2<double> const &A, vtr2<double> const &B, parameter const ¶ms)
vtr2<node> m(A.size() + 1);
m[i] = vtr<node>(B.size() + 1);
for (int i = 0; i < min((int)A.size(), params.relax + 1); i++) //include 0,0 = 0 !!
for (int i = 0; i < min((int)B.size(), params.relax + 1); i++)
if (params.isSubsequence() && calcul::lenRatio(A.size(), B.size()) < params.subsequence) {
if (A.size() < B.size())
for (size_t i = 0; i < B.size() + 1; i++)
m[0][i].value = 0;
if (B.size() < A.size())
for (size_t i = 0; i < A.size() + 1; i++)
m[i][0].value = 0;
}
int block = params.block;
//const int w = (int)(B.size() * params.w);
for (size_t i = 0; i < m.size(); i += block) //row - y
{
for (size_t j = 0; j < m[0].size(); j += block) //col - x
{
size_t iis = i + 1;
size_t jjs = j + 1;
size_t iie = iis + block < m.size() ? iis + block : m.size();
size_t jje = jjs + block < m[0].size() ? jjs + block : m[0].size();
for (size_t ii = iis; ii < iie; ii++)
{
for (size_t jj = jjs; jj < jje; jj++)
{
/*const double u = m[ii - 1][jj].value;
const double l = m[ii][jj - 1].value;
const double d = m[ii - 1][jj - 1].value;
double min = std::min({ u, l, d });*/
//_m_prefetchw(m[ii - 1][jj].value);
double min = std::min({ m[ii - 1][jj].value, m[ii][jj - 1].value, m[ii - 1][jj - 1].value });
m[ii][jj].value = calcul::distance_dtw_euklid(A[ii - 1], B[jj - 1]) + min;
auto end = get_relaxedEnds(m, params);
auto back = get_warping(m, end, params);
cout << endl << print::distanceMatrix(m);
//cout << endl << print::printPathShape(back.path, end, (int)A.size() + 1, (int)B.size() + 1);
//print::write(print::printHtmlDistanceMatrix<T>(m), "c:\\code\\data\\sc\\dm.html", false);
}
return back;
}
///sumary:
///parameters:
result_path dtw::matrix_memoized(vtr2<double> const &A, vtr2<double> const &B, parameter const ¶ms)
vtr<node2<double>> m1(B.size() + 1);
vtr<node2<double>> m2(B.size() + 1);
for (size_t i = 1; i < A.size() + 1; i++) //row y
m2[i].pa.resize(A.size() + B.size());
for (size_t j = 1; j < B.size() + 1; j++) //col x
{
const double u = m1[j].value;
const double d = m1[j - 1].value;
const double l = m2[j - 1].value;
int to = 0;
double min = std::min({ d, u, l });
int pathSize = 0;
if (min == d){
m2[j].pa = std::move(m1[j - 1].pa);
m2[j].pa[m2[j].pathSize * 2] = 0;
m2[j].pa[m2[j].pathSize * 2 + 1] = 0;
m2[j].pathSize = m1[j - 1].pathSize + 1;
else if (u < l)
to = 1;
else if (l < u)
to = 2; // L
if (m2[j - 1].pathSize < m1[j].pathSize)
to = 2;
else if (m1[j].pathSize < m2[j - 1].pathSize)
to = 1;
else
{
double coefA = A.size() / (double)i;
double coefB = B.size() / (double)j;
if (coefA > coefB)
to = 2;
else
to = 1;
m2[j].pa = m1[j].pa;
m2[j].pa[m2[j].pathSize * 2] = 0;
m2[j].pa[m2[j].pathSize * 2 + 1] = 1;
m2[j].pathSize = m1[j].pathSize + 1;
m2[j].pa = m1[j - 1].pa;
m2[j].pa[m2[j].pathSize * 2] = 1;
m2[j].pa[m2[j].pathSize * 2 + 1] = 0;
m2[j].pathSize = m2[j - 1].pathSize + 1;
m2[j].value = calcul::distance_dtw_euklid(A[i - 1], B[j - 1]) + min;
/*m1 = move(m2);
m2.resize(B.size() + 1);
m2[0].value = constant::MAX_double;*/
m1.swap(m2);
m1[0].value = constant::MAX_double;
back.start.col = 0;
back.start.row = 0;
back.end.col = (int)B.size();
back.end.row = (int)A.size();
///sumary:
///parameters:
result_path dtw::matrix_diagonal(vtr2<double> const &A, vtr2<double> const &B, parameter const ¶ms)
{
size_t simds = A[0].size() / 4;
size_t rest = A[0].size() % 4;
//double m1[A.size()] __attribute__((aligned(A.size())));
//double m2[A.size()] __attribute__((aligned(A.size())));
//double m3[A.size()] __attribute__((aligned(A.size())));
size_t shorter = std::min(A.size(), B.size());
size_t lenDiff = std::abs((int)A.size() - (int)B.size());
vtr<node2<double>> m1(shorter + 2);
vtr<node2<double>> m2(shorter + 2);
vtr<node2<double>> m3(shorter + 2);
m1[0].value = 0;
//double coef = A.size() / (double)B.size();
//__asm__(
//
//)
//size_t iEnd = shorter;
for (size_t i = 0; i < A.size(); i++)
{
size_t a = i;
size_t jEnd = i >= B.size() ? B.size() + 1 : i + 2;
for (size_t j = 1; j < jEnd; j++)
{
node2<double> l = m2[j - 1]; //l
node2<double> u = m2[j]; //u
node2<double> d = m1[j - 1]; //d
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
double min = std::min({ l.value, u.value, d.value });
//m3[j] = calcul::distance_dtw(A[a], B[j - 1]) + min;
m3[j] = calcul::distance_dtw_simd(A[a], B[j - 1], simds, rest) + min;
//cout << a - 1 << ", " << j - 1 << ", " << j - 1 << ", " << j << ", " << j - 1 << endl;
if (min == d.value)
m3[j].pathSize = d.pathSize + 1;
else if (u.value < l.value)
m3[j].pathSize = u.pathSize + 1;
else if (l.value < u.value)
m3[j].pathSize = l.pathSize + 1;
else
{
if (l.pathSize < u.pathSize)
m3[j].pathSize = l.pathSize + 1;
else if (u.pathSize < l.pathSize)
m3[j].pathSize = u.pathSize + 1;
else
{
double coefA = A.size() / (double)a;
double coefB = B.size() / (double)j;
if (coefA > coefB)
m3[j].pathSize = l.pathSize + 1;
else
m3[j].pathSize = u.pathSize + 1;
}
}
a--;
}
//cout << a - 1 << ", " << j - 1 << ", " << j - 1 << ", " << j << ", " << j - 1 << endl;
m1.swap(m2);
m2.swap(m3);
m3[0].value = constant::MAX_double;
}
rotate(m2.begin(), m2.begin() + 1, m2.end());
for (size_t i = 1; i < B.size(); i++)
{
size_t a = A.size() - 1;
size_t jEnd = A.size() - i;
// A.size() - (i - lenDiff);
if (B.size() > A.size())
{
if (i > lenDiff)
jEnd = A.size() - (i - lenDiff);
else
jEnd = A.size();
}
else
jEnd = std::min(A.size(), B.size()) - i; // -1 for decreasing size of last triangle
for (size_t j = 0; j < jEnd/* A.size() - i*/; j++)
{
node2<double> l = m2[j]; //l
node2<double> u = m2[j + 1]; //u
node2<double> d = m1[j + 1]; //d
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
//__m256d v256 = _mm256_set_pd(l.value, u.value, d.value, dMAX)
//__mm256_min
double min = std::min({ l.value, u.value, d.value });
//m3[j] = calcul::distance_dtw(A[a], B[j + i]) + min;//l,u,d
m3[j] = calcul::distance_dtw_simd(A[a], B[j + i], simds, rest) + min;
//cout << a << ", " << j << ", " << j << ", " << j + 1 << ", " << j + 1 << endl;
if (min == d.value)
m3[j].pathSize = d.pathSize + 1;
else if (u.value < l.value)
m3[j].pathSize = u.pathSize + 1;
else if (l.value < u.value)
m3[j].pathSize = l.pathSize + 1;
else
{
if (l.pathSize < u.pathSize)
m3[j].pathSize = l.pathSize + 1;
else if (u.pathSize < l.pathSize)
m3[j].pathSize = u.pathSize + 1;
else
{
double coefA = A.size() / (double)a + 1;
double coefB = B.size() / (double)j + 1;
if (coefA > coefB)
m3[j].pathSize = l.pathSize + 1;
else
m3[j].pathSize = u.pathSize + 1;
}
}
a--;
}
m1.swap(m2);
m2.swap(m3);
}
//back.path = m2[0].path;
back.scoreRaw = m2[0].value;
back.start = { 0, 0 };
back.end = { (int)A.size(), (int)B.size() };
///sumary:
///parameters:
result_path dtw::matrix_simd(vtr2<double> const &A, vtr2<double> const &B, parameter const ¶ms)
{
size_t simds = A[0].size() / 4;
size_t rest = A[0].size() % 4;
//double m1[A.size()] __attribute__((aligned(A.size())));
//double m2[A.size()] __attribute__((aligned(A.size())));
//double m3[A.size()] __attribute__((aligned(A.size())));
size_t shorter = std::min(A.size(), B.size());
size_t lenDiff = std::abs((int)A.size() - (int)B.size());
vtr<node2<double>> m1(shorter + 2);
vtr<node2<double>> m2(shorter + 2);
vtr<node2<double>> m3(shorter + 2);
m1[0].value = 0;
//double coef = A.size() / (double)B.size();
for (size_t i = 0; i < A.size(); i++)
{
size_t a = i;