calcul.cpp

#include "stdafx.h"

#include "calcul.h"
#include <numeric>
#include <deque>
#include <bitset>
#include "cstruct.h"
#include "print.h"

using namespace std;

///Calculates distance between two vectors.
///@param[in] u element of the compared time series (n dimensional vector)
///@param[in] v element of the compared time series (n dimensional vector)
///@return Euclidean distance
double calcul::distanceDtwEuklid(vtr<double> const &u, vtr<double> const &v)
{
	double sum = 0;

	for (size_t i = 0; i < u.size(); i++)
	{		
		sum += pow(u[i] - v[i], 2);
	}

	return sum;
}

///Calculates distance between two vectors.
///@param[in] u element of the compared time series (n dimensional vector)
///@param[in] v element of the compared time series (n dimensional vector)
///@return Manhattan distance
double calcul::distanceDtwManhattan(vtr<double> const &u, vtr<double> const &v)
{
	double sum = 0;

	for (size_t i = 0; i < u.size(); i++)
	{
		sum += abs(u[i] - v[i]);
	}

	return sum;
}

///Calculates chroma distance between two 12d vectors (specific data music data format needed).
///@param[in] u element of the compared time series (12d vector)
///@param[in] v element of the compared time series (12d vector)
///@param[in] threshold used for filtering insignificant tones in the input vectors
///@return Euclidean distance
double calcul::distanceDtwCsiChroma(vtr<double> const &u, vtr<double> const &v, double threshold /*def value in header*/)
{
	vtr<int> u_filter(u.size());
	vtr<int> v_filter(v.size());

	//binary filtering
	for (size_t i = 0; i < u.size(); i++)
	{
		if (u[i] > threshold)
			u_filter[i] = 1;

		if (v[i] > threshold)
			v_filter[i] = 1;
	}

	auto getRoot = [&scale = cstruct::scaleChord](vtr<double> const &vec) {
		int idx = 0;
		int min = constant::MAX_int;
		for (int i = 0; i < (int)scale.size(); i++) 
		{
			int sum = 0;
			for (size_t j = 0; j < vec.size(); j++) 
			{
				sum += (int)std::abs(vec[j] - scale[i][j]);			//max value == 1 -> pow 2 -> abs
			}

			if (sum < min) 
			{
				min = sum;
				idx = i;
			}
		}
		return idx;
	};

	auto idxU = getRoot(u);
	auto idxV = getRoot(v);

	//calculate result (zero u/v if minimal scale vectors == 1)
	double result = 0;
	for (size_t i = 0; i < u.size(); i++)
	{
		double tmpU;
		if (cstruct::scaleChord[idxU][i] == 1) 
		{
			tmpU = 0;
		}
		else 
		{
			tmpU = u[i];
		}

		double tmpV;
		if (cstruct::scaleChord[idxV][i] == 1)
		{
			tmpV = 0;
		}
		else 
		{
			tmpV = v[i];
		}

		//result += std::abs(tmpU - tmpV);
		result += std::pow(tmpU - tmpV, 2);
	}

	return result;
}

