#include <iostream>
#include <vector>
#include <omp.h>
#include <stdint.h>
#include <cstdlib>
#include <algorithm>
#include <ctime>
#include <chrono>
#include <fstream>
#include <iostream>


uint64_t n = 30000000;
uint64_t c = 15;
uint64_t m = 60;
uint64_t N = 1;
uint64_t proc_num = 4;

using namespace std;

vector< vector<double> > U(n);
vector< vector<double> > C(c);
vector< vector<uint8_t> > M(n);

static inline uint64_t rdtsc(void) {
    union {
        uint64_t tsc_64;
        struct {
            uint32_t lo_32;
            uint32_t hi_32;
        };
    } tsc;

    asm volatile("rdtsc" :
                    "=a" (tsc.lo_32),
                    "=d" (tsc.hi_32));
    return tsc.tsc_64;
}

void initZerosDoubles(vector< vector<double> > &V, int size_1, int size_2)
{
    for (uint64_t i=0; i < size_1; ++i)
    {
        V[i] = vector<double>(size_2, 0.0);
    }
}

void initZerosInt8(vector< vector<uint8_t> > &V, int size_1, int size_2)
{
    for (uint64_t i=0; i < size_1; ++i)
    {
        V[i] = vector<uint8_t>(size_2, 0.0);
    }
}

double f(double fMin, double fMax)
{
    double res = (double)rand() / RAND_MAX;
    return fMin + res * (fMax - fMin);
}
double fRand()
{
    return f(0, 1);
}

void generateMatrix()
{
    for (uint64_t i=0; i < U.size(); ++i)
    {
        vector<double> v(m);
        generate(v.begin(), v.end(), fRand);

        U[i] = v;
    }

}

void generateCenter()
{
    vector<uint64_t> v(U.size());
    iota(v.begin(), v.end(), 0);
    random_shuffle(v.begin(), v.end());
    for(uint64_t i=0; i < c; ++i)
    {
        C[i] = U[v[i]];
    }
}

double l2_norm(vector<double> const& u1, vector<double> const& u2, int size)
{
    double accum = 0.;
    for (uint64_t i=0; i < size; ++i)
    {
        accum += (u2[i] - u1[i]) * (u2[i] - u1[i]);
    }

    return accum;
}

void sum(vector<double> & u1, vector<double> const& u2, int size)
{
    for (uint64_t i=0; i < size; ++i)
    {
        u1[i] += u2[i];
    }
}

double mainCycle()
{

    double J;
    cout << "Starting.." << endl;
    auto start_ts = chrono::steady_clock::now();
    for (int numCycles=0; numCycles < N; ++numCycles)
    {
        J = 0;
        initZerosInt8(M, n, c);
        vector< vector<double> > C_new[proc_num];
        vector<uint64_t> c_newCount[proc_num];
        for (uint64_t i=0; i < proc_num; ++i)
        {
            C_new[i] = vector< vector<double> >(c);
            initZerosDoubles(C_new[i], c, m);
            c_newCount[i] = vector<uint64_t>(c, 0);
        }

        start_ts = chrono::steady_clock::now();
        #pragma omp parallel for num_threads(proc_num)
        for (uint64_t i=0; i < n; ++i)
        {
            double d_min = l2_norm(U[i], C[0], m);
            uint64_t j_min = 0;
            for (uint64_t j=0; j < c; ++j)
            {
                double d = l2_norm(U[i], C[j], m);
                if (d < d_min)
                {
                    d_min = d;
                    j_min = j;
                }

            }
            J += d_min;
            M[i][j_min] = 1;

            sum(C_new[i % proc_num][j_min], U[i], m);
            c_newCount[i % proc_num][j_min] += 1;
        }

        for (uint64_t i=0; i < c; ++i)
        {
            for (uint64_t j=0; j < m; ++j)
            {
                double c_new = 0;
                uint64_t c_new_count = 0;
                for (uint64_t p=0; p < proc_num; ++p)
                {
                    c_new += C_new[p][i][j];
                    c_new_count += c_newCount[p][i];
                }
                C[i][j] = c_new / c_new_count;
            }
        }

    }
    auto finish_ts = chrono::steady_clock::now();
    auto duration_ts = chrono::duration<double> (finish_ts - start_ts).count();
   //for (uint64_t i=0; i < c; ++i)
   //{
   //    for (uint64_t j=0; j < m; ++j)
   //    {
   //        cout << C[i][j] << " ";
   //    }
   //    cout << endl;
   //}


    double flops = 3*n*(m*c+1) / duration_ts  / 1000000000.0;
    cout << "n_obj: " << n << ", n_class: " << c << ", n_feat: " << m << ", n_proc: " << proc_num << ", n_feat: " << duration_ts << " s, GFLOPS: " << flops  << endl;

    return flops;
}

int main()
{
    std::srand(unsigned(std::time(0)));
    generateMatrix();


    ofstream outfile;
    outfile.open("res.txt", ios::out | ios::trunc );
    outfile << "n= " << n << ", c= " << c << endl;

    double cur_flops;

    for (proc_num = 1; proc_num <= 48; proc_num +=2)
    {
        for (m = 2; m <= 60; m+=4)
        {
            generateCenter();
            cur_flops = mainCycle();
            //m << ", " << proc_num << ", " << duration_ts << ", " << cur_flops << endl;
            outfile << cur_flops << ", ";
        }
        outfile << endl;
    }
}