#include // std::cout #include // std::max #include #include #include #include #include #include #include #define INF 1e20 // Pseudo-infinite number for this code #define dist(x,y) ((x-y)*(x-y)) using namespace std; /// This function is from the UCRDTW code /// /// Calculate Dynamic Time Wrapping distance /// A,B: data and query, respectively /// cb : cummulative bound used for early abandoning /// r : size of Sakoe-Chiba warpping band double dtw(double* A, double* B, double *cb, int m, int r, double best_so_far) { double *cost; double *cost_prev; double *cost_tmp; int i, j, k; double x, y, z, min_cost; /// Instead of using matrix of size O(m^2) or O(mr), we will reuse two arrays of size O(r). cost = (double*) calloc(2 * r + 1, sizeof(double)); cost_prev = (double*) calloc(2 * r + 1, sizeof(double)); for (k = 0; k < 2 * r + 1; k++) { cost[k] = INF; cost_prev[k] = INF; } for (i = 0; i < m; i++) { k = max(0, r - i); min_cost = INF; for (j = max(0, i - r); j <= min(m - 1, i + r); j++, k++) { /// Initialize all row and column if ((i == 0) && (j == 0)) { cost[k] = dist(A[0], B[0]); min_cost = cost[k]; continue; } if ((j - 1 < 0) || (k - 1 < 0)) { y = INF; } else { y = cost[k - 1]; } if ((i - 1 < 0) || (k + 1 > 2 * r)) { x = INF; } else { x = cost_prev[k + 1]; } if ((i - 1 < 0) || (j - 1 < 0)) { z = INF; } else { z = cost_prev[k]; } /// Classic DTW calculation cost[k] = min( min( x, y) , z) + dist(A[i], B[j]); /// Find minimum cost in row for early abandoning (possibly to use column instead of row). if (cost[k] < min_cost) { min_cost = cost[k]; } } /// Move current array to previous array. cost_tmp = cost; cost = cost_prev; cost_prev = cost_tmp; } k--; /// the DTW distance is in the last cell in the matrix of size O(m^2) or at the middle of our array. double final_dtw = cost_prev[k]; free(cost); free(cost_prev); return final_dtw; } /// Sorting function for the 2nd element of two vectors bool sortbysecasc(const pair &a, const pair &b) { return a.second &a, const pair &b) { return a.second>b.second; } /// Get random value from uniform distribution double uniform(double a, double b) { return (double) std::rand() / (RAND_MAX) * (b - a) + a; } /// Algorithm for calculating similarity of (multivariate) time series data int lsh(double *** D, int M, int T, int Dim, double ** Q, int L, int K, double r, double a, double sd, int * &candidates, double * &distances, double **** &hashFunctions, double ** weights, int &nrOfCandidates) { clock_t begin_time = clock(); std::random_device rd{}; std::mt19937 gen{rd()}; std::normal_distribution<> distribution0{0,1}; srand(time(NULL)); int S = 2; //3 double delta = r * 0.75; //0.07 => 20.000 int threshold = T * 0.22; //0.8 => 18/67 /// Create hash functions int value; for (int l=0;l T - 1) continue; currentMax = std::max(Q_U[t + s][d], currentMax); currentMin = std::min(Q_L[t + s][d], currentMin); } Q_U[t][d] = currentMax; Q_L[t][d] = currentMin; currentMax = -10000; currentMin = 10000; } } double temp1 = 0; double temp2 = 0; double temp3 = 0; /// Create query hashes double *** hashQ = (double***)malloc(L*sizeof(double**)) ; // Initialize Data double *** hashQ_U = (double***)malloc(L*sizeof(double**)) ; // Initialize Data double *** hashQ_L = (double***)malloc(L*sizeof(double**)) ; // Initialize Data for (int l=0;l v; v.reserve(M); bool group_collision; int hash_collision; int projection_collision; for (int m=0;m=threshold) { hash_collision++; } else break; } if (hash_collision==K) { group_collision = true; } } if (group_collision) { v.emplace_back(m); } } clock_t end_time = clock(); double time_spent = (double)(end_time - begin_time) / CLOCKS_PER_SEC; std::cout << "Hashing done in: " << time_spent << " seconds" << std::endl; std::cout << "Number of candidates pruned: " << 100 - (100 * v.size() / M) << "%" << std::endl; /// Estimate ranking using L2/DTW ratios vector> sim; sim.reserve(v.size()); double lb; double ub; double lb_it; double ub_it; double q_lb = 0; // double q_ub = sqrt(T * pow(r, 2.0)) * K * L; double q_ub = 0; for (int l=0;l> dtwsim; dtwsim.reserve(sim.size()); double distance = 0; for (int m=0;m