dtw.cpp

#include "stdafx.h"

#include "dtw.h"
#include "calcul.h"
#include "print.h"
#include "help.h"
#include "draw.h"
#include "tree.h"

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
*/
vtr<double> dtw::main(input_method const &input, input_info const &info, parameter const &params)
{
	vtr<double> result;

	if (params.segmented)
		result = main_segment(input, info, params);
	else
		result = main_pair(input, info, params);

	return result;
}

/*
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
*/
vtr<double> dtw::main_pair(input_method const &input, input_info const &info, parameter const &params)
{
	if ((input.A.size() * input.B.size()) / 131'072 > 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;
		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] = 1 - result.score[params.scoreType - 1];

		return result.score;
	}
	else
		return result.score;
}

/*
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
*/
vtr<double> 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;

		return result[0][0].score;
	}
	else
		return result[0][0].score;
}

/*
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)
{  	
	//distanceMatrix<double> dm(A, B, params);

	DISTANCE d;
	d.classic = calcul::distance_dtw_euklid;

	if (params.distance == 1)
		d.classic = calcul::distance_dtw_euklid;
	else if (params.distance == 2)
		d.classic = calcul::distance_dtw_manhattan;
	else if (params.distance == 3)
		d.csi_chroma = calcul::distance_dtw_csi_chroma;
	else if (params.distance == 4)
		d.csi_chord = calcul::distance_dtw_csi_chord;
	else if (params.distance == 5)
		d.classic = calcul::distance_dtw_euklid_mean;

	result_path back;
	result_dtw result;

	/*if(params.matrixDataType == "double")
	{*/
		if (params.isMemoization()) 
		{
			//auto begin = chrono::steady_clock::now();
			back = matrix_memoized<double>(input.A, input.B, params);
			/*cout << print::printElapsed("v2: ", chrono::duration_cast<chrono::microseconds>(chrono::steady_clock::now() - begin).count()) << ", ";

			begin = chrono::steady_clock::now();
			back = dtw::matrix<double>(A, B, params);
			cout << print::printElapsed("v1: ", chrono::duration_cast<chrono::microseconds>(chrono::steady_clock::now() - begin).count()) << endl;

			if (back2.path != back.path)
			{
				cout << "diff: " << back2.pathSize - (int)back.path.size() << ", " << back.scoreRaw << " v2: " << back2.scoreRaw << endl;
				cout << back2.path << endl;
				cout << back.path << endl;
			}
			if (back2.scoreRaw != back.scoreRaw)
				cout << "diff: " << back2.pathSize - (int)back.path.size() << ", " << back.scoreRaw << " v2: " << back2.scoreRaw << " ~~!!!~~" << endl;*/
		
		}
		else if(params.block > 0)
			back = dtw::matrix_tiled<double>(input.A, input.B, params);
		else if(params.isSimd())
			back = dtw::matrix_simd<double>(input.A, input.B, params);
		else
			result = dtw::alignment<double>(input, info, d, params);
	//}
	/*else if(params.matrixDataType == "int")
		result = dtw::alignment<int>(input, info, d, params);
	else if (params.matrixDataType == "float")
		result = dtw::alignment<float>(input, info, d, params)*/;

		//auto noacc = createMatrix_noAccumulation(A, B, 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
