#include <iostream>     // std::cout
#include <algorithm>    // std::max
#include <fstream>
#include <string.h>
#include <tuple>
#include <stdlib.h>
#include <vector>
#include <math.h>
#include <random>

#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<int,double> &a,
                   const pair<int,double> &b)
{
    return a.second<b.second;
}

/// Sorting function for the 2nd element of two vectors
bool sortbysecdesc(const pair<int,double> &a,
                   const pair<int,double> &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<L;l++)
    {
        for (int k=0;k<K;k++)
        {
            for (int d=0;d<Dim;d++)
            {
                value = distribution0(gen);
                for (int t=0;t<T;t++)
                    {
                    if (weights == NULL)
                    {
                        hashFunctions[l][k][t][d] = value;
                    } else {
                        hashFunctions[l][k][t][d] = value;// * weights[t][d];
                    }
                }
            }
        }
    }

    /// Calculate ratios for DTW approximation
    double ** ratios = (double **)malloc(L*sizeof(double*));
    for (int l=0;l<L;l++)
    {
        ratios[l] = (double *)malloc(K*sizeof(double));
        for (int k=0;k<K;k++)
        {
            ratios[l][k] = uniform(a-sd, a+sd);
        }
    }

    /// Create hashSET
    double temp0 = 0;
    double **** hashSET = (double****)malloc(M*sizeof(double***)) ; // Initialize DataD
    for (int m=0;m<M;m++)
    {
        hashSET[m] = (double***)malloc(L*sizeof(double**)) ; // Initialize Data
        for (int l=0;l<L;l++)
        {
            hashSET[m][l] = (double**)malloc(K*sizeof(double*)) ; // Initialize Data
            for (int k=0;k<K;k++)
            {
                hashSET[m][l][k] = (double*)malloc(T*sizeof(double)) ; // Initialize Data
                for (int t=0;t<T;t++)
                {
                    temp0 = 0;
                    for (int d=0;d<Dim;d++) {
                        temp0 += hashFunctions[l][k][t][d] * D[m][t][d];
                    }
                    hashSET[m][l][k][t] = temp0;
                }
            }
        }
    }

    /// Create upper and lower bound on query
    double ** Q_U = (double**)malloc((T)*sizeof(double*));
    double ** Q_L = (double**)malloc((T)*sizeof(double*));
    for (int t=0;t<T;t++)
    {
        Q_U[t] = (double*)malloc((Dim)*sizeof(double));
        Q_L[t] = (double*)malloc((Dim)*sizeof(double));
    }
    double currentMax = -10000;
    double currentMin = 10000;
    for (int t = 0; t < T; t++) {
        for (int d = 0; d < Dim; d++) {
            for (int s = -S; s <= S; s++) {
                if (t + s < 0 || t + s > 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<L;l++)
    {
        hashQ[l] = (double**)malloc(K*sizeof(double*)) ; // Initialize Data
        hashQ_U[l] = (double**)malloc(K*sizeof(double*)) ; // Initialize Data
        hashQ_L[l] = (double**)malloc(K*sizeof(double*)) ; // Initialize Data

        for (int k=0;k<K;k++)
        {
            hashQ[l][k] = (double*)malloc(T*sizeof(double)) ; // Initialize Data
            hashQ_U[l][k] = (double*)malloc(T*sizeof(double)) ; // Initialize Data
            hashQ_L[l][k] = (double*)malloc(T*sizeof(double)) ; // Initialize Data
            for (int t=0;t<T;t++)
            {
                temp1 = 0;
                temp2 = 0;
                temp3 = 0;
                for (int d=0;d<Dim;d++)
                {
                    temp1 += hashFunctions[l][k][t][d] * Q[t][d];
                    temp2 += hashFunctions[l][k][t][d] * Q_U[t][d];
                    temp3 += hashFunctions[l][k][t][d] * Q_L[t][d];
                }
                hashQ[l][k][t] = temp1;
                hashQ_U[l][k][t] = temp2;
                hashQ_L[l][k][t] = temp3;
            }
        }
    }

    /// Generate candidates
    vector<int> v;
    v.reserve(M);
    bool group_collision;
    int hash_collision;
    int projection_collision;
    for (int m=0;m<M;m++)
    {
        group_collision = false;
        for (int l=0;l<L;l++)
        {
            hash_collision = 0;
            for (int k=0;k<K;k++)
            {
                projection_collision = 0;
                for (int t=0;t<T;t++)
                {
                    if (std::abs(hashQ[l][k][t]-hashSET[m][l][k][t]) <= delta/2 ||
                    std::abs(hashQ_U[l][k][t]-hashSET[m][l][k][t]) <= delta/2 ||
                    std::abs(hashQ_L[l][k][t]-hashSET[m][l][k][t]) <= delta/2
                            )
                    {
                        projection_collision++;
                    }
                }
                if (projection_collision>=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<pair<int, double>> 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<L;l++)
    {
        for (int k=0;k<K;k++)
        {
            q_ub += ratios[l][k] * sqrt(T * pow(r, 2.0));
            q_ub += ratios[l][k] * sqrt(T * pow(r, 2.0));
            q_ub += ratios[l][k] * sqrt(T * pow(r, 2.0));
        }
    }
    for (int m=0;m<v.size();m++)
    {
        lb = 0;
        ub = 0;
        for (int l=0;l<L;l++)
        {
            for (int k=0;k<K;k++)
            {
                lb_it = 0;
                ub_it = 0;
                for (int t=0;t<T;t++)
                {
                    lb_it += pow(((2 * std::abs(hashQ[l][k][t]-hashSET[v[m]][l][k][t])) / delta), 2) * pow(r,2);
                    ub_it += pow(((2 * std::abs(hashQ[l][k][t]-hashSET[v[m]][l][k][t])) / delta) + 1, 2) * pow(r,2);
                }
                lb += ratios[l][k] * sqrt(lb_it);
                ub += ratios[l][k] * sqrt(ub_it);
//                lb += sqrt(lb_it);
//                ub += sqrt(ub_it);
            }
        }
        sim.emplace_back(v[m], (q_ub - lb)/(ub - q_lb));
    }
    sort(sim.begin(), sim.end(), sortbysecdesc);
    end_time = clock();
    time_spent = (double)(end_time - begin_time) / CLOCKS_PER_SEC;
    std::cout << "Approximated ranking done in: " << time_spent << " seconds" << std::endl;


    /// Calculate DTW distance on compressed MST data
    vector<pair<int, double>> dtwsim;
    dtwsim.reserve(sim.size());

    double distance = 0;
    for (int m=0;m<sim.size();m++)
    {
        distance = 0;
        for (int l=0;l<L;l++)
        {
            for (int k = 0; k < K; k++)
            {
                distance += dtw(hashSET[sim[m].first][l][k], hashQ[l][k], NULL, T, 0.05 * 120, INF);
            }
        }
        dtwsim.emplace_back(sim[m].first, distance);
    }
    end_time = clock();
    time_spent = (double)(end_time - begin_time) / CLOCKS_PER_SEC;
    std::cout << "DTW sim done in: " << time_spent << std::endl;

    /// Sort and return distances
    sort(dtwsim.begin(), dtwsim.end(), sortbysecasc);
    int count = 0;
    for (int m=0;m<sim.size();m++)
    {
        for (int n=0;n<50;n++)
        {
            if (sim[m].first == dtwsim[n].first) {
                count++;
            }
        }
    }
    candidates = (int*)malloc(dtwsim.size()*sizeof(int));
    distances = (double*)malloc(dtwsim.size()*sizeof(double));
    for (int i=0; i<dtwsim.size(); i++)
    {
        candidates[i] = dtwsim[i].first;
        distances[i] = dtwsim[i].second;
    }
    nrOfCandidates = dtwsim.size();

    /// Clean up
    v.clear();
    sim.clear();
    dtwsim.clear();
    for (int t=0;t<T;t++)
    {
        free(Q_U[t]);
        free(Q_L[t]);
    }
    free(Q_U);
    free(Q_L);
    for (int m=0;m<M;m++)
    {
        for (int l=0;l<L;l++)
        {
            for (int k=0;k<K;k++)
            {
                free(hashSET[m][l][k]);
            }
            free(hashSET[m][l]);
        }
        free(hashSET[m]);
    }
    free(hashSET);
    for (int l=0;l<L;l++)
    {
        for (int k=0;k<K;k++)
        {
            free(hashQ[l][k]);
            free(hashQ_L[l][k]);
            free(hashQ_U[l][k]);
        }
        free(ratios[l]);
        free(hashQ[l]);
        free(hashQ_U[l]);
        free(hashQ_L[l]);
    }
    free(ratios);
    free(hashQ);
    free(hashQ_U);
    free(hashQ_L);
    return 0;
}

int main() {
    /// Run this function for debugging
    srand((unsigned int)time(NULL));
    clock_t begin_time = clock();
    FILE *data;
    int M = 125621;
    int K = 2; //ln((ln(0.5))/ln(1-exp(-2*(0.1)^2*120)))/ln((1-exp(-2*(0.1)^2*120))/(0.5))
    int L = 6; //ln(0.05)/(ln(1-(1-exp(-2(0.1)^2*120))^4))0.0653330009864134
    int T = 120;
    int Dim = 1;

    /// Read data
    double *** D = (double***)malloc(M*sizeof(double**)) ; // Initialize Data
    for (int i=0;i<M;i++)
    {
        D[i] = (double**)malloc((T)*sizeof(double*));
        for (int t=0;t<T;t++)
        {
            D[i][t] = (double*)malloc((Dim)*sizeof(double));
        }
    }
    data = fopen("processed-data","r");
    if( data == NULL )
        exit(2);
    char d[M]; int i; int j = 0; float f;
    while(fscanf(data,"%[^\n]\n",d) != EOF && j< M)
    {
        char *p =  strtok(d," ");
        i = 0;
        while (p != NULL && i<T)
        {
            f = atof(p);
            for (int d=0;d<Dim;d++)
            {
                D[j][i][d] = f;
            }
            i++;
            p = strtok (NULL, " ");
        }
        j++;
    }

    /// Print time
    clock_t end_time = clock();
    double time_spent = (double)(end_time - begin_time) / CLOCKS_PER_SEC;
    std::cout << "Data read in: " << time_spent << std::endl;
    begin_time = clock();

    /// Create query
    double ** Q = (double**)malloc((T)*sizeof(double*));
    for (int t=0;t<T;t++)
    {
        Q[t] = (double*)malloc((Dim)*sizeof(double));
        for (int d=0;d<Dim;d++)
        {
            Q[t][d] = D[80503][t][d];
        }
    }

    int testt = 0;
    while (testt < 5) {
        /// Initialize output variables
        int *candidates;
        int nrOfCandidates;
        double *distances;
        double ****hashFunctions = (double ****) malloc(L * sizeof(double ***));
        for (int l = 0; l < L; l++) {
            hashFunctions[l] = (double ***) malloc(K * sizeof(double **));
            for (int k = 0; k < K; k++) {
                hashFunctions[l][k] = (double **) malloc(T * sizeof(double *));
                for (int t = 0; t < T; t++) {
                    hashFunctions[l][k][t] = (double *) malloc(Dim * sizeof(double));
                }
            }
        }
        double a = 3.1549093441462044;
        double sd = 0.514583453532817;
        double r = 0.21020966674727637;

        lsh(D, M, T, Dim, Q, L, K, r, a, sd, candidates, distances, hashFunctions, NULL, nrOfCandidates);

        cout << "Top-10 candidates: ";
        for (int i=0; i<10; i++)
        {
            cout << candidates[i] << ",";
        }
        cout << endl;

//        for (int l = 0; l < L; l++) {
//            for (int k = 0; k < K; k++) {
//                for (int t = 0; t < T; t++) {
//                    free(hashFunctions[l][k][t]);
//                }
//                free(hashFunctions[l][k]);
//            }
//            free(hashFunctions[l]);
//        }
//        free(hashFunctions);
//        cout << endl;
//        free(candidates);
//        free(distances);

        testt++;
    }


    end_time = clock();
    time_spent = (double)(end_time - begin_time) / CLOCKS_PER_SEC;
    std::cout << "Query done in: " << time_spent << std::endl;

    return 0;
}