///Calculates chord distance between two pairs of two 12d vectors.
///@param[in] u element of the compared time series (12d vector)
///@param[in] v element of the compared time series (12d vector)
///@param[in] uKey element of the compared key time series (12d vector)
///@param[in] vKey element of the compared key time series (12d vector)
///@return chord distance
double calcul::distanceDtwCsiChord(vtr<double> const &u, vtr<double> const &v, vtr<int> const &uKey, vtr<int> const &vKey)
{
	double replaceZero = false; //if true zero vtr is replaced by c dur 100010010000
	
	auto getRoot = [&replaceZero](vtr<double> const &chord) 
	{
		auto getShift = [](int idx) { return idx > 11 ? idx % 12 : idx; };
		for (int i = 0; i < (int)chord.size(); i++) 
		{
			if (chord[i]) 
			{
				if (chord[getShift(i + 7)]) 
				{
					if (chord[getShift(i + 3)]) 
						return i;
					
					if (chord[getShift(i + 4)])
						return i;
				}				
			}
		}
				
		if (chord == vtr<double>{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}) 
			replaceZero = true;
		
		/*else 
		{
			int c = 0;
			for (int i = 0; i < (int)chord.size(); i++)
				if (chord[i] > 0) c++;
			
			if (c != 3)
			{
				cout << print::vector(chord);
				return -1;
			}*/
		//}
		//chord = vtr<double>{ 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0 }; */

		return 0;
	};

	auto getIdxKey = [&scale = cstruct::scaleKey](vtr<int> const &vec) 
	{
		int idx = 0;
		int min = constant::MAX_int;
		for (int i = 0; i < (int)scale.size(); i++)
		{
			int sum = 0;
			for (size_t j = 0; j < vec.size(); j++)	
			{
				sum += (int)std::abs(vec[j] - scale[i][j]);			//max value == 1 -> pow 2 -> abs
			}

			if (sum < min)
			{
				min = sum;
				idx = i;
			}
		}
		return idx;
	};

	int idxU = getRoot(u);
	int idxV = getRoot(v);
	int idxUK = getIdxKey(uKey);
	int idxVK = getIdxKey(vKey);

	/*int shiftU = (idxU > 11) ? 3 : 4;
	if (shiftU == 3) idxU -= 12;

	int shiftV = (idxV > 11) ? 3 : 4;
	if (shiftV == 3) idxV -= 12;*/

	vtr<int> replacement{ 1,0,0,0,1,0,0,1,0,0,0,0 };
	auto getPitch = [&replaceZero, &replacement](vtr<double> const &vec, vtr<int> const &key, int idx)
	{
		vtr<bool> pitch(48);					//pitch space of 4 levels (level e is irrelevant)
		pitch[idx] = 1;							//level a root

		pitch[idx + 12] = 1;					//level b root (copy)
		pitch[12 + (idx + 7) % 12] = 1;			//level b 2.

		if (replaceZero)
		{
			for (size_t i = 0; i < 12; i++)
			{
				pitch[i + 24] = replacement[i];
				pitch[i + 36] = replacement[i];
			}
		}
		else
		{
			for (size_t i = 0; i < 12; i++)
			{ 		
				pitch[i + 24] = vec[i];				//add key to 4. level (level d) pitch spaces (entire key vector) 	
				pitch[i + 36] = vec[i];				//copy 3. level to 4.
			}
		}
			
		//for (size_t i = 0; i < 12; i++)			//copy 3. level to 4.

		for (size_t i = 0; i < 12; i++)
		{											//add key to 4. level (level d) pitch spaces (entire key vector) 
			if (key[i] == 1)
				pitch[i + 36] = key[i];
		}
		return pitch;
	};

	auto pitchU = getPitch(u, uKey, idxU/*, shiftU*/);			//pitch space of 4 levels (level e is irrelevant)
	auto pitchV = getPitch(v, vKey, idxV/*, shiftV*/);

	auto getCofDistance = [](int idx1, int idx2) {
		return cstruct::cofDistance.at(std::abs(idx1 - idx2));
	};

	auto dist1 = getCofDistance(idxU, idxV);		//distance between chord roots
	auto dist2 = getCofDistance(idxUK, idxVK);

	int dist3 = 0;									//calculation of pitch space distance
	for (size_t i = 0; i < 48; i++)
	{
		if (pitchU[i] != pitchV[i])
			dist3++;
	}

	/*int count = 0;
	for (size_t i = 0; i < pitchU.size() / 12; i++) {
		for (size_t j = 0; j < 12; j++)
		{
			if (pitchU[count++])
				cout << "1";
			else
				cout << "0";
		}
		cout << endl;
	}

	count = 0;
	for (size_t i = 0; i < pitchU.size() / 12; i++) {
		for (size_t j = 0; j < 12; j++)
		{
			if (pitchV[count++])
				cout << "1";
			else
				cout << "0";
		}
		cout << endl;
	}*/

	return dist1 + dist2 + dist3;
}

//Calculates sum of Euclidean distance between all pairs.
//@param[in] series time series
//@return sum of dtw distances
//double calcul::distance_dtw_many(vtr2<double> const &series)
//{
//	double sum = 0;
//
//	for (size_t i = 0; i < series.size(); i++)
//	{
//		for (size_t j = i + 1; j < series.size(); j++)
//		{
//			sum += distance_dtw_euklid(series[i], series[j]);
//		}
//	}
//
//	return sum;
//}