*/
template<class T>
result_dtw dtw::alignment(input_method const &input, input_info const &info, DISTANCE distance, parameter const &params)
{
	vtr2<node<T>> m;
	if (params.experiment)
	{
		m = dtw::matrix_noaccumulation<T>(input, distance, params);
		auto mMinims = dtw::get_minimums(m, params);
		for (auto &i : mMinims)
			m[i.row][i.col] = 0;
		dtw::accumulate(m, params);
	}
	else
		m = dtw::matrix<T>(input, distance, params);
	
	auto end = get_relaxedEnds(m, params);
	auto warping = get_warping(m, end, params);
	warping.scoreRaw = m[end.row][end.col].value;
	warping.wpEnd = end;
	
	result_dtw result;

	if (params.drawOut.size() > 0) {
		result.matrix_acc = draw::matrix(m, warping.path, params);
		result.matrix_noacc = draw::matrix(dtw::matrix_noaccumulation<T>(input, distance, params), warping.path, params);
		//draw::visualisation<T>(m, dtw::matrix_noaccumulation<T>(input, distance, params), input, warping.path, params.drawOut);
		//m, m2, s1, s2, wp, filep
	}
	 
	if (params.isDebugInfo())
	{
		cout << endl << print::distanceMatrix<T>(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);
	}
	
	vtr<double> resultScore;
	resultScore.push_back(warping.scoreRaw);
	resultScore.push_back(calcul::score_dtw_s2(warping.scoreRaw, warping.pathSize)); // back.path.size());
	resultScore.push_back(calcul::score_dtw_s3(warping.wpEnd.row - warping.wpStart.row, warping.wpEnd.col - warping.wpStart.col, warping.pathSize));
	resultScore.push_back(calcul::score_dtw_s4(warping.scoreRaw, calcul::score_dtw_max(input.A, input.B, warping.wpStart, warping.wpEnd)));
	resultScore.push_back(calcul::score_dtw_s5(warping.scoreRaw,
						  calcul::score_dtw_max(input.A, input.B, warping.wpStart, warping.wpEnd),
						  calcul::lenRatio(warping.wpEnd.row - warping.wpStart.row, warping.wpEnd.col - warping.wpStart.row)));

	result.score = resultScore;
	result.path = warping.path;
	
	return result;
}

/*
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
*/
template<class T>
vtr2<node<T>> dtw::matrix(input_method const &input, DISTANCE distance, parameter const &params)
{
	int lenA = (int)input.A.size();
	int lenB = (int)input.B.size();

	vtr2<node<T>> m(lenA + 1);
	for (int i = 0; i < lenA + 1; i++)
		m[i] = vtr<node<T>>(lenB + 1);

	for (int i = 0; i < min(lenA, params.relax.start + 1); i++) //include 0,0 = 0 !!
		m[i][0].value = 0;

	for (int i = 0; i < min(lenB, params.relax.end + 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 int w = (int)(lenB * params.w);
	double coef = lenB / static_cast<double>(lenA);
	for (int i = 1; i < lenA + 1; i++) //row - y
	{
		size_t start = max(1, (int)(ceil((i - 1) * coef + 0.0000000001) - w));
		size_t end   = min(lenB + 1, (int)(ceil(i * coef) + 1) + w);
		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 = static_cast<T>(distance.csi_chroma(input.A[i - 1], input.B[j - 1], 0.07) + min);
			else if (params.distance == 4)
				m[i][j].value = static_cast<T>(distance.csi_chord(input.A[i - 1], input.B[j - 1], input.A2[i - 1], input.B2[j - 1]) + min);
			else
				m[i][j].value = static_cast<T>(distance.classic(input.A[i - 1], input.B[j - 1]) + min);
		}
	}

	return m;
}

/*
sumary:
parameters:
*/
template<class T>
vtr2<node<T>> 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<T>> m(lenA + 1);
	for (int i = 0; i < lenA + 1; i++)
		m[i] = vtr<node<T>>(lenB + 1);

	for (int i = 0; i < min(lenA, params.relax.start + 1); i++)
		m[i][0].value = 0;

	for (int i = 0; i < min(lenB, params.relax.end + 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 int w = (int)(lenB * params.w);
	for (int i = 1; i < lenA + 1; i++) //row = y
	{
		const size_t start = max(1, (int)(ceil((i - 1) * (lenB / (double)lenA + 0.0000000001)) - w));
		const size_t end = min(lenB + 1, (int)(ceil(i * lenB / (double)lenA) + 1) + w);
		for (size_t j = start; j < end; j++) //col = x
		{
			if (params.distance == 3)
				m[i][j].value = static_cast<T>(distance.csi_chroma(input.A[i - 1], input.B[j - 1], 0.07));
			else if (params.distance == 4)
				m[i][j].value = static_cast<T>(distance.csi_chord(input.A[i - 1], input.B[j - 1], input.A2[i - 1], input.B2[j - 1]));
			else
				m[i][j].value = static_cast<T>(distance.classic(input.A[i - 1], input.B[j - 1]));
		}
	}

	return m;
}

template<class T>
void dtw::accumulate(vtr2<node<T>> &m, parameter const &params)
{	
	int lenA = (int)m.size();
	int lenB = (int)m[0].size();

	const int w = (int)(lenB * params.w);
	for (int i = 1; i < lenA; i++) //row = y
	{
		const size_t start = max(1, (int)(ceil((i - 1) * (lenB / (double)lenA + 0.0000000001)) - w));
		const size_t end = min(lenB, (int)(ceil(i * lenB / (double)lenA)) + w);
		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 });
		}
	}
}

