Skip to content
Snippets Groups Projects
dtw.cpp 32.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • #include "stdafx.h"
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    #include "dtw.h"
    #include "calcul.h"
    #include "print.h"
    #include "help.h"
    
    #include "draw.h"
    
    #include "tree.h"
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    using namespace std;
    
    
    #undef min
    #undef max
    
    
    ///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 &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    
    	if (params.segmented)
    		result = main_segment(input, info, params);
    
    		return main_pair(input, info, params);
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    
    ///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 &params)
    
    	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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    
    	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 &params)
    
    {
    	//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;
    
    /*
    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 &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {  	
    
    	DISTANCE d;
    	d.classic = calcul::distance_dtw_euklid;
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    	if (params.distance == 1)
    
    		d.classic = calcul::distance_dtw_euklid;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	else if (params.distance == 2)
    
    		d.classic = calcul::distance_dtw_manhattan;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	else if (params.distance == 3)
    
    		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;
    
    	result_path back;
    
    	result_dtw result;
    
    	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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    		else
    
    			result = dtw::alignment(input, info, d, params);		
    	}	
    
    	
    	return result;
    }
    
    ///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 &params)
    
    	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);
    
    	 
    	if (params.isDebugInfo())
    	{
    
    		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;
    
    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 &params)
    {
    	auto m = dtw::matrix_noaccumulation(input, distance, params);
    
    	auto mMin = dtw::get_minimums(m, params);
    
    	dtw::accumulate_mod(m, mMin, params);
    
    	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 &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	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++)
    
    		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 / 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;
    
    
    			double min = std::min({ u, l, d });
    
    			if(params.distance == 3)
    
    				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;
    
    vtr2<node> dtw::matrix_noaccumulation(input_method const &input, DISTANCE distance, parameter const &params)
    
    {
    	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]);
    
    	return m;
    
    void dtw::accumulate(vtr2<node> &m, parameter const &params)
    
    {	
    	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
    
    Martin Rusek's avatar
    Martin Rusek committed
    		{
    			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 &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	int lenA = (int)m.size();
    	int lenB = (int)m[0].size();
    
    
    	const double coef = lenB / (double)lenA;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    		for (size_t j = start; j < end; j++) //col = x
    
    			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 &params)
    
    
    	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));
    
    	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 &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	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());
    
    vtr<result_path> dtw::get_warpings(vtr2<node> const &m, vtr<coord> const &minims, parameter const &params)
    
    {
    	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;
    	}
    	
    
    void dtw::filterPaths(vtr2<node> const &m, vtr<result_path> &warpings, parameter const &params)
    
    	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);
    
    ///sumary:
    ///parameters:
    coord dtw::get_relaxedEnds(vtr2<node> const &m, parameter const &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	double min = constant::MAX_double;
    
    
    	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 (params.isSubsequence()) 
    
    		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++)
    
    		if (m[lenA][i].value <= min)
    
    			min = m[lenA][i].value;
    			coordMin.row = lenA;
    			coordMin.col = (int)i;
    
    result_path dtw::matrix_tiled(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
    
    	vtr2<node> m(A.size() + 1);
    
    Martin Rusek's avatar
    Martin Rusek committed
    	for (size_t i = 0; i < A.size() + 1; i++)
    
    		m[i] = vtr<node>(B.size() + 1);
    
    	for (int i = 0; i < min((int)A.size(), params.relax + 1); i++) //include 0,0 = 0 !!
    
    Martin Rusek's avatar
    Martin Rusek committed
    		m[i][0].value = 0;
    
    
    	for (int i = 0; i < min((int)B.size(), params.relax + 1); i++)
    
    Martin Rusek's avatar
    Martin Rusek committed
    		m[0][i].value = 0;
    
    
    	if (params.isSubsequence() && calcul::lenRatio(A.size(), B.size()) < params.subsequence) {
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    		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
    	{
    
    Martin Rusek's avatar
    Martin Rusek committed
    		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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    	back.scoreRaw = m[end.row][end.col].value;
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    	if (params.isDebugInfo())
    	{
    
    		cout << endl << print::distanceMatrix(m);
    
    Martin Rusek's avatar
    Martin Rusek committed
    		//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;
    }
    
    
    result_path dtw::matrix_memoized(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
    
    	vtr<node2<double>> m1(B.size() + 1);
    	vtr<node2<double>> m2(B.size() + 1);
    
    Martin Rusek's avatar
    Martin Rusek committed
    	m1[0].value = 0;
    		
    
    	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;
    				
    
    Martin Rusek's avatar
    Martin Rusek committed
    			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;
    
    Martin Rusek's avatar
    Martin Rusek committed
    			else if (u < l)
    				to = 1;
    			else if (l < u)
    				to = 2; // L
    
    Martin Rusek's avatar
    Martin Rusek committed
    				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;
    
    Martin Rusek's avatar
    Martin Rusek committed
    			if(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;
    
    Martin Rusek's avatar
    Martin Rusek committed
    			}
    			else if (to == 2)
    			{
    
    				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;
    
    	result_path back;
    
    	back.scoreRaw = m1[m2.size() - 1].value;
    
    	back.start.col = 0;
    	back.start.row = 0;
    	back.end.col = (int)B.size();
    	back.end.row = (int)A.size();
    
    result_path dtw::matrix_diagonal(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    	m1[0].value = 0;
    
    	//double coef = A.size() / (double)B.size();
    	//__asm__(
    	//	
    	//)
    
    	//size_t iEnd = shorter;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	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
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    			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;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	}
    
    	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
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    										//__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);
    	}
    
    
    	result_path back;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	//back.path     = m2[0].path;
    	back.scoreRaw = m2[0].value;
    
    	back.start = { 0, 0 };
    	back.end = { (int)A.size(), (int)B.size() };
    
    result_path dtw::matrix_simd(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	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);
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    	m1[0].value = 0;
    
    	//double coef = A.size() / (double)B.size();
    
    
    	//size_t iEnd = shorter;		
    
    Martin Rusek's avatar
    Martin Rusek committed
    	for (size_t i = 0; i < A.size(); i++)
    	{
    		size_t a = i;