///Calculates lcss distance between two elements of the compared time series.
///@param[in] u element of the compared time series (n dimensional vector)
///@param[in] v element of the compared time series (n dimensional vector)
///@param[in] idx index
///@return lcss distance
double calcul::distanceLcss(vtr<double> const &u, vtr<double> const &v, int idx)
{
	if (u.size() == 0 || v.size() == 0) {
		return 0;
	}

	return abs(u[idx] - v[idx]);
}

///Calculates dtw raw score s1 (last cell in the distance matrix).
///@param[in] scoreRaw raw score
///@return dtw raw score s1
double calcul::scoreDtwS1(double scoreRaw)
{
	return scoreRaw;
}

///Calculates dtw score s2 (s1 divided by warping path length)
///@param[in] scoreRaw raw score (last cell in the distance matrix)
///@param[in] pathLength length of the warping path
///@return dtw score s2
double calcul::scoreDtwS2(double scoreRaw, size_t pathLength)
{
	return sqrt(scoreRaw) / static_cast<double>(pathLength);
}

///Calculates dtw score s3.
///@param[in] lenA length of the time series A
///@param[in] lenB length of the time series B
///@param[in] pathLength length of the warping path
///@return dtw score s3
double calcul::scoreDtwS3(size_t lenA, size_t lenB, size_t pathLength)
{
	return 1 - (lenA / static_cast<double>(pathLength) + lenB / static_cast<double>(pathLength)) / 2;
}

///Calculates dtw score s4.
///@param[in] scoreRaw raw score (last cell in the distance matrix)
///@param[in] scoreRawMax maximal dtw score
///@return dtw score s4
double calcul::scoreDtwS4(double scoreRaw, double scoreRawMax)
{
	return (sqrt(scoreRaw) / sqrt(scoreRawMax));
}

///Calculates dtw score s5.
///@param[in] scoreRaw raw score (last cell in distance matrix)
///@param[in] scoreRawMax maximal dtw score
///@param[in] coefficient ratio of the compared time series lengths
///@return dtw score s5
double calcul::scoreDtwS5(double scoreRaw, double scoreRawMax, double coefficient)
{
	return (sqrt(scoreRaw) / sqrt(scoreRawMax)) * coefficient;
}

///Calculates maximal dtw score.
///@param[in] A time series A
///@param[in] B time series B
///@param[in] start coordinations of the warping path start (left upper end).
///@param[in] end coordinations of the warping path end (bottom right end)
///@return maximal raw dtw score
double calcul::scoreDtwMax(vtr2<double> const &A, vtr2<double> const &B, coord start, coord end)
{
	double min = constant::MAX_double;
	double max = constant::MIN_double;
		
	size_t sA = 0;
	size_t eA = (int)A.size();
	size_t sB = 0;
	size_t eB = (int)B.size();
 
	if (A.size() < B.size()) 
	{
		sB = start.col;
		eB = end.col;
	}

	if (B.size() < A.size()) 
	{
		sA = start.row;
		eA = end.row;
	}

	for (size_t i = sA; i < eA; i++)
	{
		double sum = 0;
		for (size_t j = 0; j < A[i].size(); j++)
		{
			sum += A[i][j];
		}

		if (min > sum) 
			min = sum;

		if (max < sum)
			max = sum;
	}
	
	for (size_t i = sB; i < eB; i++)
	{
		double sum = 0;
		for (size_t j = 0; j < B[i].size(); j++)
		{
			sum += B[i][j];
		}

		if (min > sum)
			min = sum;

		if (max < sum)
			max = sum;
	}

	return pow(max - min, 2) * std::max(eA - sA, eB - sB);
}

///Calculates input time series length ratio.
///@param[in] lenA length of the time series A
///@param[in] lenB length of the time series B
///@return length ratio
double calcul::lenRatio(size_t lenA, size_t lenB)
{
	if (lenA < lenB) 
	{
		return (double)lenA / lenB;
	}
	else if (lenB < lenA) 
	{
		return (double)lenB / lenA;
	}
	else 
	{
		return 1;
	}
}

