Skip to content
Snippets Groups Projects
calcul.cpp 8.87 KiB
Newer Older
  • Learn to ignore specific revisions
  • Martin Rusek's avatar
    Martin Rusek committed
    #pragma once
    
    #include "stdafx.h"
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    #include "calcul.h"
    #include <ctime>
    #include <cmath>
    #include <iostream>
    #include <iomanip>
    #include <algorithm>
    
    Martin Rusek's avatar
    Martin Rusek committed
    #include <numeric>
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    using namespace std;
    
    
    double calcul::distancesMany_dtw(vtr2<double> const &points)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	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(points[i], points[j]);
    
    Martin Rusek's avatar
    Martin Rusek committed
    		}
    	}
    
    	return sum;
    }
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::distance_dtw(vtr<double> const &u, vtr<double> const &v)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	double sum = 0;
    
    	for (size_t i = 0; i < u.size(); i++)
    	{		
    		sum += pow(u[i] - v[i], 2);
    
    		//sum += fastPow(u[i] - v[i], 2);
    
    		//sum += abs(u[i] - v[i]);
    
    Martin Rusek's avatar
    Martin Rusek committed
    		//su = u[i] - v[i];
    		//sum += (su + (su >> 31)) ^ (su >> 31);
    		//su = u[i] - v[i];
    		//sum += (su > 0) ? su : checked(-su);
    	}
    
    	//return (int)Math.Sqrt(sum);
    	return sum;
    }
    
    
    inline double calcul::fastPow(const double a, const double b)
    {
    	union {
    		double d;
    		int x[2];
    	} u = { a };
    	u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447);
    	u.x[0] = 0;
    	return static_cast<int>(u.d);
    }
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::distance_lcss(vtr<double> const &u, vtr<double> const &v, int idx)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    	if (u.size() == 0 || v.size() == 0)
    		return 0;
    
    	double sum = abs(u[idx] - v[idx]);
    
    	return sum;
    }
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::scorePair_dtw_s1(double ratioRaw)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    Martin Rusek's avatar
    Martin Rusek committed
    	return ratioRaw;
    }
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::scorePair_dtw_s2(double ratioRaw, size_t pathLength)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	return sqrt(ratioRaw) / static_cast<double>(pathLength);
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::scorePair_dtw_s3(size_t lenA, size_t lenB, size_t lenPath)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	return 1 - (lenA / static_cast<double>(lenPath) + lenB / static_cast<double>(lenPath)) / 2;
    
    Martin Rusek's avatar
    Martin Rusek committed
    }
    
    double calcul::scorePair_dtw_s4(double ratioRaw, double ratioRawMax)
    {
    
    	return /*1 -*/ (sqrt(ratioRaw) / sqrt(ratioRawMax));
    
    Martin Rusek's avatar
    Martin Rusek committed
    }
    
    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)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	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++)
    
    		for (size_t j = 0; j < A[i].size(); 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)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	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);
    //}
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::scoreMulti_dtw(vtr3<double> const &input, vtr3<double> const &output)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    	size_t sumA = 0;
    	for (auto s : input)
    		sumA += s.size();
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    
    	size_t sumB = 0;
    	for (auto s : output)
    		sumB += s.size();
    
    Martin Rusek's avatar
    Martin Rusek committed
    
    
    	return sumA / static_cast<double>(sumB);
    
    //double calcul::GetMRRratio(vtr2<int> const &orderMatrix, map<int, int> &clusters)
    
    Martin Rusek's avatar
    Martin Rusek committed
    //{
    //    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);
    //}
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    double calcul::scoreAverageRank(int ref, vtr<int> const &queryAnswer, clusters const &clusters)
    
    Martin Rusek's avatar
    Martin Rusek committed
    {
    
    Martin Rusek's avatar
    Martin Rusek committed
    	const int clusterId = clusters.strIDcID.at(ref).idCluster;
    	const int clusterSize = clusters.size.at(clusterId);
    	
    	int counter = 0;
    	double averageRank = 0;
    	for (size_t i = 0; i < queryAnswer.size(); i++)         //query (ref sequence) //column
    
    Martin Rusek's avatar
    Martin Rusek committed
    		if (clusters.strIDcID.at(queryAnswer[i]).idCluster == clusterId)
    
    Martin Rusek's avatar
    Martin Rusek committed
    			averageRank += (int)i+1;
    			counter++;
    			if (counter == clusterSize - 1)
    				break;
    
    Martin Rusek's avatar
    Martin Rusek committed
    	averageRank = averageRank / (clusterSize - 1);
    	return averageRank;
    }
    
    double calcul::scoreAveragePrecision(int ref, vtr<int> const &queryAnswer, clusters const &clusters)
    {	
    	double pr = 0;
    	const int clusterId = clusters.strIDcID.at(ref).idCluster;
    	const int clusterSize = clusters.size.at(clusterId);
    
    	int counter = 0;
    	for (size_t i = 0; i < queryAnswer.size(); i++)         //query (ref sequence) //column
    	{
    		if (clusters.strIDcID.at(queryAnswer[i]).idCluster == clusterId)
    		{
    			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 &top, clusters const &clusters)
    {
    	int cluster = clusters.strIDcID.at(ref).idCluster;
    	int refClusterSize = clusters.size.at(cluster);
    
    	int truePositives = 0;
    	int falsePositives = 0;
    	//int range = min(topThreshold, (int)top.size());
    	for (int i = 0; i < refClusterSize - 1; i++)
    	{
    		if (clusters.strIDcID.at(top[i]).idCluster == cluster)
    			truePositives++;
    		else
    			falsePositives++;
    	}
    
    	double precision = (double)truePositives / (truePositives + falsePositives);
    	return precision;
    }
    
    double calcul::scoreRecall(int ref, vtr<int> const &top, clusters const &clusters)
    {
    	int cluster = clusters.strIDcID.at(ref).idCluster;
    	int refClusterSize = clusters.size.at(cluster);
    
    	int truePositives = 0;
    	for (int i = 0; i < refClusterSize - 1; i++)
    	{
    		if (clusters.strIDcID.at(top[i]).idCluster == cluster)
    			truePositives++;
    	}
    
    	double recall = (double)truePositives / (refClusterSize - 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::mean_ts(vtr2<double> const &ts)
    {
    	double mean = 0;
    	for (size_t i = 0; i < ts.size(); i++) //lenght of sequence
    	{
    		mean += ts[i][0];
    	}
    
    Martin Rusek's avatar
    Martin Rusek committed
    	return mean / ts.size();
    
    Martin Rusek's avatar
    Martin Rusek committed
    //double calcul::std_ts(tseries const &ts)
    //{
    //
    //}
    
    
    Martin Rusek's avatar
    Martin Rusek committed
    //
    //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::GetLowerBoundKim(vtr2<double> const &, vtr2<double> const &)
    
    Martin Rusek's avatar
    Martin Rusek committed
    //{
    //    double output;
    //
    //    return output;
    //}
    
    //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;
    //}