calcul.cpp

#include "stdafx.h"

#include "calcul.h"
#include <ctime>
#include <numeric>
#include <deque>
#include <bitset>

using namespace std;

double calcul::distance_dtw_euklid(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;
}

double calcul::distance_dtw_euklid_mean(vtr<double> const &u, vtr<double> const &v)
{
	double sum = 0;
	for (size_t i = 0; i < u.size(); i++)
	{
		double ratio = 0;

		if (u[i] > v[i])
			ratio += 1 - v[i] / u[i];
		else if (v[i] > u[i])
			sum += 1 - u[i] / v[i];

		//sum += std::abs(1 - u[i] / (v[i] + 0.00000001));	
	}

	return sum;
}

double calcul::distance_dtw_manhattan(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;
}

double calcul::distance_dtw_simd(vtr<double> const &u, vtr<double> const &v, size_t simds, size_t rest)
{
	double sum = 0;

	//if (simds > 0) 
	//{
	//	//_mm256_zeroall();

	//	size_t iEnd = u.size() - 4 < 0 ? 0 : u.size();

	//	for (size_t i = 0; i < u.size() - rest; i += 4)
	//	{
	//		//__m256d u256 = _mm256_set_pd(u[i], u[i + 1], u[i + 2], u[i + 3]);
	//		//__m256d v256 = _mm256_set_pd(v[i], v[i + 1], v[i + 2], v[i + 3]);

	//		__m256d result = _mm256_sub_pd(_mm256_set_pd(u[i], u[i + 1], u[i + 2], u[i + 3]), _mm256_set_pd(v[i], v[i + 1], v[i + 2], v[i + 3]));
	//		result = _mm256_mul_pd(result, result);
	//		
	//		auto r = result.m256d_f64; // == double* f = (double*)&result;

	//		sum += r[0] + r[1] + r[2] + r[3];
	//	}
	//}

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

double calcul::distance_dtw_csiChroma(vtr<double> const &u, vtr<double> const &v, double treshold = 0.07)
{
	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] > treshold)
			u_filter[i] = 1;

		if (v[i] > treshold)
			v_filter[i] = 1;
	}
	
	int minU = constant::MAX_int;
	int minV = constant::MAX_int;
	int idxU = 0;
	int idxV = 0;

	//search minimal distance between filtered vectors and minor/major scale
	for (size_t i = 0; i < scaleChord.size(); i++)
	{
		int sumU = 0;
		int sumV = 0;
		for (size_t j = 0; j < u.size(); j++)
		{
			sumU += std::abs(u_filter[j] - scaleChord[i][j]);
			sumV += std::abs(v_filter[j] - scaleChord[i][j]);
			//max value == 1 -> pow 2 -> abs
		}

		if (sumU < minU)
		{
			minU = sumU;
			idxU = static_cast<int>(i);
		}

		if (sumV < minV)
		{
			minV = sumV;
			idxV = static_cast<int>(i);
		}
	}

	//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 (scaleChord[idxU][i] == 1)
			tmpU = 0;
		else
			tmpU = u[i];

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

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

	return result;
}