///Calculates ratio.
///@param[in] dividend
///@param[in] divisor
///@return ratio
double calcul::ratio(size_t dividend, size_t divisor)
{
	return dividend / (double)divisor;
}

///Calculates lcss score s1.
///@param[in] scoreRaw raw score (last cell in the distance matrix)
///@param[in] pathLength length of the warping path
///@return lcss score s1
double calcul::scoreLcssS1(double scoreRaw, size_t pathLength)
{
	return 1 - scoreRaw / static_cast<double>(pathLength);
}

///Calculates lcss score s2.
///@param[in] scoreRaw raw score (last cell in the distance matrix)
///@param[in] maxABLen length of the longer input time series
///@return lcss score s2
double calcul::scoreLcssS2(double scoreRaw, size_t maxABLen)
{
	return 1 - scoreRaw / maxABLen;
}

///Calculates lcss score s3.
///@param[in] scoreRaw raw score (last cell in the distance matrix)
///@return lcss score s3
double calcul::scoreLcssS3(double scoreRaw)
{
	return scoreRaw;
}

///Calculates average rank.
///@param[in] ref cluster id reference
///@param[in] queryAnswer query answer row (from id matrix (int))
///@param[in] clusters ground truth of the clusters
///@return average rank
double calcul::scoreAverageRank(int ref, vtr<int> const &queryAnswer, inputGroundTruth const &clusters)
{
	int clusterSize = clusters.getClusterSize(ref);

	int counter = 0;
	double averageRank = 0;
	for (size_t i = 0; i < queryAnswer.size(); i++)         //query (ref sequence) //column
	{
		if (clusters.getClusterID(queryAnswer[i]) == clusters.getClusterID(ref))
		{
			averageRank += (int)i + 1;
			counter++;
			if (counter == clusterSize - 1)
				break;
		}
	}

	averageRank = averageRank / (clusterSize - 1);
	return averageRank;
}

///Calculates average precision.
///@param[in] ref cluster id reference
///@param[in] queryAnswer query answer row (from id matrix (int))
///@param[in] clusters ground truth of the clusters
///@return average precision
double calcul::scoreAveragePrecision(int ref, vtr<int> const &queryAnswer, inputGroundTruth const &clusters)
{	
	int clusterSize = clusters.getClusterSize(ref);
		
	double pr = 0;
	int counter = 0;
	for (size_t i = 0; i < queryAnswer.size(); i++)         //query (ref sequence) //column
	{
		if (clusters.getClusterID(queryAnswer[i]) == clusters.getClusterID(ref))
		{
			pr += ++counter / (double)(i + 1);
			if (counter == clusterSize - 1)
				break;
		}
	}

	double ap = pr / (clusterSize - 1);
	return ap;
}

///Calculates map score (mean average precision).
///@param[in] averagePrecisions vector of average precisions
///@return map score
double calcul::scoreMap(vtr<double> const &averagePrecisions)
{	
	return accumulate(averagePrecisions.begin(), averagePrecisions.end(), 0.0) / averagePrecisions.size();
}

///Calculates precision.
///@param[in] ref cluster id reference
///@param[in] queryAnswer query answer row (from id matrix (int))
///@param[in] clusters ground truth of the clusters
///@return precision fraction
double calcul::scorePrecision(int ref, vtr<int> const &queryAnswer, inputGroundTruth const &clusters)
{
	int clusterSize = clusters.getClusterSize(ref);

	int truePositives = 0;
	int falsePositives = 0;
	//int range = min(topThreshold, (int)top.size());
	for (int i = 0; i < clusterSize - 1; i++)
	{
		if (clusters.getClusterID(queryAnswer[i]) == clusters.getClusterID(ref)) 
		{
			truePositives++;
		}
		else {
			falsePositives++;
		}
	}

	double precision = (double)truePositives / (truePositives + falsePositives);
	return precision;
}