template<class T>
void dtw::accumulate(vtr2<node<T>> &m, vtr parameter const &params)
{
	int lenA = (int)m.size();
	int lenB = (int)m[0].size();

	const int w = (int)(lenB * params.w);
	for (int i = 1; i < lenA; i++) //row = y
	{
		const size_t start = max(1, (int)(ceil((i - 1) * (lenB / (double)lenA + 0.0000000001)) - w));
		const size_t end = min(lenB, (int)(ceil(i * lenB / (double)lenA)) + w);
		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 });
		}
	}
}

template<class T>
vtr<coords> dtw::get_minimums(vtr2<node<T>> const &m, parameter const &params)
{
	vtr<coords> minims;

	int lenA = (int)m.size();
	int lenB = (int)m[0].size();

	const int w = (int)(lenB * params.w);
	for (int i = 1; i < lenA; i++) //row = y
	{
		const int start = max(1, (int)(ceil((i - 1) * (lenB / (double)lenA + 0.0000000001)) - w));
		const int end = min(lenB, (int)(ceil(i * lenB / (double)lenA)) + w);
		for (int j = start; j < end; j++) //col = x
		{
			T current = m[i][j].value;
			
			if (i + 1 < m.size() && j + 1 < 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(coords(i, j));
			}
		}
	}

	return minims;
}

/*
sumary:
parameters:
vtr2<node<T>> - 
i
j
params
*/
template<class T>
result_path dtw::get_warping(vtr2<node<T>> const &m, coords coord, parameter const &params)
{
	result_path wp;
	double sizeA = coord.row;
	double sizeB = coord.col;

	while (coord.row > 0 && coord.col > 0)
	{
		double u = m[coord.row - 1][coord.col].value;
		double l = m[coord.row][coord.col - 1].value;
		double d = m[coord.row - 1][coord.col - 1].value;

		if (min({ d, u, l }) == d)
		{
			wp.path = "M" + wp.path;
			coord.row--;
			coord.col--;
		}
		else
		{
			if (l < u)
			{
				wp.path = "L" + wp.path;
				coord.col--;

			}
			else if (u < l)
			{
				wp.path = "U" + wp.path;
				coord.row--;
			}
			else
			{
				if(sizeA / coord.row > sizeB / coord.col)
				{
					wp.path = "L" + wp.path;
					coord.col--;
				}
				else
				{
					wp.path = "U" + wp.path;
					coord.row--;
				}
			}
		}
	}

	if (!params.isSubsequence() || (params.isSubsequence() && m.size() < m[0].size()))
		while (coord.row > params.relax.start)
		{
			wp.path = "U" + wp.path;
			coord.row--;

		}

	if (!params.isSubsequence() || (params.isSubsequence() && m[0].size() < m.size()))
		while (coord.col > params.relax.end)
		{
			wp.path = "L" + wp.path;
			coord.col--;
		}

	wp.wpStart.row = coord.row;
	wp.wpStart.col = coord.col;
	wp.pathSize = (int)wp.path.size();

	return wp;
}

template<class T>
vtr<result_path> dtw::get_warpings(vtr2<node<T>> const &m, vtr<result_path> const &paths, parameter const &params)
{
	vtr<result_path> outPaths;

	for(auto &i : paths)
		outPaths.push_back(get_warping(m, i.wpStart.row, i.wpStart.col, params))

	return paths;
}

/*
sumary:
parameters:
*/
template<class T>
coords dtw::get_relaxedEnds(vtr2<node<T>> const &m, parameter const &params)
{
	double min = constant::MAX_double;
	coords coordMin;

	int lenA = (int)m.size() - 1;
	int lenB = (int)m[0].size() - 1;

	int startA = (lenA - params.relax.start) < 0 ? 0 : lenA - params.relax.start;
	int startB = (lenB - params.relax.end) < 0 ? 0 : lenB - params.relax.end;

	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;
		}
	}

	return coordMin;
}