double calcul::distance_dtw_csiChord(vtr<double> const &u, vtr<double> const &v, vtr<int> const &uKey, vtr<int> const &vKey)
{
	int minU = constant::MAX_int;
	int minV = constant::MAX_int;
	int idxU = 0;
	int idxV = 0;

	//search minimal distance between filtered vectors and minor/major scale
	for (size_t i = 0; i < scaleChord.size(); i++)
	{
		int sumU = 0;
		int sumV = 0;
		for (size_t j = 0; j < u.size(); j++)
		{			
			sumU += (int)std::abs(u[j] - scaleChord[i][j]);			//max value == 1 -> pow 2 -> abs
			sumV += (int)std::abs(v[j] - scaleChord[i][j]);
		}

		if (sumU < minU)
		{
			minU = sumU;
			idxU = static_cast<int>(i);
		}

		if (sumV < minV)
		{
			minV = sumV;
			idxV = static_cast<int>(i);
		}
	}

	if (idxU > 11) idxU -= 12;
	if (idxV > 11) idxV -= 12;
	auto dist1 = calcul::distance_circleFifth(idxU, idxV);	//distance between chord roots
	auto dist2 = calcul::distance_circleFifth(uKey[uKey.size() - 1], vKey[vKey.size() - 1]); //dist key roots
	
	auto getPitch = [](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
		pitch[idx + 19 % 12] = 1;				//level b 2.
		pitch[idx + 24] = 1;					//level c root
		pitch[idx + 31 % 24] = 1;				//level c 3.

		if (idx > 11)							//level c 2.
			pitch[idx + 27 % 24] = 1;
		else
			pitch[idx + 28 % 24] = 1;

		return pitch;
	};

	auto pitchU = getPitch(idxU);					//pitch space of 4 levels (level e is irrelevant)
	auto pitchV = getPitch(idxV);
	
	//minU = int_max;							//reinit variables
	//minV = int_max;
	//idxU = 0;
	//idxV = 0;

	////search minimal distance between filtered vectors and minor/major scale
	//for (size_t i = 0; i < scaleChord.size(); i++)
	//{
	//	int sumU = 0;
	//	int sumV = 0;
	//	for (size_t j = 0; j < uKey.size(); j++)
	//	{
	//		sumU += std::abs(uKey[j] - scaleKey[i][j]);
	//		sumV += std::abs(vKey[j] - scaleKey[i][j]);
	//		//max value == 1 -> pow 2 -> abs
	//	}

	//	if (sumU < minU)
	//	{
	//		minU = sumU;
	//		idxU = static_cast<int>(i);
	//	}

	//	if (sumV < minV)
	//	{
	//		minV = sumV;
	//		idxV = static_cast<int>(i);
	//	}
	//}

	//auto dist2 = calcul::distance_circleOfFifth(idxU > 11 ? idxU - 12 : idxU, idxV > 11 ? idxV - 12 : idxV);

	for (size_t i = 0; i < 12; i++)				//add level d to pitch spaces 
	{
		pitchU[i + 36] = uKey[i];		
		pitchV[i + 36] = vKey[i];
	}

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

	return dist1 + dist2 + dist3;				
}

double calcul::distance_dtw_many(vtr2<double> const &points)
{
	double sum = 0;

	for (size_t i = 0; i < points.size(); i++)
	{
		for (size_t j = i + 1; j < points.size(); j++)
		{
			sum += distance_dtw_euklid(points[i], points[j]);
		}
	}

	return sum;
}

double calcul::distance_lcss(vtr<double> const &u, vtr<double> const &v, int idx)
{
	if (u.size() == 0 || v.size() == 0)
		return 0;

	double sum = abs(u[idx] - v[idx]);

	return sum;
}

double calcul::score_dtw_s1(double ratioRaw)
{
	return ratioRaw;
}

double calcul::score_dtw_s2(double ratioRaw, size_t pathLength)
{
	return sqrt(ratioRaw) / static_cast<double>(pathLength);
}

double calcul::score_dtw_s3(size_t lenA, size_t lenB, size_t lenPath)
{
	return 1 - (lenA / static_cast<double>(lenPath) + lenB / static_cast<double>(lenPath)) / 2;
}

double calcul::score_dtw_s4(double ratioRaw, double ratioRawMax)
{
	return (sqrt(ratioRaw) / sqrt(ratioRawMax));
}

double calcul::score_dtw_s5(double ratioRaw, double ratioRawMax, double coeficient)
{
	return (sqrt(ratioRaw) / sqrt(ratioRawMax)) * coeficient;
}

double calcul::score_dtw_max(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)		//added
			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)		//added
			min = sum;

		if (max < sum)
			max = sum;
	}

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

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

double calcul::ratio(size_t dividend, size_t divisor)
{
	return dividend / (double)divisor;
}

double calcul::score_lcss_s1(double ratioRaw, size_t pathLength)
{
	return 1 - ratioRaw / static_cast<double>(pathLength);
}

double calcul::score_lcss_s2(double ratioRaw, size_t maxABLen)
{
	return 1 - ratioRaw / maxABLen;
}

double calcul::score_lcss_s3(double ratioRaw)
{
	return ratioRaw;
}

//double calcul::getPairRatio_lcss(int lenIN, int rawscore)
//{
//	return rawscore / (lenIN / 2.0);
//}

double calcul::score_multi_dtw(vtr3<double> const &input, vtr3<double> const &output)
{
	size_t sumA = 0;
	for (auto s : input)
		sumA += s.size();

	size_t sumB = 0;
	for (auto s : output)
		sumB += s.size();

	return sumA / static_cast<double>(sumB);
}

//double calcul::GetMRRratio(vtr2<int> const &orderMatrix, map<int, int> &clusters)
//{
//    double mapRatio = 0;
//    for (size_t i = 0; i < orderMatrix.size(); i++)         //query (ref sequence) //column
//    {
//        double pr = 0;
//        double coverIdx = 0;
//        for (size_t j = 1; j < orderMatrix.size(); j++)     // row / rank
//        {
//            if (clusters[orderMatrix[j][i]] == clusters[i + 1])
//                //pr += 1.0 / j;
//                pr += ++coverIdx / j;
//        }
//        mapRatio += pr / coverIdx;
//        cout << setw(5) << pr / coverIdx << " ";
//    }
//    cout << endl;
//
//    return (orderMatrix.size() - 1);
//}

