calcul.cpp

#pragma once
#include "calcul.h"
#include <ctime>
#include <numeric>
#include <deque>

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_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_csi(vtr<double> const &u, vtr<double> const &v)
{
	double result = 0;

	return result;
}

//double calcul::distancesMany_dtw(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::scorePair_dtw_s1(double ratioRaw)
//{
//	return ratioRaw;
//}
//
//double calcul::scorePair_dtw_s2(double ratioRaw, size_t pathLength)
//{
//	return sqrt(ratioRaw) / static_cast<double>(pathLength);
//}
//
//double calcul::scorePair_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::scorePair_dtw_s4(double ratioRaw, double ratioRawMax)
//{
//	return /*1 -*/ (sqrt(ratioRaw) / sqrt(ratioRawMax));
//}
//
//double calcul::scorePair_dtw_s5(double ratioRaw, double ratioRawMax, double coeficient)
//{
//	return (/*1 -*/ sqrt(ratioRaw) / sqrt(ratioRawMax)) * coeficient;
//}

double calcul::scorePair_dtw_max(vtr2<double> const &A, vtr2<double> const &B, coords start, coords end)
{
	double min = numeric_limits<double>::max();
	double max = numeric_limits<double>::min();
		
	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::coeficient_auto(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::coeficient(size_t dividend, size_t divisor)
//{
//	return dividend / (double)divisor;
//}

//double calcul::scorePair_lcss_s1(double ratioRaw, size_t pathLength)
//{
//	return 1 -  ratioRaw / static_cast<double>(pathLength);
//}
//
//double calcul::scorePair_lcss_s2(double ratioRaw, size_t maxABLen)
//{
//	return 1 - ratioRaw / maxABLen;
//}
//
//double calcul::scorePair_lcss_s3(double ratioRaw)
//{
//	return ratioRaw;
//}

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

double calcul::scoreMulti_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::scoreAverageRank(int ref, vtr<int> const &queryAnswer, clusters const &clusters)
{
	const 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::scoreAveragePrecision(int ref, vtr<int> const &queryAnswer, clusters const &clusters)
{	
	const 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::scoreMap(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::scorePrecision(int ref, vtr<int> const &queryAnswer, clusters const &clusters)
{
	int truePositives = 0;
	int falsePositives = 0;
	//int range = min(topThreshold, (int)top.size());
	for (int i = 0; i < clusters.getClusterSize(ref) - 1; i++)
	{
		if (clusters.getClusterID(queryAnswer[i]) == clusters.getClusterID(ref))
			truePositives++;
		else
			falsePositives++;
	}

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

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

//bool calcul::isInWindow(int row, int col, float ratio, int percent)
//{
//	if (ratio < 1)
//	{
//		float r = 1 / ratio;
//		float tmp = ceil(col * r - row);
//		if (tmp >= -(percent - 1) * r && tmp < (r * percent))
//			return true;
//	}
//	else
//	{
//		float tmp = ceil(row * ratio - col);
//		if (tmp >= -(percent - 1) * ratio && tmp < (ratio * percent))
//			return true;
//	}
//	return false;
//}

double calcul::ts_mean(vtr2<double> const &ts)
{
	double mean = 0;
	for (size_t i = 0; i < ts.size(); i++) //lenght of sequence
	{
		mean += ts[i][0];
	}

	return mean / ts.size();
}

//double calcul::ts_std(tseries const &ts)
//{
//
//}

//
//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 (size_t i = 1; i < 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((int)i);
		minfifo.push_back((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 (i - maxfifo.front() >= width) maxfifo.pop_front();
		if (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::lowerBound_keogh(vtr2<double> const &A, vtr2<double> const &B, parameter const &params)
//{
//	vtr<double> U(B.size());
//	vtr<double> L(B.size());
//	
//	int w = (int)params.w;
//	double min(2 * w);
//	double max(2 * w);
//
//	for (size_t i = 1; i < w; i++)
//	{
//		if (max < A[i][0])
//			max = A[i][0];
//		if (min > A[i][0])
//			min = A[i][0];
//	}
//
//	for (size_t i = 0; i < A.size() - w; i++)
//	{
//		if (max < A[i + w][0])
//			max = A[i +	w][0];
//		if (min > A[i +	w][0])
//			min = A[i +	w][0];
//
//		U[i] = max;
//		L[i] = min;		
//	}
//
//	double result = 0;
//	for (size_t i = 0; i < B.size(); i++)
//	{
//		if (B[i][0] > U[i])
//			result += pow(B[i][0] - U[i], 2);
//		else if (B[i][0] < L[i])
//			result += pow(B[i][0] - L[i], 2);
//	}
//
//	return result;
//}

//float calcul::GetManyDistances(Vtr<string, 2> const  &S)
//{
//	float sum = 0;
//	size_t a = S.size();
//
//	for (int i = 0; i < a; i++)
//	{
//		for (int j = i + 1; j < a; j++)
//		{
//			sum += GetDistance(S[i], S[j]);
//		}
//	}
//
//	return sum;
//}

//float calcul::GetDistance(string* const &u, string* const &v)
//{
//	int uL = sizeof(u) / sizeof(u);
//	int vL = sizeof(v) / sizeof(v);
//	double sum = 0;
//	int shorterDim = (uL < vL) ? uL : vL;
//
//	//if (u == null || v == null)
//	//    return 0;
//	//else
//	for (int i = 0; i < shorterDim; i++)
//	{
//		//sum += SubstitutionMatrixManager.GetDistance(u[i].ToString() + v[i].ToString());
//	}
//
//	return (float)sum;
//}