///Calculates recall.
///@param[in] ref cluster id reference
///@param[in] queryAnswer query answer row (from id matrix (int))
///@param[in] clusters ground truth of the clusters
///@return recall fraction
double calcul::scoreRecall(int ref, vtr<int> const &queryAnswer, inputGroundTruth const &clusters)
{	
	int clusterSize = clusters.getClusterSize(ref);

	int truePositives = 0;
	for (int i = 0; i < clusterSize - 1; i++)
	{
		if (clusters.getClusterID(queryAnswer[i]) == clusters.getClusterID(ref))
			truePositives++;
	}

	double recall = (double)truePositives / (clusterSize - 1);
	return recall;
}

///Calculates Keogh's lower bound for dtw. EXPERIMENT
///@param[in] A time series A
///@param[in] B time series B
///@param[in] params parameters
///@return Keogh's lower bound distance
double calcul::lbKeogh(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
{
	int width = 1 + 2 * (int)params.w;
	int w = (int)params.w;
	deque<int> maxfifo, minfifo;
	maxfifo.emplace_back(0);
	minfifo.emplace_back(0);

	vtr<double> maxvalues(B.size());
	vtr<double> minvalues(B.size());

	for (int i = 1; i < static_cast<int>(A.size()); ++i)
	{
		if (i >= w + 1) 
		{
			maxvalues[i - w - 1] = A[maxfifo.front()][0];
			minvalues[i - w - 1] = A[minfifo.front()][0];
		}
		if (A[i][0] > A[i - 1][0]) 
		{
			maxfifo.pop_back();
			while (maxfifo.size() > 0) 
			{
				if (A[i][0] <= A[maxfifo.back()][0]) 
					break;
				
				maxfifo.pop_back();
			}
		}
		else 
		{
			minfifo.pop_back();
			while (minfifo.size() > 0) 
			{
				if (A[i][0] >= A[minfifo.back()][0])
					break;
			
				minfifo.pop_back();
			}
		}

		maxfifo.emplace_back(static_cast<int>(i));
		minfifo.emplace_back(static_cast<int>(i));
		if (i == width + maxfifo.front()) 
		{
			maxfifo.pop_front();
		}
		else if (i == width + minfifo.front()) 
		{
			minfifo.pop_front();
		}
	}

	for (size_t i = A.size(); i <= A.size() + w; ++i) 
	{
		maxvalues[i - w - 1] = A[maxfifo.front()][0];
		minvalues[i - w - 1] = A[minfifo.front()][0];
		
		if (static_cast<int>(i) - maxfifo.front() >= width)
			maxfifo.pop_front();
		
		if (static_cast<int>(i) - minfifo.front() >= width) 
			minfifo.pop_front();
	}

	double result = 0;
	for (size_t i = 0; i < A.size(); i++)
	{
		if (B[i][0] > maxvalues[i])
		{
			result += pow(B[i][0] - maxvalues[i], 2);
		}
		else if (B[i][0] < minvalues[i])
		{
			result += pow(B[i][0] - minvalues[i], 2);
		}
	}

	return result;
}

//Calculates Cosine distance.
//@param[in] u element of the compared time series
//@param[in] v element of the compared time series
//@return cosine distance
//double calcul::distance_cosine(vtr<double> const &u, vtr<double> const &v)
//{
//	double sumAB = 0;
//	for (size_t i = 0; i < u.size(); i++)
//		sumAB += u[i] * v[i];
//
//	double sumA = 0;
//	for (size_t i = 0; i < u.size(); i++)
//		sumA += pow(u[i], 2);
//
//	double sumB = 0;
//	for (size_t i = 0; i < u.size(); i++)
//		sumB += pow(v[i], 2);
//
//	return sumAB / (sqrt(sumA) * sqrt(sumB));
//}
//
//double calcul::distance_cosine2(input_point_ref const &point)
//{
//	double sumAB = 0;
//	for (size_t i = 0; i < point.p1.size(); i++)
//		sumAB += point.p1[i] * point.p2[i];
//
//	double sumA = 0;
//	for (size_t i = 0; i < point.p1.size(); i++)
//		sumA += pow(point.p1[i], 2);
//
//	double sumB = 0;
//	for (size_t i = 0; i < point.p1.size(); i++)
//		sumB += pow(point.p2[i], 2);
//
//	return sumAB / (sqrt(sumA) * sqrt(sumB));
//}

///Calculates first cell included in the warping window.
///@param[in] row row
///@param[in] win warping window width
///@param[in] coefficient input time series length ratio
///@return index of first cell in warping window
int calcul::dtwPassStart(size_t row, int win, double coefficient)
{
	return std::max(1, (int)(ceil((row - 1) * coefficient + 0.0000000001) - win));
}

///Calculates last cell included in the warping window.
///@param[in] row row
///@param[in] win warping window width
///@param[in] coefficient input time series length ratio
///@param[in] lenB length of the time series B
///@return index of last cell in warping window
int calcul::dtwPassEnd(size_t row, int win, double coefficient, int lenB)
{
	return std::min(lenB + 1, (int)(ceil(row * coefficient) + 1) + win);
	//const size_t end = min(lenB + 1, (int)(ceil(i * lenB / (double)lenA) + 1) + w);
}

///Calculates if cell is included in the flexible warping pass.
///@param[in] row row
///@param[in] col column
///@param[in] past last cell which was included
///@param[in] w flexible warping pass width
///@return TRUE if cell is in the flexible pass
///@return FALSE if cell is not in the flexible pass
bool calcul::isFlexiblePass(int row, int col, coord const &past, int w)
{
	return std::abs((past.row - row) - (past.col - col)) < w;
}

///Calculates mean value of the input vector.
///@param[in] v vector of doubles
///@return mean
double calcul::vtrMean(vtr<double> const &v)
{
	double sum = 0;

	for (auto &i : v)
	{
		sum += i;
	}

	return sum / v.size();
}

///Calculates vector of means of the 2d input vector.
///@param[in] series two dimensional vector
///@return vector of means
vtr<double> calcul::vtrMean(vtr2<double> const &series)
{
	vtr<double> average(series.size());

	for (size_t i = 0; i < series.size(); i++)
	{
		average[i] = vtrMean(series[i]);
	}

	return average;
}

///Calculates standard deviation of the input vector.
///@param[in] v vector
///@return standard deviation
double calcul::vtrStd(vtr<double> const &v)
{
	double mean = vtrMean(v);

	double sum = 0;
	for (auto &i : v)
	{
		sum += pow(i - mean, 2);
	}

	return sqrt(sum / v.size());
}

///Finds maximal value contained in the input vector.
///@tparam data type of searched vector
///@param[in] input vector
///@return maximal value
template <typename T>
T calcul::vtrMax(vtr<T> const &input)
{
	T max = (T)constant::MIN_double;

	for (auto &i : input)
	{
		if (i > max)
			max = i;
	}

	return max;
}
template int calcul::vtrMax(vtr<int> const &input);
template double calcul::vtrMax<double>(vtr<double> const &input);

///Finds maximal value contained in the 2d input vector.
///@tparam data type of searched vector
///@param[in] input vector
///@return maximal value
template <typename T>
T calcul::vtrMax(vtr2<T> const &input)
{
	double max = constant::MIN_double;

	for (auto &&i : input)
	{
		for (auto &&j : i)
		{
			if (j > max)
				max = j;
		}
	}

	return max;
}
template double calcul::vtrMax<double>(vtr2<double> const &input);

///Finds maximal value contained in the 3d input vector.
///@param[in] input vector
///@return maximal value
template <typename T>
T calcul::vtrMax(vtr3<T> const &input)
{
	double max = constant::MIN_double;

	for (auto const &i : input) 
	{
		auto tmp = vtr_findMax(i);

		if (tmp > max)
			max = tmp;
	}

	return max;
}

///Finds minimal value contained in the input vector.
///@param[in] input vector
///@return minimal value
template <typename T>
T calcul::vtrMin(vtr<T> const &input)
{
	T min = (T)constant::MAX_double;

	for (auto &i : input)
	{
		if (i < min)
			min = i;
	}

	return min;
}
template double calcul::vtrMin<double>(vtr<double> const &input);

///Finds minimal value contained in the input vector.
///@param[in] input vector
///@return minimal value
template <typename T>
T calcul::vtrMin(vtr2<T> const &input)
{
	double min = constant::MAX_double;

	for (auto const &i : input) 
	{
		for (auto const &j : i)
		{
			if (j < min)
				min = j;
		}
	}

	return min;
}
template double calcul::vtrMin<double>(vtr2<double> const &input);