double calcul::score_averageRank(int ref, vtr<int> const &queryAnswer, input_clusters 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;
}

double calcul::score_averagePrecision(int ref, vtr<int> const &queryAnswer, input_clusters 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;
}

double calcul::score_map(vtr<double> const &averagePrecisions)
{	
	//vtr<double> ratios;
	//double mapRatio = 0;
	//for (size_t i = 0; i < orderMatrix.size(); i++)         //query (ref sequence) //column
	//{
	//	double pr = 0;
	//	double coverIdx = 0;
	//	for (size_t j = 1; j < orderMatrix.size(); j++)     // row / rank
	//	{
	//		if (clusters.strIDcID.at(orderMatrix[j][i]).idCluster == clusters.strIDcID.at((int)i + 1).idCluster)
	//			pr += ++coverIdx / j;
	//	}
	//	double tmpPrecision = 0;
	//	if (pr == 0 && coverIdx == 0)
	//		tmpPrecision = 0;
	//	else
	//		tmpPrecision += pr / coverIdx;

	//		mapRatio += tmpPrecision;
	//	
	//	ratios.push_back(tmpPrecision);
	//}

	//ratios.push_back(mapRatio / (orderMatrix.size()));	

	return accumulate(averagePrecisions.begin(), averagePrecisions.end(), 0.0) / averagePrecisions.size();
}

double calcul::score_precision(int ref, vtr<int> const &queryAnswer, input_clusters 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;
}

double calcul::score_recall(int ref, vtr<int> const &queryAnswer, input_clusters const &clusters)
{	
	int clusterSize = 0;

	//if (clusters.size.size() < ref)
	//	return 0;
	//else
		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;
}

//
//double GetMultiRatiolcss(Vtr<string, 3> input, vector<int> raws)
//{
//    size_t sumA = 0;
//    for (auto s: input)
//        sumA += s.size();
//
//    int sumB = 0;
//    for (auto s: raws)
//        sumB += s;
//
//    return sumA / (double)sumB;
//}

double calcul::lb_keogh(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.push_back(0);
	minfifo.push_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]) { //overshoot
			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.push_back(static_cast<int>(i));
		minfifo.push_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;
}

double calcul::distance_circleFifth(int idx1, int idx2)
{
	double distance = std::abs(idx1 - idx2);

	if (distance > 6)
		distance = 6 - (distance - 6);

	return distance;
}

int calcul::dtw_wpassStart(size_t i, int w, double coeficient)
{
	return std::max(1, (int)(ceil((i - 1) * coeficient + 0.0000000001) - w));
}

int calcul::dtw_wpassEnd(size_t i, int w, double coeficient, int lenB)
{
	return std::min(lenB + 1, (int)(ceil(i * coeficient) + 1) + w);
}

bool calcul::dtw_isFlexiblePass(int i, int j, int iPast, int jPast, int w, int d)
{
	//int kd = k - d;
//	return std::abs((i * k - i * kd) - (j * k - j * kd)) < w;
	return std::abs((i - iPast) - (j - jPast)) < w;
}

double calcul::vtr_mean(vtr<double> const &v)
{
	double sum = 0;
	for (auto &i : v)
		sum += i;

	return sum / v.size();
}

double calcul::vtr_std(vtr<double> const &v)
{
	double mean = vtr_mean(v);

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

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

template <typename T>
T calcul::vtr_max(vtr<T> const &input)
{
	double max = constant::MIN_double;

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

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

template <typename T>
T calcul::vtr_max(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::vtr_max<double>(vtr2<double> const &input);

template <typename T>
T calcul::vtr_max(vtr3<T> const &input)
{
	double max = constant::MIN_double;

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

		if (tmp > max)
			max = tmp;
	}

	return max;
}

template <typename T>
T calcul::vtr_min(vtr2<T> const &input)
{
	double min = constant::MAX_double;

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

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

vtr<double> calcul::vtr_mean(vtr2<double> const &serie)
{
	vtr<double> average(serie.size());
	for (size_t i = 0; i < serie.size(); i++) {
		double avg = 0;
		for (size_t j = 0; j < serie[0].size(); j++) {
			avg += serie[i][j];
		}
		average[i] = avg / serie[0].size();
	}

	return average;
}