/*
sumary:
parameters:
*/
template<class T>
result_path dtw::matrix_tiled(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
{
	vtr2<node<T>> m(A.size() + 1);
	for (size_t i = 0; i < A.size() + 1; i++)
		m[i] = vtr<node<T>>(B.size() + 1);

	for (int i = 0; i < min((int)A.size(), params.relax.start + 1); i++) //include 0,0 = 0 !!
		m[i][0].value = 0;

	for (int i = 0; i < min((int)B.size(), params.relax.end + 1); i++)
		m[0][i].value = 0;

	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
	{
		//size_t start = max(1, (int)(ceil((i - 1) * (static_cast<double>(B.size()) / (double)A.size() + 0.0000000001)) - w));
		//size_t end = min((int)B.size() + 1, (int)(ceil(i * B.size() / (double)A.size()) + 1) + w);
		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 = static_cast<T>(calcul::distance_dtw_euklid(A[ii - 1], B[jj - 1]) + min);
				}
			}			
		}
	}

	auto end = get_relaxedEnds(m, params);
	auto back = get_warping(m, end, params);
	back.scoreRaw = m[end.row][end.col].value;
	back.wpEnd = end;

	if (params.isDebugInfo())
	{
		cout << endl << print::distanceMatrix<T>(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:
*/
template<class T>
result_path dtw::matrix_memoized(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
{
	vtr<node2<T>> m1(B.size() + 1);
	vtr<node2<T>> m2(B.size() + 1);
	
	//double coef = A.size() / (double)B.size();
	
	m1[0].value = 0;
	//m1[0].path.reserve(A.size() + B.size());
	
	string path = "";
	path.reserve(A.size() + B.size());
		
	//const int w = (int)(B.size() * params.w);
	for (size_t i = 1; i < A.size() + 1; i++) //row - y
	{
		//path = "";
		for (size_t j = 1; j < B.size() + 1; j++) //row - y
		{
			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){
				//path = std::move(m1[j - 1].path);
				//path.push_back('M');
				pathSize = m1[j - 1].pathSize + 1;
			}
			else if (u < l)
				to = 1;
			else if (l < u)
				to = 2; // L
			else 
			{
				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;
				}
			}

			if(to == 1)
			{
				//path.append(m1[j].path + "U");
				//path.push_back('U');
				pathSize = m1[j].pathSize + 1;
			}
			else if (to == 2)
			{
				//path.append(m2[j - 1].path + "L");
				//path.push_back('L');
				pathSize = m2[j - 1].pathSize + 1;
			}

			//m2[j].path = std::move(path);
			m2[j].pathSize = pathSize;
			m2[j].value = static_cast<T>(calcul::distance_dtw_euklid(A[i - 1], B[j - 1])) + min;
		}
		m1[0].value = std::numeric_limits<T>::max();
		m2.swap(m1);
	}

	result_path back;
	//back.path     = m1[m2.size() - 1].path;
	back.pathSize = m1[m2.size() - 1].pathSize;
	back.scoreRaw = m1[m2.size() - 1].value;
	back.wpStart.col = 0;
	back.wpStart.row = 0;
	back.wpEnd.col = (int)B.size();
	back.wpEnd.row = (int)A.size();
	
	return back;
}

/*
sumary:
parameters:
*/
template<class T>
result_path dtw::matrix_diagonal(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
{
	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<T>> m1(shorter + 2);
	vtr<node2<T>> m2(shorter + 2);
	vtr<node2<T>> 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<T> l = m2[j - 1];	//l
			node2<T> u = m2[j];		//u
			node2<T> d = m1[j - 1];	//d

			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<T> l = m2[j];			//l
			node2<T> u = m2[j + 1];		//u
			node2<T> d = m1[j + 1];		//d

										//__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;
	//back.path     = m2[0].path;
	back.pathSize = m2[0].pathSize;
	back.scoreRaw = m2[0].value;
	back.wpStart = { 0, 0 };
	back.wpEnd = { (int)A.size(), (int)B.size() };

	return back;
}

/*
sumary:
parameters:
*/
template<class T>
result_path dtw::matrix_simd(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
{
	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<T>> m1(shorter + 2);
	vtr<node2<T>> m2(shorter + 2);
	vtr<node2<T>> m3(shorter + 2);

	m1[0].value = 0;

	//double coef = A.size() / (double)B.size();

	//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<T> l = m2[j - 1];	//l
			node2<T> u = m2[j];		//u
			node2<T> d = m1[j - 1];	//d

			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<T> l = m2[j];			//l
			node2<T> u = m2[j + 1];		//u
			node2<T> d = m1[j + 1];		//d

			//__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)